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);
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;
280 Int_t numberIterationsX,
281 Int_t numberIterationsY,
282 Int_t numberIterationsZ,
679 Int_t i, j,
x,
y, z, sampling, q1, q2, q3;
680 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;
681 if (ssizex <= 0 || ssizey <= 0 || ssizez <= 0)
682 return "Wrong parameters";
683 if (numberIterationsX < 1 || numberIterationsY < 1 || numberIterationsZ < 1)
684 return "Width of Clipping Window Must Be Positive";
685 if (ssizex < 2 * numberIterationsX + 1 || ssizey < 2 * numberIterationsY + 1 || ssizey < 2 * numberIterationsZ + 1)
686 return (
"Too Large Clipping Window");
688 for(i=0;i<ssizex;i++){
689 working_space[i] =
new Double_t*[ssizey];
690 for(j=0;j<ssizey;j++)
691 working_space[i][j]=
new Double_t[ssizez];
697 for (i = 1; i <= sampling; i++) {
699 for (z = q3; z < ssizez - q3; z++) {
700 for (y = q2; y < ssizey - q2; y++) {
701 for (x = q1; x < ssizex - q1; x++) {
702 a = spectrum[
x][
y][z];
703 p1 = spectrum[x + q1][y + q2][z - q3];
704 p2 = spectrum[x - q1][y + q2][z - q3];
705 p3 = spectrum[x + q1][y - q2][z - q3];
706 p4 = spectrum[x - q1][y - q2][z - q3];
707 p5 = spectrum[x + q1][y + q2][z + q3];
708 p6 = spectrum[x - q1][y + q2][z + q3];
709 p7 = spectrum[x + q1][y - q2][z + q3];
710 p8 = spectrum[x - q1][y - q2][z + q3];
711 s1 = spectrum[x + q1][
y ][z - q3];
712 s2 = spectrum[
x ][y + q2][z - q3];
713 s3 = spectrum[x - q1][
y ][z - q3];
714 s4 = spectrum[
x ][y - q2][z - q3];
715 s5 = spectrum[x + q1][
y ][z + q3];
716 s6 = spectrum[
x ][y + q2][z + q3];
717 s7 = spectrum[x - q1][
y ][z + q3];
718 s8 = spectrum[
x ][y - q2][z + q3];
719 s9 = spectrum[x - q1][y + q2][z ];
720 s10 = spectrum[x - q1][y - q2][z ];
721 s11 = spectrum[x + q1][y + q2][z ];
722 s12 = spectrum[x + q1][y - q2][z ];
723 r1 = spectrum[
x ][
y ][z - q3];
724 r2 = spectrum[
x ][
y ][z + q3];
725 r3 = spectrum[x - q1][
y ][z ];
726 r4 = spectrum[x + q1][
y ][z ];
727 r5 = spectrum[
x ][y + q2][z ];
728 r6 = spectrum[
x ][y - q2][z ];
765 s1 = s1 - (p1 +
p3) / 2.0;
766 s2 = s2 - (p1 +
p2) / 2.0;
767 s3 = s3 - (p2 + p4) / 2.0;
768 s4 = s4 - (p3 + p4) / 2.0;
769 s5 = s5 - (p5 + p7) / 2.0;
770 s6 = s6 - (p5 + p6) / 2.0;
771 s7 = s7 - (p6 + p8) / 2.0;
772 s8 = s8 - (p7 + p8) / 2.0;
773 s9 = s9 - (p2 + p6) / 2.0;
774 s10 = s10 - (p4 + p8) / 2.0;
775 s11 = s11 - (p1 + p5) / 2.0;
776 s12 = s12 - (p3 + p7) / 2.0;
777 b = (s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0;
780 b = (s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0;
783 b = (s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0;
786 b = (s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0;
789 b = (s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0;
792 b = (s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0;
795 r1 = r1 - ((s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0);
796 r2 = r2 - ((s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0);
797 r3 = r3 - ((s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0);
798 r4 = r4 - ((s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0);
799 r5 = r5 - ((s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0);
800 r6 = r6 - ((s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0);
801 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;
804 working_space[
x][
y][z] =
a;
808 for (z = q3; z < ssizez - q3; z++) {
809 for (y = q2; y < ssizey - q2; y++) {
810 for (x = q1; x < ssizex - q1; x++) {
811 spectrum[
x][
y][z] = working_space[
x][
y][z];
819 for (i = 1; i <= sampling; i++) {
821 for (z = q3; z < ssizez - q3; z++) {
822 for (y = q2; y < ssizey - q2; y++) {
823 for (x = q1; x < ssizex - q1; x++) {
824 a = spectrum[
x][
y][z];
825 p1 = spectrum[x + q1][y + q2][z - q3];
826 p2 = spectrum[x - q1][y + q2][z - q3];
827 p3 = spectrum[x + q1][y - q2][z - q3];
828 p4 = spectrum[x - q1][y - q2][z - q3];
829 p5 = spectrum[x + q1][y + q2][z + q3];
830 p6 = spectrum[x - q1][y + q2][z + q3];
831 p7 = spectrum[x + q1][y - q2][z + q3];
832 p8 = spectrum[x - q1][y - q2][z + q3];
833 s1 = spectrum[x + q1][
y ][z - q3];
834 s2 = spectrum[
x ][y + q2][z - q3];
835 s3 = spectrum[x - q1][
y ][z - q3];
836 s4 = spectrum[
x ][y - q2][z - q3];
837 s5 = spectrum[x + q1][
y ][z + q3];
838 s6 = spectrum[
x ][y + q2][z + q3];
839 s7 = spectrum[x - q1][
y ][z + q3];
840 s8 = spectrum[
x ][y - q2][z + q3];
841 s9 = spectrum[x - q1][y + q2][z ];
842 s10 = spectrum[x - q1][y - q2][z ];
843 s11 = spectrum[x + q1][y + q2][z ];
844 s12 = spectrum[x + q1][y - q2][z ];
845 r1 = spectrum[
x ][
y ][z - q3];
846 r2 = spectrum[
x ][
y ][z + q3];
847 r3 = spectrum[x - q1][
y ][z ];
848 r4 = spectrum[x + q1][
y ][z ];
849 r5 = spectrum[
x ][y + q2][z ];
850 r6 = spectrum[
x ][y - q2][z ];
851 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;
852 c = -(s1 + s2 + s3 + s4 + s5 + s6 + s7 + s8 + s9 + s10 + s11 + s12) / 4 + (r1 + r2 + r3 + r4 + r5 + r6) / 2;
853 d = -(p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8) / 8 + (s1 + s2 + s3 + s4 + s5 + s6 + s7 + s8 + s9 + s10 + s11 + s12) / 12;
854 if(b < a && b >= 0 && c >=0 && d >= 0)
856 working_space[
x][
y][z] =
a;
860 for (z = q3; z < ssizez - q3; z++) {
861 for (y = q2; y < ssizey - q2; y++) {
862 for (x = q1; x < ssizex - q1; x++) {
863 spectrum[
x][
y][z] = working_space[
x][
y][z];
873 for (i = sampling; i >= 1; i--) {
875 for (z = q3; z < ssizez - q3; z++) {
876 for (y = q2; y < ssizey - q2; y++) {
877 for (x = q1; x < ssizex - q1; x++) {
878 a = spectrum[
x][
y][z];
879 p1 = spectrum[x + q1][y + q2][z - q3];
880 p2 = spectrum[x - q1][y + q2][z - q3];
881 p3 = spectrum[x + q1][y - q2][z - q3];
882 p4 = spectrum[x - q1][y - q2][z - q3];
883 p5 = spectrum[x + q1][y + q2][z + q3];
884 p6 = spectrum[x - q1][y + q2][z + q3];
885 p7 = spectrum[x + q1][y - q2][z + q3];
886 p8 = spectrum[x - q1][y - q2][z + q3];
887 s1 = spectrum[x + q1][
y ][z - q3];
888 s2 = spectrum[
x ][y + q2][z - q3];
889 s3 = spectrum[x - q1][
y ][z - q3];
890 s4 = spectrum[
x ][y - q2][z - q3];
891 s5 = spectrum[x + q1][
y ][z + q3];
892 s6 = spectrum[
x ][y + q2][z + q3];
893 s7 = spectrum[x - q1][
y ][z + q3];
894 s8 = spectrum[
x ][y - q2][z + q3];
895 s9 = spectrum[x - q1][y + q2][z ];
896 s10 = spectrum[x - q1][y - q2][z ];
897 s11 = spectrum[x + q1][y + q2][z ];
898 s12 = spectrum[x + q1][y - q2][z ];
899 r1 = spectrum[
x ][
y ][z - q3];
900 r2 = spectrum[
x ][
y ][z + q3];
901 r3 = spectrum[x - q1][
y ][z ];
902 r4 = spectrum[x + q1][
y ][z ];
903 r5 = spectrum[
x ][y + q2][z ];
904 r6 = spectrum[
x ][y - q2][z ];
941 s1 = s1 - (p1 +
p3) / 2.0;
942 s2 = s2 - (p1 +
p2) / 2.0;
943 s3 = s3 - (p2 + p4) / 2.0;
944 s4 = s4 - (p3 + p4) / 2.0;
945 s5 = s5 - (p5 + p7) / 2.0;
946 s6 = s6 - (p5 + p6) / 2.0;
947 s7 = s7 - (p6 + p8) / 2.0;
948 s8 = s8 - (p7 + p8) / 2.0;
949 s9 = s9 - (p2 + p6) / 2.0;
950 s10 = s10 - (p4 + p8) / 2.0;
951 s11 = s11 - (p1 + p5) / 2.0;
952 s12 = s12 - (p3 + p7) / 2.0;
953 b = (s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0;
956 b = (s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0;
959 b = (s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0;
962 b = (s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0;
965 b = (s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0;
968 b = (s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0;
971 r1 = r1 - ((s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0);
972 r2 = r2 - ((s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0);
973 r3 = r3 - ((s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0);
974 r4 = r4 - ((s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0);
975 r5 = r5 - ((s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0);
976 r6 = r6 - ((s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0);
977 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;
980 working_space[
x][
y][z] =
a;
984 for (z = q3; z < ssizez - q3; z++) {
985 for (y = q2; y < ssizey - q2; y++) {
986 for (x = q1; x < ssizex - q1; x++) {
987 spectrum[
x][
y][z] = working_space[
x][
y][z];
995 for (i = sampling; i >= 1; i--) {
997 for (z = q3; z < ssizez - q3; z++) {
998 for (y = q2; y < ssizey - q2; y++) {
999 for (x = q1; x < ssizex - q1; x++) {
1000 a = spectrum[
x][
y][z];
1001 p1 = spectrum[x + q1][y + q2][z - q3];
1002 p2 = spectrum[x - q1][y + q2][z - q3];
1003 p3 = spectrum[x + q1][y - q2][z - q3];
1004 p4 = spectrum[x - q1][y - q2][z - q3];
1005 p5 = spectrum[x + q1][y + q2][z + q3];
1006 p6 = spectrum[x - q1][y + q2][z + q3];
1007 p7 = spectrum[x + q1][y - q2][z + q3];
1008 p8 = spectrum[x - q1][y - q2][z + q3];
1009 s1 = spectrum[x + q1][
y ][z - q3];
1010 s2 = spectrum[
x ][y + q2][z - q3];
1011 s3 = spectrum[x - q1][
y ][z - q3];
1012 s4 = spectrum[
x ][y - q2][z - q3];
1013 s5 = spectrum[x + q1][
y ][z + q3];
1014 s6 = spectrum[
x ][y + q2][z + q3];
1015 s7 = spectrum[x - q1][
y ][z + q3];
1016 s8 = spectrum[
x ][y - q2][z + q3];
1017 s9 = spectrum[x - q1][y + q2][z ];
1018 s10 = spectrum[x - q1][y - q2][z ];
1019 s11 = spectrum[x + q1][y + q2][z ];
1020 s12 = spectrum[x + q1][y - q2][z ];
1021 r1 = spectrum[
x ][
y ][z - q3];
1022 r2 = spectrum[
x ][
y ][z + q3];
1023 r3 = spectrum[x - q1][
y ][z ];
1024 r4 = spectrum[x + q1][
y ][z ];
1025 r5 = spectrum[
x ][y + q2][z ];
1026 r6 = spectrum[
x ][y - q2][z ];
1027 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;
1028 c = -(s1 + s2 + s3 + s4 + s5 + s6 + s7 + s8 + s9 + s10 + s11 + s12) / 4+(r1 + r2 + r3 + r4 + r5 + r6) / 2;
1029 d = -(p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8)/8 + (s1 + s2 + s3 + s4 + s5 + s6 + s7 + s8 + s9 + s10 + s11 + s12) / 12;
1030 if(b < a && b >= 0 && c >=0 && d >= 0)
1032 working_space[
x][
y][z] =
a;
1036 for (z = q3; z < ssizez - q3; z++) {
1037 for (y = q2; y < ssizey - q2; y++) {
1038 for (x = q1; x < ssizex - q1; x++) {
1039 spectrum[
x][
y][z] = working_space[
x][
y][z];
1046 for(i = 0;i < ssizex; i++){
1047 for(j = 0;j < ssizey; j++)
1048 delete[] working_space[i][j];
1049 delete[] working_space[i];
1051 delete[] working_space;
1310 Double_t nom,nip,nim,sp,sm,spx,smx,spy,smy,spz,smz,plocha=0;
1312 return "Averaging Window must be positive";
1314 for(i = 0;i < ssizex; i++){
1315 working_space[i] =
new Double_t*[ssizey];
1316 for(j = 0;j < ssizey; j++)
1317 working_space[i][j] =
new Double_t[ssizez];
1325 for(i = 0,maxch = 0;i < ssizex; i++){
1326 for(j = 0;j < ssizey; j++){
1327 for(k = 0;k < ssizez; k++){
1328 working_space[i][j][k] = 0;
1329 if(maxch < source[i][j][k])
1330 maxch = source[i][j][k];
1331 plocha += source[i][j][k];
1336 delete [] working_space;
1341 working_space[
xmin][
ymin][zmin] = 1;
1342 for(i = xmin;i <
xmax; i++){
1343 nip = source[i][
ymin][zmin] / maxch;
1344 nim = source[i + 1][
ymin][zmin] / maxch;
1346 for(l = 1;l <= averWindow; l++){
1348 a = source[
xmax][
ymin][zmin] / maxch;
1351 a = source[i +
l][
ymin][zmin] / maxch;
1363 if(i - l + 1 < xmin)
1364 a = source[
xmin][
ymin][zmin] / maxch;
1367 a = source[i - l + 1][
ymin][zmin] / maxch;
1381 a = working_space[i + 1][
ymin][zmin] = a * working_space[i][
ymin][zmin];
1384 for(i = ymin;i <
ymax; i++){
1385 nip = source[
xmin][i][zmin] / maxch;
1386 nim = source[
xmin][i + 1][zmin] / maxch;
1388 for(l = 1;l <= averWindow; l++){
1390 a = source[
xmin][
ymax][zmin] / maxch;
1393 a = source[
xmin][i +
l][zmin] / maxch;
1405 if(i - l + 1 < ymin)
1406 a = source[
xmin][
ymin][zmin] / maxch;
1409 a = source[
xmin][i - l + 1][zmin] / maxch;
1423 a = working_space[
xmin][i + 1][zmin] = a * working_space[
xmin][i][zmin];
1426 for(i = zmin;i < zmax; i++){
1427 nip = source[
xmin][
ymin][i] / maxch;
1428 nim = source[
xmin][
ymin][i + 1] / maxch;
1430 for(l = 1;l <= averWindow; l++){
1432 a = source[
xmin][
ymin][zmax] / maxch;
1447 if(i - l + 1 < zmin)
1448 a = source[
xmin][
ymin][zmin] / maxch;
1451 a = source[
xmin][
ymin][i - l + 1] / maxch;
1467 for(i = xmin;i <
xmax; i++){
1468 for(j = ymin;j <
ymax; j++){
1469 nip = source[i][j + 1][zmin] / maxch;
1470 nim = source[i + 1][j + 1][zmin] / maxch;
1472 for(l = 1;l <= averWindow; l++){
1474 a = source[
xmax][j][zmin] / maxch;
1477 a = source[i +
l][j][zmin] / maxch;
1489 if(i - l + 1 < xmin)
1490 a = source[
xmin][j][zmin] / maxch;
1493 a = source[i - l + 1][j][zmin] / maxch;
1507 nip = source[i + 1][j][zmin] / maxch;
1508 nim = source[i + 1][j + 1][zmin] / maxch;
1509 for(l = 1;l <= averWindow; l++){
1511 a = source[i][
ymax][zmin] / maxch;
1514 a = source[i][j +
l][zmin] / maxch;
1526 if(j - l + 1 < ymin)
1527 a = source[i][
ymin][zmin] / maxch;
1530 a = source[i][j - l + 1][zmin] / maxch;
1543 a = (spx * working_space[i][j + 1][zmin] + spy * working_space[i + 1][j][zmin]) / (smx + smy);
1544 working_space[i + 1][j + 1][zmin] =
a;
1548 for(i = xmin;i <
xmax; i++){
1549 for(j = zmin;j < zmax; j++){
1550 nip = source[i][
ymin][j + 1] / maxch;
1551 nim = source[i + 1][
ymin][j + 1] / maxch;
1553 for(l = 1;l <= averWindow; l++){
1558 a = source[i +
l][
ymin][j] / maxch;
1570 if(i - l + 1 < xmin)
1574 a = source[i - l + 1][
ymin][j] / maxch;
1588 nip = source[i + 1][
ymin][j] / maxch;
1589 nim = source[i + 1][
ymin][j + 1] / maxch;
1590 for(l = 1;l <= averWindow; l++){
1592 a = source[i][
ymin][zmax] / maxch;
1595 a = source[i][
ymin][j +
l] / maxch;
1607 if(j - l + 1 < zmin)
1608 a = source[i][
ymin][zmin] / maxch;
1611 a = source[i][
ymin][j - l + 1] / maxch;
1624 a = (spx * working_space[i][
ymin][j + 1] + spz * working_space[i + 1][
ymin][j]) / (smx + smz);
1625 working_space[i + 1][
ymin][j + 1] =
a;
1629 for(i = ymin;i <
ymax; i++){
1630 for(j = zmin;j < zmax; j++){
1631 nip = source[
xmin][i][j + 1] / maxch;
1632 nim = source[
xmin][i + 1][j + 1] / maxch;
1634 for(l = 1;l <= averWindow; l++){
1639 a = source[
xmin][i +
l][j] / maxch;
1651 if(i - l + 1 < ymin)
1655 a = source[
xmin][i - l + 1][j] / maxch;
1669 nip = source[
xmin][i + 1][j] / maxch;
1670 nim = source[
xmin][i + 1][j + 1] / maxch;
1671 for(l = 1;l <= averWindow; l++){
1673 a = source[
xmin][i][zmax] / maxch;
1676 a = source[
xmin][i][j +
l] / maxch;
1688 if(j - l + 1 < zmin)
1689 a = source[
xmin][i][zmin] / maxch;
1692 a = source[
xmin][i][j - l + 1] / maxch;
1705 a = (spy * working_space[
xmin][i][j + 1] + spz * working_space[
xmin][i + 1][j]) / (smy + smz);
1706 working_space[
xmin][i + 1][j + 1] =
a;
1710 for(i = xmin;i <
xmax; i++){
1711 for(j = ymin;j <
ymax; j++){
1712 for(k = zmin;k < zmax; k++){
1713 nip = source[i][j + 1][k + 1] / maxch;
1714 nim = source[i + 1][j + 1][k + 1] / maxch;
1716 for(l = 1;l <= averWindow; l++){
1718 a = source[
xmax][j][k] / maxch;
1721 a = source[i +
l][j][k] / maxch;
1733 if(i - l + 1 < xmin)
1734 a = source[
xmin][j][k] / maxch;
1737 a = source[i - l + 1][j][k] / maxch;
1751 nip = source[i + 1][j][k + 1] / maxch;
1752 nim = source[i + 1][j + 1][k + 1] / maxch;
1753 for(l = 1;l <= averWindow; l++){
1755 a = source[i][
ymax][k] / maxch;
1758 a = source[i][j +
l][k] / maxch;
1770 if(j - l + 1 < ymin)
1771 a = source[i][
ymin][k] / maxch;
1774 a = source[i][j - l + 1][k] / maxch;
1788 nip = source[i + 1][j + 1][k] / maxch;
1789 nim = source[i + 1][j + 1][k + 1] / maxch;
1790 for(l = 1;l <= averWindow; l++){
1792 a = source[i][j][zmax] / maxch;
1795 a = source[i][j][k +
l] / maxch;
1807 if(j - l + 1 < ymin)
1808 a = source[i][j][zmin] / maxch;
1811 a = source[i][j][k - l + 1] / maxch;
1824 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);
1825 working_space[i + 1][j + 1][k + 1] =
a;
1830 for(i = xmin;i <=
xmax; i++){
1831 for(j = ymin;j <=
ymax; j++){
1832 for(k = zmin;k <= zmax; k++){
1833 working_space[i][j][k] = working_space[i][j][k] / nom;
1837 for(i = 0;i < ssizex; i++){
1838 for(j = 0;j < ssizey; j++){
1839 for(k = 0;k < ssizez; k++){
1840 source[i][j][k] = plocha * working_space[i][j][k];
1844 for(i = 0;i < ssizex; i++){
1845 for(j = 0;j < ssizey; j++)
1846 delete[] working_space[i][j];
1847 delete[] working_space[i];
1849 delete[] working_space;
1875 Int_t numberIterations,
1876 Int_t numberRepetitions,
2363 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;
2364 Double_t lda, ldb, ldc, area, maximum = 0;
2365 if (ssizex <= 0 || ssizey <= 0 || ssizez <= 0)
2366 return "Wrong parameters";
2367 if (numberIterations <= 0)
2368 return "Number of iterations must be positive";
2369 if (numberRepetitions <= 0)
2370 return "Number of repetitions must be positive";
2372 for(i=0;i<ssizex;i++){
2373 working_space[i]=
new Double_t* [ssizey];
2374 for(j=0;j<ssizey;j++)
2375 working_space[i][j]=
new Double_t [5*ssizez];
2378 lhx = -1, lhy = -1, lhz = -1;
2379 for (i = 0; i < ssizex; i++) {
2380 for (j = 0; j < ssizey; j++) {
2381 for (k = 0; k < ssizez; k++) {
2382 lda = resp[i][j][k];
2391 working_space[i][j][k] = lda;
2393 if (lda > maximum) {
2395 positx = i, posity = j, positz = k;
2400 if (lhx == -1 || lhy == -1 || lhz == -1) {
2401 delete [] working_space;
2402 return (
"Zero response data");
2406 for (i3 = 0; i3 < ssizez; i3++) {
2407 for (i2 = 0; i2 < ssizey; i2++) {
2408 for (i1 = 0; i1 < ssizex; i1++) {
2410 for (j3 = 0; j3 <= (lhz - 1); j3++) {
2411 for (j2 = 0; j2 <= (lhy - 1); j2++) {
2412 for (j1 = 0; j1 <= (lhx - 1); j1++) {
2413 k3 = i3 + j3, k2 = i2 + j2, k1 = i1 + j1;
2414 if (k3 >= 0 && k3 < ssizez && k2 >= 0 && k2 < ssizey && k1 >= 0 && k1 < ssizex) {
2415 lda = working_space[j1][j2][j3];
2416 ldb = source[k1][k2][k3];
2417 ldc = ldc + lda * ldb;
2422 working_space[i1][i2][i3 + ssizez] = ldc;
2428 i1min = -(lhx - 1), i1max = lhx - 1;
2429 i2min = -(lhy - 1), i2max = lhy - 1;
2430 i3min = -(lhz - 1), i3max = lhz - 1;
2431 for (i3 = i3min; i3 <= i3max; i3++) {
2432 for (i2 = i2min; i2 <= i2max; i2++) {
2433 for (i1 = i1min; i1 <= i1max; i1++) {
2438 j3max = lhz - 1 - i3;
2439 if (j3max > lhz - 1)
2441 for (j3 = j3min; j3 <= j3max; j3++) {
2445 j2max = lhy - 1 - i2;
2446 if (j2max > lhy - 1)
2448 for (j2 = j2min; j2 <= j2max; j2++) {
2452 j1max = lhx - 1 - i1;
2453 if (j1max > lhx - 1)
2455 for (j1 = j1min; j1 <= j1max; j1++) {
2456 lda = working_space[j1][j2][j3];
2457 if (i1 + j1 < ssizex && i2 + j2 < ssizey)
2458 ldb = working_space[i1 + j1][i2 + j2][i3 + j3];
2461 ldc = ldc + lda * ldb;
2465 working_space[i1 - i1min][i2 - i2min][i3 - i3min + 2 * ssizez ] = ldc;
2471 for (i3 = 0; i3 < ssizez; i3++) {
2472 for (i2 = 0; i2 < ssizey; i2++) {
2473 for (i1 = 0; i1 < ssizex; i1++) {
2474 working_space[i1][i2][i3 + 3 * ssizez] = 1;
2475 working_space[i1][i2][i3 + 4 * ssizez] = 0;
2481 for (repet = 0; repet < numberRepetitions; repet++) {
2483 for (i = 0; i < ssizex; i++) {
2484 for (j = 0; j < ssizey; j++) {
2485 for (k = 0; k < ssizez; k++) {
2486 working_space[i][j][k + 3 * ssizez] =
TMath::Power(working_space[i][j][k + 3 * ssizez],boost);
2491 for (lindex = 0; lindex < numberIterations; lindex++) {
2492 for (i3 = 0; i3 < ssizez; i3++) {
2493 for (i2 = 0; i2 < ssizey; i2++) {
2494 for (i1 = 0; i1 < ssizex; i1++) {
2497 if (j3min > lhz - 1)
2500 j3max = ssizez - i3 - 1;
2501 if (j3max > lhz - 1)
2504 if (j2min > lhy - 1)
2507 j2max = ssizey - i2 - 1;
2508 if (j2max > lhy - 1)
2511 if (j1min > lhx - 1)
2514 j1max = ssizex - i1 - 1;
2515 if (j1max > lhx - 1)
2517 for (j3 = j3min; j3 <= j3max; j3++) {
2518 for (j2 = j2min; j2 <= j2max; j2++) {
2519 for (j1 = j1min; j1 <= j1max; j1++) {
2520 ldc = working_space[j1 - i1min][j2 - i2min][j3 - i3min + 2 * ssizez];
2521 lda = working_space[i1 + j1][i2 + j2][i3 + j3 + 3 * ssizez];
2522 ldb = ldb + lda * ldc;
2526 lda = working_space[i1][i2][i3 + 3 * ssizez];
2527 ldc = working_space[i1][i2][i3 + 1 * ssizez];
2528 if (ldc * lda != 0 && ldb != 0) {
2529 lda = lda * ldc / ldb;
2534 working_space[i1][i2][i3 + 4 * ssizez] = lda;
2538 for (i3 = 0; i3 < ssizez; i3++) {
2539 for (i2 = 0; i2 < ssizey; i2++) {
2540 for (i1 = 0; i1 < ssizex; i1++)
2541 working_space[i1][i2][i3 + 3 * ssizez] = working_space[i1][i2][i3 + 4 * ssizez];
2546 for (i = 0; i < ssizex; i++) {
2547 for (j = 0; j < ssizey; j++){
2548 for (k = 0; k < ssizez; k++)
2549 source[(i + positx) % ssizex][(j + posity) % ssizey][(k + positz) % ssizez] = area * working_space[i][j][k + 3 * ssizez];
2552 delete [] working_space;
2907 Int_t number_of_iterations = (
Int_t)(4 * sigma + 0.5);
2910 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;
2912 Double_t a,b,maxch,plocha = 0,plocha_markov = 0;
2913 Double_t nom,nip,nim,sp,sm,spx,spy,smx,smy,spz,smz;
2914 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;
2917 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;
2919 Error(
"SearchHighRes",
"Invalid sigma, must be greater than or equal to 1");
2923 if(threshold<=0||threshold>=100){
2924 Error(
"SearchHighRes",
"Invalid threshold, must be positive and less than 100");
2928 j = (
Int_t)(pocet_sigma*sigma+0.5);
2930 Error(
"SearchHighRes",
"Too large sigma");
2934 if (markov ==
true) {
2935 if (averWindow <= 0) {
2936 Error(
"SearchHighRes",
"Averanging window must be positive");
2941 if(backgroundRemove ==
true){
2942 if(sizex_ext < 2 * number_of_iterations + 1 || sizey_ext < 2 * number_of_iterations + 1 || sizez_ext < 2 * number_of_iterations + 1){
2943 Error(
"SearchHighRes",
"Too large clipping window");
2948 i = (
Int_t)(4 * sigma + 0.5);
2951 for(j = 0;j < ssizex + i; j++){
2952 working_space[j] =
new Double_t* [ssizey + i];
2953 for(k = 0;k < ssizey + i; k++)
2954 working_space[j][k] =
new Double_t [5 * (ssizez + i)];
2956 for(k = 0;k < sizez_ext; k++){
2957 for(j = 0;j < sizey_ext; j++){
2958 for(i = 0;i < sizex_ext; i++){
2962 working_space[i][j][k + sizez_ext] = source[0][0][0];
2964 else if(k >= ssizez + shift)
2965 working_space[i][j][k + sizez_ext] = source[0][0][ssizez - 1];
2968 working_space[i][j][k + sizez_ext] = source[0][0][k - shift];
2971 else if(j >= ssizey + shift){
2973 working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][0];
2975 else if(k >= ssizez + shift)
2976 working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][ssizez - 1];
2979 working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][k - shift];
2984 working_space[i][j][k + sizez_ext] = source[0][j - shift][0];
2986 else if(k >= ssizez + shift)
2987 working_space[i][j][k + sizez_ext] = source[0][j - shift][ssizez - 1];
2990 working_space[i][j][k + sizez_ext] = source[0][j - shift][k - shift];
2994 else if(i >= ssizex + shift){
2997 working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][0];
2999 else if(k >= ssizez + shift)
3000 working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][ssizez - 1];
3003 working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][k - shift];
3006 else if(j >= ssizey + shift){
3008 working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][0];
3010 else if(k >= ssizez + shift)
3011 working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][ssizez - 1];
3014 working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][k - shift];
3019 working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][0];
3021 else if(k >= ssizez + shift)
3022 working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][ssizez - 1];
3025 working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][k - shift];
3032 working_space[i][j][k + sizez_ext] = source[i - shift][0][0];
3034 else if(k >= ssizez + shift)
3035 working_space[i][j][k + sizez_ext] = source[i - shift][0][ssizez - 1];
3038 working_space[i][j][k + sizez_ext] = source[i - shift][0][k - shift];
3041 else if(j >= ssizey + shift){
3043 working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][0];
3045 else if(k >= ssizez + shift)
3046 working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][ssizez - 1];
3049 working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][k - shift];
3054 working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][0];
3056 else if(k >= ssizez + shift)
3057 working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][ssizez - 1];
3060 working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][k - shift];
3066 if(backgroundRemove ==
true){
3067 for(i = 1;i <= number_of_iterations; i++){
3068 for(z = i;z < sizez_ext - i; z++){
3069 for(y = i;y < sizey_ext - i; y++){
3070 for(x = i;x < sizex_ext - i; x++){
3071 a = working_space[
x][
y][z + sizez_ext];
3072 p1 = working_space[x + i][y + i][z - i + sizez_ext];
3073 p2 = working_space[x - i][y + i][z - i + sizez_ext];
3074 p3 = working_space[x + i][y - i][z - i + sizez_ext];
3075 p4 = working_space[x - i][y - i][z - i + sizez_ext];
3076 p5 = working_space[x + i][y + i][z + i + sizez_ext];
3077 p6 = working_space[x - i][y + i][z + i + sizez_ext];
3078 p7 = working_space[x + i][y - i][z + i + sizez_ext];
3079 p8 = working_space[x - i][y - i][z + i + sizez_ext];
3080 s1 = working_space[x + i][
y ][z - i + sizez_ext];
3081 s2 = working_space[
x ][y + i][z - i + sizez_ext];
3082 s3 = working_space[x - i][
y ][z - i + sizez_ext];
3083 s4 = working_space[
x ][y - i][z - i + sizez_ext];
3084 s5 = working_space[x + i][
y ][z + i + sizez_ext];
3085 s6 = working_space[
x ][y + i][z + i + sizez_ext];
3086 s7 = working_space[x - i][
y ][z + i + sizez_ext];
3087 s8 = working_space[
x ][y - i][z + i + sizez_ext];
3088 s9 = working_space[x - i][y + i][z + sizez_ext];
3089 s10 = working_space[x - i][y - i][z +sizez_ext];
3090 s11 = working_space[x + i][y + i][z +sizez_ext];
3091 s12 = working_space[x + i][y - i][z +sizez_ext];
3092 r1 = working_space[
x ][
y ][z - i + sizez_ext];
3093 r2 = working_space[
x ][
y ][z + i + sizez_ext];
3094 r3 = working_space[x - i][
y ][z + sizez_ext];
3095 r4 = working_space[x + i][
y ][z + sizez_ext];
3096 r5 = working_space[
x ][y + i][z + sizez_ext];
3097 r6 = working_space[
x ][y - i][z + sizez_ext];
3098 b = (p1 +
p3) / 2.0;
3102 b = (p1 +
p2) / 2.0;
3106 b = (p2 + p4) / 2.0;
3110 b = (p3 + p4) / 2.0;
3114 b = (p5 + p7) / 2.0;
3118 b = (p5 + p6) / 2.0;
3122 b = (p6 + p8) / 2.0;
3126 b = (p7 + p8) / 2.0;
3130 b = (p2 + p6) / 2.0;
3134 b = (p4 + p8) / 2.0;
3138 b = (p1 + p5) / 2.0;
3142 b = (p3 + p7) / 2.0;
3146 s1 = s1 - (p1 +
p3) / 2.0;
3147 s2 = s2 - (p1 +
p2) / 2.0;
3148 s3 = s3 - (p2 + p4) / 2.0;
3149 s4 = s4 - (p3 + p4) / 2.0;
3150 s5 = s5 - (p5 + p7) / 2.0;
3151 s6 = s6 - (p5 + p6) / 2.0;
3152 s7 = s7 - (p6 + p8) / 2.0;
3153 s8 = s8 - (p7 + p8) / 2.0;
3154 s9 = s9 - (p2 + p6) / 2.0;
3155 s10 = s10 - (p4 + p8) / 2.0;
3156 s11 = s11 - (p1 + p5) / 2.0;
3157 s12 = s12 - (p3 + p7) / 2.0;
3158 b = (s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0;
3162 b = (s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0;
3166 b = (s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0;
3170 b = (s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0;
3174 b = (s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0;
3178 b = (s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0;
3182 r1 = r1 - ((s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0);
3183 r2 = r2 - ((s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0);
3184 r3 = r3 - ((s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0);
3185 r4 = r4 - ((s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0);
3186 r5 = r5 - ((s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0);
3187 r6 = r6 - ((s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0);
3188 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;
3192 working_space[
x][
y][z] =
a;
3196 for(z = i;z < sizez_ext - i; z++){
3197 for(y = i;y < sizey_ext - i; y++){
3198 for(x = i;x < sizex_ext - i; x++){
3199 working_space[
x][
y][z + sizez_ext] = working_space[
x][
y][z];
3204 for(k = 0;k < sizez_ext; k++){
3205 for(j = 0;j < sizey_ext; j++){
3206 for(i = 0;i < sizex_ext; i++){
3210 working_space[i][j][k + sizez_ext] = source[0][0][0] - working_space[i][j][k + sizez_ext];
3212 else if(k >= ssizez + shift)
3213 working_space[i][j][k + sizez_ext] = source[0][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
3216 working_space[i][j][k + sizez_ext] = source[0][0][k - shift] - working_space[i][j][k + sizez_ext];
3219 else if(j >= ssizey + shift){
3221 working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
3223 else if(k >= ssizez + shift)
3224 working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
3227 working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
3232 working_space[i][j][k + sizez_ext] = source[0][j - shift][0] - working_space[i][j][k + sizez_ext];
3234 else if(k >= ssizez + shift)
3235 working_space[i][j][k + sizez_ext] = source[0][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
3238 working_space[i][j][k + sizez_ext] = source[0][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
3242 else if(i >= ssizex + shift){
3245 working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][0] - working_space[i][j][k + sizez_ext];
3247 else if(k >= ssizez + shift)
3248 working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
3251 working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][k - shift] - working_space[i][j][k + sizez_ext];
3254 else if(j >= ssizey + shift){
3256 working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
3258 else if(k >= ssizez + shift)
3259 working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
3262 working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
3267 working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][0] - working_space[i][j][k + sizez_ext];
3269 else if(k >= ssizez + shift)
3270 working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
3273 working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
3280 working_space[i][j][k + sizez_ext] = source[i - shift][0][0] - working_space[i][j][k + sizez_ext];
3282 else if(k >= ssizez + shift)
3283 working_space[i][j][k + sizez_ext] = source[i - shift][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
3286 working_space[i][j][k + sizez_ext] = source[i - shift][0][k - shift] - working_space[i][j][k + sizez_ext];
3289 else if(j >= ssizey + shift){
3291 working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
3293 else if(k >= ssizez + shift)
3294 working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
3297 working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
3302 working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][0] - working_space[i][j][k + sizez_ext];
3304 else if(k >= ssizez + shift)
3305 working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
3308 working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
3317 for(i = 0;i < sizex_ext; i++){
3318 for(j = 0;j < sizey_ext; j++){
3319 for(k = 0;k < sizez_ext; k++){
3320 working_space[i][j][k + 2 * sizez_ext] = working_space[i][j][sizez_ext + k];
3321 plocha_markov = plocha_markov + working_space[i][j][sizez_ext + k];
3326 xmax = sizex_ext - 1;
3328 ymax = sizey_ext - 1;
3330 zmax = sizez_ext - 1;
3331 for(i = 0,maxch = 0;i < sizex_ext; i++){
3332 for(j = 0;j < sizey_ext;j++){
3333 for(k = 0;k < sizez_ext;k++){
3334 working_space[i][j][k] = 0;
3335 if(maxch < working_space[i][j][k + 2 * sizez_ext])
3336 maxch = working_space[i][j][k + 2 * sizez_ext];
3338 plocha += working_space[i][j][k + 2 * sizez_ext];
3343 delete [] working_space;
3347 working_space[
xmin][
ymin][zmin] = 1;
3348 for(i = xmin;i <
xmax; i++){
3349 nip = working_space[i][
ymin][zmin + 2 * sizez_ext] / maxch;
3350 nim = working_space[i + 1][
ymin][zmin + 2 * sizez_ext] / maxch;
3352 for(l = 1;l <= averWindow; l++){
3354 a = working_space[xmax][ymin][zmin + 2 * sizez_ext] / maxch;
3357 a = working_space[i +
l][
ymin][zmin + 2 * sizez_ext] / maxch;
3369 if(i - l + 1 < xmin)
3370 a = working_space[
xmin][
ymin][zmin + 2 * sizez_ext] / maxch;
3373 a = working_space[i - l + 1][
ymin][zmin + 2 * sizez_ext] / maxch;
3387 a = working_space[i + 1][
ymin][zmin] = a * working_space[i][
ymin][zmin];
3390 for(i = ymin;i <
ymax; i++){
3391 nip = working_space[
xmin][i][zmin + 2 * sizez_ext] / maxch;
3392 nim = working_space[
xmin][i + 1][zmin + 2 * sizez_ext] / maxch;
3394 for(l = 1;l <= averWindow; l++){
3396 a = working_space[
xmin][
ymax][zmin + 2 * sizez_ext] / maxch;
3399 a = working_space[
xmin][i +
l][zmin + 2 * sizez_ext] / maxch;
3411 if(i - l + 1 < ymin)
3412 a = working_space[
xmin][
ymin][zmin + 2 * sizez_ext] / maxch;
3415 a = working_space[
xmin][i - l + 1][zmin + 2 * sizez_ext] / maxch;
3429 a = working_space[
xmin][i + 1][zmin] = a * working_space[
xmin][i][zmin];
3432 for(i = zmin;i < zmax;i++){
3433 nip = working_space[
xmin][
ymin][i + 2 * sizez_ext] / maxch;
3434 nim = working_space[
xmin][
ymin][i + 1 + 2 * sizez_ext] / maxch;
3436 for(l = 1;l <= averWindow;l++){
3438 a = working_space[
xmin][
ymin][zmax + 2 * sizez_ext] / maxch;
3441 a = working_space[
xmin][
ymin][i + l + 2 * sizez_ext] / maxch;
3453 if(i - l + 1 < zmin)
3454 a = working_space[
xmin][
ymin][zmin + 2 * sizez_ext] / maxch;
3457 a = working_space[
xmin][
ymin][i - l + 1 + 2 * sizez_ext] / maxch;
3474 for(i = xmin;i <
xmax; i++){
3475 for(j = ymin;j <
ymax; j++){
3476 nip = working_space[i][j + 1][zmin + 2 * sizez_ext] / maxch;
3477 nim = working_space[i + 1][j + 1][zmin + 2 * sizez_ext] / maxch;
3479 for(l = 1;l <= averWindow; l++){
3481 a = working_space[
xmax][j][zmin + 2 * sizez_ext] / maxch;
3484 a = working_space[i +
l][j][zmin + 2 * sizez_ext] / maxch;
3496 if(i - l + 1 < xmin)
3497 a = working_space[
xmin][j][zmin + 2 * sizez_ext] / maxch;
3500 a = working_space[i - l + 1][j][zmin + 2 * sizez_ext] / maxch;
3514 nip = working_space[i + 1][j][zmin + 2 * sizez_ext] / maxch;
3515 nim = working_space[i + 1][j + 1][zmin + 2 * sizez_ext] / maxch;
3516 for(l = 1;l <= averWindow; l++){
3518 a = working_space[i][
ymax][zmin + 2 * sizez_ext] / maxch;
3521 a = working_space[i][j +
l][zmin + 2 * sizez_ext] / maxch;
3533 if(j - l + 1 < ymin)
3534 a = working_space[i][
ymin][zmin + 2 * sizez_ext] / maxch;
3537 a = working_space[i][j - l + 1][zmin + 2 * sizez_ext] / maxch;
3550 a = (spx * working_space[i][j + 1][zmin] + spy * working_space[i + 1][j][zmin]) / (smx + smy);
3551 working_space[i + 1][j + 1][zmin] =
a;
3555 for(i = xmin;i <
xmax;i++){
3556 for(j = zmin;j < zmax;j++){
3557 nip = working_space[i][
ymin][j + 1 + 2 * sizez_ext] / maxch;
3558 nim = working_space[i + 1][
ymin][j + 1 + 2 * sizez_ext] / maxch;
3560 for(l = 1;l <= averWindow; l++){
3562 a = working_space[
xmax][
ymin][j + 2 * sizez_ext] / maxch;
3565 a = working_space[i +
l][
ymin][j + 2 * sizez_ext] / maxch;
3577 if(i - l + 1 < xmin)
3578 a = working_space[
xmin][
ymin][j + 2 * sizez_ext] / maxch;
3581 a = working_space[i - l + 1][
ymin][j + 2 * sizez_ext] / maxch;
3595 nip = working_space[i + 1][
ymin][j + 2 * sizez_ext] / maxch;
3596 nim = working_space[i + 1][
ymin][j + 1 + 2 * sizez_ext] / maxch;
3597 for(l = 1;l <= averWindow; l++){
3599 a = working_space[i][
ymin][zmax + 2 * sizez_ext] / maxch;
3602 a = working_space[i][
ymin][j + l + 2 * sizez_ext] / maxch;
3614 if(j - l + 1 < zmin)
3615 a = working_space[i][
ymin][zmin + 2 * sizez_ext] / maxch;
3618 a = working_space[i][
ymin][j - l + 1 + 2 * sizez_ext] / maxch;
3631 a = (spx * working_space[i][
ymin][j + 1] + spz * working_space[i + 1][
ymin][j]) / (smx + smz);
3632 working_space[i + 1][
ymin][j + 1] =
a;
3636 for(i = ymin;i <
ymax;i++){
3637 for(j = zmin;j < zmax;j++){
3638 nip = working_space[
xmin][i][j + 1 + 2 * sizez_ext] / maxch;
3639 nim = working_space[
xmin][i + 1][j + 1 + 2 * sizez_ext] / maxch;
3641 for(l = 1;l <= averWindow; l++){
3643 a = working_space[
xmin][
ymax][j + 2 * sizez_ext] / maxch;
3646 a = working_space[
xmin][i +
l][j + 2 * sizez_ext] / maxch;
3658 if(i - l + 1 < ymin)
3659 a = working_space[
xmin][
ymin][j + 2 * sizez_ext] / maxch;
3662 a = working_space[
xmin][i - l + 1][j + 2 * sizez_ext] / maxch;
3676 nip = working_space[
xmin][i + 1][j + 2 * sizez_ext] / maxch;
3677 nim = working_space[
xmin][i + 1][j + 1 + 2 * sizez_ext] / maxch;
3678 for(l = 1;l <= averWindow; l++){
3680 a = working_space[
xmin][i][zmax + 2 * sizez_ext] / maxch;
3683 a = working_space[
xmin][i][j + l + 2 * sizez_ext] / maxch;
3695 if(j - l + 1 < zmin)
3696 a = working_space[
xmin][i][zmin + 2 * sizez_ext] / maxch;
3699 a = working_space[
xmin][i][j - l + 1 + 2 * sizez_ext] / maxch;
3712 a = (spy * working_space[
xmin][i][j + 1] + spz * working_space[
xmin][i + 1][j]) / (smy + smz);
3713 working_space[
xmin][i + 1][j + 1] =
a;
3717 for(i = xmin;i <
xmax; i++){
3718 for(j = ymin;j <
ymax; j++){
3719 for(k = zmin;k < zmax; k++){
3720 nip = working_space[i][j + 1][k + 1 + 2 * sizez_ext] / maxch;
3721 nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
3723 for(l = 1;l <= averWindow; l++){
3725 a = working_space[
xmax][j][k + 2 * sizez_ext] / maxch;
3728 a = working_space[i +
l][j][k + 2 * sizez_ext] / maxch;
3740 if(i - l + 1 < xmin)
3741 a = working_space[
xmin][j][k + 2 * sizez_ext] / maxch;
3744 a = working_space[i - l + 1][j][k + 2 * sizez_ext] / maxch;
3758 nip = working_space[i + 1][j][k + 1 + 2 * sizez_ext] / maxch;
3759 nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
3760 for(l = 1;l <= averWindow; l++){
3762 a = working_space[i][
ymax][k + 2 * sizez_ext] / maxch;
3765 a = working_space[i][j +
l][k + 2 * sizez_ext] / maxch;
3777 if(j - l + 1 < ymin)
3778 a = working_space[i][
ymin][k + 2 * sizez_ext] / maxch;
3781 a = working_space[i][j - l + 1][k + 2 * sizez_ext] / maxch;
3795 nip = working_space[i + 1][j + 1][k + 2 * sizez_ext] / maxch;
3796 nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
3797 for(l = 1;l <= averWindow; l++ ){
3799 a = working_space[i][j][zmax + 2 * sizez_ext] / maxch;
3802 a = working_space[i][j][k + l + 2 * sizez_ext] / maxch;
3814 if(j - l + 1 < ymin)
3815 a = working_space[i][j][zmin + 2 * sizez_ext] / maxch;
3818 a = working_space[i][j][k - l + 1 + 2 * sizez_ext] / maxch;
3831 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);
3832 working_space[i + 1][j + 1][k + 1] =
a;
3838 for(i = xmin;i <=
xmax; i++){
3839 for(j = ymin;j <=
ymax; j++){
3840 for(k = zmin;k <= zmax; k++){
3841 working_space[i][j][k] = working_space[i][j][k] / nom;
3842 a+=working_space[i][j][k];
3846 for(i = 0;i < sizex_ext; i++){
3847 for(j = 0;j < sizey_ext; j++){
3848 for(k = 0;k < sizez_ext; k++){
3849 working_space[i][j][k + sizez_ext] = working_space[i][j][k] * plocha_markov /
a;
3856 lhx = -1,lhy = -1,lhz = -1;
3857 positx = 0,posity = 0,positz = 0;
3860 for(i = 0;i < sizex_ext; i++){
3861 for(j = 0;j < sizey_ext; j++){
3862 for(k = 0;k < sizez_ext; k++){
3866 lda = (lda * lda + ldb * ldb + ldc * ldc) / (2 * sigma * sigma);
3879 working_space[i][j][k] = lda;
3883 positx = i,posity = j,positz = k;
3889 for(i = 0;i < sizex_ext; i++){
3890 for(j = 0;j < sizey_ext; j++){
3891 for(k = 0;k < sizez_ext; k++){
3892 working_space[i][j][k + 2 * sizez_ext] =
TMath::Abs(working_space[i][j][k + sizez_ext]);
3897 for (i3 = 0; i3 < sizez_ext; i3++) {
3898 for (i2 = 0; i2 < sizey_ext; i2++) {
3899 for (i1 = 0; i1 < sizex_ext; i1++) {
3901 for (j3 = 0; j3 <= (lhz - 1); j3++) {
3902 for (j2 = 0; j2 <= (lhy - 1); j2++) {
3903 for (j1 = 0; j1 <= (lhx - 1); j1++) {
3904 k3 = i3 + j3, k2 = i2 + j2, k1 = i1 + j1;
3905 if (k3 >= 0 && k3 < sizez_ext && k2 >= 0 && k2 < sizey_ext && k1 >= 0 && k1 < sizex_ext) {
3906 lda = working_space[j1][j2][j3];
3907 ldb = working_space[k1][k2][k3+2*sizez_ext];
3908 ldc = ldc + lda * ldb;
3913 working_space[i1][i2][i3 + sizez_ext] = ldc;
3918 i1min = -(lhx - 1), i1max = lhx - 1;
3919 i2min = -(lhy - 1), i2max = lhy - 1;
3920 i3min = -(lhz - 1), i3max = lhz - 1;
3921 for (i3 = i3min; i3 <= i3max; i3++) {
3922 for (i2 = i2min; i2 <= i2max; i2++) {
3923 for (i1 = i1min; i1 <= i1max; i1++) {
3929 j3max = lhz - 1 - i3;
3930 if (j3max > lhz - 1)
3933 for (j3 = j3min; j3 <= j3max; j3++) {
3938 j2max = lhy - 1 - i2;
3939 if (j2max > lhy - 1)
3942 for (j2 = j2min; j2 <= j2max; j2++) {
3947 j1max = lhx - 1 - i1;
3948 if (j1max > lhx - 1)
3951 for (j1 = j1min; j1 <= j1max; j1++) {
3952 lda = working_space[j1][j2][j3];
3953 if (i1 + j1 < sizex_ext && i2 + j2 < sizey_ext)
3954 ldb = working_space[i1 + j1][i2 + j2][i3 + j3];
3959 ldc = ldc + lda * ldb;
3963 working_space[i1 - i1min][i2 - i2min][i3 - i3min + 2 * sizez_ext ] = ldc;
3968 for (i3 = 0; i3 < sizez_ext; i3++) {
3969 for (i2 = 0; i2 < sizey_ext; i2++) {
3970 for (i1 = 0; i1 < sizex_ext; i1++) {
3971 working_space[i1][i2][i3 + 3 * sizez_ext] = 1;
3972 working_space[i1][i2][i3 + 4 * sizez_ext] = 0;
3978 for (lindex=0;lindex<deconIterations;lindex++){
3979 for (i3 = 0; i3 < sizez_ext; i3++) {
3980 for (i2 = 0; i2 < sizey_ext; i2++) {
3981 for (i1 = 0; i1 < sizex_ext; i1++) {
3982 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){
3985 if (j3min > lhz - 1)
3989 j3max = sizez_ext - i3 - 1;
3990 if (j3max > lhz - 1)
3994 if (j2min > lhy - 1)
3998 j2max = sizey_ext - i2 - 1;
3999 if (j2max > lhy - 1)
4003 if (j1min > lhx - 1)
4007 j1max = sizex_ext - i1 - 1;
4008 if (j1max > lhx - 1)
4011 for (j3 = j3min; j3 <= j3max; j3++) {
4012 for (j2 = j2min; j2 <= j2max; j2++) {
4013 for (j1 = j1min; j1 <= j1max; j1++) {
4014 ldc = working_space[j1 - i1min][j2 - i2min][j3 - i3min + 2 * sizez_ext];
4015 lda = working_space[i1 + j1][i2 + j2][i3 + j3 + 3 * sizez_ext];
4016 ldb = ldb + lda * ldc;
4020 lda = working_space[i1][i2][i3 + 3 * sizez_ext];
4021 ldc = working_space[i1][i2][i3 + 1 * sizez_ext];
4022 if (ldc * lda != 0 && ldb != 0) {
4023 lda = lda * ldc / ldb;
4028 working_space[i1][i2][i3 + 4 * sizez_ext] = lda;
4033 for (i3 = 0; i3 < sizez_ext; i3++) {
4034 for (i2 = 0; i2 < sizey_ext; i2++) {
4035 for (i1 = 0; i1 < sizex_ext; i1++)
4036 working_space[i1][i2][i3 + 3 * sizez_ext] = working_space[i1][i2][i3 + 4 * sizez_ext];
4042 for(i = 0;i < sizex_ext; i++){
4043 for(j = 0;j < sizey_ext; j++){
4044 for(k = 0;k < sizez_ext; k++){
4045 working_space[(i + positx) % sizex_ext][(j + posity) % sizey_ext][(k + positz) % sizez_ext] = area * working_space[i][j][k + 3 * sizez_ext];
4046 if(maximum < working_space[(i + positx) % sizex_ext][(j + posity) % sizey_ext][(k + positz) % sizez_ext])
4047 maximum = working_space[(i + positx) % sizex_ext][(j + posity) % sizey_ext][(k + positz) % sizez_ext];
4052 for(i = 1;i < sizex_ext - 1; i++){
4053 for(j = 1;j < sizey_ext - 1; j++){
4054 for(l = 1;l < sizez_ext - 1; l++){
4055 a = working_space[i][j][
l];
4056 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]){
4057 if(i >= shift && i < ssizex + shift && j >= shift && j < ssizey + shift && l >= shift && l < ssizez + shift){
4058 if(working_space[i][j][l] > threshold * maximum / 100.0){
4060 for(k = i - 1,a = 0,b = 0;k <= i + 1; k++){
4061 a += (
Double_t)(k - shift) * working_space[k][j][
l];
4062 b += working_space[k][j][
l];
4065 for(k = j - 1,a = 0,b = 0;k <= j + 1; k++){
4066 a += (
Double_t)(k - shift) * working_space[i][k][
l];
4067 b += working_space[i][k][
l];
4070 for(k = l - 1,a = 0,b = 0;k <= l + 1; k++){
4071 a += (
Double_t)(k - shift) * working_space[i][j][k];
4072 b += working_space[i][j][k];
4083 for(i = 0;i < ssizex; i++){
4084 for(j = 0;j < ssizey; j++){
4085 for(k = 0;k < ssizez; k++){
4086 dest[i][j][k] = working_space[i + shift][j + shift][k + shift];
4090 k = (
Int_t)(4 * sigma + 0.5);
4092 for(i = 0;i < ssizex + k; i++){
4093 for(j = 0;j < ssizey + k; j++)
4094 delete[] working_space[i][j];
4095 delete[] working_space[i];
4097 delete[] working_space;
4131 Int_t i,j,k,
l,li,lj,lk,lmin,lmax,
xmin,
xmax,
ymin,
ymax,zmin,zmax;
4132 Double_t maxch,plocha = 0,plocha_markov = 0;
4133 Double_t nom,nip,nim,sp,sm,spx,spy,smx,smy,spz,smz;
4134 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;
4137 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;
4139 Int_t number_of_iterations=(
Int_t)(4 * sigma + 0.5);
4140 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;
4143 Error(
"SearchFast",
"Invalid sigma, must be greater than or equal to 1");
4147 if(threshold<=0||threshold>=100){
4148 Error(
"SearchFast",
"Invalid threshold, must be positive and less than 100");
4152 j = (
Int_t)(pocet_sigma*sigma+0.5);
4154 Error(
"SearchFast",
"Too large sigma");
4158 if (markov ==
true) {
4159 if (averWindow <= 0) {
4160 Error(
"SearchFast",
"Averanging window must be positive");
4165 if(sizex_ext < 2 * number_of_iterations + 1 || sizey_ext < 2 * number_of_iterations + 1 || sizez_ext < 2 * number_of_iterations + 1){
4166 Error(
"SearchFast",
"Too large clipping window");
4170 i = (
Int_t)(4 * sigma + 0.5);
4173 for(j = 0;j < ssizex + i; j++){
4174 working_space[j] =
new Double_t* [ssizey + i];
4175 for(k = 0;k < ssizey + i; k++)
4176 working_space[j][k] =
new Double_t [4 * (ssizez + i)];
4179 for(k = 0;k < sizez_ext; k++){
4180 for(j = 0;j < sizey_ext; j++){
4181 for(i = 0;i < sizex_ext; i++){
4185 working_space[i][j][k + sizez_ext] = source[0][0][0];
4187 else if(k >= ssizez + shift)
4188 working_space[i][j][k + sizez_ext] = source[0][0][ssizez - 1];
4191 working_space[i][j][k + sizez_ext] = source[0][0][k - shift];
4194 else if(j >= ssizey + shift){
4196 working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][0];
4198 else if(k >= ssizez + shift)
4199 working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][ssizez - 1];
4202 working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][k - shift];
4207 working_space[i][j][k + sizez_ext] = source[0][j - shift][0];
4209 else if(k >= ssizez + shift)
4210 working_space[i][j][k + sizez_ext] = source[0][j - shift][ssizez - 1];
4213 working_space[i][j][k + sizez_ext] = source[0][j - shift][k - shift];
4217 else if(i >= ssizex + shift){
4220 working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][0];
4222 else if(k >= ssizez + shift)
4223 working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][ssizez - 1];
4226 working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][k - shift];
4229 else if(j >= ssizey + shift){
4231 working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][0];
4233 else if(k >= ssizez + shift)
4234 working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][ssizez - 1];
4237 working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][k - shift];
4242 working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][0];
4244 else if(k >= ssizez + shift)
4245 working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][ssizez - 1];
4248 working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][k - shift];
4255 working_space[i][j][k + sizez_ext] = source[i - shift][0][0];
4257 else if(k >= ssizez + shift)
4258 working_space[i][j][k + sizez_ext] = source[i - shift][0][ssizez - 1];
4261 working_space[i][j][k + sizez_ext] = source[i - shift][0][k - shift];
4264 else if(j >= ssizey + shift){
4266 working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][0];
4268 else if(k >= ssizez + shift)
4269 working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][ssizez - 1];
4272 working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][k - shift];
4277 working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][0];
4279 else if(k >= ssizez + shift)
4280 working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][ssizez - 1];
4283 working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][k - shift];
4289 for(i = 1;i <= number_of_iterations; i++){
4290 for(z = i;z < sizez_ext - i; z++){
4291 for(y = i;y < sizey_ext - i; y++){
4292 for(x = i;x < sizex_ext - i; x++){
4293 a = working_space[
x][
y][z + sizez_ext];
4294 p1 = working_space[x + i][y + i][z - i + sizez_ext];
4295 p2 = working_space[x - i][y + i][z - i + sizez_ext];
4296 p3 = working_space[x + i][y - i][z - i + sizez_ext];
4297 p4 = working_space[x - i][y - i][z - i + sizez_ext];
4298 p5 = working_space[x + i][y + i][z + i + sizez_ext];
4299 p6 = working_space[x - i][y + i][z + i + sizez_ext];
4300 p7 = working_space[x + i][y - i][z + i + sizez_ext];
4301 p8 = working_space[x - i][y - i][z + i + sizez_ext];
4302 s1 = working_space[x + i][
y ][z - i + sizez_ext];
4303 s2 = working_space[
x ][y + i][z - i + sizez_ext];
4304 s3 = working_space[x - i][
y ][z - i + sizez_ext];
4305 s4 = working_space[
x ][y - i][z - i + sizez_ext];
4306 s5 = working_space[x + i][
y ][z + i + sizez_ext];
4307 s6 = working_space[
x ][y + i][z + i + sizez_ext];
4308 s7 = working_space[x - i][
y ][z + i + sizez_ext];
4309 s8 = working_space[
x ][y - i][z + i + sizez_ext];
4310 s9 = working_space[x - i][y + i][z + sizez_ext];
4311 s10 = working_space[x - i][y - i][z +sizez_ext];
4312 s11 = working_space[x + i][y + i][z +sizez_ext];
4313 s12 = working_space[x + i][y - i][z +sizez_ext];
4314 r1 = working_space[
x ][
y ][z - i + sizez_ext];
4315 r2 = working_space[
x ][
y ][z + i + sizez_ext];
4316 r3 = working_space[x - i][
y ][z + sizez_ext];
4317 r4 = working_space[x + i][
y ][z + sizez_ext];
4318 r5 = working_space[
x ][y + i][z + sizez_ext];
4319 r6 = working_space[
x ][y - i][z + sizez_ext];
4320 b = (p1 +
p3) / 2.0;
4324 b = (p1 +
p2) / 2.0;
4328 b = (p2 + p4) / 2.0;
4332 b = (p3 + p4) / 2.0;
4336 b = (p5 + p7) / 2.0;
4340 b = (p5 + p6) / 2.0;
4344 b = (p6 + p8) / 2.0;
4348 b = (p7 + p8) / 2.0;
4352 b = (p2 + p6) / 2.0;
4356 b = (p4 + p8) / 2.0;
4360 b = (p1 + p5) / 2.0;
4364 b = (p3 + p7) / 2.0;
4368 s1 = s1 - (p1 +
p3) / 2.0;
4369 s2 = s2 - (p1 +
p2) / 2.0;
4370 s3 = s3 - (p2 + p4) / 2.0;
4371 s4 = s4 - (p3 + p4) / 2.0;
4372 s5 = s5 - (p5 + p7) / 2.0;
4373 s6 = s6 - (p5 + p6) / 2.0;
4374 s7 = s7 - (p6 + p8) / 2.0;
4375 s8 = s8 - (p7 + p8) / 2.0;
4376 s9 = s9 - (p2 + p6) / 2.0;
4377 s10 = s10 - (p4 + p8) / 2.0;
4378 s11 = s11 - (p1 + p5) / 2.0;
4379 s12 = s12 - (p3 + p7) / 2.0;
4380 b = (s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0;
4384 b = (s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0;
4388 b = (s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0;
4392 b = (s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0;
4396 b = (s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0;
4400 b = (s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0;
4404 r1 = r1 - ((s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0);
4405 r2 = r2 - ((s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0);
4406 r3 = r3 - ((s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0);
4407 r4 = r4 - ((s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0);
4408 r5 = r5 - ((s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0);
4409 r6 = r6 - ((s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0);
4410 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;
4414 working_space[
x][
y][z] =
a;
4418 for(z = i;z < sizez_ext - i; z++){
4419 for(y = i;y < sizey_ext - i; y++){
4420 for(x = i;x < sizex_ext - i; x++){
4421 working_space[
x][
y][z + sizez_ext] = working_space[
x][
y][z];
4426 for(k = 0;k < sizez_ext; k++){
4427 for(j = 0;j < sizey_ext; j++){
4428 for(i = 0;i < sizex_ext; i++){
4432 working_space[i][j][k + 3 * sizez_ext] = source[0][0][0] - working_space[i][j][k + sizez_ext];
4434 else if(k >= ssizez + shift)
4435 working_space[i][j][k + 3 * sizez_ext] = source[0][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
4438 working_space[i][j][k + 3 * sizez_ext] = source[0][0][k - shift] - working_space[i][j][k + sizez_ext];
4441 else if(j >= ssizey + shift){
4443 working_space[i][j][k + 3 * sizez_ext] = source[0][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
4445 else if(k >= ssizez + shift)
4446 working_space[i][j][k + 3 * sizez_ext] = source[0][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
4449 working_space[i][j][k + 3 * sizez_ext] = source[0][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
4454 working_space[i][j][k + 3 * sizez_ext] = source[0][j - shift][0] - working_space[i][j][k + sizez_ext];
4456 else if(k >= ssizez + shift)
4457 working_space[i][j][k + 3 * sizez_ext] = source[0][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
4460 working_space[i][j][k + 3 * sizez_ext] = source[0][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
4464 else if(i >= ssizex + shift){
4467 working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][0][0] - working_space[i][j][k + sizez_ext];
4469 else if(k >= ssizez + shift)
4470 working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
4473 working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][0][k - shift] - working_space[i][j][k + sizez_ext];
4476 else if(j >= ssizey + shift){
4478 working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
4480 else if(k >= ssizez + shift)
4481 working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
4484 working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
4489 working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][j - shift][0] - working_space[i][j][k + sizez_ext];
4491 else if(k >= ssizez + shift)
4492 working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
4495 working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
4502 working_space[i][j][k + 3 * sizez_ext] = source[i - shift][0][0] - working_space[i][j][k + sizez_ext];
4504 else if(k >= ssizez + shift)
4505 working_space[i][j][k + 3 * sizez_ext] = source[i - shift][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
4508 working_space[i][j][k + 3 * sizez_ext] = source[i - shift][0][k - shift] - working_space[i][j][k + sizez_ext];
4511 else if(j >= ssizey + shift){
4513 working_space[i][j][k + 3 * sizez_ext] = source[i - shift][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
4515 else if(k >= ssizez + shift)
4516 working_space[i][j][k + 3 * sizez_ext] = source[i - shift][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
4519 working_space[i][j][k + 3 * sizez_ext] = source[i - shift][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
4524 working_space[i][j][k + 3 * sizez_ext] = source[i - shift][j - shift][0] - working_space[i][j][k + sizez_ext];
4526 else if(k >= ssizez + shift)
4527 working_space[i][j][k + 3 * sizez_ext] = source[i - shift][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
4530 working_space[i][j][k + 3 * sizez_ext] = source[i - shift][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
4537 for(i = 0;i < sizex_ext; i++){
4538 for(j = 0;j < sizey_ext; j++){
4539 for(k = 0;k < sizez_ext; k++){
4540 if(i >= shift && i < ssizex + shift && j >= shift && j < ssizey + shift && k >= shift && k < ssizez + shift){
4541 working_space[i][j][k + 2 * sizez_ext] = source[i - shift][j - shift][k - shift];
4542 plocha_markov = plocha_markov + source[i - shift][j - shift][k - shift];
4545 working_space[i][j][k + 2 * sizez_ext] = 0;
4552 xmax = sizex_ext - 1;
4554 ymax = sizey_ext - 1;
4556 zmax = sizez_ext - 1;
4557 for(i = 0,maxch = 0;i < sizex_ext; i++){
4558 for(j = 0;j < sizey_ext;j++){
4559 for(k = 0;k < sizez_ext;k++){
4560 working_space[i][j][k] = 0;
4561 if(maxch < working_space[i][j][k + 2 * sizez_ext])
4562 maxch = working_space[i][j][k + 2 * sizez_ext];
4564 plocha += working_space[i][j][k + 2 * sizez_ext];
4569 delete [] working_space;
4574 working_space[
xmin][
ymin][zmin] = 1;
4575 for(i = xmin;i <
xmax; i++){
4576 nip = working_space[i][
ymin][zmin + 2 * sizez_ext] / maxch;
4577 nim = working_space[i + 1][
ymin][zmin + 2 * sizez_ext] / maxch;
4579 for(l = 1;l <= averWindow; l++){
4581 a = working_space[xmax][ymin][zmin + 2 * sizez_ext] / maxch;
4584 a = working_space[i +
l][
ymin][zmin + 2 * sizez_ext] / maxch;
4596 if(i - l + 1 < xmin)
4597 a = working_space[
xmin][
ymin][zmin + 2 * sizez_ext] / maxch;
4600 a = working_space[i - l + 1][
ymin][zmin + 2 * sizez_ext] / maxch;
4614 a = working_space[i + 1][
ymin][zmin] = a * working_space[i][
ymin][zmin];
4617 for(i = ymin;i <
ymax; i++){
4618 nip = working_space[
xmin][i][zmin + 2 * sizez_ext] / maxch;
4619 nim = working_space[
xmin][i + 1][zmin + 2 * sizez_ext] / maxch;
4621 for(l = 1;l <= averWindow; l++){
4623 a = working_space[
xmin][
ymax][zmin + 2 * sizez_ext] / maxch;
4626 a = working_space[
xmin][i +
l][zmin + 2 * sizez_ext] / maxch;
4638 if(i - l + 1 < ymin)
4639 a = working_space[
xmin][
ymin][zmin + 2 * sizez_ext] / maxch;
4642 a = working_space[
xmin][i - l + 1][zmin + 2 * sizez_ext] / maxch;
4656 a = working_space[
xmin][i + 1][zmin] = a * working_space[
xmin][i][zmin];
4659 for(i = zmin;i < zmax;i++){
4660 nip = working_space[
xmin][
ymin][i + 2 * sizez_ext] / maxch;
4661 nim = working_space[
xmin][
ymin][i + 1 + 2 * sizez_ext] / maxch;
4663 for(l = 1;l <= averWindow;l++){
4665 a = working_space[
xmin][
ymin][zmax + 2 * sizez_ext] / maxch;
4668 a = working_space[
xmin][
ymin][i + l + 2 * sizez_ext] / maxch;
4680 if(i - l + 1 < zmin)
4681 a = working_space[
xmin][
ymin][zmin + 2 * sizez_ext] / maxch;
4684 a = working_space[
xmin][
ymin][i - l + 1 + 2 * sizez_ext] / maxch;
4701 for(i = xmin;i <
xmax; i++){
4702 for(j = ymin;j <
ymax; j++){
4703 nip = working_space[i][j + 1][zmin + 2 * sizez_ext] / maxch;
4704 nim = working_space[i + 1][j + 1][zmin + 2 * sizez_ext] / maxch;
4706 for(l = 1;l <= averWindow; l++){
4708 a = working_space[
xmax][j][zmin + 2 * sizez_ext] / maxch;
4711 a = working_space[i +
l][j][zmin + 2 * sizez_ext] / maxch;
4723 if(i - l + 1 < xmin)
4724 a = working_space[
xmin][j][zmin + 2 * sizez_ext] / maxch;
4727 a = working_space[i - l + 1][j][zmin + 2 * sizez_ext] / maxch;
4741 nip = working_space[i + 1][j][zmin + 2 * sizez_ext] / maxch;
4742 nim = working_space[i + 1][j + 1][zmin + 2 * sizez_ext] / maxch;
4743 for(l = 1;l <= averWindow; l++){
4745 a = working_space[i][
ymax][zmin + 2 * sizez_ext] / maxch;
4748 a = working_space[i][j +
l][zmin + 2 * sizez_ext] / maxch;
4760 if(j - l + 1 < ymin)
4761 a = working_space[i][
ymin][zmin + 2 * sizez_ext] / maxch;
4764 a = working_space[i][j - l + 1][zmin + 2 * sizez_ext] / maxch;
4777 a = (spx * working_space[i][j + 1][zmin] + spy * working_space[i + 1][j][zmin]) / (smx + smy);
4778 working_space[i + 1][j + 1][zmin] =
a;
4782 for(i = xmin;i <
xmax;i++){
4783 for(j = zmin;j < zmax;j++){
4784 nip = working_space[i][
ymin][j + 1 + 2 * sizez_ext] / maxch;
4785 nim = working_space[i + 1][
ymin][j + 1 + 2 * sizez_ext] / maxch;
4787 for(l = 1;l <= averWindow; l++){
4789 a = working_space[
xmax][
ymin][j + 2 * sizez_ext] / maxch;
4792 a = working_space[i +
l][
ymin][j + 2 * sizez_ext] / maxch;
4804 if(i - l + 1 < xmin)
4805 a = working_space[
xmin][
ymin][j + 2 * sizez_ext] / maxch;
4808 a = working_space[i - l + 1][
ymin][j + 2 * sizez_ext] / maxch;
4822 nip = working_space[i + 1][
ymin][j + 2 * sizez_ext] / maxch;
4823 nim = working_space[i + 1][
ymin][j + 1 + 2 * sizez_ext] / maxch;
4824 for(l = 1;l <= averWindow; l++){
4826 a = working_space[i][
ymin][zmax + 2 * sizez_ext] / maxch;
4829 a = working_space[i][
ymin][j + l + 2 * sizez_ext] / maxch;
4841 if(j - l + 1 < zmin)
4842 a = working_space[i][
ymin][zmin + 2 * sizez_ext] / maxch;
4845 a = working_space[i][
ymin][j - l + 1 + 2 * sizez_ext] / maxch;
4858 a = (spx * working_space[i][
ymin][j + 1] + spz * working_space[i + 1][
ymin][j]) / (smx + smz);
4859 working_space[i + 1][
ymin][j + 1] =
a;
4863 for(i = ymin;i <
ymax;i++){
4864 for(j = zmin;j < zmax;j++){
4865 nip = working_space[
xmin][i][j + 1 + 2 * sizez_ext] / maxch;
4866 nim = working_space[
xmin][i + 1][j + 1 + 2 * sizez_ext] / maxch;
4868 for(l = 1;l <= averWindow; l++){
4870 a = working_space[
xmin][
ymax][j + 2 * sizez_ext] / maxch;
4873 a = working_space[
xmin][i +
l][j + 2 * sizez_ext] / maxch;
4885 if(i - l + 1 < ymin)
4886 a = working_space[
xmin][
ymin][j + 2 * sizez_ext] / maxch;
4889 a = working_space[
xmin][i - l + 1][j + 2 * sizez_ext] / maxch;
4903 nip = working_space[
xmin][i + 1][j + 2 * sizez_ext] / maxch;
4904 nim = working_space[
xmin][i + 1][j + 1 + 2 * sizez_ext] / maxch;
4905 for(l = 1;l <= averWindow; l++){
4907 a = working_space[
xmin][i][zmax + 2 * sizez_ext] / maxch;
4910 a = working_space[
xmin][i][j + l + 2 * sizez_ext] / maxch;
4922 if(j - l + 1 < zmin)
4923 a = working_space[
xmin][i][zmin + 2 * sizez_ext] / maxch;
4926 a = working_space[
xmin][i][j - l + 1 + 2 * sizez_ext] / maxch;
4939 a = (spy * working_space[
xmin][i][j + 1] + spz * working_space[
xmin][i + 1][j]) / (smy + smz);
4940 working_space[
xmin][i + 1][j + 1] =
a;
4944 for(i = xmin;i <
xmax; i++){
4945 for(j = ymin;j <
ymax; j++){
4946 for(k = zmin;k < zmax; k++){
4947 nip = working_space[i][j + 1][k + 1 + 2 * sizez_ext] / maxch;
4948 nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
4950 for(l = 1;l <= averWindow; l++){
4952 a = working_space[
xmax][j][k + 2 * sizez_ext] / maxch;
4955 a = working_space[i +
l][j][k + 2 * sizez_ext] / maxch;
4967 if(i - l + 1 < xmin)
4968 a = working_space[
xmin][j][k + 2 * sizez_ext] / maxch;
4971 a = working_space[i - l + 1][j][k + 2 * sizez_ext] / maxch;
4985 nip = working_space[i + 1][j][k + 1 + 2 * sizez_ext] / maxch;
4986 nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
4987 for(l = 1;l <= averWindow; l++){
4989 a = working_space[i][
ymax][k + 2 * sizez_ext] / maxch;
4992 a = working_space[i][j +
l][k + 2 * sizez_ext] / maxch;
5004 if(j - l + 1 < ymin)
5005 a = working_space[i][
ymin][k + 2 * sizez_ext] / maxch;
5008 a = working_space[i][j - l + 1][k + 2 * sizez_ext] / maxch;
5022 nip = working_space[i + 1][j + 1][k + 2 * sizez_ext] / maxch;
5023 nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
5024 for(l = 1;l <= averWindow; l++ ){
5026 a = working_space[i][j][zmax + 2 * sizez_ext] / maxch;
5029 a = working_space[i][j][k + l + 2 * sizez_ext] / maxch;
5041 if(j - l + 1 < ymin)
5042 a = working_space[i][j][zmin + 2 * sizez_ext] / maxch;
5045 a = working_space[i][j][k - l + 1 + 2 * sizez_ext] / maxch;
5058 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);
5059 working_space[i + 1][j + 1][k + 1] =
a;
5065 for(i = xmin;i <=
xmax; i++){
5066 for(j = ymin;j <=
ymax; j++){
5067 for(k = zmin;k <= zmax; k++){
5068 working_space[i][j][k] = working_space[i][j][k] / nom;
5069 a+=working_space[i][j][k];
5073 for(i = 0;i < sizex_ext; i++){
5074 for(j = 0;j < sizey_ext; j++){
5075 for(k = 0;k < sizez_ext; k++){
5076 working_space[i][j][k + 2 * sizez_ext] = working_space[i][j][k] * plocha_markov /
a;
5083 for(k = 0;k < ssizez; k++){
5084 for(j = 0;j < ssizey; j++){
5085 for(i = 0;i < ssizex; i++){
5086 working_space[i][j][k] = 0;
5087 working_space[i][j][sizez_ext + k] = 0;
5088 if(working_space[i][j][k + 3 * sizez_ext] > maximum)
5089 maximum=working_space[i][j][k+3*sizez_ext];
5096 j = (
Int_t)(pocet_sigma * sigma + 0.5);
5097 for(i = -j;i <= j; i++){
5100 b = 2.0 * sigma * sigma;
5105 s = s - sigma * sigma;
5106 s = s / (sigma * sigma * sigma * sigma);
5108 c[PEAK_WINDOW / 2 + i] = s;
5115 c[i] = c[i] / norma;
5117 a = pocet_sigma * sigma + 0.5;
5120 zmax = sizez_ext - i - 1;
5122 ymax = sizey_ext - i - 1;
5124 xmax = sizex_ext - i - 1;
5125 lmin = PEAK_WINDOW / 2 - i;
5126 lmax = PEAK_WINDOW / 2 + i;
5127 for(i = xmin;i <=
xmax; i++){
5128 for(j = ymin;j <=
ymax; j++){
5129 for(k = zmin;k <= zmax; k++){
5131 for(li = lmin;li <= lmax; li++){
5132 for(lj = lmin;lj <= lmax; lj++){
5133 for(lk = lmin;lk <= lmax; lk++){
5134 a = working_space[i + li - PEAK_WINDOW / 2][j + lj - PEAK_WINDOW / 2][k + lk - PEAK_WINDOW / 2 + 2 * sizez_ext];
5135 b = c[li] * c[lj] * c[lk];
5141 working_space[i][j][k] = s;
5142 working_space[i][j][k + sizez_ext] =
TMath::Sqrt(f);
5146 for(x = xmin;x <=
xmax; x++){
5147 for(y = ymin + 1;y <
ymax; y++){
5148 for(z = zmin + 1;z < zmax; z++){
5149 val = working_space[
x][
y][z];
5150 val1 = working_space[x - 1][y - 1][z - 1];
5151 val2 = working_space[
x ][y - 1][z - 1];
5152 val3 = working_space[x + 1][y - 1][z - 1];
5153 val4 = working_space[x - 1][
y ][z - 1];
5154 val5 = working_space[
x ][
y ][z - 1];
5155 val6 = working_space[x + 1][
y ][z - 1];
5156 val7 = working_space[x - 1][y + 1][z - 1];
5157 val8 = working_space[
x ][y + 1][z - 1];
5158 val9 = working_space[x + 1][y + 1][z - 1];
5159 val10 = working_space[x - 1][y - 1][z ];
5160 val11 = working_space[
x ][y - 1][z ];
5161 val12 = working_space[x + 1][y - 1][z ];
5162 val13 = working_space[x - 1][
y ][z ];
5163 val14 = working_space[x + 1][
y ][z ];
5164 val15 = working_space[x - 1][y + 1][z ];
5165 val16 = working_space[
x ][y + 1][z ];
5166 val17 = working_space[x + 1][y + 1][z ];
5167 val18 = working_space[x - 1][y - 1][z + 1];
5168 val19 = working_space[
x ][y - 1][z + 1];
5169 val20 = working_space[x + 1][y - 1][z + 1];
5170 val21 = working_space[x - 1][
y ][z + 1];
5171 val22 = working_space[
x ][
y ][z + 1];
5172 val23 = working_space[x + 1][
y ][z + 1];
5173 val24 = working_space[x - 1][y + 1][z + 1];
5174 val25 = working_space[
x ][y + 1][z + 1];
5175 val26 = working_space[x + 1][y + 1][z + 1];
5176 f = -s_f_ratio_peaks * working_space[
x][
y][z + sizez_ext];
5177 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){
5179 for(li = lmin;li <= lmax; li++){
5180 a = working_space[x + li - PEAK_WINDOW / 2][
y][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 + li - PEAK_WINDOW / 2][z + 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 a = working_space[
x][
y][z + li - PEAK_WINDOW / 2 + 2 * sizez_ext];
5198 f += a * c[li] * c[li];
5200 f = -s_f_ratio_peaks *
sqrt(f);
5203 for(li = lmin;li <= lmax; li++){
5204 for(lj = lmin;lj <= lmax; lj++){
5205 a = working_space[x + li - PEAK_WINDOW / 2][y + lj - PEAK_WINDOW / 2][z + 2 * sizez_ext];
5206 s += a * c[li] * c[lj];
5207 f += a * c[li] * c[li] * c[lj] * c[lj];
5210 f = s_f_ratio_peaks *
sqrt(f);
5213 for(li = lmin;li <= lmax; li++){
5214 for(lj = lmin;lj <= lmax; lj++){
5215 a = working_space[x + li - PEAK_WINDOW / 2][
y][z + lj - PEAK_WINDOW / 2 + 2 * sizez_ext];
5216 s += a * c[li] * c[lj];
5217 f += a * c[li] * c[li] * c[lj] * c[lj];
5220 f = s_f_ratio_peaks *
sqrt(f);
5223 for(li = lmin;li <= lmax; li++){
5224 for(lj=lmin;lj<=lmax;lj++){
5225 a = working_space[
x][y + li - PEAK_WINDOW / 2][z + lj - PEAK_WINDOW / 2 + 2 * sizez_ext];
5226 s += a * c[li] * c[lj];
5227 f += a * c[li] * c[li] * c[lj] * c[lj];
5230 f = s_f_ratio_peaks *
sqrt(f);
5232 if(x >= shift && x < ssizex + shift && y >= shift && y < ssizey + shift && z >= shift && z < ssizez + shift){
5233 if(working_space[x][y][z + 3 * sizez_ext] > threshold * maximum / 100.0){
5235 for(k = x - 1,a = 0,b = 0;k <= x + 1; k++){
5236 a += (
Double_t)(k - shift) * working_space[k][
y][z];
5237 b += working_space[k][
y][z];
5240 for(k = y - 1,a = 0,b = 0;k <= y + 1; k++){
5241 a += (
Double_t)(k - shift) * working_space[
x][k][z];
5242 b += working_space[
x][k][z];
5245 for(k = z - 1,a = 0,b = 0;k <= z + 1; k++){
5246 a += (
Double_t)(k - shift) * working_space[
x][
y][k];
5247 b += working_space[
x][
y][k];
5264 for(i = 0;i < ssizex; i++){
5265 for(j = 0;j < ssizey; j++){
5266 for(k = 0;k < ssizez; k++){
5267 val = -working_space[i + shift][j + shift][k + shift];
5270 dest[i][j][k] = val;
5274 k = (
Int_t)(4 * sigma + 0.5);
5276 for(i = 0;i < ssizex + k; i++){
5277 for(j = 0;j < ssizey + k; j++)
5278 delete[] working_space[i][j];
5279 delete[] working_space[i];
5281 delete[] working_space;
void SetResolution(Double_t resolution=1)
resolution: determines resolution of the neighboring peaks default value is 1 correspond to 3 sigma d...
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
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)
virtual Int_t GetDimension() const
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...
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)
ClassImp(TSpectrum3) TSpectrum3
Constructor.
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
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 void Error(const char *method, const char *msgfmt,...) const
Issue error message.
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 const char * GetName() const
Returns name of object.
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)
virtual void Print(Option_t *option="") const
Print the array of positions.
ClassImp(TMCParticle) void TMCParticle printf(": p=(%7.3f,%7.3f,%9.3f) ;", fPx, fPy, fPz)
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
#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 //...
unsigned int r2[N_CITIES]
virtual ~TSpectrum3()
Destructor.
ROOT 6.06/09 - Reference Guide Generated on Sun Dec 10 2017 22:25:14 using Doxygen 1.8.9.
