63 #define PEAK_WINDOW 1024 134 Error(
"Background",
"function not yet implemented: h=%s, iter=%d, option=%sn" 135 , h->
GetName(), number_of_iterations, option);
144 printf(
"\nNumber of positions = %d\n",
fNPeaks);
184 if (dimension != 3) {
185 Error(
"Search",
"Must be a 3-d histogram");
192 Int_t i, j, k, binx,biny,binz, npeaks;
195 for (i = 0; i < sizex; i++) {
198 for (j = 0; j < sizey; j++) {
201 for (k = 0; k < sizez; k++)
206 npeaks =
SearchHighRes((
const Double_t***)source, dest, sizex, sizey, sizez, sigma, 100*threshold,
kTRUE, 3,
kFALSE, 3);
211 for (i = 0; i < npeaks; i++) {
219 for (i = 0; i < sizex; i++) {
220 for (j = 0; j < sizey; j++){
221 delete [] source[i][j];
222 delete [] dest[i][j];
230 if (strstr(option,
"goff"))
232 if (!npeaks)
return 0;
279 Int_t numberIterationsX,
280 Int_t numberIterationsY,
281 Int_t numberIterationsZ,
677 Int_t i, j,
x,
y, z, sampling, q1, q2, q3;
678 Double_t a, b,
c, d,
p1,
p2,
p3, p4, p5, p6, p7, p8, s1, s2, s3, s4, s5, s6, s7, s8, s9, s10, s11, s12,
r1,
r2,
r3, r4, r5, r6;
679 if (ssizex <= 0 || ssizey <= 0 || ssizez <= 0)
680 return "Wrong parameters";
681 if (numberIterationsX < 1 || numberIterationsY < 1 || numberIterationsZ < 1)
682 return "Width of Clipping Window Must Be Positive";
683 if (ssizex < 2 * numberIterationsX + 1 || ssizey < 2 * numberIterationsY + 1 || ssizey < 2 * numberIterationsZ + 1)
684 return (
"Too Large Clipping Window");
686 for(i=0;i<ssizex;i++){
687 working_space[i] =
new Double_t*[ssizey];
688 for(j=0;j<ssizey;j++)
689 working_space[i][j]=
new Double_t[ssizez];
695 for (i = 1; i <= sampling; i++) {
697 for (z = q3; z < ssizez - q3; z++) {
698 for (y = q2; y < ssizey - q2; y++) {
699 for (x = q1; x < ssizex - q1; x++) {
700 a = spectrum[
x][
y][z];
701 p1 = spectrum[x + q1][y + q2][z - q3];
702 p2 = spectrum[x - q1][y + q2][z - q3];
703 p3 = spectrum[x + q1][y - q2][z - q3];
704 p4 = spectrum[x - q1][y - q2][z - q3];
705 p5 = spectrum[x + q1][y + q2][z + q3];
706 p6 = spectrum[x - q1][y + q2][z + q3];
707 p7 = spectrum[x + q1][y - q2][z + q3];
708 p8 = spectrum[x - q1][y - q2][z + q3];
709 s1 = spectrum[x + q1][
y ][z - q3];
710 s2 = spectrum[
x ][y + q2][z - q3];
711 s3 = spectrum[x - q1][
y ][z - q3];
712 s4 = spectrum[
x ][y - q2][z - q3];
713 s5 = spectrum[x + q1][
y ][z + q3];
714 s6 = spectrum[
x ][y + q2][z + q3];
715 s7 = spectrum[x - q1][
y ][z + q3];
716 s8 = spectrum[
x ][y - q2][z + q3];
717 s9 = spectrum[x - q1][y + q2][z ];
718 s10 = spectrum[x - q1][y - q2][z ];
719 s11 = spectrum[x + q1][y + q2][z ];
720 s12 = spectrum[x + q1][y - q2][z ];
721 r1 = spectrum[
x ][
y ][z - q3];
722 r2 = spectrum[
x ][
y ][z + q3];
723 r3 = spectrum[x - q1][
y ][z ];
724 r4 = spectrum[x + q1][
y ][z ];
725 r5 = spectrum[
x ][y + q2][z ];
726 r6 = spectrum[
x ][y - q2][z ];
763 s1 = s1 - (p1 +
p3) / 2.0;
764 s2 = s2 - (p1 +
p2) / 2.0;
765 s3 = s3 - (p2 + p4) / 2.0;
766 s4 = s4 - (p3 + p4) / 2.0;
767 s5 = s5 - (p5 + p7) / 2.0;
768 s6 = s6 - (p5 + p6) / 2.0;
769 s7 = s7 - (p6 + p8) / 2.0;
770 s8 = s8 - (p7 + p8) / 2.0;
771 s9 = s9 - (p2 + p6) / 2.0;
772 s10 = s10 - (p4 + p8) / 2.0;
773 s11 = s11 - (p1 + p5) / 2.0;
774 s12 = s12 - (p3 + p7) / 2.0;
775 b = (s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0;
778 b = (s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0;
781 b = (s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0;
784 b = (s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0;
787 b = (s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0;
790 b = (s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0;
793 r1 = r1 - ((s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0);
794 r2 = r2 - ((s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0);
795 r3 = r3 - ((s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0);
796 r4 = r4 - ((s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0);
797 r5 = r5 - ((s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0);
798 r6 = r6 - ((s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0);
799 b = (r1 +
r2) / 2.0 + (r3 + r4) / 2.0 + (r5 + r6) / 2.0 + (s1 + s3 + s5 + s7) / 4.0 + (s2 + s4 + s6 + s8) / 4.0 + (s9 + s10 + s11 + s12) / 4.0 + (p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8) / 8.0;
802 working_space[
x][
y][z] =
a;
806 for (z = q3; z < ssizez - q3; z++) {
807 for (y = q2; y < ssizey - q2; y++) {
808 for (x = q1; x < ssizex - q1; x++) {
809 spectrum[
x][
y][z] = working_space[
x][
y][z];
817 for (i = 1; i <= sampling; i++) {
819 for (z = q3; z < ssizez - q3; z++) {
820 for (y = q2; y < ssizey - q2; y++) {
821 for (x = q1; x < ssizex - q1; x++) {
822 a = spectrum[
x][
y][z];
823 p1 = spectrum[x + q1][y + q2][z - q3];
824 p2 = spectrum[x - q1][y + q2][z - q3];
825 p3 = spectrum[x + q1][y - q2][z - q3];
826 p4 = spectrum[x - q1][y - q2][z - q3];
827 p5 = spectrum[x + q1][y + q2][z + q3];
828 p6 = spectrum[x - q1][y + q2][z + q3];
829 p7 = spectrum[x + q1][y - q2][z + q3];
830 p8 = spectrum[x - q1][y - q2][z + q3];
831 s1 = spectrum[x + q1][
y ][z - q3];
832 s2 = spectrum[
x ][y + q2][z - q3];
833 s3 = spectrum[x - q1][
y ][z - q3];
834 s4 = spectrum[
x ][y - q2][z - q3];
835 s5 = spectrum[x + q1][
y ][z + q3];
836 s6 = spectrum[
x ][y + q2][z + q3];
837 s7 = spectrum[x - q1][
y ][z + q3];
838 s8 = spectrum[
x ][y - q2][z + q3];
839 s9 = spectrum[x - q1][y + q2][z ];
840 s10 = spectrum[x - q1][y - q2][z ];
841 s11 = spectrum[x + q1][y + q2][z ];
842 s12 = spectrum[x + q1][y - q2][z ];
843 r1 = spectrum[
x ][
y ][z - q3];
844 r2 = spectrum[
x ][
y ][z + q3];
845 r3 = spectrum[x - q1][
y ][z ];
846 r4 = spectrum[x + q1][
y ][z ];
847 r5 = spectrum[
x ][y + q2][z ];
848 r6 = spectrum[
x ][y - q2][z ];
849 b=(p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8) / 8 - (s1 + s2 + s3 + s4 + s5 + s6 + s7 + s8 + s9 + s10 + s11 + s12) / 4 + (r1 + r2 + r3 + r4 + r5 + r6) / 2;
850 c = -(s1 + s2 + s3 + s4 + s5 + s6 + s7 + s8 + s9 + s10 + s11 + s12) / 4 + (r1 + r2 + r3 + r4 + r5 + r6) / 2;
851 d = -(p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8) / 8 + (s1 + s2 + s3 + s4 + s5 + s6 + s7 + s8 + s9 + s10 + s11 + s12) / 12;
852 if(b < a && b >= 0 && c >=0 && d >= 0)
854 working_space[
x][
y][z] =
a;
858 for (z = q3; z < ssizez - q3; z++) {
859 for (y = q2; y < ssizey - q2; y++) {
860 for (x = q1; x < ssizex - q1; x++) {
861 spectrum[
x][
y][z] = working_space[
x][
y][z];
871 for (i = sampling; i >= 1; i--) {
873 for (z = q3; z < ssizez - q3; z++) {
874 for (y = q2; y < ssizey - q2; y++) {
875 for (x = q1; x < ssizex - q1; x++) {
876 a = spectrum[
x][
y][z];
877 p1 = spectrum[x + q1][y + q2][z - q3];
878 p2 = spectrum[x - q1][y + q2][z - q3];
879 p3 = spectrum[x + q1][y - q2][z - q3];
880 p4 = spectrum[x - q1][y - q2][z - q3];
881 p5 = spectrum[x + q1][y + q2][z + q3];
882 p6 = spectrum[x - q1][y + q2][z + q3];
883 p7 = spectrum[x + q1][y - q2][z + q3];
884 p8 = spectrum[x - q1][y - q2][z + q3];
885 s1 = spectrum[x + q1][
y ][z - q3];
886 s2 = spectrum[
x ][y + q2][z - q3];
887 s3 = spectrum[x - q1][
y ][z - q3];
888 s4 = spectrum[
x ][y - q2][z - q3];
889 s5 = spectrum[x + q1][
y ][z + q3];
890 s6 = spectrum[
x ][y + q2][z + q3];
891 s7 = spectrum[x - q1][
y ][z + q3];
892 s8 = spectrum[
x ][y - q2][z + q3];
893 s9 = spectrum[x - q1][y + q2][z ];
894 s10 = spectrum[x - q1][y - q2][z ];
895 s11 = spectrum[x + q1][y + q2][z ];
896 s12 = spectrum[x + q1][y - q2][z ];
897 r1 = spectrum[
x ][
y ][z - q3];
898 r2 = spectrum[
x ][
y ][z + q3];
899 r3 = spectrum[x - q1][
y ][z ];
900 r4 = spectrum[x + q1][
y ][z ];
901 r5 = spectrum[
x ][y + q2][z ];
902 r6 = spectrum[
x ][y - q2][z ];
939 s1 = s1 - (p1 +
p3) / 2.0;
940 s2 = s2 - (p1 +
p2) / 2.0;
941 s3 = s3 - (p2 + p4) / 2.0;
942 s4 = s4 - (p3 + p4) / 2.0;
943 s5 = s5 - (p5 + p7) / 2.0;
944 s6 = s6 - (p5 + p6) / 2.0;
945 s7 = s7 - (p6 + p8) / 2.0;
946 s8 = s8 - (p7 + p8) / 2.0;
947 s9 = s9 - (p2 + p6) / 2.0;
948 s10 = s10 - (p4 + p8) / 2.0;
949 s11 = s11 - (p1 + p5) / 2.0;
950 s12 = s12 - (p3 + p7) / 2.0;
951 b = (s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0;
954 b = (s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0;
957 b = (s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0;
960 b = (s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0;
963 b = (s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0;
966 b = (s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0;
969 r1 = r1 - ((s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0);
970 r2 = r2 - ((s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0);
971 r3 = r3 - ((s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0);
972 r4 = r4 - ((s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0);
973 r5 = r5 - ((s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0);
974 r6 = r6 - ((s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0);
975 b = (r1 +
r2) / 2.0 + (r3 + r4) / 2.0 + (r5 + r6) / 2.0 + (s1 + s3 + s5 + s7) / 4.0 + (s2 + s4 + s6 + s8) / 4.0 + (s9 + s10 + s11 + s12) / 4.0 + (p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8) / 8.0;
978 working_space[
x][
y][z] =
a;
982 for (z = q3; z < ssizez - q3; z++) {
983 for (y = q2; y < ssizey - q2; y++) {
984 for (x = q1; x < ssizex - q1; x++) {
985 spectrum[
x][
y][z] = working_space[
x][
y][z];
993 for (i = sampling; i >= 1; i--) {
995 for (z = q3; z < ssizez - q3; z++) {
996 for (y = q2; y < ssizey - q2; y++) {
997 for (x = q1; x < ssizex - q1; x++) {
998 a = spectrum[
x][
y][z];
999 p1 = spectrum[x + q1][y + q2][z - q3];
1000 p2 = spectrum[x - q1][y + q2][z - q3];
1001 p3 = spectrum[x + q1][y - q2][z - q3];
1002 p4 = spectrum[x - q1][y - q2][z - q3];
1003 p5 = spectrum[x + q1][y + q2][z + q3];
1004 p6 = spectrum[x - q1][y + q2][z + q3];
1005 p7 = spectrum[x + q1][y - q2][z + q3];
1006 p8 = spectrum[x - q1][y - q2][z + q3];
1007 s1 = spectrum[x + q1][
y ][z - q3];
1008 s2 = spectrum[
x ][y + q2][z - q3];
1009 s3 = spectrum[x - q1][
y ][z - q3];
1010 s4 = spectrum[
x ][y - q2][z - q3];
1011 s5 = spectrum[x + q1][
y ][z + q3];
1012 s6 = spectrum[
x ][y + q2][z + q3];
1013 s7 = spectrum[x - q1][
y ][z + q3];
1014 s8 = spectrum[
x ][y - q2][z + q3];
1015 s9 = spectrum[x - q1][y + q2][z ];
1016 s10 = spectrum[x - q1][y - q2][z ];
1017 s11 = spectrum[x + q1][y + q2][z ];
1018 s12 = spectrum[x + q1][y - q2][z ];
1019 r1 = spectrum[
x ][
y ][z - q3];
1020 r2 = spectrum[
x ][
y ][z + q3];
1021 r3 = spectrum[x - q1][
y ][z ];
1022 r4 = spectrum[x + q1][
y ][z ];
1023 r5 = spectrum[
x ][y + q2][z ];
1024 r6 = spectrum[
x ][y - q2][z ];
1025 b = (p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8) / 8 - (s1 + s2 + s3 + s4 + s5 + s6 + s7 + s8 + s9 + s10 + s11 + s12) / 4 + (r1 + r2 + r3 + r4 + r5 + r6) / 2;
1026 c = -(s1 + s2 + s3 + s4 + s5 + s6 + s7 + s8 + s9 + s10 + s11 + s12) / 4+(r1 + r2 + r3 + r4 + r5 + r6) / 2;
1027 d = -(p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8)/8 + (s1 + s2 + s3 + s4 + s5 + s6 + s7 + s8 + s9 + s10 + s11 + s12) / 12;
1028 if(b < a && b >= 0 && c >=0 && d >= 0)
1030 working_space[
x][
y][z] =
a;
1034 for (z = q3; z < ssizez - q3; z++) {
1035 for (y = q2; y < ssizey - q2; y++) {
1036 for (x = q1; x < ssizex - q1; x++) {
1037 spectrum[
x][
y][z] = working_space[
x][
y][z];
1044 for(i = 0;i < ssizex; i++){
1045 for(j = 0;j < ssizey; j++)
1046 delete[] working_space[i][j];
1047 delete[] working_space[i];
1049 delete[] working_space;
1306 Double_t nom,nip,nim,sp,sm,spx,smx,spy,smy,spz,smz,plocha=0;
1308 return "Averaging Window must be positive";
1310 for(i = 0;i < ssizex; i++){
1311 working_space[i] =
new Double_t*[ssizey];
1312 for(j = 0;j < ssizey; j++)
1313 working_space[i][j] =
new Double_t[ssizez];
1321 for(i = 0,maxch = 0;i < ssizex; i++){
1322 for(j = 0;j < ssizey; j++){
1323 for(k = 0;k < ssizez; k++){
1324 working_space[i][j][k] = 0;
1325 if(maxch < source[i][j][k])
1326 maxch = source[i][j][k];
1327 plocha += source[i][j][k];
1332 delete [] working_space;
1337 working_space[
xmin][
ymin][zmin] = 1;
1338 for(i = xmin;i <
xmax; i++){
1339 nip = source[i][
ymin][zmin] / maxch;
1340 nim = source[i + 1][
ymin][zmin] / maxch;
1342 for(l = 1;l <= averWindow; l++){
1344 a = source[
xmax][
ymin][zmin] / maxch;
1347 a = source[i +
l][
ymin][zmin] / maxch;
1359 if(i - l + 1 < xmin)
1360 a = source[
xmin][
ymin][zmin] / maxch;
1363 a = source[i - l + 1][
ymin][zmin] / maxch;
1377 a = working_space[i + 1][
ymin][zmin] = a * working_space[i][
ymin][zmin];
1380 for(i = ymin;i <
ymax; i++){
1381 nip = source[
xmin][i][zmin] / maxch;
1382 nim = source[
xmin][i + 1][zmin] / maxch;
1384 for(l = 1;l <= averWindow; l++){
1386 a = source[
xmin][
ymax][zmin] / maxch;
1389 a = source[
xmin][i +
l][zmin] / maxch;
1401 if(i - l + 1 < ymin)
1402 a = source[
xmin][
ymin][zmin] / maxch;
1405 a = source[
xmin][i - l + 1][zmin] / maxch;
1419 a = working_space[
xmin][i + 1][zmin] = a * working_space[
xmin][i][zmin];
1422 for(i = zmin;i < zmax; i++){
1423 nip = source[
xmin][
ymin][i] / maxch;
1424 nim = source[
xmin][
ymin][i + 1] / maxch;
1426 for(l = 1;l <= averWindow; l++){
1428 a = source[
xmin][
ymin][zmax] / maxch;
1443 if(i - l + 1 < zmin)
1444 a = source[
xmin][
ymin][zmin] / maxch;
1447 a = source[
xmin][
ymin][i - l + 1] / maxch;
1463 for(i = xmin;i <
xmax; i++){
1464 for(j = ymin;j <
ymax; j++){
1465 nip = source[i][j + 1][zmin] / maxch;
1466 nim = source[i + 1][j + 1][zmin] / maxch;
1468 for(l = 1;l <= averWindow; l++){
1470 a = source[
xmax][j][zmin] / maxch;
1473 a = source[i +
l][j][zmin] / maxch;
1485 if(i - l + 1 < xmin)
1486 a = source[
xmin][j][zmin] / maxch;
1489 a = source[i - l + 1][j][zmin] / maxch;
1503 nip = source[i + 1][j][zmin] / maxch;
1504 nim = source[i + 1][j + 1][zmin] / maxch;
1505 for(l = 1;l <= averWindow; l++){
1507 a = source[i][
ymax][zmin] / maxch;
1510 a = source[i][j +
l][zmin] / maxch;
1522 if(j - l + 1 < ymin)
1523 a = source[i][
ymin][zmin] / maxch;
1526 a = source[i][j - l + 1][zmin] / maxch;
1539 a = (spx * working_space[i][j + 1][zmin] + spy * working_space[i + 1][j][zmin]) / (smx + smy);
1540 working_space[i + 1][j + 1][zmin] =
a;
1544 for(i = xmin;i <
xmax; i++){
1545 for(j = zmin;j < zmax; j++){
1546 nip = source[i][
ymin][j + 1] / maxch;
1547 nim = source[i + 1][
ymin][j + 1] / maxch;
1549 for(l = 1;l <= averWindow; l++){
1554 a = source[i +
l][
ymin][j] / maxch;
1566 if(i - l + 1 < xmin)
1570 a = source[i - l + 1][
ymin][j] / maxch;
1584 nip = source[i + 1][
ymin][j] / maxch;
1585 nim = source[i + 1][
ymin][j + 1] / maxch;
1586 for(l = 1;l <= averWindow; l++){
1588 a = source[i][
ymin][zmax] / maxch;
1591 a = source[i][
ymin][j +
l] / maxch;
1603 if(j - l + 1 < zmin)
1604 a = source[i][
ymin][zmin] / maxch;
1607 a = source[i][
ymin][j - l + 1] / maxch;
1620 a = (spx * working_space[i][
ymin][j + 1] + spz * working_space[i + 1][
ymin][j]) / (smx + smz);
1621 working_space[i + 1][
ymin][j + 1] =
a;
1625 for(i = ymin;i <
ymax; i++){
1626 for(j = zmin;j < zmax; j++){
1627 nip = source[
xmin][i][j + 1] / maxch;
1628 nim = source[
xmin][i + 1][j + 1] / maxch;
1630 for(l = 1;l <= averWindow; l++){
1635 a = source[
xmin][i +
l][j] / maxch;
1647 if(i - l + 1 < ymin)
1651 a = source[
xmin][i - l + 1][j] / maxch;
1665 nip = source[
xmin][i + 1][j] / maxch;
1666 nim = source[
xmin][i + 1][j + 1] / maxch;
1667 for(l = 1;l <= averWindow; l++){
1669 a = source[
xmin][i][zmax] / maxch;
1672 a = source[
xmin][i][j +
l] / maxch;
1684 if(j - l + 1 < zmin)
1685 a = source[
xmin][i][zmin] / maxch;
1688 a = source[
xmin][i][j - l + 1] / maxch;
1701 a = (spy * working_space[
xmin][i][j + 1] + spz * working_space[
xmin][i + 1][j]) / (smy + smz);
1702 working_space[
xmin][i + 1][j + 1] =
a;
1706 for(i = xmin;i <
xmax; i++){
1707 for(j = ymin;j <
ymax; j++){
1708 for(k = zmin;k < zmax; k++){
1709 nip = source[i][j + 1][k + 1] / maxch;
1710 nim = source[i + 1][j + 1][k + 1] / maxch;
1712 for(l = 1;l <= averWindow; l++){
1714 a = source[
xmax][j][k] / maxch;
1717 a = source[i +
l][j][k] / maxch;
1729 if(i - l + 1 < xmin)
1730 a = source[
xmin][j][k] / maxch;
1733 a = source[i - l + 1][j][k] / maxch;
1747 nip = source[i + 1][j][k + 1] / maxch;
1748 nim = source[i + 1][j + 1][k + 1] / maxch;
1749 for(l = 1;l <= averWindow; l++){
1751 a = source[i][
ymax][k] / maxch;
1754 a = source[i][j +
l][k] / maxch;
1766 if(j - l + 1 < ymin)
1767 a = source[i][
ymin][k] / maxch;
1770 a = source[i][j - l + 1][k] / maxch;
1784 nip = source[i + 1][j + 1][k] / maxch;
1785 nim = source[i + 1][j + 1][k + 1] / maxch;
1786 for(l = 1;l <= averWindow; l++){
1788 a = source[i][j][zmax] / maxch;
1791 a = source[i][j][k +
l] / maxch;
1803 if(j - l + 1 < ymin)
1804 a = source[i][j][zmin] / maxch;
1807 a = source[i][j][k - l + 1] / maxch;
1820 a = (spx * working_space[i][j + 1][k + 1] + spy * working_space[i + 1][j][k + 1] + spz * working_space[i + 1][j + 1][k]) / (smx + smy + smz);
1821 working_space[i + 1][j + 1][k + 1] =
a;
1826 for(i = xmin;i <=
xmax; i++){
1827 for(j = ymin;j <=
ymax; j++){
1828 for(k = zmin;k <= zmax; k++){
1829 working_space[i][j][k] = working_space[i][j][k] / nom;
1833 for(i = 0;i < ssizex; i++){
1834 for(j = 0;j < ssizey; j++){
1835 for(k = 0;k < ssizez; k++){
1836 source[i][j][k] = plocha * working_space[i][j][k];
1840 for(i = 0;i < ssizex; i++){
1841 for(j = 0;j < ssizey; j++)
1842 delete[] working_space[i][j];
1843 delete[] working_space[i];
1845 delete[] working_space;
1870 Int_t numberIterations,
1871 Int_t numberRepetitions,
2357 Int_t i, j, k, lhx, lhy, lhz, i1, i2, i3, j1, j2, j3, k1, k2, k3, lindex, i1min, i1max, i2min, i2max, i3min, i3max, j1min, j1max, j2min, j2max, j3min, j3max, positx = 0, posity = 0, positz = 0, repet;
2358 Double_t lda, ldb, ldc, area, maximum = 0;
2359 if (ssizex <= 0 || ssizey <= 0 || ssizez <= 0)
2360 return "Wrong parameters";
2361 if (numberIterations <= 0)
2362 return "Number of iterations must be positive";
2363 if (numberRepetitions <= 0)
2364 return "Number of repetitions must be positive";
2366 for(i=0;i<ssizex;i++){
2367 working_space[i]=
new Double_t* [ssizey];
2368 for(j=0;j<ssizey;j++)
2369 working_space[i][j]=
new Double_t [5*ssizez];
2372 lhx = -1, lhy = -1, lhz = -1;
2373 for (i = 0; i < ssizex; i++) {
2374 for (j = 0; j < ssizey; j++) {
2375 for (k = 0; k < ssizez; k++) {
2376 lda = resp[i][j][k];
2385 working_space[i][j][k] = lda;
2387 if (lda > maximum) {
2389 positx = i, posity = j, positz = k;
2394 if (lhx == -1 || lhy == -1 || lhz == -1) {
2395 delete [] working_space;
2396 return (
"Zero response data");
2400 for (i3 = 0; i3 < ssizez; i3++) {
2401 for (i2 = 0; i2 < ssizey; i2++) {
2402 for (i1 = 0; i1 < ssizex; i1++) {
2404 for (j3 = 0; j3 <= (lhz - 1); j3++) {
2405 for (j2 = 0; j2 <= (lhy - 1); j2++) {
2406 for (j1 = 0; j1 <= (lhx - 1); j1++) {
2407 k3 = i3 + j3, k2 = i2 + j2, k1 = i1 + j1;
2408 if (k3 >= 0 && k3 < ssizez && k2 >= 0 && k2 < ssizey && k1 >= 0 && k1 < ssizex) {
2409 lda = working_space[j1][j2][j3];
2410 ldb = source[k1][k2][k3];
2411 ldc = ldc + lda * ldb;
2416 working_space[i1][i2][i3 + ssizez] = ldc;
2422 i1min = -(lhx - 1), i1max = lhx - 1;
2423 i2min = -(lhy - 1), i2max = lhy - 1;
2424 i3min = -(lhz - 1), i3max = lhz - 1;
2425 for (i3 = i3min; i3 <= i3max; i3++) {
2426 for (i2 = i2min; i2 <= i2max; i2++) {
2427 for (i1 = i1min; i1 <= i1max; i1++) {
2432 j3max = lhz - 1 - i3;
2433 if (j3max > lhz - 1)
2435 for (j3 = j3min; j3 <= j3max; j3++) {
2439 j2max = lhy - 1 - i2;
2440 if (j2max > lhy - 1)
2442 for (j2 = j2min; j2 <= j2max; j2++) {
2446 j1max = lhx - 1 - i1;
2447 if (j1max > lhx - 1)
2449 for (j1 = j1min; j1 <= j1max; j1++) {
2450 lda = working_space[j1][j2][j3];
2451 if (i1 + j1 < ssizex && i2 + j2 < ssizey)
2452 ldb = working_space[i1 + j1][i2 + j2][i3 + j3];
2455 ldc = ldc + lda * ldb;
2459 working_space[i1 - i1min][i2 - i2min][i3 - i3min + 2 * ssizez ] = ldc;
2465 for (i3 = 0; i3 < ssizez; i3++) {
2466 for (i2 = 0; i2 < ssizey; i2++) {
2467 for (i1 = 0; i1 < ssizex; i1++) {
2468 working_space[i1][i2][i3 + 3 * ssizez] = 1;
2469 working_space[i1][i2][i3 + 4 * ssizez] = 0;
2475 for (repet = 0; repet < numberRepetitions; repet++) {
2477 for (i = 0; i < ssizex; i++) {
2478 for (j = 0; j < ssizey; j++) {
2479 for (k = 0; k < ssizez; k++) {
2480 working_space[i][j][k + 3 * ssizez] =
TMath::Power(working_space[i][j][k + 3 * ssizez],boost);
2485 for (lindex = 0; lindex < numberIterations; lindex++) {
2486 for (i3 = 0; i3 < ssizez; i3++) {
2487 for (i2 = 0; i2 < ssizey; i2++) {
2488 for (i1 = 0; i1 < ssizex; i1++) {
2491 if (j3min > lhz - 1)
2494 j3max = ssizez - i3 - 1;
2495 if (j3max > lhz - 1)
2498 if (j2min > lhy - 1)
2501 j2max = ssizey - i2 - 1;
2502 if (j2max > lhy - 1)
2505 if (j1min > lhx - 1)
2508 j1max = ssizex - i1 - 1;
2509 if (j1max > lhx - 1)
2511 for (j3 = j3min; j3 <= j3max; j3++) {
2512 for (j2 = j2min; j2 <= j2max; j2++) {
2513 for (j1 = j1min; j1 <= j1max; j1++) {
2514 ldc = working_space[j1 - i1min][j2 - i2min][j3 - i3min + 2 * ssizez];
2515 lda = working_space[i1 + j1][i2 + j2][i3 + j3 + 3 * ssizez];
2516 ldb = ldb + lda * ldc;
2520 lda = working_space[i1][i2][i3 + 3 * ssizez];
2521 ldc = working_space[i1][i2][i3 + 1 * ssizez];
2522 if (ldc * lda != 0 && ldb != 0) {
2523 lda = lda * ldc / ldb;
2528 working_space[i1][i2][i3 + 4 * ssizez] = lda;
2532 for (i3 = 0; i3 < ssizez; i3++) {
2533 for (i2 = 0; i2 < ssizey; i2++) {
2534 for (i1 = 0; i1 < ssizex; i1++)
2535 working_space[i1][i2][i3 + 3 * ssizez] = working_space[i1][i2][i3 + 4 * ssizez];
2540 for (i = 0; i < ssizex; i++) {
2541 for (j = 0; j < ssizey; j++){
2542 for (k = 0; k < ssizez; k++)
2543 source[(i + positx) % ssizex][(j + posity) % ssizey][(k + positz) % ssizez] = area * working_space[i][j][k + 3 * ssizez];
2546 delete [] working_space;
2899 Int_t number_of_iterations = (
Int_t)(4 * sigma + 0.5);
2902 Int_t xmin,
xmax,
l,peak_index = 0,sizex_ext=ssizex + 4 * number_of_iterations,sizey_ext = ssizey + 4 * number_of_iterations,sizez_ext = ssizez + 4 * number_of_iterations,shift = 2 * number_of_iterations;
2904 Double_t a,b,maxch,plocha = 0,plocha_markov = 0;
2905 Double_t nom,nip,nim,sp,sm,spx,spy,smx,smy,spz,smz;
2906 Double_t p1,
p2,
p3,p4,p5,p6,p7,p8,s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,
r1,
r2,
r3,r4,r5,r6;
2909 Int_t lhx,lhy,lhz,i1,i2,i3,j1,j2,j3,k1,k2,k3,i1min,i1max,i2min,i2max,i3min,i3max,j1min,j1max,j2min,j2max,j3min,j3max,positx,posity,positz;
2911 Error(
"SearchHighRes",
"Invalid sigma, must be greater than or equal to 1");
2915 if(threshold<=0||threshold>=100){
2916 Error(
"SearchHighRes",
"Invalid threshold, must be positive and less than 100");
2920 j = (
Int_t)(pocet_sigma*sigma+0.5);
2922 Error(
"SearchHighRes",
"Too large sigma");
2926 if (markov ==
true) {
2927 if (averWindow <= 0) {
2928 Error(
"SearchHighRes",
"Averanging window must be positive");
2933 if(backgroundRemove ==
true){
2934 if(sizex_ext < 2 * number_of_iterations + 1 || sizey_ext < 2 * number_of_iterations + 1 || sizez_ext < 2 * number_of_iterations + 1){
2935 Error(
"SearchHighRes",
"Too large clipping window");
2940 i = (
Int_t)(4 * sigma + 0.5);
2943 for(j = 0;j < ssizex + i; j++){
2944 working_space[j] =
new Double_t* [ssizey + i];
2945 for(k = 0;k < ssizey + i; k++)
2946 working_space[j][k] =
new Double_t [5 * (ssizez + i)];
2948 for(k = 0;k < sizez_ext; k++){
2949 for(j = 0;j < sizey_ext; j++){
2950 for(i = 0;i < sizex_ext; i++){
2954 working_space[i][j][k + sizez_ext] = source[0][0][0];
2956 else if(k >= ssizez + shift)
2957 working_space[i][j][k + sizez_ext] = source[0][0][ssizez - 1];
2960 working_space[i][j][k + sizez_ext] = source[0][0][k - shift];
2963 else if(j >= ssizey + shift){
2965 working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][0];
2967 else if(k >= ssizez + shift)
2968 working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][ssizez - 1];
2971 working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][k - shift];
2976 working_space[i][j][k + sizez_ext] = source[0][j - shift][0];
2978 else if(k >= ssizez + shift)
2979 working_space[i][j][k + sizez_ext] = source[0][j - shift][ssizez - 1];
2982 working_space[i][j][k + sizez_ext] = source[0][j - shift][k - shift];
2986 else if(i >= ssizex + shift){
2989 working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][0];
2991 else if(k >= ssizez + shift)
2992 working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][ssizez - 1];
2995 working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][k - shift];
2998 else if(j >= ssizey + shift){
3000 working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][0];
3002 else if(k >= ssizez + shift)
3003 working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][ssizez - 1];
3006 working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][k - shift];
3011 working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][0];
3013 else if(k >= ssizez + shift)
3014 working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][ssizez - 1];
3017 working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][k - shift];
3024 working_space[i][j][k + sizez_ext] = source[i - shift][0][0];
3026 else if(k >= ssizez + shift)
3027 working_space[i][j][k + sizez_ext] = source[i - shift][0][ssizez - 1];
3030 working_space[i][j][k + sizez_ext] = source[i - shift][0][k - shift];
3033 else if(j >= ssizey + shift){
3035 working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][0];
3037 else if(k >= ssizez + shift)
3038 working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][ssizez - 1];
3041 working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][k - shift];
3046 working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][0];
3048 else if(k >= ssizez + shift)
3049 working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][ssizez - 1];
3052 working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][k - shift];
3058 if(backgroundRemove ==
true){
3059 for(i = 1;i <= number_of_iterations; i++){
3060 for(z = i;z < sizez_ext - i; z++){
3061 for(y = i;y < sizey_ext - i; y++){
3062 for(x = i;x < sizex_ext - i; x++){
3063 a = working_space[
x][
y][z + sizez_ext];
3064 p1 = working_space[x + i][y + i][z - i + sizez_ext];
3065 p2 = working_space[x - i][y + i][z - i + sizez_ext];
3066 p3 = working_space[x + i][y - i][z - i + sizez_ext];
3067 p4 = working_space[x - i][y - i][z - i + sizez_ext];
3068 p5 = working_space[x + i][y + i][z + i + sizez_ext];
3069 p6 = working_space[x - i][y + i][z + i + sizez_ext];
3070 p7 = working_space[x + i][y - i][z + i + sizez_ext];
3071 p8 = working_space[x - i][y - i][z + i + sizez_ext];
3072 s1 = working_space[x + i][
y ][z - i + sizez_ext];
3073 s2 = working_space[
x ][y + i][z - i + sizez_ext];
3074 s3 = working_space[x - i][
y ][z - i + sizez_ext];
3075 s4 = working_space[
x ][y - i][z - i + sizez_ext];
3076 s5 = working_space[x + i][
y ][z + i + sizez_ext];
3077 s6 = working_space[
x ][y + i][z + i + sizez_ext];
3078 s7 = working_space[x - i][
y ][z + i + sizez_ext];
3079 s8 = working_space[
x ][y - i][z + i + sizez_ext];
3080 s9 = working_space[x - i][y + i][z + sizez_ext];
3081 s10 = working_space[x - i][y - i][z +sizez_ext];
3082 s11 = working_space[x + i][y + i][z +sizez_ext];
3083 s12 = working_space[x + i][y - i][z +sizez_ext];
3084 r1 = working_space[
x ][
y ][z - i + sizez_ext];
3085 r2 = working_space[
x ][
y ][z + i + sizez_ext];
3086 r3 = working_space[x - i][
y ][z + sizez_ext];
3087 r4 = working_space[x + i][
y ][z + sizez_ext];
3088 r5 = working_space[
x ][y + i][z + sizez_ext];
3089 r6 = working_space[
x ][y - i][z + sizez_ext];
3090 b = (p1 +
p3) / 2.0;
3094 b = (p1 +
p2) / 2.0;
3098 b = (p2 + p4) / 2.0;
3102 b = (p3 + p4) / 2.0;
3106 b = (p5 + p7) / 2.0;
3110 b = (p5 + p6) / 2.0;
3114 b = (p6 + p8) / 2.0;
3118 b = (p7 + p8) / 2.0;
3122 b = (p2 + p6) / 2.0;
3126 b = (p4 + p8) / 2.0;
3130 b = (p1 + p5) / 2.0;
3134 b = (p3 + p7) / 2.0;
3138 s1 = s1 - (p1 +
p3) / 2.0;
3139 s2 = s2 - (p1 +
p2) / 2.0;
3140 s3 = s3 - (p2 + p4) / 2.0;
3141 s4 = s4 - (p3 + p4) / 2.0;
3142 s5 = s5 - (p5 + p7) / 2.0;
3143 s6 = s6 - (p5 + p6) / 2.0;
3144 s7 = s7 - (p6 + p8) / 2.0;
3145 s8 = s8 - (p7 + p8) / 2.0;
3146 s9 = s9 - (p2 + p6) / 2.0;
3147 s10 = s10 - (p4 + p8) / 2.0;
3148 s11 = s11 - (p1 + p5) / 2.0;
3149 s12 = s12 - (p3 + p7) / 2.0;
3150 b = (s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0;
3154 b = (s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0;
3158 b = (s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0;
3162 b = (s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0;
3166 b = (s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0;
3170 b = (s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0;
3174 r1 = r1 - ((s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0);
3175 r2 = r2 - ((s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0);
3176 r3 = r3 - ((s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0);
3177 r4 = r4 - ((s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0);
3178 r5 = r5 - ((s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0);
3179 r6 = r6 - ((s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0);
3180 b = (r1 +
r2) / 2.0 + (r3 + r4) / 2.0 + (r5 + r6) / 2.0 + (s1 + s3 + s5 + s7) / 4.0 + (s2 + s4 + s6 + s8) / 4.0 + (s9 + s10 + s11 + s12) / 4.0 + (p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8) / 8.0;
3184 working_space[
x][
y][z] =
a;
3188 for(z = i;z < sizez_ext - i; z++){
3189 for(y = i;y < sizey_ext - i; y++){
3190 for(x = i;x < sizex_ext - i; x++){
3191 working_space[
x][
y][z + sizez_ext] = working_space[
x][
y][z];
3196 for(k = 0;k < sizez_ext; k++){
3197 for(j = 0;j < sizey_ext; j++){
3198 for(i = 0;i < sizex_ext; i++){
3202 working_space[i][j][k + sizez_ext] = source[0][0][0] - working_space[i][j][k + sizez_ext];
3204 else if(k >= ssizez + shift)
3205 working_space[i][j][k + sizez_ext] = source[0][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
3208 working_space[i][j][k + sizez_ext] = source[0][0][k - shift] - working_space[i][j][k + sizez_ext];
3211 else if(j >= ssizey + shift){
3213 working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
3215 else if(k >= ssizez + shift)
3216 working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
3219 working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
3224 working_space[i][j][k + sizez_ext] = source[0][j - shift][0] - working_space[i][j][k + sizez_ext];
3226 else if(k >= ssizez + shift)
3227 working_space[i][j][k + sizez_ext] = source[0][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
3230 working_space[i][j][k + sizez_ext] = source[0][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
3234 else if(i >= ssizex + shift){
3237 working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][0] - working_space[i][j][k + sizez_ext];
3239 else if(k >= ssizez + shift)
3240 working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
3243 working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][k - shift] - working_space[i][j][k + sizez_ext];
3246 else if(j >= ssizey + shift){
3248 working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
3250 else if(k >= ssizez + shift)
3251 working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
3254 working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
3259 working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][0] - working_space[i][j][k + sizez_ext];
3261 else if(k >= ssizez + shift)
3262 working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
3265 working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
3272 working_space[i][j][k + sizez_ext] = source[i - shift][0][0] - working_space[i][j][k + sizez_ext];
3274 else if(k >= ssizez + shift)
3275 working_space[i][j][k + sizez_ext] = source[i - shift][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
3278 working_space[i][j][k + sizez_ext] = source[i - shift][0][k - shift] - working_space[i][j][k + sizez_ext];
3281 else if(j >= ssizey + shift){
3283 working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
3285 else if(k >= ssizez + shift)
3286 working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
3289 working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
3294 working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][0] - working_space[i][j][k + sizez_ext];
3296 else if(k >= ssizez + shift)
3297 working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
3300 working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
3309 for(i = 0;i < sizex_ext; i++){
3310 for(j = 0;j < sizey_ext; j++){
3311 for(k = 0;k < sizez_ext; k++){
3312 working_space[i][j][k + 2 * sizez_ext] = working_space[i][j][sizez_ext + k];
3313 plocha_markov = plocha_markov + working_space[i][j][sizez_ext + k];
3318 xmax = sizex_ext - 1;
3320 ymax = sizey_ext - 1;
3322 zmax = sizez_ext - 1;
3323 for(i = 0,maxch = 0;i < sizex_ext; i++){
3324 for(j = 0;j < sizey_ext;j++){
3325 for(k = 0;k < sizez_ext;k++){
3326 working_space[i][j][k] = 0;
3327 if(maxch < working_space[i][j][k + 2 * sizez_ext])
3328 maxch = working_space[i][j][k + 2 * sizez_ext];
3330 plocha += working_space[i][j][k + 2 * sizez_ext];
3335 delete [] working_space;
3339 working_space[
xmin][
ymin][zmin] = 1;
3340 for(i = xmin;i <
xmax; i++){
3341 nip = working_space[i][
ymin][zmin + 2 * sizez_ext] / maxch;
3342 nim = working_space[i + 1][
ymin][zmin + 2 * sizez_ext] / maxch;
3344 for(l = 1;l <= averWindow; l++){
3346 a = working_space[xmax][ymin][zmin + 2 * sizez_ext] / maxch;
3349 a = working_space[i +
l][
ymin][zmin + 2 * sizez_ext] / maxch;
3361 if(i - l + 1 < xmin)
3362 a = working_space[
xmin][
ymin][zmin + 2 * sizez_ext] / maxch;
3365 a = working_space[i - l + 1][
ymin][zmin + 2 * sizez_ext] / maxch;
3379 a = working_space[i + 1][
ymin][zmin] = a * working_space[i][
ymin][zmin];
3382 for(i = ymin;i <
ymax; i++){
3383 nip = working_space[
xmin][i][zmin + 2 * sizez_ext] / maxch;
3384 nim = working_space[
xmin][i + 1][zmin + 2 * sizez_ext] / maxch;
3386 for(l = 1;l <= averWindow; l++){
3388 a = working_space[
xmin][
ymax][zmin + 2 * sizez_ext] / maxch;
3391 a = working_space[
xmin][i +
l][zmin + 2 * sizez_ext] / maxch;
3403 if(i - l + 1 < ymin)
3404 a = working_space[
xmin][
ymin][zmin + 2 * sizez_ext] / maxch;
3407 a = working_space[
xmin][i - l + 1][zmin + 2 * sizez_ext] / maxch;
3421 a = working_space[
xmin][i + 1][zmin] = a * working_space[
xmin][i][zmin];
3424 for(i = zmin;i < zmax;i++){
3425 nip = working_space[
xmin][
ymin][i + 2 * sizez_ext] / maxch;
3426 nim = working_space[
xmin][
ymin][i + 1 + 2 * sizez_ext] / maxch;
3428 for(l = 1;l <= averWindow;l++){
3430 a = working_space[
xmin][
ymin][zmax + 2 * sizez_ext] / maxch;
3433 a = working_space[
xmin][
ymin][i + l + 2 * sizez_ext] / maxch;
3445 if(i - l + 1 < zmin)
3446 a = working_space[
xmin][
ymin][zmin + 2 * sizez_ext] / maxch;
3449 a = working_space[
xmin][
ymin][i - l + 1 + 2 * sizez_ext] / maxch;
3466 for(i = xmin;i <
xmax; i++){
3467 for(j = ymin;j <
ymax; j++){
3468 nip = working_space[i][j + 1][zmin + 2 * sizez_ext] / maxch;
3469 nim = working_space[i + 1][j + 1][zmin + 2 * sizez_ext] / maxch;
3471 for(l = 1;l <= averWindow; l++){
3473 a = working_space[
xmax][j][zmin + 2 * sizez_ext] / maxch;
3476 a = working_space[i +
l][j][zmin + 2 * sizez_ext] / maxch;
3488 if(i - l + 1 < xmin)
3489 a = working_space[
xmin][j][zmin + 2 * sizez_ext] / maxch;
3492 a = working_space[i - l + 1][j][zmin + 2 * sizez_ext] / maxch;
3506 nip = working_space[i + 1][j][zmin + 2 * sizez_ext] / maxch;
3507 nim = working_space[i + 1][j + 1][zmin + 2 * sizez_ext] / maxch;
3508 for(l = 1;l <= averWindow; l++){
3510 a = working_space[i][
ymax][zmin + 2 * sizez_ext] / maxch;
3513 a = working_space[i][j +
l][zmin + 2 * sizez_ext] / maxch;
3525 if(j - l + 1 < ymin)
3526 a = working_space[i][
ymin][zmin + 2 * sizez_ext] / maxch;
3529 a = working_space[i][j - l + 1][zmin + 2 * sizez_ext] / maxch;
3542 a = (spx * working_space[i][j + 1][zmin] + spy * working_space[i + 1][j][zmin]) / (smx + smy);
3543 working_space[i + 1][j + 1][zmin] =
a;
3547 for(i = xmin;i <
xmax;i++){
3548 for(j = zmin;j < zmax;j++){
3549 nip = working_space[i][
ymin][j + 1 + 2 * sizez_ext] / maxch;
3550 nim = working_space[i + 1][
ymin][j + 1 + 2 * sizez_ext] / maxch;
3552 for(l = 1;l <= averWindow; l++){
3554 a = working_space[
xmax][
ymin][j + 2 * sizez_ext] / maxch;
3557 a = working_space[i +
l][
ymin][j + 2 * sizez_ext] / maxch;
3569 if(i - l + 1 < xmin)
3570 a = working_space[
xmin][
ymin][j + 2 * sizez_ext] / maxch;
3573 a = working_space[i - l + 1][
ymin][j + 2 * sizez_ext] / maxch;
3587 nip = working_space[i + 1][
ymin][j + 2 * sizez_ext] / maxch;
3588 nim = working_space[i + 1][
ymin][j + 1 + 2 * sizez_ext] / maxch;
3589 for(l = 1;l <= averWindow; l++){
3591 a = working_space[i][
ymin][zmax + 2 * sizez_ext] / maxch;
3594 a = working_space[i][
ymin][j + l + 2 * sizez_ext] / maxch;
3606 if(j - l + 1 < zmin)
3607 a = working_space[i][
ymin][zmin + 2 * sizez_ext] / maxch;
3610 a = working_space[i][
ymin][j - l + 1 + 2 * sizez_ext] / maxch;
3623 a = (spx * working_space[i][
ymin][j + 1] + spz * working_space[i + 1][
ymin][j]) / (smx + smz);
3624 working_space[i + 1][
ymin][j + 1] =
a;
3628 for(i = ymin;i <
ymax;i++){
3629 for(j = zmin;j < zmax;j++){
3630 nip = working_space[
xmin][i][j + 1 + 2 * sizez_ext] / maxch;
3631 nim = working_space[
xmin][i + 1][j + 1 + 2 * sizez_ext] / maxch;
3633 for(l = 1;l <= averWindow; l++){
3635 a = working_space[
xmin][
ymax][j + 2 * sizez_ext] / maxch;
3638 a = working_space[
xmin][i +
l][j + 2 * sizez_ext] / maxch;
3650 if(i - l + 1 < ymin)
3651 a = working_space[
xmin][
ymin][j + 2 * sizez_ext] / maxch;
3654 a = working_space[
xmin][i - l + 1][j + 2 * sizez_ext] / maxch;
3668 nip = working_space[
xmin][i + 1][j + 2 * sizez_ext] / maxch;
3669 nim = working_space[
xmin][i + 1][j + 1 + 2 * sizez_ext] / maxch;
3670 for(l = 1;l <= averWindow; l++){
3672 a = working_space[
xmin][i][zmax + 2 * sizez_ext] / maxch;
3675 a = working_space[
xmin][i][j + l + 2 * sizez_ext] / maxch;
3687 if(j - l + 1 < zmin)
3688 a = working_space[
xmin][i][zmin + 2 * sizez_ext] / maxch;
3691 a = working_space[
xmin][i][j - l + 1 + 2 * sizez_ext] / maxch;
3704 a = (spy * working_space[
xmin][i][j + 1] + spz * working_space[
xmin][i + 1][j]) / (smy + smz);
3705 working_space[
xmin][i + 1][j + 1] =
a;
3709 for(i = xmin;i <
xmax; i++){
3710 for(j = ymin;j <
ymax; j++){
3711 for(k = zmin;k < zmax; k++){
3712 nip = working_space[i][j + 1][k + 1 + 2 * sizez_ext] / maxch;
3713 nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
3715 for(l = 1;l <= averWindow; l++){
3717 a = working_space[
xmax][j][k + 2 * sizez_ext] / maxch;
3720 a = working_space[i +
l][j][k + 2 * sizez_ext] / maxch;
3732 if(i - l + 1 < xmin)
3733 a = working_space[
xmin][j][k + 2 * sizez_ext] / maxch;
3736 a = working_space[i - l + 1][j][k + 2 * sizez_ext] / maxch;
3750 nip = working_space[i + 1][j][k + 1 + 2 * sizez_ext] / maxch;
3751 nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
3752 for(l = 1;l <= averWindow; l++){
3754 a = working_space[i][
ymax][k + 2 * sizez_ext] / maxch;
3757 a = working_space[i][j +
l][k + 2 * sizez_ext] / maxch;
3769 if(j - l + 1 < ymin)
3770 a = working_space[i][
ymin][k + 2 * sizez_ext] / maxch;
3773 a = working_space[i][j - l + 1][k + 2 * sizez_ext] / maxch;
3787 nip = working_space[i + 1][j + 1][k + 2 * sizez_ext] / maxch;
3788 nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
3789 for(l = 1;l <= averWindow; l++ ){
3791 a = working_space[i][j][zmax + 2 * sizez_ext] / maxch;
3794 a = working_space[i][j][k + l + 2 * sizez_ext] / maxch;
3806 if(j - l + 1 < ymin)
3807 a = working_space[i][j][zmin + 2 * sizez_ext] / maxch;
3810 a = working_space[i][j][k - l + 1 + 2 * sizez_ext] / maxch;
3823 a = (spx * working_space[i][j + 1][k + 1] + spy * working_space[i + 1][j][k + 1] + spz * working_space[i + 1][j + 1][k]) / (smx + smy + smz);
3824 working_space[i + 1][j + 1][k + 1] =
a;
3830 for(i = xmin;i <=
xmax; i++){
3831 for(j = ymin;j <=
ymax; j++){
3832 for(k = zmin;k <= zmax; k++){
3833 working_space[i][j][k] = working_space[i][j][k] / nom;
3834 a+=working_space[i][j][k];
3838 for(i = 0;i < sizex_ext; i++){
3839 for(j = 0;j < sizey_ext; j++){
3840 for(k = 0;k < sizez_ext; k++){
3841 working_space[i][j][k + sizez_ext] = working_space[i][j][k] * plocha_markov /
a;
3848 lhx = -1,lhy = -1,lhz = -1;
3849 positx = 0,posity = 0,positz = 0;
3852 for(i = 0;i < sizex_ext; i++){
3853 for(j = 0;j < sizey_ext; j++){
3854 for(k = 0;k < sizez_ext; k++){
3858 lda = (lda * lda + ldb * ldb + ldc * ldc) / (2 * sigma * sigma);
3871 working_space[i][j][k] = lda;
3875 positx = i,posity = j,positz = k;
3881 for(i = 0;i < sizex_ext; i++){
3882 for(j = 0;j < sizey_ext; j++){
3883 for(k = 0;k < sizez_ext; k++){
3884 working_space[i][j][k + 2 * sizez_ext] =
TMath::Abs(working_space[i][j][k + sizez_ext]);
3889 for (i3 = 0; i3 < sizez_ext; i3++) {
3890 for (i2 = 0; i2 < sizey_ext; i2++) {
3891 for (i1 = 0; i1 < sizex_ext; i1++) {
3893 for (j3 = 0; j3 <= (lhz - 1); j3++) {
3894 for (j2 = 0; j2 <= (lhy - 1); j2++) {
3895 for (j1 = 0; j1 <= (lhx - 1); j1++) {
3896 k3 = i3 + j3, k2 = i2 + j2, k1 = i1 + j1;
3897 if (k3 >= 0 && k3 < sizez_ext && k2 >= 0 && k2 < sizey_ext && k1 >= 0 && k1 < sizex_ext) {
3898 lda = working_space[j1][j2][j3];
3899 ldb = working_space[k1][k2][k3+2*sizez_ext];
3900 ldc = ldc + lda * ldb;
3905 working_space[i1][i2][i3 + sizez_ext] = ldc;
3910 i1min = -(lhx - 1), i1max = lhx - 1;
3911 i2min = -(lhy - 1), i2max = lhy - 1;
3912 i3min = -(lhz - 1), i3max = lhz - 1;
3913 for (i3 = i3min; i3 <= i3max; i3++) {
3914 for (i2 = i2min; i2 <= i2max; i2++) {
3915 for (i1 = i1min; i1 <= i1max; i1++) {
3921 j3max = lhz - 1 - i3;
3922 if (j3max > lhz - 1)
3925 for (j3 = j3min; j3 <= j3max; j3++) {
3930 j2max = lhy - 1 - i2;
3931 if (j2max > lhy - 1)
3934 for (j2 = j2min; j2 <= j2max; j2++) {
3939 j1max = lhx - 1 - i1;
3940 if (j1max > lhx - 1)
3943 for (j1 = j1min; j1 <= j1max; j1++) {
3944 lda = working_space[j1][j2][j3];
3945 if (i1 + j1 < sizex_ext && i2 + j2 < sizey_ext)
3946 ldb = working_space[i1 + j1][i2 + j2][i3 + j3];
3951 ldc = ldc + lda * ldb;
3955 working_space[i1 - i1min][i2 - i2min][i3 - i3min + 2 * sizez_ext ] = ldc;
3960 for (i3 = 0; i3 < sizez_ext; i3++) {
3961 for (i2 = 0; i2 < sizey_ext; i2++) {
3962 for (i1 = 0; i1 < sizex_ext; i1++) {
3963 working_space[i1][i2][i3 + 3 * sizez_ext] = 1;
3964 working_space[i1][i2][i3 + 4 * sizez_ext] = 0;
3970 for (lindex=0;lindex<deconIterations;lindex++){
3971 for (i3 = 0; i3 < sizez_ext; i3++) {
3972 for (i2 = 0; i2 < sizey_ext; i2++) {
3973 for (i1 = 0; i1 < sizex_ext; i1++) {
3974 if (
TMath::Abs(working_space[i1][i2][i3 + 3 * sizez_ext])>1e-6 &&
TMath::Abs(working_space[i1][i2][i3 + 1 * sizez_ext])>1e-6){
3977 if (j3min > lhz - 1)
3981 j3max = sizez_ext - i3 - 1;
3982 if (j3max > lhz - 1)
3986 if (j2min > lhy - 1)
3990 j2max = sizey_ext - i2 - 1;
3991 if (j2max > lhy - 1)
3995 if (j1min > lhx - 1)
3999 j1max = sizex_ext - i1 - 1;
4000 if (j1max > lhx - 1)
4003 for (j3 = j3min; j3 <= j3max; j3++) {
4004 for (j2 = j2min; j2 <= j2max; j2++) {
4005 for (j1 = j1min; j1 <= j1max; j1++) {
4006 ldc = working_space[j1 - i1min][j2 - i2min][j3 - i3min + 2 * sizez_ext];
4007 lda = working_space[i1 + j1][i2 + j2][i3 + j3 + 3 * sizez_ext];
4008 ldb = ldb + lda * ldc;
4012 lda = working_space[i1][i2][i3 + 3 * sizez_ext];
4013 ldc = working_space[i1][i2][i3 + 1 * sizez_ext];
4014 if (ldc * lda != 0 && ldb != 0) {
4015 lda = lda * ldc / ldb;
4020 working_space[i1][i2][i3 + 4 * sizez_ext] = lda;
4025 for (i3 = 0; i3 < sizez_ext; i3++) {
4026 for (i2 = 0; i2 < sizey_ext; i2++) {
4027 for (i1 = 0; i1 < sizex_ext; i1++)
4028 working_space[i1][i2][i3 + 3 * sizez_ext] = working_space[i1][i2][i3 + 4 * sizez_ext];
4034 for(i = 0;i < sizex_ext; i++){
4035 for(j = 0;j < sizey_ext; j++){
4036 for(k = 0;k < sizez_ext; k++){
4037 working_space[(i + positx) % sizex_ext][(j + posity) % sizey_ext][(k + positz) % sizez_ext] = area * working_space[i][j][k + 3 * sizez_ext];
4038 if(maximum < working_space[(i + positx) % sizex_ext][(j + posity) % sizey_ext][(k + positz) % sizez_ext])
4039 maximum = working_space[(i + positx) % sizex_ext][(j + posity) % sizey_ext][(k + positz) % sizez_ext];
4044 for(i = 1;i < sizex_ext - 1; i++){
4045 for(j = 1;j < sizey_ext - 1; j++){
4046 for(l = 1;l < sizez_ext - 1; l++){
4047 a = working_space[i][j][
l];
4048 if(a > working_space[i][j][l - 1] && a > working_space[i - 1][j][l - 1] && a > working_space[i - 1][j - 1][l - 1] && a > working_space[i][j - 1][l - 1] && a > working_space[i + 1][j - 1][l - 1] && a > working_space[i + 1][j][l - 1] && a > working_space[i + 1][j + 1][l - 1] && a > working_space[i][j + 1][l - 1] && a > working_space[i - 1][j + 1][l - 1] && a > working_space[i - 1][j][l] && a > working_space[i - 1][j - 1][l] && a > working_space[i][j - 1][l] && a > working_space[i + 1][j - 1][l] && a > working_space[i + 1][j][l] && a > working_space[i + 1][j + 1][l] && a > working_space[i][j + 1][l] && a > working_space[i - 1][j + 1][l] && a > working_space[i][j][l + 1] && a > working_space[i - 1][j][l + 1] && a > working_space[i - 1][j - 1][l + 1] && a > working_space[i][j - 1][l + 1] && a > working_space[i + 1][j - 1][l + 1] && a > working_space[i + 1][j][l + 1] && a > working_space[i + 1][j + 1][l + 1] && a > working_space[i][j + 1][l + 1] && a > working_space[i - 1][j + 1][l + 1]){
4049 if(i >= shift && i < ssizex + shift && j >= shift && j < ssizey + shift && l >= shift && l < ssizez + shift){
4050 if(working_space[i][j][l] > threshold * maximum / 100.0){
4052 for(k = i - 1,a = 0,b = 0;k <= i + 1; k++){
4053 a += (
Double_t)(k - shift) * working_space[k][j][
l];
4054 b += working_space[k][j][
l];
4057 for(k = j - 1,a = 0,b = 0;k <= j + 1; k++){
4058 a += (
Double_t)(k - shift) * working_space[i][k][
l];
4059 b += working_space[i][k][
l];
4062 for(k = l - 1,a = 0,b = 0;k <= l + 1; k++){
4063 a += (
Double_t)(k - shift) * working_space[i][j][k];
4064 b += working_space[i][j][k];
4075 for(i = 0;i < ssizex; i++){
4076 for(j = 0;j < ssizey; j++){
4077 for(k = 0;k < ssizez; k++){
4078 dest[i][j][k] = working_space[i + shift][j + shift][k + shift];
4082 k = (
Int_t)(4 * sigma + 0.5);
4084 for(i = 0;i < ssizex + k; i++){
4085 for(j = 0;j < ssizey + k; j++)
4086 delete[] working_space[i][j];
4087 delete[] working_space[i];
4089 delete[] working_space;
4123 Int_t i,j,k,
l,li,lj,lk,lmin,lmax,
xmin,
xmax,
ymin,
ymax,zmin,zmax;
4124 Double_t maxch,plocha = 0,plocha_markov = 0;
4125 Double_t nom,nip,nim,sp,sm,spx,spy,smx,smy,spz,smz;
4126 Double_t norma,val,val1,val2,val3,val4,val5,val6,val7,val8,val9,val10,val11,val12,val13,val14,val15,val16,val17,val18,val19,val20,val21,val22,val23,val24,val25,val26;
4129 Double_t p1,
p2,
p3,p4,p5,p6,p7,p8,s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,
r1,
r2,
r3,r4,r5,r6;
4131 Int_t number_of_iterations=(
Int_t)(4 * sigma + 0.5);
4132 Int_t sizex_ext=ssizex + 4 * number_of_iterations,sizey_ext = ssizey + 4 * number_of_iterations,sizez_ext = ssizez + 4 * number_of_iterations,shift = 2 * number_of_iterations;
4135 Error(
"SearchFast",
"Invalid sigma, must be greater than or equal to 1");
4139 if(threshold<=0||threshold>=100){
4140 Error(
"SearchFast",
"Invalid threshold, must be positive and less than 100");
4144 j = (
Int_t)(pocet_sigma*sigma+0.5);
4146 Error(
"SearchFast",
"Too large sigma");
4150 if (markov ==
true) {
4151 if (averWindow <= 0) {
4152 Error(
"SearchFast",
"Averanging window must be positive");
4157 if(sizex_ext < 2 * number_of_iterations + 1 || sizey_ext < 2 * number_of_iterations + 1 || sizez_ext < 2 * number_of_iterations + 1){
4158 Error(
"SearchFast",
"Too large clipping window");
4162 i = (
Int_t)(4 * sigma + 0.5);
4165 for(j = 0;j < ssizex + i; j++){
4166 working_space[j] =
new Double_t* [ssizey + i];
4167 for(k = 0;k < ssizey + i; k++)
4168 working_space[j][k] =
new Double_t [4 * (ssizez + i)];
4171 for(k = 0;k < sizez_ext; k++){
4172 for(j = 0;j < sizey_ext; j++){
4173 for(i = 0;i < sizex_ext; i++){
4177 working_space[i][j][k + sizez_ext] = source[0][0][0];
4179 else if(k >= ssizez + shift)
4180 working_space[i][j][k + sizez_ext] = source[0][0][ssizez - 1];
4183 working_space[i][j][k + sizez_ext] = source[0][0][k - shift];
4186 else if(j >= ssizey + shift){
4188 working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][0];
4190 else if(k >= ssizez + shift)
4191 working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][ssizez - 1];
4194 working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][k - shift];
4199 working_space[i][j][k + sizez_ext] = source[0][j - shift][0];
4201 else if(k >= ssizez + shift)
4202 working_space[i][j][k + sizez_ext] = source[0][j - shift][ssizez - 1];
4205 working_space[i][j][k + sizez_ext] = source[0][j - shift][k - shift];
4209 else if(i >= ssizex + shift){
4212 working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][0];
4214 else if(k >= ssizez + shift)
4215 working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][ssizez - 1];
4218 working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][k - shift];
4221 else if(j >= ssizey + shift){
4223 working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][0];
4225 else if(k >= ssizez + shift)
4226 working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][ssizez - 1];
4229 working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][k - shift];
4234 working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][0];
4236 else if(k >= ssizez + shift)
4237 working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][ssizez - 1];
4240 working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][k - shift];
4247 working_space[i][j][k + sizez_ext] = source[i - shift][0][0];
4249 else if(k >= ssizez + shift)
4250 working_space[i][j][k + sizez_ext] = source[i - shift][0][ssizez - 1];
4253 working_space[i][j][k + sizez_ext] = source[i - shift][0][k - shift];
4256 else if(j >= ssizey + shift){
4258 working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][0];
4260 else if(k >= ssizez + shift)
4261 working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][ssizez - 1];
4264 working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][k - shift];
4269 working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][0];
4271 else if(k >= ssizez + shift)
4272 working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][ssizez - 1];
4275 working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][k - shift];
4281 for(i = 1;i <= number_of_iterations; i++){
4282 for(z = i;z < sizez_ext - i; z++){
4283 for(y = i;y < sizey_ext - i; y++){
4284 for(x = i;x < sizex_ext - i; x++){
4285 a = working_space[
x][
y][z + sizez_ext];
4286 p1 = working_space[x + i][y + i][z - i + sizez_ext];
4287 p2 = working_space[x - i][y + i][z - i + sizez_ext];
4288 p3 = working_space[x + i][y - i][z - i + sizez_ext];
4289 p4 = working_space[x - i][y - i][z - i + sizez_ext];
4290 p5 = working_space[x + i][y + i][z + i + sizez_ext];
4291 p6 = working_space[x - i][y + i][z + i + sizez_ext];
4292 p7 = working_space[x + i][y - i][z + i + sizez_ext];
4293 p8 = working_space[x - i][y - i][z + i + sizez_ext];
4294 s1 = working_space[x + i][
y ][z - i + sizez_ext];
4295 s2 = working_space[
x ][y + i][z - i + sizez_ext];
4296 s3 = working_space[x - i][
y ][z - i + sizez_ext];
4297 s4 = working_space[
x ][y - i][z - i + sizez_ext];
4298 s5 = working_space[x + i][
y ][z + i + sizez_ext];
4299 s6 = working_space[
x ][y + i][z + i + sizez_ext];
4300 s7 = working_space[x - i][
y ][z + i + sizez_ext];
4301 s8 = working_space[
x ][y - i][z + i + sizez_ext];
4302 s9 = working_space[x - i][y + i][z + sizez_ext];
4303 s10 = working_space[x - i][y - i][z +sizez_ext];
4304 s11 = working_space[x + i][y + i][z +sizez_ext];
4305 s12 = working_space[x + i][y - i][z +sizez_ext];
4306 r1 = working_space[
x ][
y ][z - i + sizez_ext];
4307 r2 = working_space[
x ][
y ][z + i + sizez_ext];
4308 r3 = working_space[x - i][
y ][z + sizez_ext];
4309 r4 = working_space[x + i][
y ][z + sizez_ext];
4310 r5 = working_space[
x ][y + i][z + sizez_ext];
4311 r6 = working_space[
x ][y - i][z + sizez_ext];
4312 b = (p1 +
p3) / 2.0;
4316 b = (p1 +
p2) / 2.0;
4320 b = (p2 + p4) / 2.0;
4324 b = (p3 + p4) / 2.0;
4328 b = (p5 + p7) / 2.0;
4332 b = (p5 + p6) / 2.0;
4336 b = (p6 + p8) / 2.0;
4340 b = (p7 + p8) / 2.0;
4344 b = (p2 + p6) / 2.0;
4348 b = (p4 + p8) / 2.0;
4352 b = (p1 + p5) / 2.0;
4356 b = (p3 + p7) / 2.0;
4360 s1 = s1 - (p1 +
p3) / 2.0;
4361 s2 = s2 - (p1 +
p2) / 2.0;
4362 s3 = s3 - (p2 + p4) / 2.0;
4363 s4 = s4 - (p3 + p4) / 2.0;
4364 s5 = s5 - (p5 + p7) / 2.0;
4365 s6 = s6 - (p5 + p6) / 2.0;
4366 s7 = s7 - (p6 + p8) / 2.0;
4367 s8 = s8 - (p7 + p8) / 2.0;
4368 s9 = s9 - (p2 + p6) / 2.0;
4369 s10 = s10 - (p4 + p8) / 2.0;
4370 s11 = s11 - (p1 + p5) / 2.0;
4371 s12 = s12 - (p3 + p7) / 2.0;
4372 b = (s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0;
4376 b = (s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0;
4380 b = (s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0;
4384 b = (s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0;
4388 b = (s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0;
4392 b = (s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0;
4396 r1 = r1 - ((s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0);
4397 r2 = r2 - ((s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0);
4398 r3 = r3 - ((s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0);
4399 r4 = r4 - ((s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0);
4400 r5 = r5 - ((s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0);
4401 r6 = r6 - ((s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0);
4402 b = (r1 +
r2) / 2.0 + (r3 + r4) / 2.0 + (r5 + r6) / 2.0 + (s1 + s3 + s5 + s7) / 4.0 + (s2 + s4 + s6 + s8) / 4.0 + (s9 + s10 + s11 + s12) / 4.0 + (p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8) / 8.0;
4406 working_space[
x][
y][z] =
a;
4410 for(z = i;z < sizez_ext - i; z++){
4411 for(y = i;y < sizey_ext - i; y++){
4412 for(x = i;x < sizex_ext - i; x++){
4413 working_space[
x][
y][z + sizez_ext] = working_space[
x][
y][z];
4418 for(k = 0;k < sizez_ext; k++){
4419 for(j = 0;j < sizey_ext; j++){
4420 for(i = 0;i < sizex_ext; i++){
4424 working_space[i][j][k + 3 * sizez_ext] = source[0][0][0] - working_space[i][j][k + sizez_ext];
4426 else if(k >= ssizez + shift)
4427 working_space[i][j][k + 3 * sizez_ext] = source[0][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
4430 working_space[i][j][k + 3 * sizez_ext] = source[0][0][k - shift] - working_space[i][j][k + sizez_ext];
4433 else if(j >= ssizey + shift){
4435 working_space[i][j][k + 3 * sizez_ext] = source[0][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
4437 else if(k >= ssizez + shift)
4438 working_space[i][j][k + 3 * sizez_ext] = source[0][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
4441 working_space[i][j][k + 3 * sizez_ext] = source[0][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
4446 working_space[i][j][k + 3 * sizez_ext] = source[0][j - shift][0] - working_space[i][j][k + sizez_ext];
4448 else if(k >= ssizez + shift)
4449 working_space[i][j][k + 3 * sizez_ext] = source[0][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
4452 working_space[i][j][k + 3 * sizez_ext] = source[0][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
4456 else if(i >= ssizex + shift){
4459 working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][0][0] - working_space[i][j][k + sizez_ext];
4461 else if(k >= ssizez + shift)
4462 working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
4465 working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][0][k - shift] - working_space[i][j][k + sizez_ext];
4468 else if(j >= ssizey + shift){
4470 working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
4472 else if(k >= ssizez + shift)
4473 working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
4476 working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
4481 working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][j - shift][0] - working_space[i][j][k + sizez_ext];
4483 else if(k >= ssizez + shift)
4484 working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
4487 working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
4494 working_space[i][j][k + 3 * sizez_ext] = source[i - shift][0][0] - working_space[i][j][k + sizez_ext];
4496 else if(k >= ssizez + shift)
4497 working_space[i][j][k + 3 * sizez_ext] = source[i - shift][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
4500 working_space[i][j][k + 3 * sizez_ext] = source[i - shift][0][k - shift] - working_space[i][j][k + sizez_ext];
4503 else if(j >= ssizey + shift){
4505 working_space[i][j][k + 3 * sizez_ext] = source[i - shift][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
4507 else if(k >= ssizez + shift)
4508 working_space[i][j][k + 3 * sizez_ext] = source[i - shift][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
4511 working_space[i][j][k + 3 * sizez_ext] = source[i - shift][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
4516 working_space[i][j][k + 3 * sizez_ext] = source[i - shift][j - shift][0] - working_space[i][j][k + sizez_ext];
4518 else if(k >= ssizez + shift)
4519 working_space[i][j][k + 3 * sizez_ext] = source[i - shift][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
4522 working_space[i][j][k + 3 * sizez_ext] = source[i - shift][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
4529 for(i = 0;i < sizex_ext; i++){
4530 for(j = 0;j < sizey_ext; j++){
4531 for(k = 0;k < sizez_ext; k++){
4532 if(i >= shift && i < ssizex + shift && j >= shift && j < ssizey + shift && k >= shift && k < ssizez + shift){
4533 working_space[i][j][k + 2 * sizez_ext] = source[i - shift][j - shift][k - shift];
4534 plocha_markov = plocha_markov + source[i - shift][j - shift][k - shift];
4537 working_space[i][j][k + 2 * sizez_ext] = 0;
4544 xmax = sizex_ext - 1;
4546 ymax = sizey_ext - 1;
4548 zmax = sizez_ext - 1;
4549 for(i = 0,maxch = 0;i < sizex_ext; i++){
4550 for(j = 0;j < sizey_ext;j++){
4551 for(k = 0;k < sizez_ext;k++){
4552 working_space[i][j][k] = 0;
4553 if(maxch < working_space[i][j][k + 2 * sizez_ext])
4554 maxch = working_space[i][j][k + 2 * sizez_ext];
4556 plocha += working_space[i][j][k + 2 * sizez_ext];
4561 delete [] working_space;
4566 working_space[
xmin][
ymin][zmin] = 1;
4567 for(i = xmin;i <
xmax; i++){
4568 nip = working_space[i][
ymin][zmin + 2 * sizez_ext] / maxch;
4569 nim = working_space[i + 1][
ymin][zmin + 2 * sizez_ext] / maxch;
4571 for(l = 1;l <= averWindow; l++){
4573 a = working_space[xmax][ymin][zmin + 2 * sizez_ext] / maxch;
4576 a = working_space[i +
l][
ymin][zmin + 2 * sizez_ext] / maxch;
4588 if(i - l + 1 < xmin)
4589 a = working_space[
xmin][
ymin][zmin + 2 * sizez_ext] / maxch;
4592 a = working_space[i - l + 1][
ymin][zmin + 2 * sizez_ext] / maxch;
4606 a = working_space[i + 1][
ymin][zmin] = a * working_space[i][
ymin][zmin];
4609 for(i = ymin;i <
ymax; i++){
4610 nip = working_space[
xmin][i][zmin + 2 * sizez_ext] / maxch;
4611 nim = working_space[
xmin][i + 1][zmin + 2 * sizez_ext] / maxch;
4613 for(l = 1;l <= averWindow; l++){
4615 a = working_space[
xmin][
ymax][zmin + 2 * sizez_ext] / maxch;
4618 a = working_space[
xmin][i +
l][zmin + 2 * sizez_ext] / maxch;
4630 if(i - l + 1 < ymin)
4631 a = working_space[
xmin][
ymin][zmin + 2 * sizez_ext] / maxch;
4634 a = working_space[
xmin][i - l + 1][zmin + 2 * sizez_ext] / maxch;
4648 a = working_space[
xmin][i + 1][zmin] = a * working_space[
xmin][i][zmin];
4651 for(i = zmin;i < zmax;i++){
4652 nip = working_space[
xmin][
ymin][i + 2 * sizez_ext] / maxch;
4653 nim = working_space[
xmin][
ymin][i + 1 + 2 * sizez_ext] / maxch;
4655 for(l = 1;l <= averWindow;l++){
4657 a = working_space[
xmin][
ymin][zmax + 2 * sizez_ext] / maxch;
4660 a = working_space[
xmin][
ymin][i + l + 2 * sizez_ext] / maxch;
4672 if(i - l + 1 < zmin)
4673 a = working_space[
xmin][
ymin][zmin + 2 * sizez_ext] / maxch;
4676 a = working_space[
xmin][
ymin][i - l + 1 + 2 * sizez_ext] / maxch;
4693 for(i = xmin;i <
xmax; i++){
4694 for(j = ymin;j <
ymax; j++){
4695 nip = working_space[i][j + 1][zmin + 2 * sizez_ext] / maxch;
4696 nim = working_space[i + 1][j + 1][zmin + 2 * sizez_ext] / maxch;
4698 for(l = 1;l <= averWindow; l++){
4700 a = working_space[
xmax][j][zmin + 2 * sizez_ext] / maxch;
4703 a = working_space[i +
l][j][zmin + 2 * sizez_ext] / maxch;
4715 if(i - l + 1 < xmin)
4716 a = working_space[
xmin][j][zmin + 2 * sizez_ext] / maxch;
4719 a = working_space[i - l + 1][j][zmin + 2 * sizez_ext] / maxch;
4733 nip = working_space[i + 1][j][zmin + 2 * sizez_ext] / maxch;
4734 nim = working_space[i + 1][j + 1][zmin + 2 * sizez_ext] / maxch;
4735 for(l = 1;l <= averWindow; l++){
4737 a = working_space[i][
ymax][zmin + 2 * sizez_ext] / maxch;
4740 a = working_space[i][j +
l][zmin + 2 * sizez_ext] / maxch;
4752 if(j - l + 1 < ymin)
4753 a = working_space[i][
ymin][zmin + 2 * sizez_ext] / maxch;
4756 a = working_space[i][j - l + 1][zmin + 2 * sizez_ext] / maxch;
4769 a = (spx * working_space[i][j + 1][zmin] + spy * working_space[i + 1][j][zmin]) / (smx + smy);
4770 working_space[i + 1][j + 1][zmin] =
a;
4774 for(i = xmin;i <
xmax;i++){
4775 for(j = zmin;j < zmax;j++){
4776 nip = working_space[i][
ymin][j + 1 + 2 * sizez_ext] / maxch;
4777 nim = working_space[i + 1][
ymin][j + 1 + 2 * sizez_ext] / maxch;
4779 for(l = 1;l <= averWindow; l++){
4781 a = working_space[
xmax][
ymin][j + 2 * sizez_ext] / maxch;
4784 a = working_space[i +
l][
ymin][j + 2 * sizez_ext] / maxch;
4796 if(i - l + 1 < xmin)
4797 a = working_space[
xmin][
ymin][j + 2 * sizez_ext] / maxch;
4800 a = working_space[i - l + 1][
ymin][j + 2 * sizez_ext] / maxch;
4814 nip = working_space[i + 1][
ymin][j + 2 * sizez_ext] / maxch;
4815 nim = working_space[i + 1][
ymin][j + 1 + 2 * sizez_ext] / maxch;
4816 for(l = 1;l <= averWindow; l++){
4818 a = working_space[i][
ymin][zmax + 2 * sizez_ext] / maxch;
4821 a = working_space[i][
ymin][j + l + 2 * sizez_ext] / maxch;
4833 if(j - l + 1 < zmin)
4834 a = working_space[i][
ymin][zmin + 2 * sizez_ext] / maxch;
4837 a = working_space[i][
ymin][j - l + 1 + 2 * sizez_ext] / maxch;
4850 a = (spx * working_space[i][
ymin][j + 1] + spz * working_space[i + 1][
ymin][j]) / (smx + smz);
4851 working_space[i + 1][
ymin][j + 1] =
a;
4855 for(i = ymin;i <
ymax;i++){
4856 for(j = zmin;j < zmax;j++){
4857 nip = working_space[
xmin][i][j + 1 + 2 * sizez_ext] / maxch;
4858 nim = working_space[
xmin][i + 1][j + 1 + 2 * sizez_ext] / maxch;
4860 for(l = 1;l <= averWindow; l++){
4862 a = working_space[
xmin][
ymax][j + 2 * sizez_ext] / maxch;
4865 a = working_space[
xmin][i +
l][j + 2 * sizez_ext] / maxch;
4877 if(i - l + 1 < ymin)
4878 a = working_space[
xmin][
ymin][j + 2 * sizez_ext] / maxch;
4881 a = working_space[
xmin][i - l + 1][j + 2 * sizez_ext] / maxch;
4895 nip = working_space[
xmin][i + 1][j + 2 * sizez_ext] / maxch;
4896 nim = working_space[
xmin][i + 1][j + 1 + 2 * sizez_ext] / maxch;
4897 for(l = 1;l <= averWindow; l++){
4899 a = working_space[
xmin][i][zmax + 2 * sizez_ext] / maxch;
4902 a = working_space[
xmin][i][j + l + 2 * sizez_ext] / maxch;
4914 if(j - l + 1 < zmin)
4915 a = working_space[
xmin][i][zmin + 2 * sizez_ext] / maxch;
4918 a = working_space[
xmin][i][j - l + 1 + 2 * sizez_ext] / maxch;
4931 a = (spy * working_space[
xmin][i][j + 1] + spz * working_space[
xmin][i + 1][j]) / (smy + smz);
4932 working_space[
xmin][i + 1][j + 1] =
a;
4936 for(i = xmin;i <
xmax; i++){
4937 for(j = ymin;j <
ymax; j++){
4938 for(k = zmin;k < zmax; k++){
4939 nip = working_space[i][j + 1][k + 1 + 2 * sizez_ext] / maxch;
4940 nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
4942 for(l = 1;l <= averWindow; l++){
4944 a = working_space[
xmax][j][k + 2 * sizez_ext] / maxch;
4947 a = working_space[i +
l][j][k + 2 * sizez_ext] / maxch;
4959 if(i - l + 1 < xmin)
4960 a = working_space[
xmin][j][k + 2 * sizez_ext] / maxch;
4963 a = working_space[i - l + 1][j][k + 2 * sizez_ext] / maxch;
4977 nip = working_space[i + 1][j][k + 1 + 2 * sizez_ext] / maxch;
4978 nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
4979 for(l = 1;l <= averWindow; l++){
4981 a = working_space[i][
ymax][k + 2 * sizez_ext] / maxch;
4984 a = working_space[i][j +
l][k + 2 * sizez_ext] / maxch;
4996 if(j - l + 1 < ymin)
4997 a = working_space[i][
ymin][k + 2 * sizez_ext] / maxch;
5000 a = working_space[i][j - l + 1][k + 2 * sizez_ext] / maxch;
5014 nip = working_space[i + 1][j + 1][k + 2 * sizez_ext] / maxch;
5015 nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
5016 for(l = 1;l <= averWindow; l++ ){
5018 a = working_space[i][j][zmax + 2 * sizez_ext] / maxch;
5021 a = working_space[i][j][k + l + 2 * sizez_ext] / maxch;
5033 if(j - l + 1 < ymin)
5034 a = working_space[i][j][zmin + 2 * sizez_ext] / maxch;
5037 a = working_space[i][j][k - l + 1 + 2 * sizez_ext] / maxch;
5050 a = (spx * working_space[i][j + 1][k + 1] + spy * working_space[i + 1][j][k + 1] + spz * working_space[i + 1][j + 1][k]) / (smx + smy + smz);
5051 working_space[i + 1][j + 1][k + 1] =
a;
5057 for(i = xmin;i <=
xmax; i++){
5058 for(j = ymin;j <=
ymax; j++){
5059 for(k = zmin;k <= zmax; k++){
5060 working_space[i][j][k] = working_space[i][j][k] / nom;
5061 a+=working_space[i][j][k];
5065 for(i = 0;i < sizex_ext; i++){
5066 for(j = 0;j < sizey_ext; j++){
5067 for(k = 0;k < sizez_ext; k++){
5068 working_space[i][j][k + 2 * sizez_ext] = working_space[i][j][k] * plocha_markov /
a;
5075 for(k = 0;k < ssizez; k++){
5076 for(j = 0;j < ssizey; j++){
5077 for(i = 0;i < ssizex; i++){
5078 working_space[i][j][k] = 0;
5079 working_space[i][j][sizez_ext + k] = 0;
5080 if(working_space[i][j][k + 3 * sizez_ext] > maximum)
5081 maximum=working_space[i][j][k+3*sizez_ext];
5088 j = (
Int_t)(pocet_sigma * sigma + 0.5);
5089 for(i = -j;i <= j; i++){
5092 b = 2.0 * sigma * sigma;
5097 s = s - sigma * sigma;
5098 s = s / (sigma * sigma * sigma * sigma);
5100 c[PEAK_WINDOW / 2 + i] = s;
5107 c[i] = c[i] / norma;
5109 a = pocet_sigma * sigma + 0.5;
5112 zmax = sizez_ext - i - 1;
5114 ymax = sizey_ext - i - 1;
5116 xmax = sizex_ext - i - 1;
5117 lmin = PEAK_WINDOW / 2 - i;
5118 lmax = PEAK_WINDOW / 2 + i;
5119 for(i = xmin;i <=
xmax; i++){
5120 for(j = ymin;j <=
ymax; j++){
5121 for(k = zmin;k <= zmax; k++){
5123 for(li = lmin;li <= lmax; li++){
5124 for(lj = lmin;lj <= lmax; lj++){
5125 for(lk = lmin;lk <= lmax; lk++){
5126 a = working_space[i + li - PEAK_WINDOW / 2][j + lj - PEAK_WINDOW / 2][k + lk - PEAK_WINDOW / 2 + 2 * sizez_ext];
5127 b = c[li] * c[lj] * c[lk];
5133 working_space[i][j][k] = s;
5134 working_space[i][j][k + sizez_ext] =
TMath::Sqrt(f);
5138 for(x = xmin;x <=
xmax; x++){
5139 for(y = ymin + 1;y <
ymax; y++){
5140 for(z = zmin + 1;z < zmax; z++){
5141 val = working_space[
x][
y][z];
5142 val1 = working_space[x - 1][y - 1][z - 1];
5143 val2 = working_space[
x ][y - 1][z - 1];
5144 val3 = working_space[x + 1][y - 1][z - 1];
5145 val4 = working_space[x - 1][
y ][z - 1];
5146 val5 = working_space[
x ][
y ][z - 1];
5147 val6 = working_space[x + 1][
y ][z - 1];
5148 val7 = working_space[x - 1][y + 1][z - 1];
5149 val8 = working_space[
x ][y + 1][z - 1];
5150 val9 = working_space[x + 1][y + 1][z - 1];
5151 val10 = working_space[x - 1][y - 1][z ];
5152 val11 = working_space[
x ][y - 1][z ];
5153 val12 = working_space[x + 1][y - 1][z ];
5154 val13 = working_space[x - 1][
y ][z ];
5155 val14 = working_space[x + 1][
y ][z ];
5156 val15 = working_space[x - 1][y + 1][z ];
5157 val16 = working_space[
x ][y + 1][z ];
5158 val17 = working_space[x + 1][y + 1][z ];
5159 val18 = working_space[x - 1][y - 1][z + 1];
5160 val19 = working_space[
x ][y - 1][z + 1];
5161 val20 = working_space[x + 1][y - 1][z + 1];
5162 val21 = working_space[x - 1][
y ][z + 1];
5163 val22 = working_space[
x ][
y ][z + 1];
5164 val23 = working_space[x + 1][
y ][z + 1];
5165 val24 = working_space[x - 1][y + 1][z + 1];
5166 val25 = working_space[
x ][y + 1][z + 1];
5167 val26 = working_space[x + 1][y + 1][z + 1];
5168 f = -s_f_ratio_peaks * working_space[
x][
y][z + sizez_ext];
5169 if(val < f && val < val1 && val < val2 && val < val3 && val < val4 && val < val5 && val < val6 && val < val7 && val < val8 && val < val9 && val < val10 && val < val11 && val < val12 && val < val13 && val < val14 && val < val15 && val < val16 && val < val17 && val < val18 && val < val19 && val < val20 && val < val21 && val < val22 && val < val23 && val < val24 && val < val25 && val < val26){
5171 for(li = lmin;li <= lmax; li++){
5172 a = working_space[x + li - PEAK_WINDOW / 2][
y][z + 2 * sizez_ext];
5174 f += a * c[li] * c[li];
5176 f = -s_f_ratio_peaks *
sqrt(f);
5179 for(li = lmin;li <= lmax; li++){
5180 a = working_space[
x][y + li - PEAK_WINDOW / 2][z + 2 * sizez_ext];
5182 f += a * c[li] * c[li];
5184 f = -s_f_ratio_peaks *
sqrt(f);
5187 for(li = lmin;li <= lmax; li++){
5188 a = working_space[
x][
y][z + li - PEAK_WINDOW / 2 + 2 * sizez_ext];
5190 f += a * c[li] * c[li];
5192 f = -s_f_ratio_peaks *
sqrt(f);
5195 for(li = lmin;li <= lmax; li++){
5196 for(lj = lmin;lj <= lmax; lj++){
5197 a = working_space[x + li - PEAK_WINDOW / 2][y + lj - PEAK_WINDOW / 2][z + 2 * sizez_ext];
5198 s += a * c[li] * c[lj];
5199 f += a * c[li] * c[li] * c[lj] * c[lj];
5202 f = s_f_ratio_peaks *
sqrt(f);
5205 for(li = lmin;li <= lmax; li++){
5206 for(lj = lmin;lj <= lmax; lj++){
5207 a = working_space[x + li - PEAK_WINDOW / 2][
y][z + lj - PEAK_WINDOW / 2 + 2 * sizez_ext];
5208 s += a * c[li] * c[lj];
5209 f += a * c[li] * c[li] * c[lj] * c[lj];
5212 f = s_f_ratio_peaks *
sqrt(f);
5215 for(li = lmin;li <= lmax; li++){
5216 for(lj=lmin;lj<=lmax;lj++){
5217 a = working_space[
x][y + li - PEAK_WINDOW / 2][z + lj - PEAK_WINDOW / 2 + 2 * sizez_ext];
5218 s += a * c[li] * c[lj];
5219 f += a * c[li] * c[li] * c[lj] * c[lj];
5222 f = s_f_ratio_peaks *
sqrt(f);
5224 if(x >= shift && x < ssizex + shift && y >= shift && y < ssizey + shift && z >= shift && z < ssizez + shift){
5225 if(working_space[x][y][z + 3 * sizez_ext] > threshold * maximum / 100.0){
5227 for(k = x - 1,a = 0,b = 0;k <= x + 1; k++){
5228 a += (
Double_t)(k - shift) * working_space[k][
y][z];
5229 b += working_space[k][
y][z];
5232 for(k = y - 1,a = 0,b = 0;k <= y + 1; k++){
5233 a += (
Double_t)(k - shift) * working_space[
x][k][z];
5234 b += working_space[
x][k][z];
5237 for(k = z - 1,a = 0,b = 0;k <= z + 1; k++){
5238 a += (
Double_t)(k - shift) * working_space[
x][
y][k];
5239 b += working_space[
x][
y][k];
5256 for(i = 0;i < ssizex; i++){
5257 for(j = 0;j < ssizey; j++){
5258 for(k = 0;k < ssizez; k++){
5259 val = -working_space[i + shift][j + shift][k + shift];
5262 dest[i][j][k] = val;
5266 k = (
Int_t)(4 * sigma + 0.5);
5268 for(i = 0;i < ssizex + k; i++){
5269 for(j = 0;j < ssizey + k; j++)
5270 delete[] working_space[i][j];
5271 delete[] working_space[i];
5273 delete[] working_space;
virtual const char * GetName() const
Returns name of object.
void SetResolution(Double_t resolution=1)
resolution: determines resolution of the neighboring peaks default value is 1 correspond to 3 sigma d...
static double p3(double t, double a, double b, double c, double d)
Int_t SearchFast(const Double_t ***source, Double_t ***dest, Int_t ssizex, Int_t ssizey, Int_t ssizez, Double_t sigma, Double_t threshold, Bool_t markov, Int_t averWindow)
const char * SmoothMarkov(Double_t ***source, Int_t ssizex, Int_t ssizey, Int_t ssizez, Int_t averWindow)
THREE-DIMENSIONAL MARKOV SPECTRUM SMOOTHING FUNCTION // // This function calculates smoothed spectrum...
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
const char * Deconvolution(Double_t ***source, const Double_t ***resp, Int_t ssizex, Int_t ssizey, Int_t ssizez, Int_t numberIterations, Int_t numberRepetitions, Double_t boost)
THREE-DIMENSIONAL DECONVOLUTION FUNCTION // This function calculates deconvolution from source spectr...
Short_t Min(Short_t a, Short_t b)
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
virtual Int_t GetDimension() const
The TNamed class is the base class for all named ROOT classes.
unsigned int r3[N_CITIES]
static double p2(double t, double a, double b, double c)
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
virtual const char * Background(const TH1 *hist, Int_t niter, Option_t *option="goff")
ONE-DIMENSIONAL BACKGROUND ESTIMATION FUNCTION // This function calculates background spectrum from s...
unsigned int r1[N_CITIES]
LVector boost(const LVector &v, const BoostVector &b)
Boost a generic Lorentz Vector class using a generic 3D Vector class describing the boost The only re...
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Advanced 3-dimentional spectra processing functions.
static double p1(double t, double a, double b)
Int_t SearchHighRes(const Double_t ***source, Double_t ***dest, Int_t ssizex, Int_t ssizey, Int_t ssizez, Double_t sigma, Double_t threshold, Bool_t backgroundRemove, Int_t deconIterations, Bool_t markov, Int_t averWindow)
#define dest(otri, vertexptr)
Short_t Max(Short_t a, Short_t b)
Double_t Sqrt(Double_t x)
virtual Int_t Search(const TH1 *hist, Double_t sigma=2, Option_t *option="goff", Double_t threshold=0.05)
ONE-DIMENSIONAL PEAK SEARCH FUNCTION // This function searches for peaks in source spectrum in hin //...
virtual void Print(Option_t *option="") const
Print the array of positions.
unsigned int r2[N_CITIES]
virtual ~TSpectrum3()
Destructor.
ROOT 6.06/08 - Reference Guide Generated on Thu Sep 29 2016 21:20:57 using Doxygen 1.8.12.