40 #define PEAK_WINDOW 1024
67 :
TNamed(
"Spectrum",
"Miroslav Morhac peak finder")
151 if (h == 0)
return 0;
154 Error(
"Search",
"Only implemented for 1-d histograms");
181 Int_t size = last-first+1;
184 for (i = 0; i < size; i++) source[i] = h->
GetBinContent(i + first);
187 Background(source,size,numberIterations, direction, filterOrder,smoothing,
188 smoothWindow,compton);
193 char *
hbname =
new char[nch+20];
194 snprintf(hbname,nch+20,
"%s_background",h->
GetName());
197 hb->GetListOfFunctions()->
Delete();
199 for (i=0; i< size; i++) hb->SetBinContent(i+first,source[i]);
200 hb->SetEntries(size);
204 if (
gPad)
delete gPad->GetPrimitive(hbname);
265 if (hin == 0)
return 0;
268 Error(
"Search",
"Only implemented for 1-d and 2-d histograms");
271 if (threshold <=0 || threshold >= 1) {
272 Warning(
"Search",
"threshold must 0<threshold<1, threshold=0.05 assumed");
292 if (dimension == 1) {
295 Int_t size = last-first+1;
296 Int_t i, bin, npeaks;
299 for (i = 0; i < size; i++) source[i] = hin->
GetBinContent(i + first);
302 if (sigma < 1) sigma = 1;
303 if (sigma > 8) sigma = 8;
305 npeaks =
SearchHighRes(source, dest, size, sigma, 100*threshold,
308 for (i = 0; i < npeaks; i++) {
318 if (!npeaks)
return 0;
774 int numberIterations,
775 int direction,
int filterOrder,
776 bool smoothing,
int smoothWindow,
779 int i, j,
w, bw, b1, b2, priz;
780 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;
782 return "Wrong Parameters";
783 if (numberIterations < 1)
784 return "Width of Clipping Window Must Be Positive";
785 if (ssize < 2 * numberIterations + 1)
786 return "Too Large Clipping Window";
788 return "Incorrect width of smoothing window";
790 for (i = 0; i < ssize; i++){
791 working_space[i] = spectrum[i];
792 working_space[i + ssize] = spectrum[i];
794 bw=(smoothWindow-1)/2;
798 i = numberIterations;
801 for (j = i; j < ssize - i; j++) {
803 a = working_space[ssize + j];
804 b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
807 working_space[j] =
a;
810 else if (smoothing ==
kTRUE){
811 a = working_space[ssize + j];
814 for (w = j - bw; w <= j + bw; w++){
815 if ( w >= 0 && w < ssize){
816 av += working_space[ssize +
w];
823 for (w = j - i - bw; w <= j - i + bw; w++){
824 if ( w >= 0 && w < ssize){
825 b += working_space[ssize +
w];
832 for (w = j + i - bw; w <= j + i + bw; w++){
833 if ( w >= 0 && w < ssize){
834 c += working_space[ssize +
w];
845 for (j = i; j < ssize - i; j++)
846 working_space[ssize + j] = working_space[j];
856 for (j = i; j < ssize - i; j++) {
858 a = working_space[ssize + j];
859 b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
862 c -= working_space[ssize + j - (
Int_t) (2 * ai)] / 6;
863 c += 4 * working_space[ssize + j - (
Int_t) ai] / 6;
864 c += 4 * working_space[ssize + j + (
Int_t) ai] / 6;
865 c -= working_space[ssize + j + (
Int_t) (2 * ai)] / 6;
870 working_space[j] =
a;
873 else if (smoothing ==
kTRUE){
874 a = working_space[ssize + j];
877 for (w = j - bw; w <= j + bw; w++){
878 if ( w >= 0 && w < ssize){
879 av += working_space[ssize +
w];
886 for (w = j - i - bw; w <= j - i + bw; w++){
887 if ( w >= 0 && w < ssize){
888 b += working_space[ssize +
w];
895 for (w = j + i - bw; w <= j + i + bw; w++){
896 if ( w >= 0 && w < ssize){
897 c += working_space[ssize +
w];
905 for (w = j - (
Int_t)(2 * ai) - bw; w <= j - (
Int_t)(2 * ai) + bw; w++){
906 if (w >= 0 && w < ssize){
907 b4 += working_space[ssize +
w];
913 for (w = j - (
Int_t)ai - bw; w <= j - (
Int_t)ai + bw; w++){
914 if (w >= 0 && w < ssize){
915 c4 += working_space[ssize +
w];
921 for (w = j + (
Int_t)ai - bw; w <= j + (
Int_t)ai + bw; w++){
922 if (w >= 0 && w < ssize){
923 d4 += working_space[ssize +
w];
929 for (w = j + (
Int_t)(2 * ai) - bw; w <= j + (
Int_t)(2 * ai) + bw; w++){
930 if (w >= 0 && w < ssize){
931 e4 += working_space[ssize +
w];
936 b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6;
944 for (j = i; j < ssize - i; j++)
945 working_space[ssize + j] = working_space[j];
955 for (j = i; j < ssize - i; j++) {
957 a = working_space[ssize + j];
958 b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
961 c -= working_space[ssize + j - (
Int_t) (2 * ai)] / 6;
962 c += 4 * working_space[ssize + j - (
Int_t) ai] / 6;
963 c += 4 * working_space[ssize + j + (
Int_t) ai] / 6;
964 c -= working_space[ssize + j + (
Int_t) (2 * ai)] / 6;
967 d += working_space[ssize + j - (
Int_t) (3 * ai)] / 20;
968 d -= 6 * working_space[ssize + j - (
Int_t) (2 * ai)] / 20;
969 d += 15 * working_space[ssize + j - (
Int_t) ai] / 20;
970 d += 15 * working_space[ssize + j + (
Int_t) ai] / 20;
971 d -= 6 * working_space[ssize + j + (
Int_t) (2 * ai)] / 20;
972 d += working_space[ssize + j + (
Int_t) (3 * ai)] / 20;
979 working_space[j] =
a;
982 else if (smoothing ==
kTRUE){
983 a = working_space[ssize + j];
986 for (w = j - bw; w <= j + bw; w++){
987 if ( w >= 0 && w < ssize){
988 av += working_space[ssize +
w];
995 for (w = j - i - bw; w <= j - i + bw; w++){
996 if ( w >= 0 && w < ssize){
997 b += working_space[ssize +
w];
1004 for (w = j + i - bw; w <= j + i + bw; w++){
1005 if ( w >= 0 && w < ssize){
1006 c += working_space[ssize +
w];
1014 for (w = j - (
Int_t)(2 * ai) - bw; w <= j - (
Int_t)(2 * ai) + bw; w++){
1015 if (w >= 0 && w < ssize){
1016 b4 += working_space[ssize +
w];
1022 for (w = j - (
Int_t)ai - bw; w <= j - (
Int_t)ai + bw; w++){
1023 if (w >= 0 && w < ssize){
1024 c4 += working_space[ssize +
w];
1030 for (w = j + (
Int_t)ai - bw; w <= j + (
Int_t)ai + bw; w++){
1031 if (w >= 0 && w < ssize){
1032 d4 += working_space[ssize +
w];
1038 for (w = j + (
Int_t)(2 * ai) - bw; w <= j + (
Int_t)(2 * ai) + bw; w++){
1039 if (w >= 0 && w < ssize){
1040 e4 += working_space[ssize +
w];
1045 b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6;
1048 for (w = j - (
Int_t)(3 * ai) - bw; w <= j - (
Int_t)(3 * ai) + bw; w++){
1049 if (w >= 0 && w < ssize){
1050 b6 += working_space[ssize +
w];
1056 for (w = j - (
Int_t)(2 * ai) - bw; w <= j - (
Int_t)(2 * ai) + bw; w++){
1057 if (w >= 0 && w < ssize){
1058 c6 += working_space[ssize +
w];
1064 for (w = j - (
Int_t)ai - bw; w <= j - (
Int_t)ai + bw; w++){
1065 if (w >= 0 && w < ssize){
1066 d6 += working_space[ssize +
w];
1072 for (w = j + (
Int_t)ai - bw; w <= j + (
Int_t)ai + bw; w++){
1073 if (w >= 0 && w < ssize){
1074 e6 += working_space[ssize +
w];
1080 for (w = j + (
Int_t)(2 * ai) - bw; w <= j + (
Int_t)(2 * ai) + bw; w++){
1081 if (w >= 0 && w < ssize){
1082 f6 += working_space[ssize +
w];
1088 for (w = j + (
Int_t)(3 * ai) - bw; w <= j + (
Int_t)(3 * ai) + bw; w++){
1089 if (w >= 0 && w < ssize){
1090 g6 += working_space[ssize +
w];
1095 b6 = (b6 - 6 * c6 + 15 * d6 + 15 * e6 - 6 * f6 + g6) / 20;
1102 working_space[j]=av;
1105 for (j = i; j < ssize - i; j++)
1106 working_space[ssize + j] = working_space[j];
1116 for (j = i; j < ssize - i; j++) {
1117 if (smoothing ==
kFALSE){
1118 a = working_space[ssize + j];
1119 b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
1122 c -= working_space[ssize + j - (
Int_t) (2 * ai)] / 6;
1123 c += 4 * working_space[ssize + j - (
Int_t) ai] / 6;
1124 c += 4 * working_space[ssize + j + (
Int_t) ai] / 6;
1125 c -= working_space[ssize + j + (
Int_t) (2 * ai)] / 6;
1128 d += working_space[ssize + j - (
Int_t) (3 * ai)] / 20;
1129 d -= 6 * working_space[ssize + j - (
Int_t) (2 * ai)] / 20;
1130 d += 15 * working_space[ssize + j - (
Int_t) ai] / 20;
1131 d += 15 * working_space[ssize + j + (
Int_t) ai] / 20;
1132 d -= 6 * working_space[ssize + j + (
Int_t) (2 * ai)] / 20;
1133 d += working_space[ssize + j + (
Int_t) (3 * ai)] / 20;
1136 e -= working_space[ssize + j - (
Int_t) (4 * ai)] / 70;
1137 e += 8 * working_space[ssize + j - (
Int_t) (3 * ai)] / 70;
1138 e -= 28 * working_space[ssize + j - (
Int_t) (2 * ai)] / 70;
1139 e += 56 * working_space[ssize + j - (
Int_t) ai] / 70;
1140 e += 56 * working_space[ssize + j + (
Int_t) ai] / 70;
1141 e -= 28 * working_space[ssize + j + (
Int_t) (2 * ai)] / 70;
1142 e += 8 * working_space[ssize + j + (
Int_t) (3 * ai)] / 70;
1143 e -= working_space[ssize + j + (
Int_t) (4 * ai)] / 70;
1152 working_space[j] =
a;
1155 else if (smoothing ==
kTRUE){
1156 a = working_space[ssize + j];
1159 for (w = j - bw; w <= j + bw; w++){
1160 if ( w >= 0 && w < ssize){
1161 av += working_space[ssize +
w];
1168 for (w = j - i - bw; w <= j - i + bw; w++){
1169 if ( w >= 0 && w < ssize){
1170 b += working_space[ssize +
w];
1177 for (w = j + i - bw; w <= j + i + bw; w++){
1178 if ( w >= 0 && w < ssize){
1179 c += working_space[ssize +
w];
1187 for (w = j - (
Int_t)(2 * ai) - bw; w <= j - (
Int_t)(2 * ai) + bw; w++){
1188 if (w >= 0 && w < ssize){
1189 b4 += working_space[ssize +
w];
1195 for (w = j - (
Int_t)ai - bw; w <= j - (
Int_t)ai + bw; w++){
1196 if (w >= 0 && w < ssize){
1197 c4 += working_space[ssize +
w];
1203 for (w = j + (
Int_t)ai - bw; w <= j + (
Int_t)ai + bw; w++){
1204 if (w >= 0 && w < ssize){
1205 d4 += working_space[ssize +
w];
1211 for (w = j + (
Int_t)(2 * ai) - bw; w <= j + (
Int_t)(2 * ai) + bw; w++){
1212 if (w >= 0 && w < ssize){
1213 e4 += working_space[ssize +
w];
1218 b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6;
1221 for (w = j - (
Int_t)(3 * ai) - bw; w <= j - (
Int_t)(3 * ai) + bw; w++){
1222 if (w >= 0 && w < ssize){
1223 b6 += working_space[ssize +
w];
1229 for (w = j - (
Int_t)(2 * ai) - bw; w <= j - (
Int_t)(2 * ai) + bw; w++){
1230 if (w >= 0 && w < ssize){
1231 c6 += working_space[ssize +
w];
1237 for (w = j - (
Int_t)ai - bw; w <= j - (
Int_t)ai + bw; w++){
1238 if (w >= 0 && w < ssize){
1239 d6 += working_space[ssize +
w];
1245 for (w = j + (
Int_t)ai - bw; w <= j + (
Int_t)ai + bw; w++){
1246 if (w >= 0 && w < ssize){
1247 e6 += working_space[ssize +
w];
1253 for (w = j + (
Int_t)(2 * ai) - bw; w <= j + (
Int_t)(2 * ai) + bw; w++){
1254 if (w >= 0 && w < ssize){
1255 f6 += working_space[ssize +
w];
1261 for (w = j + (
Int_t)(3 * ai) - bw; w <= j + (
Int_t)(3 * ai) + bw; w++){
1262 if (w >= 0 && w < ssize){
1263 g6 += working_space[ssize +
w];
1268 b6 = (b6 - 6 * c6 + 15 * d6 + 15 * e6 - 6 * f6 + g6) / 20;
1271 for (w = j - (
Int_t)(4 * ai) - bw; w <= j - (
Int_t)(4 * ai) + bw; w++){
1272 if (w >= 0 && w < ssize){
1273 b8 += working_space[ssize +
w];
1279 for (w = j - (
Int_t)(3 * ai) - bw; w <= j - (
Int_t)(3 * ai) + bw; w++){
1280 if (w >= 0 && w < ssize){
1281 c8 += working_space[ssize +
w];
1287 for (w = j - (
Int_t)(2 * ai) - bw; w <= j - (
Int_t)(2 * ai) + bw; w++){
1288 if (w >= 0 && w < ssize){
1289 d8 += working_space[ssize +
w];
1295 for (w = j - (
Int_t)ai - bw; w <= j - (
Int_t)ai + bw; w++){
1296 if (w >= 0 && w < ssize){
1297 e8 += working_space[ssize +
w];
1303 for (w = j + (
Int_t)ai - bw; w <= j + (
Int_t)ai + bw; w++){
1304 if (w >= 0 && w < ssize){
1305 f8 += working_space[ssize +
w];
1311 for (w = j + (
Int_t)(2 * ai) - bw; w <= j + (
Int_t)(2 * ai) + bw; w++){
1312 if (w >= 0 && w < ssize){
1313 g8 += working_space[ssize +
w];
1319 for (w = j + (
Int_t)(3 * ai) - bw; w <= j + (
Int_t)(3 * ai) + bw; w++){
1320 if (w >= 0 && w < ssize){
1321 h8 += working_space[ssize +
w];
1327 for (w = j + (
Int_t)(4 * ai) - bw; w <= j + (
Int_t)(4 * ai) + bw; w++){
1328 if (w >= 0 && w < ssize){
1329 i8 += working_space[ssize +
w];
1334 b8 = ( -b8 + 8 * c8 - 28 * d8 + 56 * e8 - 56 * f8 - 28 * g8 + 8 * h8 - i8)/70;
1343 working_space[j]=av;
1346 for (j = i; j < ssize - i; j++)
1347 working_space[ssize + j] = working_space[j];
1355 if (compton ==
kTRUE) {
1356 for (i = 0, b2 = 0; i < ssize; i++){
1358 a = working_space[i], b = spectrum[i];
1364 yb1 = working_space[b1];
1365 for (b2 = b1 + 1, c = 0, priz = 0; priz == 0 && b2 < ssize; b2++){
1366 a = working_space[b2], b = spectrum[b2];
1375 yb2 = working_space[b2];
1377 for (j = b1, c = 0; j <= b2; j++){
1382 c = (yb2 - yb1) / c;
1383 for (j = b1, d = 0; j <= b2 && j < ssize; j++){
1387 working_space[ssize + j] =
a;
1393 for (j = b2, c = 0; j >= b1; j--){
1398 c = (yb1 - yb2) / c;
1399 for (j = b2, d = 0;j >= b1 && j >= 0; j--){
1403 working_space[ssize + j] =
a;
1412 for (j = 0; j < ssize; j++)
1413 spectrum[j] = working_space[ssize + j];
1414 delete[]working_space;
1488 Double_t nom, nip, nim, sp, sm, area = 0;
1490 return "Averaging Window must be positive";
1492 xmin = 0,xmax = ssize - 1;
1493 for(i = 0, maxch = 0; i < ssize; i++){
1495 if(maxch < source[i])
1501 delete [] working_space;
1506 working_space[
xmin] = 1;
1507 for(i = xmin; i <
xmax; i++){
1508 nip = source[i] / maxch;
1509 nim = source[i + 1] / maxch;
1511 for(l = 1; l <= averWindow; l++){
1513 a = source[
xmax] / maxch;
1516 a = source[i +
l] / maxch;
1526 if((i - l + 1) < xmin)
1527 a = source[
xmin] / maxch;
1530 a = source[i - l + 1] / maxch;
1541 a = working_space[i + 1] = working_space[i] *
a;
1544 for(i = xmin; i <=
xmax; i++){
1545 working_space[i] = working_space[i] / nom;
1547 for(i = 0; i < ssize; i++)
1548 source[i] = working_space[i] * area;
1549 delete[]working_space;
1854 int ssize,
int numberIterations,
1858 return "Wrong Parameters";
1860 if (numberRepetitions <= 0)
1861 return "Wrong Parameters";
1866 int i, j, k, lindex, posit, lh_gold,
l, repet;
1867 Double_t lda, ldb, ldc, area, maximum;
1874 for (i = 0; i < ssize; i++) {
1878 working_space[i] = lda;
1880 if (lda > maximum) {
1885 if (lh_gold == -1) {
1886 delete [] working_space;
1887 return "ZERO RESPONSE VECTOR";
1891 for (i = 0; i < ssize; i++)
1892 working_space[2 * ssize + i] = source[i];
1895 for (i = 0; i < ssize; i++){
1897 for (j = 0; j < ssize; j++){
1898 ldb = working_space[j];
1901 ldc = working_space[k];
1902 lda = lda + ldb * ldc;
1905 working_space[ssize + i] = lda;
1907 for (k = 0; k < ssize; k++){
1910 ldb = working_space[
l];
1911 ldc = working_space[2 * ssize + k];
1912 lda = lda + ldb * ldc;
1915 working_space[3 * ssize + i]=lda;
1919 for (i = 0; i < ssize; i++){
1920 working_space[2 * ssize + i] = working_space[3 * ssize + i];
1924 for (i = 0; i < ssize; i++)
1925 working_space[i] = 1;
1928 for (repet = 0; repet < numberRepetitions; repet++) {
1930 for (i = 0; i < ssize; i++)
1931 working_space[i] =
TMath::Power(working_space[i], boost);
1933 for (lindex = 0; lindex < numberIterations; lindex++) {
1934 for (i = 0; i < ssize; i++) {
1935 if (working_space[2 * ssize + i] > 0.000001
1936 && working_space[i] > 0.000001) {
1938 for (j = 0; j < lh_gold; j++) {
1939 ldb = working_space[j + ssize];
1944 ldc = working_space[k];
1947 ldc += working_space[k];
1951 ldc = working_space[i];
1952 lda = lda + ldb * ldc;
1954 ldb = working_space[2 * ssize + i];
1960 ldb = working_space[i];
1962 working_space[3 * ssize + i] = lda;
1965 for (i = 0; i < ssize; i++)
1966 working_space[i] = working_space[3 * ssize + i];
1971 for (i = 0; i < ssize; i++) {
1972 lda = working_space[i];
1975 working_space[ssize + j] = lda;
1979 for (i = 0; i < ssize; i++)
1980 source[i] = area * working_space[ssize + i];
1981 delete[]working_space;
2134 int ssize,
int numberIterations,
2138 return "Wrong Parameters";
2140 if (numberRepetitions <= 0)
2141 return "Wrong Parameters";
2146 int i, j, k, lindex, posit, lh_gold, repet, kmin, kmax;
2153 for (i = 0; i < ssize; i++) {
2157 working_space[ssize + i] = lda;
2158 if (lda > maximum) {
2163 if (lh_gold == -1) {
2164 delete [] working_space;
2165 return "ZERO RESPONSE VECTOR";
2169 for (i = 0; i < ssize; i++)
2170 working_space[2 * ssize + i] = source[i];
2173 for (i = 0; i < ssize; i++){
2174 if (i <= ssize - lh_gold)
2175 working_space[i] = 1;
2178 working_space[i] = 0;
2182 for (repet = 0; repet < numberRepetitions; repet++) {
2184 for (i = 0; i < ssize; i++)
2185 working_space[i] =
TMath::Power(working_space[i], boost);
2187 for (lindex = 0; lindex < numberIterations; lindex++) {
2188 for (i = 0; i <= ssize - lh_gold; i++){
2190 if (working_space[i] > 0){
2191 for (j = i; j < i + lh_gold; j++){
2192 ldb = working_space[2 * ssize + j];
2196 if (kmax > lh_gold - 1)
2198 kmin = j + lh_gold - ssize;
2202 for (k = kmax; k >= kmin; k--){
2203 ldc += working_space[ssize + k] * working_space[j - k];
2211 ldb = ldb * working_space[ssize + j - i];
2215 lda = lda * working_space[i];
2217 working_space[3 * ssize + i] = lda;
2219 for (i = 0; i < ssize; i++)
2220 working_space[i] = working_space[3 * ssize + i];
2225 for (i = 0; i < ssize; i++) {
2226 lda = working_space[i];
2229 working_space[ssize + j] = lda;
2233 for (i = 0; i < ssize; i++)
2234 source[i] = working_space[ssize + i];
2235 delete[]working_space;
2360 int ssizex,
int ssizey,
2361 int numberIterations,
2364 int i, j, k, lindex, lhx = 0, repet;
2366 if (ssizex <= 0 || ssizey <= 0)
2367 return "Wrong Parameters";
2368 if (ssizex < ssizey)
2369 return "Sizex must be greater than sizey)";
2370 if (numberIterations <= 0)
2371 return "Number of iterations must be positive";
2373 new Double_t[ssizex * ssizey + 2 * ssizey * ssizey + 4 * ssizex];
2376 for (j = 0; j < ssizey && lhx != -1; j++) {
2379 for (i = 0; i < ssizex; i++) {
2380 lda = respMatrix[j][i];
2384 working_space[j * ssizex + i] = lda;
2388 for (i = 0; i < ssizex; i++)
2389 working_space[j * ssizex + i] /= area;
2393 delete [] working_space;
2394 return (
"ZERO COLUMN IN RESPONSE MATRIX");
2398 for (i = 0; i < ssizex; i++)
2399 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i] =
2403 for (i = 0; i < ssizey; i++) {
2404 for (j = 0; j < ssizey; j++) {
2406 for (k = 0; k < ssizex; k++) {
2407 ldb = working_space[ssizex * i + k];
2408 ldc = working_space[ssizex * j + k];
2409 lda = lda + ldb * ldc;
2411 working_space[ssizex * ssizey + ssizey * i + j] = lda;
2414 for (k = 0; k < ssizex; k++) {
2415 ldb = working_space[ssizex * i + k];
2417 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex +
2419 lda = lda + ldb * ldc;
2421 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i] =
2426 for (i = 0; i < ssizey; i++)
2427 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i] =
2428 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i];
2431 for (i = 0; i < ssizey; i++) {
2432 for (j = 0; j < ssizey; j++) {
2434 for (k = 0; k < ssizey; k++) {
2435 ldb = working_space[ssizex * ssizey + ssizey * i + k];
2436 ldc = working_space[ssizex * ssizey + ssizey * j + k];
2437 lda = lda + ldb * ldc;
2439 working_space[ssizex * ssizey + ssizey * ssizey + ssizey * i + j] =
2443 for (k = 0; k < ssizey; k++) {
2444 ldb = working_space[ssizex * ssizey + ssizey * i + k];
2446 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex +
2448 lda = lda + ldb * ldc;
2450 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i] =
2455 for (i = 0; i < ssizey; i++)
2456 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i] =
2457 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 + i] = 1;
2464 for (repet = 0; repet < numberRepetitions; repet++) {
2466 for (i = 0; i < ssizey; i++)
2467 working_space[ssizex * ssizey + 2 * ssizey * ssizey + i] =
TMath::Power(working_space[ssizex * ssizey + 2 * ssizey * ssizey + i], boost);
2469 for (lindex = 0; lindex < numberIterations; lindex++) {
2470 for (i = 0; i < ssizey; i++) {
2472 for (j = 0; j < ssizey; j++) {
2474 working_space[ssizex * ssizey + ssizey * ssizey + ssizey * i + j];
2475 ldc = working_space[ssizex * ssizey + 2 * ssizey * ssizey + j];
2476 lda = lda + ldb * ldc;
2479 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i];
2486 ldb = working_space[ssizex * ssizey + 2 * ssizey * ssizey + i];
2488 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i] = lda;
2490 for (i = 0; i < ssizey; i++)
2491 working_space[ssizex * ssizey + 2 * ssizey * ssizey + i] =
2492 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i];
2497 for (i = 0; i < ssizex; i++) {
2499 source[i] = working_space[ssizex * ssizey + 2 * ssizey * ssizey + i];
2504 delete[]working_space;
2732 bool backgroundRemove,
int deconIterations,
2733 bool markov,
int averWindow)
2735 int i, j, numberIterations = (
Int_t)(7 * sigma + 0.5);
2737 int k, lindex, posit, imin, imax, jmin, jmax, lh_gold, priz;
2738 Double_t lda, ldb, ldc, area, maximum, maximum_decon;
2739 int xmin,
xmax,
l, peak_index = 0, size_ext = ssize + 2 * numberIterations, shift = numberIterations, bw = 2,
w;
2741 Double_t nom, nip, nim, sp, sm, plocha = 0;
2742 Double_t m0low=0,m1low=0,m2low=0,l0low=0,l1low=0,detlow,av,men;
2744 Error(
"SearchHighRes",
"Invalid sigma, must be greater than or equal to 1");
2748 if(threshold<=0 || threshold>=100){
2749 Error(
"SearchHighRes",
"Invalid threshold, must be positive and less than 100");
2753 j = (
Int_t) (5.0 * sigma + 0.5);
2755 Error(
"SearchHighRes",
"Too large sigma");
2759 if (markov ==
true) {
2760 if (averWindow <= 0) {
2761 Error(
"SearchHighRes",
"Averaging window must be positive");
2766 if(backgroundRemove ==
true){
2767 if(ssize < 2 * numberIterations + 1){
2768 Error(
"SearchHighRes",
"Too large clipping window");
2773 k = int(2 * sigma+0.5);
2775 for(i = 0;i < k;i++){
2776 a = i,b = source[i];
2777 m0low += 1,m1low +=
a,m2low += a *
a,l0low += b,l1low += a * b;
2779 detlow = m0low * m2low - m1low * m1low;
2781 l1low = (-l0low * m1low + l1low * m0low) / detlow;
2793 i = (
Int_t)(7 * sigma + 0.5);
2796 for (j=0;j<7 * (ssize + i);j++) working_space[j] = 0;
2797 for(i = 0; i < size_ext; i++){
2800 working_space[i + size_ext] = source[0] + l1low *
a;
2801 if(working_space[i + size_ext] < 0)
2802 working_space[i + size_ext]=0;
2805 else if(i >= ssize + shift){
2806 a = i - (ssize - 1 + shift);
2807 working_space[i + size_ext] = source[ssize - 1];
2808 if(working_space[i + size_ext] < 0)
2809 working_space[i + size_ext]=0;
2813 working_space[i + size_ext] = source[i - shift];
2816 if(backgroundRemove ==
true){
2817 for(i = 1; i <= numberIterations; i++){
2818 for(j = i; j < size_ext - i; j++){
2819 if(markov ==
false){
2820 a = working_space[size_ext + j];
2821 b = (working_space[size_ext + j - i] + working_space[size_ext + j + i]) / 2.0;
2829 a = working_space[size_ext + j];
2832 for (
w = j - bw;
w <= j + bw;
w++){
2833 if (
w >= 0 &&
w < size_ext){
2834 av += working_space[size_ext +
w];
2841 for (
w = j - i - bw;
w <= j - i + bw;
w++){
2842 if (
w >= 0 &&
w < size_ext){
2843 b += working_space[size_ext +
w];
2850 for (
w = j + i - bw;
w <= j + i + bw;
w++){
2851 if (
w >= 0 &&
w < size_ext){
2852 c += working_space[size_ext +
w];
2860 working_space[j]=av;
2863 for(j = i; j < size_ext - i; j++)
2864 working_space[size_ext + j] = working_space[j];
2866 for(j = 0;j < size_ext; j++){
2869 b = source[0] + l1low *
a;
2871 working_space[size_ext + j] = b - working_space[size_ext + j];
2874 else if(j >= ssize + shift){
2875 a = j - (ssize - 1 + shift);
2876 b = source[ssize - 1];
2878 working_space[size_ext + j] = b - working_space[size_ext + j];
2882 working_space[size_ext + j] = source[j - shift] - working_space[size_ext + j];
2885 for(j = 0;j < size_ext; j++){
2886 if(working_space[size_ext + j] < 0) working_space[size_ext + j] = 0;
2890 for(i = 0; i < size_ext; i++){
2891 working_space[i + 6*size_ext] = working_space[i + size_ext];
2895 for(j = 0; j < size_ext; j++)
2896 working_space[2 * size_ext + j] = working_space[size_ext + j];
2897 xmin = 0,xmax = size_ext - 1;
2898 for(i = 0, maxch = 0; i < size_ext; i++){
2899 working_space[i] = 0;
2900 if(maxch < working_space[2 * size_ext + i])
2901 maxch = working_space[2 * size_ext + i];
2902 plocha += working_space[2 * size_ext + i];
2905 delete [] working_space;
2910 working_space[
xmin] = 1;
2911 for(i = xmin; i <
xmax; i++){
2912 nip = working_space[2 * size_ext + i] / maxch;
2913 nim = working_space[2 * size_ext + i + 1] / maxch;
2915 for(l = 1; l <= averWindow; l++){
2917 a = working_space[2 * size_ext + xmax] / maxch;
2920 a = working_space[2 * size_ext + i +
l] / maxch;
2932 if((i - l + 1) < xmin)
2933 a = working_space[2 * size_ext +
xmin] / maxch;
2936 a = working_space[2 * size_ext + i - l + 1] / maxch;
2950 a = working_space[i + 1] = working_space[i] *
a;
2953 for(i = xmin; i <=
xmax; i++){
2954 working_space[i] = working_space[i] / nom;
2956 for(j = 0; j < size_ext; j++)
2957 working_space[size_ext + j] = working_space[j] * plocha;
2958 for(j = 0; j < size_ext; j++){
2959 working_space[2 * size_ext + j] = working_space[size_ext + j];
2961 if(backgroundRemove ==
true){
2962 for(i = 1; i <= numberIterations; i++){
2963 for(j = i; j < size_ext - i; j++){
2964 a = working_space[size_ext + j];
2965 b = (working_space[size_ext + j - i] + working_space[size_ext + j + i]) / 2.0;
2968 working_space[j] =
a;
2970 for(j = i; j < size_ext - i; j++)
2971 working_space[size_ext + j] = working_space[j];
2973 for(j = 0; j < size_ext; j++){
2974 working_space[size_ext + j] = working_space[2 * size_ext + j] - working_space[size_ext + j];
2984 for(i = 0; i < size_ext; i++){
2986 lda = lda * lda / (2 * sigma *
sigma);
2992 working_space[i] = lda;
3000 for(i = 0; i < size_ext; i++)
3001 working_space[2 * size_ext + i] =
TMath::Abs(working_space[size_ext + i]);
3008 for(i = imin; i <= imax; i++){
3013 jmax = lh_gold - 1 - i;
3014 if(jmax > (lh_gold - 1))
3017 for(j = jmin;j <= jmax; j++){
3018 ldb = working_space[j];
3019 ldc = working_space[i + j];
3020 lda = lda + ldb * ldc;
3022 working_space[size_ext + i - imin] = lda;
3026 imin = -i,imax = size_ext + i - 1;
3027 for(i = imin; i <= imax; i++){
3029 for(j = 0; j <= (lh_gold - 1); j++){
3030 ldb = working_space[j];
3032 if(k >= 0 && k < size_ext){
3033 ldc = working_space[2 * size_ext + k];
3034 lda = lda + ldb * ldc;
3038 working_space[4 * size_ext + i - imin] = lda;
3041 for(i = imin; i <= imax; i++)
3042 working_space[2 * size_ext + i - imin] = working_space[4 * size_ext + i - imin];
3044 for(i = 0; i < size_ext; i++)
3045 working_space[i] = 1;
3047 for(lindex = 0; lindex < deconIterations; lindex++){
3048 for(i = 0; i < size_ext; i++){
3049 if(
TMath::Abs(working_space[2 * size_ext + i]) > 0.00001 &&
TMath::Abs(working_space[i]) > 0.00001){
3057 if(jmax > (size_ext - 1 - i))
3060 for(j = jmin; j <= jmax; j++){
3061 ldb = working_space[j + lh_gold - 1 + size_ext];
3062 ldc = working_space[i + j];
3063 lda = lda + ldb * ldc;
3065 ldb = working_space[2 * size_ext + i];
3072 ldb = working_space[i];
3074 working_space[3 * size_ext + i] = lda;
3077 for(i = 0; i < size_ext; i++){
3078 working_space[i] = working_space[3 * size_ext + i];
3082 for(i=0;i<size_ext;i++){
3083 lda = working_space[i];
3086 working_space[size_ext + j] = lda;
3089 maximum = 0, maximum_decon = 0;
3091 for(i = 0; i < size_ext - j; i++){
3092 if(i >= shift && i < ssize + shift){
3093 working_space[i] = area * working_space[size_ext + i + j];
3094 if(maximum_decon < working_space[i])
3095 maximum_decon = working_space[i];
3096 if(maximum < working_space[6 * size_ext + i])
3097 maximum = working_space[6 * size_ext + i];
3101 working_space[i] = 0;
3109 for(i = 1; i < size_ext - 1; i++){
3110 if(working_space[i] > working_space[i - 1] && working_space[i] > working_space[i + 1]){
3111 if(i >= shift && i < ssize + shift){
3112 if(working_space[i] > lda*maximum_decon && working_space[6 * size_ext + i] > threshold * maximum / 100.0){
3113 for(j = i - 1, a = 0, b = 0; j <= i + 1; j++){
3114 a += (
Double_t)(j - shift) * working_space[j];
3115 b += working_space[j];
3123 if(peak_index == 0){
3129 for(j = 0, priz = 0; j < peak_index && priz == 0; j++){
3130 if(working_space[6 * size_ext + shift + (
Int_t)
a] > working_space[6 * size_ext + shift + (
Int_t)
fPositionX[j]])
3140 for(k = peak_index; k >= j; k--){
3155 for(i = 0; i < ssize; i++) destVector[i] = working_space[i + shift];
3156 delete [] working_space;
3159 Warning(
"SearchHighRes",
"Peak buffer full");
3169 bool backgroundRemove,
int deconIterations,
3170 bool markov,
int averWindow)
3174 return SearchHighRes(source,destVector,ssize,sigma,threshold,backgroundRemove,
3175 deconIterations,markov,averWindow);
3185 return s.
Search(hist,sigma,option,threshold);
Int_t GetFirst() const
Return first bin on the axis i.e.
TList * GetListOfFunctions() const
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 GetBinContent(Int_t bin) const
Return content of bin number bin.
TString & ReplaceAll(const TString &s1, const TString &s2)
virtual Int_t GetDimension() const
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.
ClassImp(TSpectrum) TSpectrum
Constructor.
void ToLower()
Change string to lower-case.
virtual TObject * FindObject(const char *name) const
Find an object in this list using its name.
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
const char * Data() 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)
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
virtual void Delete(Option_t *option="")
Delete this object.
virtual Double_t GetBinCenter(Int_t bin) const
return bin center for 1D historam Better to use h1.GetXaxis().GetBinCenter(bin)
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.
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
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 const char * GetName() const
Returns name of object.
virtual void SetMarkerStyle(Style_t mstyle=1)
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)
A PolyMarker is defined by an array on N points in a 2-D space.
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)
ClassImp(TMCParticle) void TMCParticle printf(": p=(%7.3f,%7.3f,%9.3f) ;", fPx, fPy, fPz)
Int_t fMaxPeaks
Maximum number of peaks to be found.
void SetResolution(Double_t resolution=1)
resolution: determines resolution of the neighbouring peaks default value is 1 correspond to 3 sigma ...
Advanced Spectra Processing.
Int_t GetLast() const
Return last bin on the axis i.e.
Double_t * fPosition
[fNPeaks] array of current peak positions
virtual void Print(Option_t *option="") const
Print the array of positions.
virtual void Add(TObject *obj)
#define dest(otri, vertexptr)
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
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.
Int_t fNPeaks
number of peaks found
Double_t Sqrt(Double_t x)
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
resolution of the neighboring peaks
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.