44 #define PEAK_WINDOW 1024 71 :
TNamed(
"Spectrum",
"Miroslav Morhac peak finder")
155 if (h == 0)
return 0;
158 Error(
"Search",
"Only implemented for 1-d histograms");
185 Int_t size = last-first+1;
188 for (i = 0; i < size; i++) source[i] = h->
GetBinContent(i + first);
191 Background(source,size,numberIterations, direction, filterOrder,smoothing,
192 smoothWindow,compton);
197 char *
hbname =
new char[nch+20];
201 hb->GetListOfFunctions()->
Delete();
203 for (i=0; i< size; i++) hb->SetBinContent(i+first,source[i]);
204 hb->SetEntries(size);
208 if (
gPad)
delete gPad->GetPrimitive(hbname);
221 printf(
"\nNumber of positions = %d\n",
fNPeaks);
269 if (hin == 0)
return 0;
272 Error(
"Search",
"Only implemented for 1-d and 2-d histograms");
275 if (threshold <=0 || threshold >= 1) {
276 Warning(
"Search",
"threshold must 0<threshold<1, threshold=0.05 assumed");
296 if (dimension == 1) {
299 Int_t size = last-first+1;
300 Int_t i, bin, npeaks;
303 for (i = 0; i < size; i++) source[i] = hin->
GetBinContent(i + first);
306 if (sigma < 1) sigma = 1;
307 if (sigma > 8) sigma = 8;
309 npeaks =
SearchHighRes(source, dest, size, sigma, 100*threshold,
312 for (i = 0; i < npeaks; i++) {
322 if (!npeaks)
return 0;
779 int numberIterations,
780 int direction,
int filterOrder,
781 bool smoothing,
int smoothWindow,
784 int i, j, w, bw, b1, b2, priz;
785 Double_t a,
b,
c, d,
e, yb1, yb2, ai, av, men, b4, c4, d4, e4, b6, c6, d6, e6, f6, g6, b8, c8, d8, e8, f8, g8, h8, i8;
787 return "Wrong Parameters";
788 if (numberIterations < 1)
789 return "Width of Clipping Window Must Be Positive";
790 if (ssize < 2 * numberIterations + 1)
791 return "Too Large Clipping Window";
793 return "Incorrect width of smoothing window";
795 for (i = 0; i < ssize; i++){
796 working_space[i] = spectrum[i];
797 working_space[i + ssize] = spectrum[i];
799 bw=(smoothWindow-1)/2;
803 i = numberIterations;
806 for (j = i; j < ssize - i; j++) {
808 a = working_space[ssize + j];
809 b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
812 working_space[j] =
a;
815 else if (smoothing ==
kTRUE){
816 a = working_space[ssize + j];
819 for (w = j - bw; w <= j + bw; w++){
820 if ( w >= 0 && w < ssize){
821 av += working_space[ssize + w];
828 for (w = j - i - bw; w <= j - i + bw; w++){
829 if ( w >= 0 && w < ssize){
830 b += working_space[ssize + w];
837 for (w = j + i - bw; w <= j + i + bw; w++){
838 if ( w >= 0 && w < ssize){
839 c += working_space[ssize + w];
850 for (j = i; j < ssize - i; j++)
851 working_space[ssize + j] = working_space[j];
861 for (j = i; j < ssize - i; j++) {
863 a = working_space[ssize + j];
864 b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
867 c -= working_space[ssize + j - (
Int_t) (2 * ai)] / 6;
868 c += 4 * working_space[ssize + j - (
Int_t) ai] / 6;
869 c += 4 * working_space[ssize + j + (
Int_t) ai] / 6;
870 c -= working_space[ssize + j + (
Int_t) (2 * ai)] / 6;
875 working_space[j] =
a;
878 else if (smoothing ==
kTRUE){
879 a = working_space[ssize + j];
882 for (w = j - bw; w <= j + bw; w++){
883 if ( w >= 0 && w < ssize){
884 av += working_space[ssize + w];
891 for (w = j - i - bw; w <= j - i + bw; w++){
892 if ( w >= 0 && w < ssize){
893 b += working_space[ssize + w];
900 for (w = j + i - bw; w <= j + i + bw; w++){
901 if ( w >= 0 && w < ssize){
902 c += working_space[ssize + w];
910 for (w = j - (
Int_t)(2 * ai) - bw; w <= j - (
Int_t)(2 * ai) + bw; w++){
911 if (w >= 0 && w < ssize){
912 b4 += working_space[ssize + w];
918 for (w = j - (
Int_t)ai - bw; w <= j - (
Int_t)ai + bw; w++){
919 if (w >= 0 && w < ssize){
920 c4 += working_space[ssize + w];
926 for (w = j + (
Int_t)ai - bw; w <= j + (
Int_t)ai + bw; w++){
927 if (w >= 0 && w < ssize){
928 d4 += working_space[ssize + w];
934 for (w = j + (
Int_t)(2 * ai) - bw; w <= j + (
Int_t)(2 * ai) + bw; w++){
935 if (w >= 0 && w < ssize){
936 e4 += working_space[ssize + w];
941 b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6;
949 for (j = i; j < ssize - i; j++)
950 working_space[ssize + j] = working_space[j];
960 for (j = i; j < ssize - i; j++) {
962 a = working_space[ssize + j];
963 b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
966 c -= working_space[ssize + j - (
Int_t) (2 * ai)] / 6;
967 c += 4 * working_space[ssize + j - (
Int_t) ai] / 6;
968 c += 4 * working_space[ssize + j + (
Int_t) ai] / 6;
969 c -= working_space[ssize + j + (
Int_t) (2 * ai)] / 6;
972 d += working_space[ssize + j - (
Int_t) (3 * ai)] / 20;
973 d -= 6 * working_space[ssize + j - (
Int_t) (2 * ai)] / 20;
974 d += 15 * working_space[ssize + j - (
Int_t) ai] / 20;
975 d += 15 * working_space[ssize + j + (
Int_t) ai] / 20;
976 d -= 6 * working_space[ssize + j + (
Int_t) (2 * ai)] / 20;
977 d += working_space[ssize + j + (
Int_t) (3 * ai)] / 20;
984 working_space[j] =
a;
987 else if (smoothing ==
kTRUE){
988 a = working_space[ssize + j];
991 for (w = j - bw; w <= j + bw; w++){
992 if ( w >= 0 && w < ssize){
993 av += working_space[ssize + w];
1000 for (w = j - i - bw; w <= j - i + bw; w++){
1001 if ( w >= 0 && w < ssize){
1002 b += working_space[ssize + w];
1009 for (w = j + i - bw; w <= j + i + bw; w++){
1010 if ( w >= 0 && w < ssize){
1011 c += working_space[ssize + w];
1019 for (w = j - (
Int_t)(2 * ai) - bw; w <= j - (
Int_t)(2 * ai) + bw; w++){
1020 if (w >= 0 && w < ssize){
1021 b4 += working_space[ssize + w];
1027 for (w = j - (
Int_t)ai - bw; w <= j - (
Int_t)ai + bw; w++){
1028 if (w >= 0 && w < ssize){
1029 c4 += working_space[ssize + w];
1035 for (w = j + (
Int_t)ai - bw; w <= j + (
Int_t)ai + bw; w++){
1036 if (w >= 0 && w < ssize){
1037 d4 += working_space[ssize + w];
1043 for (w = j + (
Int_t)(2 * ai) - bw; w <= j + (
Int_t)(2 * ai) + bw; w++){
1044 if (w >= 0 && w < ssize){
1045 e4 += working_space[ssize + w];
1050 b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6;
1053 for (w = j - (
Int_t)(3 * ai) - bw; w <= j - (
Int_t)(3 * ai) + bw; w++){
1054 if (w >= 0 && w < ssize){
1055 b6 += working_space[ssize + w];
1061 for (w = j - (
Int_t)(2 * ai) - bw; w <= j - (
Int_t)(2 * ai) + bw; w++){
1062 if (w >= 0 && w < ssize){
1063 c6 += working_space[ssize + w];
1069 for (w = j - (
Int_t)ai - bw; w <= j - (
Int_t)ai + bw; w++){
1070 if (w >= 0 && w < ssize){
1071 d6 += working_space[ssize + w];
1077 for (w = j + (
Int_t)ai - bw; w <= j + (
Int_t)ai + bw; w++){
1078 if (w >= 0 && w < ssize){
1079 e6 += working_space[ssize + w];
1085 for (w = j + (
Int_t)(2 * ai) - bw; w <= j + (
Int_t)(2 * ai) + bw; w++){
1086 if (w >= 0 && w < ssize){
1087 f6 += working_space[ssize + w];
1093 for (w = j + (
Int_t)(3 * ai) - bw; w <= j + (
Int_t)(3 * ai) + bw; w++){
1094 if (w >= 0 && w < ssize){
1095 g6 += working_space[ssize + w];
1100 b6 = (b6 - 6 * c6 + 15 * d6 + 15 * e6 - 6 * f6 + g6) / 20;
1107 working_space[j]=av;
1110 for (j = i; j < ssize - i; j++)
1111 working_space[ssize + j] = working_space[j];
1121 for (j = i; j < ssize - i; j++) {
1122 if (smoothing ==
kFALSE){
1123 a = working_space[ssize + j];
1124 b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
1127 c -= working_space[ssize + j - (
Int_t) (2 * ai)] / 6;
1128 c += 4 * working_space[ssize + j - (
Int_t) ai] / 6;
1129 c += 4 * working_space[ssize + j + (
Int_t) ai] / 6;
1130 c -= working_space[ssize + j + (
Int_t) (2 * ai)] / 6;
1133 d += working_space[ssize + j - (
Int_t) (3 * ai)] / 20;
1134 d -= 6 * working_space[ssize + j - (
Int_t) (2 * ai)] / 20;
1135 d += 15 * working_space[ssize + j - (
Int_t) ai] / 20;
1136 d += 15 * working_space[ssize + j + (
Int_t) ai] / 20;
1137 d -= 6 * working_space[ssize + j + (
Int_t) (2 * ai)] / 20;
1138 d += working_space[ssize + j + (
Int_t) (3 * ai)] / 20;
1141 e -= working_space[ssize + j - (
Int_t) (4 * ai)] / 70;
1142 e += 8 * working_space[ssize + j - (
Int_t) (3 * ai)] / 70;
1143 e -= 28 * working_space[ssize + j - (
Int_t) (2 * ai)] / 70;
1144 e += 56 * working_space[ssize + j - (
Int_t) ai] / 70;
1145 e += 56 * working_space[ssize + j + (
Int_t) ai] / 70;
1146 e -= 28 * working_space[ssize + j + (
Int_t) (2 * ai)] / 70;
1147 e += 8 * working_space[ssize + j + (
Int_t) (3 * ai)] / 70;
1148 e -= working_space[ssize + j + (
Int_t) (4 * ai)] / 70;
1157 working_space[j] =
a;
1160 else if (smoothing ==
kTRUE){
1161 a = working_space[ssize + j];
1164 for (w = j - bw; w <= j + bw; w++){
1165 if ( w >= 0 && w < ssize){
1166 av += working_space[ssize + w];
1173 for (w = j - i - bw; w <= j - i + bw; w++){
1174 if ( w >= 0 && w < ssize){
1175 b += working_space[ssize + w];
1182 for (w = j + i - bw; w <= j + i + bw; w++){
1183 if ( w >= 0 && w < ssize){
1184 c += working_space[ssize + w];
1192 for (w = j - (
Int_t)(2 * ai) - bw; w <= j - (
Int_t)(2 * ai) + bw; w++){
1193 if (w >= 0 && w < ssize){
1194 b4 += working_space[ssize + w];
1200 for (w = j - (
Int_t)ai - bw; w <= j - (
Int_t)ai + bw; w++){
1201 if (w >= 0 && w < ssize){
1202 c4 += working_space[ssize + w];
1208 for (w = j + (
Int_t)ai - bw; w <= j + (
Int_t)ai + bw; w++){
1209 if (w >= 0 && w < ssize){
1210 d4 += working_space[ssize + w];
1216 for (w = j + (
Int_t)(2 * ai) - bw; w <= j + (
Int_t)(2 * ai) + bw; w++){
1217 if (w >= 0 && w < ssize){
1218 e4 += working_space[ssize + w];
1223 b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6;
1226 for (w = j - (
Int_t)(3 * ai) - bw; w <= j - (
Int_t)(3 * ai) + bw; w++){
1227 if (w >= 0 && w < ssize){
1228 b6 += working_space[ssize + w];
1234 for (w = j - (
Int_t)(2 * ai) - bw; w <= j - (
Int_t)(2 * ai) + bw; w++){
1235 if (w >= 0 && w < ssize){
1236 c6 += working_space[ssize + w];
1242 for (w = j - (
Int_t)ai - bw; w <= j - (
Int_t)ai + bw; w++){
1243 if (w >= 0 && w < ssize){
1244 d6 += working_space[ssize + w];
1250 for (w = j + (
Int_t)ai - bw; w <= j + (
Int_t)ai + bw; w++){
1251 if (w >= 0 && w < ssize){
1252 e6 += working_space[ssize + w];
1258 for (w = j + (
Int_t)(2 * ai) - bw; w <= j + (
Int_t)(2 * ai) + bw; w++){
1259 if (w >= 0 && w < ssize){
1260 f6 += working_space[ssize + w];
1266 for (w = j + (
Int_t)(3 * ai) - bw; w <= j + (
Int_t)(3 * ai) + bw; w++){
1267 if (w >= 0 && w < ssize){
1268 g6 += working_space[ssize + w];
1273 b6 = (b6 - 6 * c6 + 15 * d6 + 15 * e6 - 6 * f6 + g6) / 20;
1276 for (w = j - (
Int_t)(4 * ai) - bw; w <= j - (
Int_t)(4 * ai) + bw; w++){
1277 if (w >= 0 && w < ssize){
1278 b8 += working_space[ssize + w];
1284 for (w = j - (
Int_t)(3 * ai) - bw; w <= j - (
Int_t)(3 * ai) + bw; w++){
1285 if (w >= 0 && w < ssize){
1286 c8 += working_space[ssize + w];
1292 for (w = j - (
Int_t)(2 * ai) - bw; w <= j - (
Int_t)(2 * ai) + bw; w++){
1293 if (w >= 0 && w < ssize){
1294 d8 += working_space[ssize + w];
1300 for (w = j - (
Int_t)ai - bw; w <= j - (
Int_t)ai + bw; w++){
1301 if (w >= 0 && w < ssize){
1302 e8 += working_space[ssize + w];
1308 for (w = j + (
Int_t)ai - bw; w <= j + (
Int_t)ai + bw; w++){
1309 if (w >= 0 && w < ssize){
1310 f8 += working_space[ssize + w];
1316 for (w = j + (
Int_t)(2 * ai) - bw; w <= j + (
Int_t)(2 * ai) + bw; w++){
1317 if (w >= 0 && w < ssize){
1318 g8 += working_space[ssize + w];
1324 for (w = j + (
Int_t)(3 * ai) - bw; w <= j + (
Int_t)(3 * ai) + bw; w++){
1325 if (w >= 0 && w < ssize){
1326 h8 += working_space[ssize + w];
1332 for (w = j + (
Int_t)(4 * ai) - bw; w <= j + (
Int_t)(4 * ai) + bw; w++){
1333 if (w >= 0 && w < ssize){
1334 i8 += working_space[ssize + w];
1339 b8 = ( -b8 + 8 * c8 - 28 * d8 + 56 * e8 - 56 * f8 - 28 * g8 + 8 * h8 - i8)/70;
1348 working_space[j]=av;
1351 for (j = i; j < ssize - i; j++)
1352 working_space[ssize + j] = working_space[j];
1360 if (compton ==
kTRUE) {
1361 for (i = 0, b2 = 0; i < ssize; i++){
1363 a = working_space[i], b = spectrum[i];
1369 yb1 = working_space[b1];
1370 for (b2 = b1 + 1, c = 0, priz = 0; priz == 0 && b2 < ssize; b2++){
1371 a = working_space[b2], b = spectrum[b2];
1380 yb2 = working_space[b2];
1382 for (j = b1, c = 0; j <= b2; j++){
1387 c = (yb2 - yb1) / c;
1388 for (j = b1, d = 0; j <= b2 && j < ssize; j++){
1392 working_space[ssize + j] =
a;
1398 for (j = b2, c = 0; j >= b1; j--){
1403 c = (yb1 - yb2) / c;
1404 for (j = b2, d = 0;j >= b1 && j >= 0; j--){
1408 working_space[ssize + j] =
a;
1417 for (j = 0; j < ssize; j++)
1418 spectrum[j] = working_space[ssize + j];
1419 delete[]working_space;
1493 Double_t nom, nip, nim, sp, sm, area = 0;
1495 return "Averaging Window must be positive";
1497 xmin = 0,xmax = ssize - 1;
1498 for(i = 0, maxch = 0; i < ssize; i++){
1500 if(maxch < source[i])
1506 delete [] working_space;
1511 working_space[
xmin] = 1;
1512 for(i = xmin; i <
xmax; i++){
1513 nip = source[i] / maxch;
1514 nim = source[i + 1] / maxch;
1516 for(l = 1; l <= averWindow; l++){
1518 a = source[
xmax] / maxch;
1521 a = source[i +
l] / maxch;
1531 if((i - l + 1) < xmin)
1532 a = source[
xmin] / maxch;
1535 a = source[i - l + 1] / maxch;
1546 a = working_space[i + 1] = working_space[i] *
a;
1549 for(i = xmin; i <=
xmax; i++){
1550 working_space[i] = working_space[i] / nom;
1552 for(i = 0; i < ssize; i++)
1553 source[i] = working_space[i] * area;
1554 delete[]working_space;
1859 int ssize,
int numberIterations,
1863 return "Wrong Parameters";
1865 if (numberRepetitions <= 0)
1866 return "Wrong Parameters";
1871 int i, j, k, lindex, posit, lh_gold,
l, repet;
1872 Double_t lda, ldb, ldc, area, maximum;
1879 for (i = 0; i < ssize; i++) {
1883 working_space[i] = lda;
1885 if (lda > maximum) {
1890 if (lh_gold == -1) {
1891 delete [] working_space;
1892 return "ZERO RESPONSE VECTOR";
1896 for (i = 0; i < ssize; i++)
1897 working_space[2 * ssize + i] = source[i];
1900 for (i = 0; i < ssize; i++){
1902 for (j = 0; j < ssize; j++){
1903 ldb = working_space[j];
1906 ldc = working_space[k];
1907 lda = lda + ldb * ldc;
1910 working_space[ssize + i] = lda;
1912 for (k = 0; k < ssize; k++){
1915 ldb = working_space[
l];
1916 ldc = working_space[2 * ssize + k];
1917 lda = lda + ldb * ldc;
1920 working_space[3 * ssize + i]=lda;
1924 for (i = 0; i < ssize; i++){
1925 working_space[2 * ssize + i] = working_space[3 * ssize + i];
1929 for (i = 0; i < ssize; i++)
1930 working_space[i] = 1;
1933 for (repet = 0; repet < numberRepetitions; repet++) {
1935 for (i = 0; i < ssize; i++)
1936 working_space[i] =
TMath::Power(working_space[i], boost);
1938 for (lindex = 0; lindex < numberIterations; lindex++) {
1939 for (i = 0; i < ssize; i++) {
1940 if (working_space[2 * ssize + i] > 0.000001
1941 && working_space[i] > 0.000001) {
1943 for (j = 0; j < lh_gold; j++) {
1944 ldb = working_space[j + ssize];
1949 ldc = working_space[k];
1952 ldc += working_space[k];
1956 ldc = working_space[i];
1957 lda = lda + ldb * ldc;
1959 ldb = working_space[2 * ssize + i];
1965 ldb = working_space[i];
1967 working_space[3 * ssize + i] = lda;
1970 for (i = 0; i < ssize; i++)
1971 working_space[i] = working_space[3 * ssize + i];
1976 for (i = 0; i < ssize; i++) {
1977 lda = working_space[i];
1980 working_space[ssize + j] = lda;
1984 for (i = 0; i < ssize; i++)
1985 source[i] = area * working_space[ssize + i];
1986 delete[]working_space;
2139 int ssize,
int numberIterations,
2143 return "Wrong Parameters";
2145 if (numberRepetitions <= 0)
2146 return "Wrong Parameters";
2151 int i, j, k, lindex, posit, lh_gold, repet, kmin, kmax;
2158 for (i = 0; i < ssize; i++) {
2162 working_space[ssize + i] = lda;
2163 if (lda > maximum) {
2168 if (lh_gold == -1) {
2169 delete [] working_space;
2170 return "ZERO RESPONSE VECTOR";
2174 for (i = 0; i < ssize; i++)
2175 working_space[2 * ssize + i] = source[i];
2178 for (i = 0; i < ssize; i++){
2179 if (i <= ssize - lh_gold)
2180 working_space[i] = 1;
2183 working_space[i] = 0;
2187 for (repet = 0; repet < numberRepetitions; repet++) {
2189 for (i = 0; i < ssize; i++)
2190 working_space[i] =
TMath::Power(working_space[i], boost);
2192 for (lindex = 0; lindex < numberIterations; lindex++) {
2193 for (i = 0; i <= ssize - lh_gold; i++){
2195 if (working_space[i] > 0){
2196 for (j = i; j < i + lh_gold; j++){
2197 ldb = working_space[2 * ssize + j];
2201 if (kmax > lh_gold - 1)
2203 kmin = j + lh_gold - ssize;
2207 for (k = kmax; k >= kmin; k--){
2208 ldc += working_space[ssize + k] * working_space[j - k];
2216 ldb = ldb * working_space[ssize + j - i];
2220 lda = lda * working_space[i];
2222 working_space[3 * ssize + i] = lda;
2224 for (i = 0; i < ssize; i++)
2225 working_space[i] = working_space[3 * ssize + i];
2230 for (i = 0; i < ssize; i++) {
2231 lda = working_space[i];
2234 working_space[ssize + j] = lda;
2238 for (i = 0; i < ssize; i++)
2239 source[i] = working_space[ssize + i];
2240 delete[]working_space;
2365 int ssizex,
int ssizey,
2366 int numberIterations,
2369 int i, j, k, lindex, lhx = 0, repet;
2371 if (ssizex <= 0 || ssizey <= 0)
2372 return "Wrong Parameters";
2373 if (ssizex < ssizey)
2374 return "Sizex must be greater than sizey)";
2375 if (numberIterations <= 0)
2376 return "Number of iterations must be positive";
2378 new Double_t[ssizex * ssizey + 2 * ssizey * ssizey + 4 * ssizex];
2381 for (j = 0; j < ssizey && lhx != -1; j++) {
2384 for (i = 0; i < ssizex; i++) {
2385 lda = respMatrix[j][i];
2389 working_space[j * ssizex + i] = lda;
2393 for (i = 0; i < ssizex; i++)
2394 working_space[j * ssizex + i] /= area;
2398 delete [] working_space;
2399 return (
"ZERO COLUMN IN RESPONSE MATRIX");
2403 for (i = 0; i < ssizex; i++)
2404 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i] =
2408 for (i = 0; i < ssizey; i++) {
2409 for (j = 0; j < ssizey; j++) {
2411 for (k = 0; k < ssizex; k++) {
2412 ldb = working_space[ssizex * i + k];
2413 ldc = working_space[ssizex * j + k];
2414 lda = lda + ldb * ldc;
2416 working_space[ssizex * ssizey + ssizey * i + j] = lda;
2419 for (k = 0; k < ssizex; k++) {
2420 ldb = working_space[ssizex * i + k];
2422 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex +
2424 lda = lda + ldb * ldc;
2426 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i] =
2431 for (i = 0; i < ssizey; i++)
2432 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i] =
2433 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i];
2436 for (i = 0; i < ssizey; i++) {
2437 for (j = 0; j < ssizey; j++) {
2439 for (k = 0; k < ssizey; k++) {
2440 ldb = working_space[ssizex * ssizey + ssizey * i + k];
2441 ldc = working_space[ssizex * ssizey + ssizey * j + k];
2442 lda = lda + ldb * ldc;
2444 working_space[ssizex * ssizey + ssizey * ssizey + ssizey * i + j] =
2448 for (k = 0; k < ssizey; k++) {
2449 ldb = working_space[ssizex * ssizey + ssizey * i + k];
2451 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex +
2453 lda = lda + ldb * ldc;
2455 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i] =
2460 for (i = 0; i < ssizey; i++)
2461 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i] =
2462 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i];
2465 for (i = 0; i < ssizey; i++)
2466 working_space[ssizex * ssizey + 2 * ssizey * ssizey + i] = 1;
2469 for (repet = 0; repet < numberRepetitions; repet++) {
2471 for (i = 0; i < ssizey; i++)
2472 working_space[ssizex * ssizey + 2 * ssizey * ssizey + i] =
TMath::Power(working_space[ssizex * ssizey + 2 * ssizey * ssizey + i], boost);
2474 for (lindex = 0; lindex < numberIterations; lindex++) {
2475 for (i = 0; i < ssizey; i++) {
2477 for (j = 0; j < ssizey; j++) {
2479 working_space[ssizex * ssizey + ssizey * ssizey + ssizey * i + j];
2480 ldc = working_space[ssizex * ssizey + 2 * ssizey * ssizey + j];
2481 lda = lda + ldb * ldc;
2484 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i];
2491 ldb = working_space[ssizex * ssizey + 2 * ssizey * ssizey + i];
2493 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i] = lda;
2495 for (i = 0; i < ssizey; i++)
2496 working_space[ssizex * ssizey + 2 * ssizey * ssizey + i] =
2497 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i];
2502 for (i = 0; i < ssizex; i++) {
2504 source[i] = working_space[ssizex * ssizey + 2 * ssizey * ssizey + i];
2509 delete[]working_space;
2737 bool backgroundRemove,
int deconIterations,
2738 bool markov,
int averWindow)
2740 int i, j, numberIterations = (
Int_t)(7 * sigma + 0.5);
2742 int k, lindex, posit, imin, imax, jmin, jmax, lh_gold, priz;
2743 Double_t lda, ldb, ldc, area, maximum, maximum_decon;
2744 int xmin,
xmax,
l, peak_index = 0, size_ext = ssize + 2 * numberIterations, shift = numberIterations, bw = 2, w;
2746 Double_t nom, nip, nim, sp, sm, plocha = 0;
2747 Double_t m0low=0,m1low=0,m2low=0,l0low=0,l1low=0,detlow,av,men;
2749 Error(
"SearchHighRes",
"Invalid sigma, must be greater than or equal to 1");
2753 if(threshold<=0 || threshold>=100){
2754 Error(
"SearchHighRes",
"Invalid threshold, must be positive and less than 100");
2758 j = (
Int_t) (5.0 * sigma + 0.5);
2760 Error(
"SearchHighRes",
"Too large sigma");
2764 if (markov ==
true) {
2765 if (averWindow <= 0) {
2766 Error(
"SearchHighRes",
"Averaging window must be positive");
2771 if(backgroundRemove ==
true){
2772 if(ssize < 2 * numberIterations + 1){
2773 Error(
"SearchHighRes",
"Too large clipping window");
2778 k = int(2 * sigma+0.5);
2780 for(i = 0;i < k;i++){
2781 a = i,b = source[i];
2782 m0low += 1,m1low +=
a,m2low += a *
a,l0low +=
b,l1low += a *
b;
2784 detlow = m0low * m2low - m1low * m1low;
2786 l1low = (-l0low * m1low + l1low * m0low) / detlow;
2798 i = (
Int_t)(7 * sigma + 0.5);
2801 for (j=0;j<7 * (ssize + i);j++) working_space[j] = 0;
2802 for(i = 0; i < size_ext; i++){
2805 working_space[i + size_ext] = source[0] + l1low *
a;
2806 if(working_space[i + size_ext] < 0)
2807 working_space[i + size_ext]=0;
2810 else if(i >= ssize + shift){
2811 a = i - (ssize - 1 + shift);
2812 working_space[i + size_ext] = source[ssize - 1];
2813 if(working_space[i + size_ext] < 0)
2814 working_space[i + size_ext]=0;
2818 working_space[i + size_ext] = source[i - shift];
2821 if(backgroundRemove ==
true){
2822 for(i = 1; i <= numberIterations; i++){
2823 for(j = i; j < size_ext - i; j++){
2824 if(markov ==
false){
2825 a = working_space[size_ext + j];
2826 b = (working_space[size_ext + j - i] + working_space[size_ext + j + i]) / 2.0;
2834 a = working_space[size_ext + j];
2837 for (w = j - bw; w <= j + bw; w++){
2838 if ( w >= 0 && w < size_ext){
2839 av += working_space[size_ext + w];
2846 for (w = j - i - bw; w <= j - i + bw; w++){
2847 if ( w >= 0 && w < size_ext){
2848 b += working_space[size_ext + w];
2855 for (w = j + i - bw; w <= j + i + bw; w++){
2856 if ( w >= 0 && w < size_ext){
2857 c += working_space[size_ext + w];
2865 working_space[j]=av;
2868 for(j = i; j < size_ext - i; j++)
2869 working_space[size_ext + j] = working_space[j];
2871 for(j = 0;j < size_ext; j++){
2874 b = source[0] + l1low *
a;
2876 working_space[size_ext + j] = b - working_space[size_ext + j];
2879 else if(j >= ssize + shift){
2880 a = j - (ssize - 1 + shift);
2881 b = source[ssize - 1];
2883 working_space[size_ext + j] = b - working_space[size_ext + j];
2887 working_space[size_ext + j] = source[j - shift] - working_space[size_ext + j];
2890 for(j = 0;j < size_ext; j++){
2891 if(working_space[size_ext + j] < 0) working_space[size_ext + j] = 0;
2895 for(i = 0; i < size_ext; i++){
2896 working_space[i + 6*size_ext] = working_space[i + size_ext];
2900 for(j = 0; j < size_ext; j++)
2901 working_space[2 * size_ext + j] = working_space[size_ext + j];
2902 xmin = 0,xmax = size_ext - 1;
2903 for(i = 0, maxch = 0; i < size_ext; i++){
2904 working_space[i] = 0;
2905 if(maxch < working_space[2 * size_ext + i])
2906 maxch = working_space[2 * size_ext + i];
2907 plocha += working_space[2 * size_ext + i];
2910 delete [] working_space;
2915 working_space[
xmin] = 1;
2916 for(i = xmin; i <
xmax; i++){
2917 nip = working_space[2 * size_ext + i] / maxch;
2918 nim = working_space[2 * size_ext + i + 1] / maxch;
2920 for(l = 1; l <= averWindow; l++){
2922 a = working_space[2 * size_ext + xmax] / maxch;
2925 a = working_space[2 * size_ext + i +
l] / maxch;
2937 if((i - l + 1) < xmin)
2938 a = working_space[2 * size_ext +
xmin] / maxch;
2941 a = working_space[2 * size_ext + i - l + 1] / maxch;
2955 a = working_space[i + 1] = working_space[i] *
a;
2958 for(i = xmin; i <=
xmax; i++){
2959 working_space[i] = working_space[i] / nom;
2961 for(j = 0; j < size_ext; j++)
2962 working_space[size_ext + j] = working_space[j] * plocha;
2963 for(j = 0; j < size_ext; j++){
2964 working_space[2 * size_ext + j] = working_space[size_ext + j];
2966 if(backgroundRemove ==
true){
2967 for(i = 1; i <= numberIterations; i++){
2968 for(j = i; j < size_ext - i; j++){
2969 a = working_space[size_ext + j];
2970 b = (working_space[size_ext + j - i] + working_space[size_ext + j + i]) / 2.0;
2973 working_space[j] =
a;
2975 for(j = i; j < size_ext - i; j++)
2976 working_space[size_ext + j] = working_space[j];
2978 for(j = 0; j < size_ext; j++){
2979 working_space[size_ext + j] = working_space[2 * size_ext + j] - working_space[size_ext + j];
2989 for(i = 0; i < size_ext; i++){
2991 lda = lda * lda / (2 * sigma *
sigma);
2997 working_space[i] = lda;
3005 for(i = 0; i < size_ext; i++)
3006 working_space[2 * size_ext + i] =
TMath::Abs(working_space[size_ext + i]);
3013 for(i = imin; i <= imax; i++){
3018 jmax = lh_gold - 1 - i;
3019 if(jmax > (lh_gold - 1))
3022 for(j = jmin;j <= jmax; j++){
3023 ldb = working_space[j];
3024 ldc = working_space[i + j];
3025 lda = lda + ldb * ldc;
3027 working_space[size_ext + i - imin] = lda;
3031 imin = -i,imax = size_ext + i - 1;
3032 for(i = imin; i <= imax; i++){
3034 for(j = 0; j <= (lh_gold - 1); j++){
3035 ldb = working_space[j];
3037 if(k >= 0 && k < size_ext){
3038 ldc = working_space[2 * size_ext + k];
3039 lda = lda + ldb * ldc;
3043 working_space[4 * size_ext + i - imin] = lda;
3046 for(i = imin; i <= imax; i++)
3047 working_space[2 * size_ext + i - imin] = working_space[4 * size_ext + i - imin];
3049 for(i = 0; i < size_ext; i++)
3050 working_space[i] = 1;
3052 for(lindex = 0; lindex < deconIterations; lindex++){
3053 for(i = 0; i < size_ext; i++){
3054 if(
TMath::Abs(working_space[2 * size_ext + i]) > 0.00001 &&
TMath::Abs(working_space[i]) > 0.00001){
3062 if(jmax > (size_ext - 1 - i))
3065 for(j = jmin; j <= jmax; j++){
3066 ldb = working_space[j + lh_gold - 1 + size_ext];
3067 ldc = working_space[i + j];
3068 lda = lda + ldb * ldc;
3070 ldb = working_space[2 * size_ext + i];
3077 ldb = working_space[i];
3079 working_space[3 * size_ext + i] = lda;
3082 for(i = 0; i < size_ext; i++){
3083 working_space[i] = working_space[3 * size_ext + i];
3087 for(i=0;i<size_ext;i++){
3088 lda = working_space[i];
3091 working_space[size_ext + j] = lda;
3094 maximum = 0, maximum_decon = 0;
3096 for(i = 0; i < size_ext - j; i++){
3097 if(i >= shift && i < ssize + shift){
3098 working_space[i] = area * working_space[size_ext + i + j];
3099 if(maximum_decon < working_space[i])
3100 maximum_decon = working_space[i];
3101 if(maximum < working_space[6 * size_ext + i])
3102 maximum = working_space[6 * size_ext + i];
3106 working_space[i] = 0;
3114 for(i = 1; i < size_ext - 1; i++){
3115 if(working_space[i] > working_space[i - 1] && working_space[i] > working_space[i + 1]){
3116 if(i >= shift && i < ssize + shift){
3117 if(working_space[i] > lda*maximum_decon && working_space[6 * size_ext + i] > threshold * maximum / 100.0){
3118 for(j = i - 1, a = 0, b = 0; j <= i + 1; j++){
3119 a += (
Double_t)(j - shift) * working_space[j];
3120 b += working_space[j];
3128 if(peak_index == 0){
3134 for(j = 0, priz = 0; j < peak_index && priz == 0; j++){
3135 if(working_space[6 * size_ext + shift + (
Int_t)
a] > working_space[6 * size_ext + shift + (
Int_t)
fPositionX[j]])
3145 for(k = peak_index; k >= j; k--){
3160 for(i = 0; i < ssize; i++) destVector[i] = working_space[i + shift];
3161 delete [] working_space;
3164 Warning(
"SearchHighRes",
"Peak buffer full");
3174 bool backgroundRemove,
int deconIterations,
3175 bool markov,
int averWindow)
3179 return SearchHighRes(source,destVector,ssize,sigma,threshold,backgroundRemove,
3180 deconIterations,markov,averWindow);
3190 return s.
Search(hist,sigma,option,threshold);
virtual const char * GetName() const
Returns name of object.
Int_t Search1HighRes(Double_t *source, Double_t *destVector, Int_t ssize, Double_t sigma, Double_t threshold, bool backgroundRemove, Int_t deconIterations, bool markov, Int_t averWindow)
Old name of SearcHighRes introduced for back compatibility.
virtual Double_t GetBinCenter(Int_t bin) const
Return bin center for 1D histogram.
Int_t GetFirst() const
Return first bin on the axis i.e.
TString & ReplaceAll(const TString &s1, const TString &s2)
static TH1 * StaticBackground(const TH1 *hist, Int_t niter=20, Option_t *option="")
Static function, interface to TSpectrum::Background.
const char * DeconvolutionRL(Double_t *source, const Double_t *response, Int_t ssize, Int_t numberIterations, Int_t numberRepetitions, Double_t boost)
One-dimensional deconvolution function.
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
void ToLower()
Change string to lower-case.
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
virtual TObject * FindObject(const char *name) const
Find an object in this list using its name.
virtual Int_t GetDimension() const
The TNamed class is the base class for all named ROOT classes.
Double_t * fPositionX
[fNPeaks] X position of peaks
Int_t SearchHighRes(Double_t *source, Double_t *destVector, Int_t ssize, Double_t sigma, Double_t threshold, bool backgroundRemove, Int_t deconIterations, bool markov, Int_t averWindow)
One-dimensional high-resolution peak search function.
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
virtual void Delete(Option_t *option="")
Delete this object.
virtual Int_t Search(const TH1 *hist, Double_t sigma=2, Option_t *option="", Double_t threshold=0.05)
One-dimensional peak search function.
virtual ~TSpectrum()
Destructor.
Int_t GetLast() const
Return last bin on the axis i.e.
virtual TObject * Remove(TObject *obj)
Remove object from the list.
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...
TH1 * fHistogram
resulting histogram
Double_t * fPositionY
[fNPeaks] Y position of peaks
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
const char * Deconvolution(Double_t *source, const Double_t *response, Int_t ssize, Int_t numberIterations, Int_t numberRepetitions, Double_t boost)
One-dimensional deconvolution function.
static void SetAverageWindow(Int_t w=3)
Static function: Set average window of searched peaks (see TSpectrum::SearchHighRes).
static Int_t fgAverageWindow
Average window of searched peaks.
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
static void SetDeconIterations(Int_t n=3)
Static function: Set max number of decon iterations in deconvolution operation (see TSpectrum::Search...
virtual void SetMarkerSize(Size_t msize=1)
Set the marker size.
A PolyMarker is defined by an array on N points in a 2-D space.
virtual void Print(Option_t *option="") const
Print the array of positions.
virtual TH1 * Background(const TH1 *hist, Int_t niter=20, Option_t *option="")
One-dimensional background estimation function.
static Int_t fgIterations
Maximum number of decon iterations (default=3)
Int_t fMaxPeaks
Maximum number of peaks to be found.
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
void SetResolution(Double_t resolution=1)
NOT USED resolution: determines resolution of the neighbouring peaks default value is 1 correspond to...
Advanced Spectra Processing.
Double_t * fPosition
[fNPeaks] array of current peak positions
virtual void Add(TObject *obj)
#define dest(otri, vertexptr)
static Int_t StaticSearch(const TH1 *hist, Double_t sigma=2, Option_t *option="goff", Double_t threshold=0.05)
Static function, interface to TSpectrum::Search.
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Int_t fNPeaks
number of peaks found
Double_t Sqrt(Double_t x)
TList * GetListOfFunctions() const
const char * SmoothMarkov(Double_t *source, Int_t ssize, Int_t averWindow)
One-dimensional markov spectrum smoothing function.
const char * Unfolding(Double_t *source, const Double_t **respMatrix, Int_t ssizex, Int_t ssizey, Int_t numberIterations, Int_t numberRepetitions, Double_t boost)
One-dimensional unfolding function.
Double_t fResolution
NOT USED resolution of the neighboring peaks
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
const char * Data() const