NCFS-Pack
A generic (astro)particle physics analysis framework
 
Loading...
Searching...
No Matches
NcSpectrum.cxx
Go to the documentation of this file.
1
2
11
13
58
59#include "NcSpectrum.h"
60
61// Added by Nick van Eijndhoven, VUB-IIHE, Brussels, May 7, 2025 07:01Z
62#include "Riostream.h"
63
64// Commented out by Nick van Eijndhoven, VUB-IIHE, Brussels, May 7, 2025 07:01Z
65// since these statements are included in the file NcSpectrum.h
66//@@@ #include "TPolyMarker.h"
67//@@@ #include "TVirtualPad.h"
68//@@@ #include "TList.h"
69//@@@ #include "TH1.h"
70//@@@ #include "TMath.h"
71
72Int_t NcSpectrum::fgIterations = 3;
73Int_t NcSpectrum::fgAverageWindow = 3;
74
75#define PEAK_WINDOW 1024
76
77// Modified by Nick van Eijndhoven, VUB-IIHE, Brussels, May 7, 2025 07:01Z
78//@@@ ClassImp(TSpectrum)
79ClassImp(NcSpectrum); // Class implementation to enable ROOT I/O
80
81
82//______________________________________________________________________________
83NcSpectrum::NcSpectrum() :TNamed("Spectrum", "Miroslav Morhac peak finder")
84{
90
91 Int_t n = 100;
92 fMaxPeaks = n;
93 fPosition = new Float_t[n];
94 fPositionX = new Float_t[n];
95 fPositionY = new Float_t[n];
96 fResolution = 1;
97 fHistogram = 0;
98 fNPeaks = 0;
99}
100
101
102//______________________________________________________________________________
103NcSpectrum::NcSpectrum(Int_t maxpositions, Float_t resolution)
104 :TNamed("Spectrum", "Miroslav Morhac peak finder")
105{
116
117 Int_t n = maxpositions;
118 if (n <= 0) n = 1;
119 fMaxPeaks = n;
120 fPosition = new Float_t[n];
121 fPositionX = new Float_t[n];
122 fPositionY = new Float_t[n];
123 fHistogram = 0;
124 fNPeaks = 0;
125 SetResolution(resolution);
126}
127
128
129//______________________________________________________________________________
130NcSpectrum::~NcSpectrum()
131{
137
138 delete [] fPosition;
139 delete [] fPositionX;
140 delete [] fPositionY;
141 delete fHistogram;
142}
143
144
145//______________________________________________________________________________
146void NcSpectrum::SetAverageWindow(Int_t w)
147{
154
155 fgAverageWindow = w;
156}
157
158
159//______________________________________________________________________________
160void NcSpectrum::SetDeconIterations(Int_t n)
161{
168
169 fgIterations = n;
170}
171
172
173//______________________________________________________________________________
174TH1* NcSpectrum::Background(const TH1 * h, int numberIterations, Option_t * option)
175{
215
216 if (h == 0) return 0;
217 Int_t dimension = h->GetDimension();
218 if (dimension > 1) {
219 Error("Search", "Only implemented for 1-d histograms");
220 return 0;
221 }
222 TString opt = option;
223 opt.ToLower();
224
225 //set options
226 Int_t direction = kBackDecreasingWindow;
227 if (opt.Contains("backincreasingwindow")) direction = kBackIncreasingWindow;
228 Int_t filterOrder = kBackOrder2;
229 if (opt.Contains("backorder4")) filterOrder = kBackOrder4;
230 if (opt.Contains("backorder6")) filterOrder = kBackOrder6;
231 if (opt.Contains("backorder8")) filterOrder = kBackOrder8;
232 Bool_t smoothing = kTRUE;
233 if (opt.Contains("nosmoothing")) smoothing = kFALSE;
234 Int_t smoothWindow = kBackSmoothing3;
235 if (opt.Contains("backsmoothing5")) smoothWindow = kBackSmoothing5;
236 if (opt.Contains("backsmoothing7")) smoothWindow = kBackSmoothing7;
237 if (opt.Contains("backsmoothing9")) smoothWindow = kBackSmoothing9;
238 if (opt.Contains("backsmoothing11")) smoothWindow = kBackSmoothing11;
239 if (opt.Contains("backsmoothing13")) smoothWindow = kBackSmoothing13;
240 if (opt.Contains("backsmoothing15")) smoothWindow = kBackSmoothing15;
241 Bool_t compton = kFALSE;
242 if (opt.Contains("compton")) compton = kTRUE;
243
244 Int_t first = h->GetXaxis()->GetFirst();
245 Int_t last = h->GetXaxis()->GetLast();
246 Int_t size = last-first+1;
247 Int_t i;
248 Float_t * source = new float[size];
249 for (i = 0; i < size; i++) source[i] = h->GetBinContent(i + first);
250
251 //find background (source is input and in output contains the background
252 Background(source,size,numberIterations, direction, filterOrder,smoothing,
253 smoothWindow,compton);
254
255 //create output histogram containing backgound
256 //only bins in the range of the input histogram are filled
257 Int_t nch = strlen(h->GetName());
258 char *hbname = new char[nch+20];
259 snprintf(hbname,nch+20,"%s_background",h->GetName());
260 TH1 *hb = (TH1*)h->Clone(hbname);
261 hb->Reset();
262 hb->GetListOfFunctions()->Delete();
263 hb->SetLineColor(2);
264 for (i=0; i< size; i++) hb->SetBinContent(i+first,source[i]);
265 hb->SetEntries(size);
266
267 //if option "same is specified, draw the result in the pad
268 if (opt.Contains("same")) {
269 if (gPad) delete gPad->GetPrimitive(hbname);
270 hb->Draw("same");
271 }
272 delete [] source;
273 delete [] hbname;
274 return hb;
275}
276
277
278//______________________________________________________________________________
279void NcSpectrum::Print(Option_t *) const
280{
286
287 printf("\nNumber of positions = %d\n",fNPeaks);
288 for (Int_t i=0;i<fNPeaks;i++) {
289 printf(" x[%d] = %g, y[%d] = %g\n",i,fPositionX[i],i,fPositionY[i]);
290 }
291}
292
293
294//______________________________________________________________________________
295Int_t NcSpectrum::Search(const TH1 * hin, Double_t sigma, Option_t * option,
296 Double_t threshold)
297{
337
338 if (hin == 0) return 0;
339 Int_t dimension = hin->GetDimension();
340 if (dimension > 2) {
341 Error("Search", "Only implemented for 1-d and 2-d histograms");
342 return 0;
343 }
344 if (threshold <=0 || threshold >= 1) {
345 Warning("Search","threshold must 0<threshold<1, threshol=0.05 assumed");
346 threshold = 0.05;
347 }
348 TString opt = option;
349 opt.ToLower();
350 Bool_t background = kTRUE;
351 if (opt.Contains("nobackground")) {
352 background = kFALSE;
353 opt.ReplaceAll("nobackground","");
354 }
355 Bool_t markov = kTRUE;
356 if (opt.Contains("nomarkov")) {
357 markov = kFALSE;
358 opt.ReplaceAll("nomarkov","");
359 }
360 Bool_t draw = kTRUE;
361 if (opt.Contains("nodraw")) {
362 draw = kFALSE;
363 opt.ReplaceAll("nodraw","");
364 }
365 if (dimension == 1) {
366 Int_t first = hin->GetXaxis()->GetFirst();
367 Int_t last = hin->GetXaxis()->GetLast();
368 Int_t size = last-first+1;
369 Int_t i, bin, npeaks;
370 Float_t * source = new float[size];
371 Float_t * dest = new float[size];
372 for (i = 0; i < size; i++) source[i] = hin->GetBinContent(i + first);
373 if (sigma < 1) {
374 sigma = size/fMaxPeaks;
375 if (sigma < 1) sigma = 1;
376 if (sigma > 8) sigma = 8;
377 }
378 npeaks = SearchHighRes(source, dest, size, sigma, 100*threshold,
379 background, fgIterations, markov, fgAverageWindow);
380
381 for (i = 0; i < npeaks; i++) {
382 bin = first + Int_t(fPositionX[i] + 0.5);
383 fPositionX[i] = hin->GetBinCenter(bin);
384 fPositionY[i] = hin->GetBinContent(bin);
385 }
386 delete [] source;
387 delete [] dest;
388
389 if (opt.Contains("goff"))
390 return npeaks;
391 if (!npeaks) return 0;
392 TPolyMarker * pm =
393 (TPolyMarker*)hin->GetListOfFunctions()->FindObject("TPolyMarker");
394 if (pm) {
395 hin->GetListOfFunctions()->Remove(pm);
396 delete pm;
397 }
398 pm = new TPolyMarker(npeaks, fPositionX, fPositionY);
399 hin->GetListOfFunctions()->Add(pm);
400 pm->SetMarkerStyle(23);
401 pm->SetMarkerColor(kRed);
402 pm->SetMarkerSize(1.3);
403 opt.ReplaceAll(" ","");
404 opt.ReplaceAll(",","");
405 if (draw)
406 ((TH1*)hin)->Draw(opt.Data());
407 return npeaks;
408 }
409 return 0;
410}
411
412
413//______________________________________________________________________________
414void NcSpectrum::SetResolution(Float_t resolution)
415{
425
426 if (resolution > 1)
427 fResolution = resolution;
428 else
429 fResolution = 1;
430}
431
432
433//______________________________________________________________________________
434const char* NcSpectrum::Background(float *spectrum, int ssize,
435 int numberIterations,
436 int direction, int filterOrder,
437 bool smoothing,int smoothWindow,
438 bool compton)
439{
874
875 int i, j, w, bw, b1, b2, priz;
876 float 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;
877 if (ssize <= 0)
878 return "Wrong Parameters";
879 if (numberIterations < 1)
880 return "Width of Clipping Window Must Be Positive";
881 if (ssize < 2 * numberIterations + 1)
882 return "Too Large Clipping Window";
883 if (smoothing == kTRUE && smoothWindow != kBackSmoothing3 && smoothWindow != kBackSmoothing5 && smoothWindow != kBackSmoothing7 && smoothWindow != kBackSmoothing9 && smoothWindow != kBackSmoothing11 && smoothWindow != kBackSmoothing13 && smoothWindow != kBackSmoothing15)
884 return "Incorrect width of smoothing window";
885 float *working_space = new float[2 * ssize];
886 for (i = 0; i < ssize; i++){
887 working_space[i] = spectrum[i];
888 working_space[i + ssize] = spectrum[i];
889 }
890 bw=(smoothWindow-1)/2;
891 if (direction == kBackIncreasingWindow)
892 i = 1;
893 else if(direction == kBackDecreasingWindow)
894 i = numberIterations;
895 if (filterOrder == kBackOrder2) {
896 do{
897 for (j = i; j < ssize - i; j++) {
898 if (smoothing == kFALSE){
899 a = working_space[ssize + j];
900 b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
901 if (b < a)
902 a = b;
903 working_space[j] = a;
904 }
905
906 else if (smoothing == kTRUE){
907 a = working_space[ssize + j];
908 av = 0;
909 men = 0;
910 for (w = j - bw; w <= j + bw; w++){
911 if ( w >= 0 && w < ssize){
912 av += working_space[ssize + w];
913 men +=1;
914 }
915 }
916 av = av / men;
917 b = 0;
918 men = 0;
919 for (w = j - i - bw; w <= j - i + bw; w++){
920 if ( w >= 0 && w < ssize){
921 b += working_space[ssize + w];
922 men +=1;
923 }
924 }
925 b = b / men;
926 c = 0;
927 men = 0;
928 for (w = j + i - bw; w <= j + i + bw; w++){
929 if ( w >= 0 && w < ssize){
930 c += working_space[ssize + w];
931 men +=1;
932 }
933 }
934 c = c / men;
935 b = (b + c) / 2;
936 if (b < a)
937 av = b;
938 working_space[j]=av;
939 }
940 }
941 for (j = i; j < ssize - i; j++)
942 working_space[ssize + j] = working_space[j];
943 if (direction == kBackIncreasingWindow)
944 i+=1;
945 else if(direction == kBackDecreasingWindow)
946 i-=1;
947 }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1));
948 }
949
950 else if (filterOrder == kBackOrder4) {
951 do{
952 for (j = i; j < ssize - i; j++) {
953 if (smoothing == kFALSE){
954 a = working_space[ssize + j];
955 b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
956 c = 0;
957 ai = i / 2;
958 c -= working_space[ssize + j - (int) (2 * ai)] / 6;
959 c += 4 * working_space[ssize + j - (int) ai] / 6;
960 c += 4 * working_space[ssize + j + (int) ai] / 6;
961 c -= working_space[ssize + j + (int) (2 * ai)] / 6;
962 if (b < c)
963 b = c;
964 if (b < a)
965 a = b;
966 working_space[j] = a;
967 }
968
969 else if (smoothing == kTRUE){
970 a = working_space[ssize + j];
971 av = 0;
972 men = 0;
973 for (w = j - bw; w <= j + bw; w++){
974 if ( w >= 0 && w < ssize){
975 av += working_space[ssize + w];
976 men +=1;
977 }
978 }
979 av = av / men;
980 b = 0;
981 men = 0;
982 for (w = j - i - bw; w <= j - i + bw; w++){
983 if ( w >= 0 && w < ssize){
984 b += working_space[ssize + w];
985 men +=1;
986 }
987 }
988 b = b / men;
989 c = 0;
990 men = 0;
991 for (w = j + i - bw; w <= j + i + bw; w++){
992 if ( w >= 0 && w < ssize){
993 c += working_space[ssize + w];
994 men +=1;
995 }
996 }
997 c = c / men;
998 b = (b + c) / 2;
999 ai = i / 2;
1000 b4 = 0, men = 0;
1001 for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){
1002 if (w >= 0 && w < ssize){
1003 b4 += working_space[ssize + w];
1004 men +=1;
1005 }
1006 }
1007 b4 = b4 / men;
1008 c4 = 0, men = 0;
1009 for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){
1010 if (w >= 0 && w < ssize){
1011 c4 += working_space[ssize + w];
1012 men +=1;
1013 }
1014 }
1015 c4 = c4 / men;
1016 d4 = 0, men = 0;
1017 for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){
1018 if (w >= 0 && w < ssize){
1019 d4 += working_space[ssize + w];
1020 men +=1;
1021 }
1022 }
1023 d4 = d4 / men;
1024 e4 = 0, men = 0;
1025 for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){
1026 if (w >= 0 && w < ssize){
1027 e4 += working_space[ssize + w];
1028 men +=1;
1029 }
1030 }
1031 e4 = e4 / men;
1032 b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6;
1033 if (b < b4)
1034 b = b4;
1035 if (b < a)
1036 av = b;
1037 working_space[j]=av;
1038 }
1039 }
1040 for (j = i; j < ssize - i; j++)
1041 working_space[ssize + j] = working_space[j];
1042 if (direction == kBackIncreasingWindow)
1043 i+=1;
1044 else if(direction == kBackDecreasingWindow)
1045 i-=1;
1046 }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1));
1047 }
1048
1049 else if (filterOrder == kBackOrder6) {
1050 do{
1051 for (j = i; j < ssize - i; j++) {
1052 if (smoothing == kFALSE){
1053 a = working_space[ssize + j];
1054 b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
1055 c = 0;
1056 ai = i / 2;
1057 c -= working_space[ssize + j - (int) (2 * ai)] / 6;
1058 c += 4 * working_space[ssize + j - (int) ai] / 6;
1059 c += 4 * working_space[ssize + j + (int) ai] / 6;
1060 c -= working_space[ssize + j + (int) (2 * ai)] / 6;
1061 d = 0;
1062 ai = i / 3;
1063 d += working_space[ssize + j - (int) (3 * ai)] / 20;
1064 d -= 6 * working_space[ssize + j - (int) (2 * ai)] / 20;
1065 d += 15 * working_space[ssize + j - (int) ai] / 20;
1066 d += 15 * working_space[ssize + j + (int) ai] / 20;
1067 d -= 6 * working_space[ssize + j + (int) (2 * ai)] / 20;
1068 d += working_space[ssize + j + (int) (3 * ai)] / 20;
1069 if (b < d)
1070 b = d;
1071 if (b < c)
1072 b = c;
1073 if (b < a)
1074 a = b;
1075 working_space[j] = a;
1076 }
1077
1078 else if (smoothing == kTRUE){
1079 a = working_space[ssize + j];
1080 av = 0;
1081 men = 0;
1082 for (w = j - bw; w <= j + bw; w++){
1083 if ( w >= 0 && w < ssize){
1084 av += working_space[ssize + w];
1085 men +=1;
1086 }
1087 }
1088 av = av / men;
1089 b = 0;
1090 men = 0;
1091 for (w = j - i - bw; w <= j - i + bw; w++){
1092 if ( w >= 0 && w < ssize){
1093 b += working_space[ssize + w];
1094 men +=1;
1095 }
1096 }
1097 b = b / men;
1098 c = 0;
1099 men = 0;
1100 for (w = j + i - bw; w <= j + i + bw; w++){
1101 if ( w >= 0 && w < ssize){
1102 c += working_space[ssize + w];
1103 men +=1;
1104 }
1105 }
1106 c = c / men;
1107 b = (b + c) / 2;
1108 ai = i / 2;
1109 b4 = 0, men = 0;
1110 for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){
1111 if (w >= 0 && w < ssize){
1112 b4 += working_space[ssize + w];
1113 men +=1;
1114 }
1115 }
1116 b4 = b4 / men;
1117 c4 = 0, men = 0;
1118 for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){
1119 if (w >= 0 && w < ssize){
1120 c4 += working_space[ssize + w];
1121 men +=1;
1122 }
1123 }
1124 c4 = c4 / men;
1125 d4 = 0, men = 0;
1126 for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){
1127 if (w >= 0 && w < ssize){
1128 d4 += working_space[ssize + w];
1129 men +=1;
1130 }
1131 }
1132 d4 = d4 / men;
1133 e4 = 0, men = 0;
1134 for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){
1135 if (w >= 0 && w < ssize){
1136 e4 += working_space[ssize + w];
1137 men +=1;
1138 }
1139 }
1140 e4 = e4 / men;
1141 b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6;
1142 ai = i / 3;
1143 b6 = 0, men = 0;
1144 for (w = j - (int)(3 * ai) - bw; w <= j - (int)(3 * ai) + bw; w++){
1145 if (w >= 0 && w < ssize){
1146 b6 += working_space[ssize + w];
1147 men +=1;
1148 }
1149 }
1150 b6 = b6 / men;
1151 c6 = 0, men = 0;
1152 for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){
1153 if (w >= 0 && w < ssize){
1154 c6 += working_space[ssize + w];
1155 men +=1;
1156 }
1157 }
1158 c6 = c6 / men;
1159 d6 = 0, men = 0;
1160 for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){
1161 if (w >= 0 && w < ssize){
1162 d6 += working_space[ssize + w];
1163 men +=1;
1164 }
1165 }
1166 d6 = d6 / men;
1167 e6 = 0, men = 0;
1168 for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){
1169 if (w >= 0 && w < ssize){
1170 e6 += working_space[ssize + w];
1171 men +=1;
1172 }
1173 }
1174 e6 = e6 / men;
1175 f6 = 0, men = 0;
1176 for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){
1177 if (w >= 0 && w < ssize){
1178 f6 += working_space[ssize + w];
1179 men +=1;
1180 }
1181 }
1182 f6 = f6 / men;
1183 g6 = 0, men = 0;
1184 for (w = j + (int)(3 * ai) - bw; w <= j + (int)(3 * ai) + bw; w++){
1185 if (w >= 0 && w < ssize){
1186 g6 += working_space[ssize + w];
1187 men +=1;
1188 }
1189 }
1190 g6 = g6 / men;
1191 b6 = (b6 - 6 * c6 + 15 * d6 + 15 * e6 - 6 * f6 + g6) / 20;
1192 if (b < b6)
1193 b = b6;
1194 if (b < b4)
1195 b = b4;
1196 if (b < a)
1197 av = b;
1198 working_space[j]=av;
1199 }
1200 }
1201 for (j = i; j < ssize - i; j++)
1202 working_space[ssize + j] = working_space[j];
1203 if (direction == kBackIncreasingWindow)
1204 i+=1;
1205 else if(direction == kBackDecreasingWindow)
1206 i-=1;
1207 }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1));
1208 }
1209
1210 else if (filterOrder == kBackOrder8) {
1211 do{
1212 for (j = i; j < ssize - i; j++) {
1213 if (smoothing == kFALSE){
1214 a = working_space[ssize + j];
1215 b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
1216 c = 0;
1217 ai = i / 2;
1218 c -= working_space[ssize + j - (int) (2 * ai)] / 6;
1219 c += 4 * working_space[ssize + j - (int) ai] / 6;
1220 c += 4 * working_space[ssize + j + (int) ai] / 6;
1221 c -= working_space[ssize + j + (int) (2 * ai)] / 6;
1222 d = 0;
1223 ai = i / 3;
1224 d += working_space[ssize + j - (int) (3 * ai)] / 20;
1225 d -= 6 * working_space[ssize + j - (int) (2 * ai)] / 20;
1226 d += 15 * working_space[ssize + j - (int) ai] / 20;
1227 d += 15 * working_space[ssize + j + (int) ai] / 20;
1228 d -= 6 * working_space[ssize + j + (int) (2 * ai)] / 20;
1229 d += working_space[ssize + j + (int) (3 * ai)] / 20;
1230 e = 0;
1231 ai = i / 4;
1232 e -= working_space[ssize + j - (int) (4 * ai)] / 70;
1233 e += 8 * working_space[ssize + j - (int) (3 * ai)] / 70;
1234 e -= 28 * working_space[ssize + j - (int) (2 * ai)] / 70;
1235 e += 56 * working_space[ssize + j - (int) ai] / 70;
1236 e += 56 * working_space[ssize + j + (int) ai] / 70;
1237 e -= 28 * working_space[ssize + j + (int) (2 * ai)] / 70;
1238 e += 8 * working_space[ssize + j + (int) (3 * ai)] / 70;
1239 e -= working_space[ssize + j + (int) (4 * ai)] / 70;
1240 if (b < e)
1241 b = e;
1242 if (b < d)
1243 b = d;
1244 if (b < c)
1245 b = c;
1246 if (b < a)
1247 a = b;
1248 working_space[j] = a;
1249 }
1250
1251 else if (smoothing == kTRUE){
1252 a = working_space[ssize + j];
1253 av = 0;
1254 men = 0;
1255 for (w = j - bw; w <= j + bw; w++){
1256 if ( w >= 0 && w < ssize){
1257 av += working_space[ssize + w];
1258 men +=1;
1259 }
1260 }
1261 av = av / men;
1262 b = 0;
1263 men = 0;
1264 for (w = j - i - bw; w <= j - i + bw; w++){
1265 if ( w >= 0 && w < ssize){
1266 b += working_space[ssize + w];
1267 men +=1;
1268 }
1269 }
1270 b = b / men;
1271 c = 0;
1272 men = 0;
1273 for (w = j + i - bw; w <= j + i + bw; w++){
1274 if ( w >= 0 && w < ssize){
1275 c += working_space[ssize + w];
1276 men +=1;
1277 }
1278 }
1279 c = c / men;
1280 b = (b + c) / 2;
1281 ai = i / 2;
1282 b4 = 0, men = 0;
1283 for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){
1284 if (w >= 0 && w < ssize){
1285 b4 += working_space[ssize + w];
1286 men +=1;
1287 }
1288 }
1289 b4 = b4 / men;
1290 c4 = 0, men = 0;
1291 for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){
1292 if (w >= 0 && w < ssize){
1293 c4 += working_space[ssize + w];
1294 men +=1;
1295 }
1296 }
1297 c4 = c4 / men;
1298 d4 = 0, men = 0;
1299 for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){
1300 if (w >= 0 && w < ssize){
1301 d4 += working_space[ssize + w];
1302 men +=1;
1303 }
1304 }
1305 d4 = d4 / men;
1306 e4 = 0, men = 0;
1307 for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){
1308 if (w >= 0 && w < ssize){
1309 e4 += working_space[ssize + w];
1310 men +=1;
1311 }
1312 }
1313 e4 = e4 / men;
1314 b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6;
1315 ai = i / 3;
1316 b6 = 0, men = 0;
1317 for (w = j - (int)(3 * ai) - bw; w <= j - (int)(3 * ai) + bw; w++){
1318 if (w >= 0 && w < ssize){
1319 b6 += working_space[ssize + w];
1320 men +=1;
1321 }
1322 }
1323 b6 = b6 / men;
1324 c6 = 0, men = 0;
1325 for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){
1326 if (w >= 0 && w < ssize){
1327 c6 += working_space[ssize + w];
1328 men +=1;
1329 }
1330 }
1331 c6 = c6 / men;
1332 d6 = 0, men = 0;
1333 for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){
1334 if (w >= 0 && w < ssize){
1335 d6 += working_space[ssize + w];
1336 men +=1;
1337 }
1338 }
1339 d6 = d6 / men;
1340 e6 = 0, men = 0;
1341 for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){
1342 if (w >= 0 && w < ssize){
1343 e6 += working_space[ssize + w];
1344 men +=1;
1345 }
1346 }
1347 e6 = e6 / men;
1348 f6 = 0, men = 0;
1349 for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){
1350 if (w >= 0 && w < ssize){
1351 f6 += working_space[ssize + w];
1352 men +=1;
1353 }
1354 }
1355 f6 = f6 / men;
1356 g6 = 0, men = 0;
1357 for (w = j + (int)(3 * ai) - bw; w <= j + (int)(3 * ai) + bw; w++){
1358 if (w >= 0 && w < ssize){
1359 g6 += working_space[ssize + w];
1360 men +=1;
1361 }
1362 }
1363 g6 = g6 / men;
1364 b6 = (b6 - 6 * c6 + 15 * d6 + 15 * e6 - 6 * f6 + g6) / 20;
1365 ai = i / 4;
1366 b8 = 0, men = 0;
1367 for (w = j - (int)(4 * ai) - bw; w <= j - (int)(4 * ai) + bw; w++){
1368 if (w >= 0 && w < ssize){
1369 b8 += working_space[ssize + w];
1370 men +=1;
1371 }
1372 }
1373 b8 = b8 / men;
1374 c8 = 0, men = 0;
1375 for (w = j - (int)(3 * ai) - bw; w <= j - (int)(3 * ai) + bw; w++){
1376 if (w >= 0 && w < ssize){
1377 c8 += working_space[ssize + w];
1378 men +=1;
1379 }
1380 }
1381 c8 = c8 / men;
1382 d8 = 0, men = 0;
1383 for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){
1384 if (w >= 0 && w < ssize){
1385 d8 += working_space[ssize + w];
1386 men +=1;
1387 }
1388 }
1389 d8 = d8 / men;
1390 e8 = 0, men = 0;
1391 for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){
1392 if (w >= 0 && w < ssize){
1393 e8 += working_space[ssize + w];
1394 men +=1;
1395 }
1396 }
1397 e8 = e8 / men;
1398 f8 = 0, men = 0;
1399 for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){
1400 if (w >= 0 && w < ssize){
1401 f8 += working_space[ssize + w];
1402 men +=1;
1403 }
1404 }
1405 f8 = f8 / men;
1406 g8 = 0, men = 0;
1407 for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){
1408 if (w >= 0 && w < ssize){
1409 g8 += working_space[ssize + w];
1410 men +=1;
1411 }
1412 }
1413 g8 = g8 / men;
1414 h8 = 0, men = 0;
1415 for (w = j + (int)(3 * ai) - bw; w <= j + (int)(3 * ai) + bw; w++){
1416 if (w >= 0 && w < ssize){
1417 h8 += working_space[ssize + w];
1418 men +=1;
1419 }
1420 }
1421 h8 = h8 / men;
1422 i8 = 0, men = 0;
1423 for (w = j + (int)(4 * ai) - bw; w <= j + (int)(4 * ai) + bw; w++){
1424 if (w >= 0 && w < ssize){
1425 i8 += working_space[ssize + w];
1426 men +=1;
1427 }
1428 }
1429 i8 = i8 / men;
1430 b8 = ( -b8 + 8 * c8 - 28 * d8 + 56 * e8 - 56 * f8 - 28 * g8 + 8 * h8 - i8)/70;
1431 if (b < b8)
1432 b = b8;
1433 if (b < b6)
1434 b = b6;
1435 if (b < b4)
1436 b = b4;
1437 if (b < a)
1438 av = b;
1439 working_space[j]=av;
1440 }
1441 }
1442 for (j = i; j < ssize - i; j++)
1443 working_space[ssize + j] = working_space[j];
1444 if (direction == kBackIncreasingWindow)
1445 i += 1;
1446 else if(direction == kBackDecreasingWindow)
1447 i -= 1;
1448 }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1));
1449 }
1450
1451 if (compton == kTRUE) {
1452 for (i = 0, b2 = 0; i < ssize; i++){
1453 b1 = b2;
1454 a = working_space[i], b = spectrum[i];
1455 j = i;
1456 if (TMath::Abs(a - b) >= 1) {
1457 b1 = i - 1;
1458 if (b1 < 0)
1459 b1 = 0;
1460 yb1 = working_space[b1];
1461 for (b2 = b1 + 1, c = 0, priz = 0; priz == 0 && b2 < ssize; b2++){
1462 a = working_space[b2], b = spectrum[b2];
1463 c = c + b - yb1;
1464 if (TMath::Abs(a - b) < 1) {
1465 priz = 1;
1466 yb2 = b;
1467 }
1468 }
1469 if (b2 == ssize)
1470 b2 -= 1;
1471 yb2 = working_space[b2];
1472 if (yb1 <= yb2){
1473 for (j = b1, c = 0; j <= b2; j++){
1474 b = spectrum[j];
1475 c = c + b - yb1;
1476 }
1477 if (c > 1){
1478 c = (yb2 - yb1) / c;
1479 for (j = b1, d = 0; j <= b2 && j < ssize; j++){
1480 b = spectrum[j];
1481 d = d + b - yb1;
1482 a = c * d + yb1;
1483 working_space[ssize + j] = a;
1484 }
1485 }
1486 }
1487
1488 else{
1489 for (j = b2, c = 0; j >= b1; j--){
1490 b = spectrum[j];
1491 c = c + b - yb2;
1492 }
1493 if (c > 1){
1494 c = (yb1 - yb2) / c;
1495 for (j = b2, d = 0;j >= b1 && j >= 0; j--){
1496 b = spectrum[j];
1497 d = d + b - yb2;
1498 a = c * d + yb2;
1499 working_space[ssize + j] = a;
1500 }
1501 }
1502 }
1503 i=b2;
1504 }
1505 }
1506 }
1507
1508 for (j = 0; j < ssize; j++)
1509 spectrum[j] = working_space[ssize + j];
1510 delete[]working_space;
1511 return 0;
1512}
1513
1514
1515//______________________________________________________________________________
1516const char* NcSpectrum::SmoothMarkov(float *source, int ssize, int averWindow)
1517{
1589
1590 int xmin, xmax, i, l;
1591 float a, b, maxch;
1592 float nom, nip, nim, sp, sm, area = 0;
1593 if(averWindow <= 0)
1594 return "Averaging Window must be positive";
1595 float *working_space = new float[ssize];
1596 xmin = 0,xmax = ssize - 1;
1597 for(i = 0, maxch = 0; i < ssize; i++){
1598 working_space[i]=0;
1599 if(maxch < source[i])
1600 maxch = source[i];
1601
1602 area += source[i];
1603 }
1604 if(maxch == 0) {
1605 delete [] working_space;
1606 return 0 ;
1607 }
1608
1609 nom = 1;
1610 working_space[xmin] = 1;
1611 for(i = xmin; i < xmax; i++){
1612 nip = source[i] / maxch;
1613 nim = source[i + 1] / maxch;
1614 sp = 0,sm = 0;
1615 for(l = 1; l <= averWindow; l++){
1616 if((i + l) > xmax)
1617 a = source[xmax] / maxch;
1618
1619 else
1620 a = source[i + l] / maxch;
1621 b = a - nip;
1622 if(a + nip <= 0)
1623 a = 1;
1624
1625 else
1626 a = TMath::Sqrt(a + nip);
1627 b = b / a;
1628 b = TMath::Exp(b);
1629 sp = sp + b;
1630 if((i - l + 1) < xmin)
1631 a = source[xmin] / maxch;
1632
1633 else
1634 a = source[i - l + 1] / maxch;
1635 b = a - nim;
1636 if(a + nim <= 0)
1637 a = 1;
1638 else
1639 a = TMath::Sqrt(a + nim);
1640 b = b / a;
1641 b = TMath::Exp(b);
1642 sm = sm + b;
1643 }
1644 a = sp / sm;
1645 a = working_space[i + 1] = working_space[i] * a;
1646 nom = nom + a;
1647 }
1648 for(i = xmin; i <= xmax; i++){
1649 working_space[i] = working_space[i] / nom;
1650 }
1651 for(i = 0; i < ssize; i++)
1652 source[i] = working_space[i] * area;
1653 delete[]working_space;
1654 return 0;
1655}
1656
1657
1658//______________________________________________________________________________
1659const char* NcSpectrum::Deconvolution(float *source, const float *response,
1660 int ssize, int numberIterations,
1661 int numberRepetitions, double boost )
1662{
1932
1933 if (ssize <= 0)
1934 return "Wrong Parameters";
1935
1936 if (numberRepetitions <= 0)
1937 return "Wrong Parameters";
1938
1939 // working_space-pointer to the working vector
1940 // (its size must be 4*ssize of source spectrum)
1941 double *working_space = new double[4 * ssize];
1942 int i, j, k, lindex, posit, lh_gold, l, repet;
1943 double lda, ldb, ldc, area, maximum;
1944 area = 0;
1945 lh_gold = -1;
1946 posit = 0;
1947 maximum = 0;
1948
1949//read response vector
1950 for (i = 0; i < ssize; i++) {
1951 lda = response[i];
1952 if (lda != 0)
1953 lh_gold = i + 1;
1954 working_space[i] = lda;
1955 area += lda;
1956 if (lda > maximum) {
1957 maximum = lda;
1958 posit = i;
1959 }
1960 }
1961 if (lh_gold == -1) {
1962 delete [] working_space;
1963 return "ZERO RESPONSE VECTOR";
1964 }
1965
1966//read source vector
1967 for (i = 0; i < ssize; i++)
1968 working_space[2 * ssize + i] = source[i];
1969
1970// create matrix at*a and vector at*y
1971 for (i = 0; i < ssize; i++){
1972 lda = 0;
1973 for (j = 0; j < ssize; j++){
1974 ldb = working_space[j];
1975 k = i + j;
1976 if (k < ssize){
1977 ldc = working_space[k];
1978 lda = lda + ldb * ldc;
1979 }
1980 }
1981 working_space[ssize + i] = lda;
1982 lda = 0;
1983 for (k = 0; k < ssize; k++){
1984 l = k - i;
1985 if (l >= 0){
1986 ldb = working_space[l];
1987 ldc = working_space[2 * ssize + k];
1988 lda = lda + ldb * ldc;
1989 }
1990 }
1991 working_space[3 * ssize + i]=lda;
1992 }
1993
1994// move vector at*y
1995 for (i = 0; i < ssize; i++){
1996 working_space[2 * ssize + i] = working_space[3 * ssize + i];
1997 }
1998
1999//initialization of resulting vector
2000 for (i = 0; i < ssize; i++)
2001 working_space[i] = 1;
2002
2003 //**START OF ITERATIONS**
2004 for (repet = 0; repet < numberRepetitions; repet++) {
2005 if (repet != 0) {
2006 for (i = 0; i < ssize; i++)
2007 working_space[i] = TMath::Power(working_space[i], boost);
2008 }
2009 for (lindex = 0; lindex < numberIterations; lindex++) {
2010 for (i = 0; i < ssize; i++) {
2011 if (working_space[2 * ssize + i] > 0.000001
2012 && working_space[i] > 0.000001) {
2013 lda = 0;
2014 for (j = 0; j < lh_gold; j++) {
2015 ldb = working_space[j + ssize];
2016 if (j != 0){
2017 k = i + j;
2018 ldc = 0;
2019 if (k < ssize)
2020 ldc = working_space[k];
2021 k = i - j;
2022 if (k >= 0)
2023 ldc += working_space[k];
2024 }
2025
2026 else
2027 ldc = working_space[i];
2028 lda = lda + ldb * ldc;
2029 }
2030 ldb = working_space[2 * ssize + i];
2031 if (lda != 0)
2032 lda = ldb / lda;
2033
2034 else
2035 lda = 0;
2036 ldb = working_space[i];
2037 lda = lda * ldb;
2038 working_space[3 * ssize + i] = lda;
2039 }
2040 }
2041 for (i = 0; i < ssize; i++)
2042 working_space[i] = working_space[3 * ssize + i];
2043 }
2044 }
2045
2046//shift resulting spectrum
2047 for (i = 0; i < ssize; i++) {
2048 lda = working_space[i];
2049 j = i + posit;
2050 j = j % ssize;
2051 working_space[ssize + j] = lda;
2052 }
2053
2054//write back resulting spectrum
2055 for (i = 0; i < ssize; i++)
2056 source[i] = area * working_space[ssize + i];
2057 delete[]working_space;
2058 return 0;
2059}
2060
2061
2062//______________________________________________________________________________
2063const char* NcSpectrum::DeconvolutionRL(float *source, const float *response,
2064 int ssize, int numberIterations,
2065 int numberRepetitions, double boost )
2066{
2207
2208 if (ssize <= 0)
2209 return "Wrong Parameters";
2210
2211 if (numberRepetitions <= 0)
2212 return "Wrong Parameters";
2213
2214 // working_space-pointer to the working vector
2215 // (its size must be 4*ssize of source spectrum)
2216 double *working_space = new double[4 * ssize];
2217 int i, j, k, lindex, posit, lh_gold, repet, kmin, kmax;
2218 double lda, ldb, ldc, maximum;
2219 lh_gold = -1;
2220 posit = 0;
2221 maximum = 0;
2222
2223//read response vector
2224 for (i = 0; i < ssize; i++) {
2225 lda = response[i];
2226 if (lda != 0)
2227 lh_gold = i + 1;
2228 working_space[ssize + i] = lda;
2229 if (lda > maximum) {
2230 maximum = lda;
2231 posit = i;
2232 }
2233 }
2234 if (lh_gold == -1) {
2235 delete [] working_space;
2236 return "ZERO RESPONSE VECTOR";
2237 }
2238
2239//read source vector
2240 for (i = 0; i < ssize; i++)
2241 working_space[2 * ssize + i] = source[i];
2242
2243//initialization of resulting vector
2244 for (i = 0; i < ssize; i++){
2245 if (i <= ssize - lh_gold)
2246 working_space[i] = 1;
2247
2248 else
2249 working_space[i] = 0;
2250
2251 }
2252 //**START OF ITERATIONS**
2253 for (repet = 0; repet < numberRepetitions; repet++) {
2254 if (repet != 0) {
2255 for (i = 0; i < ssize; i++)
2256 working_space[i] = TMath::Power(working_space[i], boost);
2257 }
2258 for (lindex = 0; lindex < numberIterations; lindex++) {
2259 for (i = 0; i <= ssize - lh_gold; i++){
2260 lda = 0;
2261 if (working_space[i] > 0){//x[i]
2262 for (j = i; j < i + lh_gold; j++){
2263 ldb = working_space[2 * ssize + j];//y[j]
2264 if (j < ssize){
2265 if (ldb > 0){//y[j]
2266 kmax = j;
2267 if (kmax > lh_gold - 1)
2268 kmax = lh_gold - 1;
2269 kmin = j + lh_gold - ssize;
2270 if (kmin < 0)
2271 kmin = 0;
2272 ldc = 0;
2273 for (k = kmax; k >= kmin; k--){
2274 ldc += working_space[ssize + k] * working_space[j - k];//h[k]*x[j-k]
2275 }
2276 if (ldc > 0)
2277 ldb = ldb / ldc;
2278
2279 else
2280 ldb = 0;
2281 }
2282 ldb = ldb * working_space[ssize + j - i];//y[j]*h[j-i]/suma(h[j][k]x[k])
2283 }
2284 lda += ldb;
2285 }
2286 lda = lda * working_space[i];
2287 }
2288 working_space[3 * ssize + i] = lda;
2289 }
2290 for (i = 0; i < ssize; i++)
2291 working_space[i] = working_space[3 * ssize + i];
2292 }
2293 }
2294
2295//shift resulting spectrum
2296 for (i = 0; i < ssize; i++) {
2297 lda = working_space[i];
2298 j = i + posit;
2299 j = j % ssize;
2300 working_space[ssize + j] = lda;
2301 }
2302
2303//write back resulting spectrum
2304 for (i = 0; i < ssize; i++)
2305 source[i] = working_space[ssize + i];
2306 delete[]working_space;
2307 return 0;
2308}
2309
2310
2311//______________________________________________________________________________
2312const char* NcSpectrum::Unfolding(float *source,
2313 const float **respMatrix,
2314 int ssizex, int ssizey,
2315 int numberIterations,
2316 int numberRepetitions, double boost)
2317{
2417
2418 int i, j, k, lindex, lhx = 0, repet;
2419 double lda, ldb, ldc, area;
2420 if (ssizex <= 0 || ssizey <= 0)
2421 return "Wrong Parameters";
2422 if (ssizex < ssizey)
2423 return "Sizex must be greater than sizey)";
2424 if (numberIterations <= 0)
2425 return "Number of iterations must be positive";
2426 double *working_space =
2427 new double[ssizex * ssizey + 2 * ssizey * ssizey + 4 * ssizex];
2428
2429/*read response matrix*/
2430 for (j = 0; j < ssizey && lhx != -1; j++) {
2431 area = 0;
2432 lhx = -1;
2433 for (i = 0; i < ssizex; i++) {
2434 lda = respMatrix[j][i];
2435 if (lda != 0) {
2436 lhx = i + 1;
2437 }
2438 working_space[j * ssizex + i] = lda;
2439 area = area + lda;
2440 }
2441 if (lhx != -1) {
2442 for (i = 0; i < ssizex; i++)
2443 working_space[j * ssizex + i] /= area;
2444 }
2445 }
2446 if (lhx == -1) {
2447 delete [] working_space;
2448 return ("ZERO COLUMN IN RESPONSE MATRIX");
2449 }
2450
2451/*read source vector*/
2452 for (i = 0; i < ssizex; i++)
2453 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i] =
2454 source[i];
2455
2456/*create matrix at*a + at*y */
2457 for (i = 0; i < ssizey; i++) {
2458 for (j = 0; j < ssizey; j++) {
2459 lda = 0;
2460 for (k = 0; k < ssizex; k++) {
2461 ldb = working_space[ssizex * i + k];
2462 ldc = working_space[ssizex * j + k];
2463 lda = lda + ldb * ldc;
2464 }
2465 working_space[ssizex * ssizey + ssizey * i + j] = lda;
2466 }
2467 lda = 0;
2468 for (k = 0; k < ssizex; k++) {
2469 ldb = working_space[ssizex * i + k];
2470 ldc =
2471 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex +
2472 k];
2473 lda = lda + ldb * ldc;
2474 }
2475 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i] =
2476 lda;
2477 }
2478
2479/*move vector at*y*/
2480 for (i = 0; i < ssizey; i++)
2481 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i] =
2482 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i];
2483
2484/*create matrix at*a*at*a + vector at*a*at*y */
2485 for (i = 0; i < ssizey; i++) {
2486 for (j = 0; j < ssizey; j++) {
2487 lda = 0;
2488 for (k = 0; k < ssizey; k++) {
2489 ldb = working_space[ssizex * ssizey + ssizey * i + k];
2490 ldc = working_space[ssizex * ssizey + ssizey * j + k];
2491 lda = lda + ldb * ldc;
2492 }
2493 working_space[ssizex * ssizey + ssizey * ssizey + ssizey * i + j] =
2494 lda;
2495 }
2496 lda = 0;
2497 for (k = 0; k < ssizey; k++) {
2498 ldb = working_space[ssizex * ssizey + ssizey * i + k];
2499 ldc =
2500 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex +
2501 k];
2502 lda = lda + ldb * ldc;
2503 }
2504 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i] =
2505 lda;
2506 }
2507
2508/*move at*a*at*y*/
2509 for (i = 0; i < ssizey; i++)
2510 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i] =
2511 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i];
2512
2513/*initialization in resulting vector */
2514 for (i = 0; i < ssizey; i++)
2515 working_space[ssizex * ssizey + 2 * ssizey * ssizey + i] = 1;
2516
2517
2518 for (repet = 0; repet < numberRepetitions; repet++) {
2519 if (repet != 0) {
2520 for (i = 0; i < ssizey; i++)
2521 working_space[ssizex * ssizey + 2 * ssizey * ssizey + i] = TMath::Power(working_space[ssizex * ssizey + 2 * ssizey * ssizey + i], boost);
2522 }
2523 for (lindex = 0; lindex < numberIterations; lindex++) {
2524 for (i = 0; i < ssizey; i++) {
2525 lda = 0;
2526 for (j = 0; j < ssizey; j++) {
2527 ldb =
2528 working_space[ssizex * ssizey + ssizey * ssizey + ssizey * i + j];
2529 ldc = working_space[ssizex * ssizey + 2 * ssizey * ssizey + j];
2530 lda = lda + ldb * ldc;
2531 }
2532 ldb =
2533 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i];
2534 if (lda != 0) {
2535 lda = ldb / lda;
2536 }
2537
2538 else
2539 lda = 0;
2540 ldb = working_space[ssizex * ssizey + 2 * ssizey * ssizey + i];
2541 lda = lda * ldb;
2542 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i] = lda;
2543 }
2544 for (i = 0; i < ssizey; i++)
2545 working_space[ssizex * ssizey + 2 * ssizey * ssizey + i] =
2546 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i];
2547 }
2548 }
2549
2550/*write back resulting spectrum*/
2551 for (i = 0; i < ssizex; i++) {
2552 if (i < ssizey)
2553 source[i] = working_space[ssizex * ssizey + 2 * ssizey * ssizey + i];
2554
2555 else
2556 source[i] = 0;
2557 }
2558 delete[]working_space;
2559 return 0;
2560}
2561
2562
2563//______________________________________________________________________________
2564Int_t NcSpectrum::SearchHighRes(float *source,float *destVector, int ssize,
2565 float sigma, double threshold,
2566 bool backgroundRemove,int deconIterations,
2567 bool markov, int averWindow)
2568{
2784
2785 int i, j, numberIterations = (int)(7 * sigma + 0.5);
2786 double a, b, c;
2787 int k, lindex, posit, imin, imax, jmin, jmax, lh_gold, priz;
2788 double lda, ldb, ldc, area, maximum, maximum_decon;
2789 int xmin, xmax, l, peak_index = 0, size_ext = ssize + 2 * numberIterations, shift = numberIterations, bw = 2, w;
2790 double maxch;
2791 double nom, nip, nim, sp, sm, plocha = 0;
2792 double m0low=0,m1low=0,m2low=0,l0low=0,l1low=0,detlow,av,men;
2793 if (sigma < 1) {
2794 Error("SearchHighRes", "Invalid sigma, must be greater than or equal to 1");
2795 return 0;
2796 }
2797
2798 if(threshold<=0 || threshold>=100){
2799 Error("SearchHighRes", "Invalid threshold, must be positive and less than 100");
2800 return 0;
2801 }
2802
2803 j = (int) (5.0 * sigma + 0.5);
2804 if (j >= PEAK_WINDOW / 2) {
2805 Error("SearchHighRes", "Too large sigma");
2806 return 0;
2807 }
2808
2809 if (markov == true) {
2810 if (averWindow <= 0) {
2811 Error("SearchHighRes", "Averanging window must be positive");
2812 return 0;
2813 }
2814 }
2815
2816 if(backgroundRemove == true){
2817 if(ssize < 2 * numberIterations + 1){
2818 Error("SearchHighRes", "Too large clipping window");
2819 return 0;
2820 }
2821 }
2822
2823 k = int(2 * sigma+0.5);
2824 if(k >= 2){
2825 for(i = 0;i < k;i++){
2826 a = i,b = source[i];
2827 m0low += 1,m1low += a,m2low += a * a,l0low += b,l1low += a * b;
2828 }
2829 detlow = m0low * m2low - m1low * m1low;
2830 if(detlow != 0)
2831 l1low = (-l0low * m1low + l1low * m0low) / detlow;
2832
2833 else
2834 l1low = 0;
2835 if(l1low > 0)
2836 l1low=0;
2837 }
2838
2839 else{
2840 l1low = 0;
2841 }
2842
2843 i = (int)(7 * sigma + 0.5);
2844 i = 2 * i;
2845 double *working_space = new double [7 * (ssize + i)];
2846 for (j=0;j<7 * (ssize + i);j++) working_space[j] = 0;
2847 for(i = 0; i < size_ext; i++){
2848 if(i < shift){
2849 a = i - shift;
2850 working_space[i + size_ext] = source[0] + l1low * a;
2851 if(working_space[i + size_ext] < 0)
2852 working_space[i + size_ext]=0;
2853 }
2854
2855 else if(i >= ssize + shift){
2856 a = i - (ssize - 1 + shift);
2857 working_space[i + size_ext] = source[ssize - 1];
2858 if(working_space[i + size_ext] < 0)
2859 working_space[i + size_ext]=0;
2860 }
2861
2862 else
2863 working_space[i + size_ext] = source[i - shift];
2864 }
2865
2866 if(backgroundRemove == true){
2867 for(i = 1; i <= numberIterations; i++){
2868 for(j = i; j < size_ext - i; j++){
2869 if(markov == false){
2870 a = working_space[size_ext + j];
2871 b = (working_space[size_ext + j - i] + working_space[size_ext + j + i]) / 2.0;
2872 if(b < a)
2873 a = b;
2874
2875 working_space[j]=a;
2876 }
2877
2878 else{
2879 a = working_space[size_ext + j];
2880 av = 0;
2881 men = 0;
2882 for (w = j - bw; w <= j + bw; w++){
2883 if ( w >= 0 && w < size_ext){
2884 av += working_space[size_ext + w];
2885 men +=1;
2886 }
2887 }
2888 av = av / men;
2889 b = 0;
2890 men = 0;
2891 for (w = j - i - bw; w <= j - i + bw; w++){
2892 if ( w >= 0 && w < size_ext){
2893 b += working_space[size_ext + w];
2894 men +=1;
2895 }
2896 }
2897 b = b / men;
2898 c = 0;
2899 men = 0;
2900 for (w = j + i - bw; w <= j + i + bw; w++){
2901 if ( w >= 0 && w < size_ext){
2902 c += working_space[size_ext + w];
2903 men +=1;
2904 }
2905 }
2906 c = c / men;
2907 b = (b + c) / 2;
2908 if (b < a)
2909 av = b;
2910 working_space[j]=av;
2911 }
2912 }
2913 for(j = i; j < size_ext - i; j++)
2914 working_space[size_ext + j] = working_space[j];
2915 }
2916 for(j = 0;j < size_ext; j++){
2917 if(j < shift){
2918 a = j - shift;
2919 b = source[0] + l1low * a;
2920 if (b < 0) b = 0;
2921 working_space[size_ext + j] = b - working_space[size_ext + j];
2922 }
2923
2924 else if(j >= ssize + shift){
2925 a = j - (ssize - 1 + shift);
2926 b = source[ssize - 1];
2927 if (b < 0) b = 0;
2928 working_space[size_ext + j] = b - working_space[size_ext + j];
2929 }
2930
2931 else{
2932 working_space[size_ext + j] = source[j - shift] - working_space[size_ext + j];
2933 }
2934 }
2935 for(j = 0;j < size_ext; j++){
2936 if(working_space[size_ext + j] < 0) working_space[size_ext + j] = 0;
2937 }
2938 }
2939
2940 for(i = 0; i < size_ext; i++){
2941 working_space[i + 6*size_ext] = working_space[i + size_ext];
2942 }
2943
2944 if(markov == true){
2945 for(j = 0; j < size_ext; j++)
2946 working_space[2 * size_ext + j] = working_space[size_ext + j];
2947 xmin = 0,xmax = size_ext - 1;
2948 for(i = 0, maxch = 0; i < size_ext; i++){
2949 working_space[i] = 0;
2950 if(maxch < working_space[2 * size_ext + i])
2951 maxch = working_space[2 * size_ext + i];
2952 plocha += working_space[2 * size_ext + i];
2953 }
2954 if(maxch == 0) {
2955 delete [] working_space;
2956 return 0;
2957 }
2958
2959 nom = 1;
2960 working_space[xmin] = 1;
2961 for(i = xmin; i < xmax; i++){
2962 nip = working_space[2 * size_ext + i] / maxch;
2963 nim = working_space[2 * size_ext + i + 1] / maxch;
2964 sp = 0,sm = 0;
2965 for(l = 1; l <= averWindow; l++){
2966 if((i + l) > xmax)
2967 a = working_space[2 * size_ext + xmax] / maxch;
2968
2969 else
2970 a = working_space[2 * size_ext + i + l] / maxch;
2971
2972 b = a - nip;
2973 if(a + nip <= 0)
2974 a=1;
2975
2976 else
2977 a = TMath::Sqrt(a + nip);
2978
2979 b = b / a;
2980 b = TMath::Exp(b);
2981 sp = sp + b;
2982 if((i - l + 1) < xmin)
2983 a = working_space[2 * size_ext + xmin] / maxch;
2984
2985 else
2986 a = working_space[2 * size_ext + i - l + 1] / maxch;
2987
2988 b = a - nim;
2989 if(a + nim <= 0)
2990 a = 1;
2991
2992 else
2993 a = TMath::Sqrt(a + nim);
2994
2995 b = b / a;
2996 b = TMath::Exp(b);
2997 sm = sm + b;
2998 }
2999 a = sp / sm;
3000 a = working_space[i + 1] = working_space[i] * a;
3001 nom = nom + a;
3002 }
3003 for(i = xmin; i <= xmax; i++){
3004 working_space[i] = working_space[i] / nom;
3005 }
3006 for(j = 0; j < size_ext; j++)
3007 working_space[size_ext + j] = working_space[j] * plocha;
3008 for(j = 0; j < size_ext; j++){
3009 working_space[2 * size_ext + j] = working_space[size_ext + j];
3010 }
3011 if(backgroundRemove == true){
3012 for(i = 1; i <= numberIterations; i++){
3013 for(j = i; j < size_ext - i; j++){
3014 a = working_space[size_ext + j];
3015 b = (working_space[size_ext + j - i] + working_space[size_ext + j + i]) / 2.0;
3016 if(b < a)
3017 a = b;
3018 working_space[j] = a;
3019 }
3020 for(j = i; j < size_ext - i; j++)
3021 working_space[size_ext + j] = working_space[j];
3022 }
3023 for(j = 0; j < size_ext; j++){
3024 working_space[size_ext + j] = working_space[2 * size_ext + j] - working_space[size_ext + j];
3025 }
3026 }
3027 }
3028//deconvolution starts
3029 area = 0;
3030 lh_gold = -1;
3031 posit = 0;
3032 maximum = 0;
3033//generate response vector
3034 for(i = 0; i < size_ext; i++){
3035 lda = (double)i - 3 * sigma;
3036 lda = lda * lda / (2 * sigma * sigma);
3037 j = (int)(1000 * TMath::Exp(-lda));
3038 lda = j;
3039 if(lda != 0)
3040 lh_gold = i + 1;
3041
3042 working_space[i] = lda;
3043 area = area + lda;
3044 if(lda > maximum){
3045 maximum = lda;
3046 posit = i;
3047 }
3048 }
3049//read source vector
3050 for(i = 0; i < size_ext; i++)
3051 working_space[2 * size_ext + i] = TMath::Abs(working_space[size_ext + i]);
3052//create matrix at*a(vector b)
3053 i = lh_gold - 1;
3054 if(i > size_ext)
3055 i = size_ext;
3056
3057 imin = -i,imax = i;
3058 for(i = imin; i <= imax; i++){
3059 lda = 0;
3060 jmin = 0;
3061 if(i < 0)
3062 jmin = -i;
3063 jmax = lh_gold - 1 - i;
3064 if(jmax > (lh_gold - 1))
3065 jmax = lh_gold - 1;
3066
3067 for(j = jmin;j <= jmax; j++){
3068 ldb = working_space[j];
3069 ldc = working_space[i + j];
3070 lda = lda + ldb * ldc;
3071 }
3072 working_space[size_ext + i - imin] = lda;
3073 }
3074//create vector p
3075 i = lh_gold - 1;
3076 imin = -i,imax = size_ext + i - 1;
3077 for(i = imin; i <= imax; i++){
3078 lda = 0;
3079 for(j = 0; j <= (lh_gold - 1); j++){
3080 ldb = working_space[j];
3081 k = i + j;
3082 if(k >= 0 && k < size_ext){
3083 ldc = working_space[2 * size_ext + k];
3084 lda = lda + ldb * ldc;
3085 }
3086
3087 }
3088 working_space[4 * size_ext + i - imin] = lda;
3089 }
3090//move vector p
3091 for(i = imin; i <= imax; i++)
3092 working_space[2 * size_ext + i - imin] = working_space[4 * size_ext + i - imin];
3093//initialization of resulting vector
3094 for(i = 0; i < size_ext; i++)
3095 working_space[i] = 1;
3096//START OF ITERATIONS
3097 for(lindex = 0; lindex < deconIterations; lindex++){
3098 for(i = 0; i < size_ext; i++){
3099 if(TMath::Abs(working_space[2 * size_ext + i]) > 0.00001 && TMath::Abs(working_space[i]) > 0.00001){
3100 lda=0;
3101 jmin = lh_gold - 1;
3102 if(jmin > i)
3103 jmin = i;
3104
3105 jmin = -jmin;
3106 jmax = lh_gold - 1;
3107 if(jmax > (size_ext - 1 - i))
3108 jmax=size_ext-1-i;
3109
3110 for(j = jmin; j <= jmax; j++){
3111 ldb = working_space[j + lh_gold - 1 + size_ext];
3112 ldc = working_space[i + j];
3113 lda = lda + ldb * ldc;
3114 }
3115 ldb = working_space[2 * size_ext + i];
3116 if(lda != 0)
3117 lda = ldb / lda;
3118
3119 else
3120 lda = 0;
3121
3122 ldb = working_space[i];
3123 lda = lda * ldb;
3124 working_space[3 * size_ext + i] = lda;
3125 }
3126 }
3127 for(i = 0; i < size_ext; i++){
3128 working_space[i] = working_space[3 * size_ext + i];
3129 }
3130 }
3131//shift resulting spectrum
3132 for(i=0;i<size_ext;i++){
3133 lda = working_space[i];
3134 j = i + posit;
3135 j = j % size_ext;
3136 working_space[size_ext + j] = lda;
3137 }
3138//write back resulting spectrum
3139 maximum = 0, maximum_decon = 0;
3140 j = lh_gold - 1;
3141 for(i = 0; i < size_ext - j; i++){
3142 if(i >= shift && i < ssize + shift){
3143 working_space[i] = area * working_space[size_ext + i + j];
3144 if(maximum_decon < working_space[i])
3145 maximum_decon = working_space[i];
3146 if(maximum < working_space[6 * size_ext + i])
3147 maximum = working_space[6 * size_ext + i];
3148 }
3149
3150 else
3151 working_space[i] = 0;
3152 }
3153 lda=1;
3154 if(lda>threshold)
3155 lda=threshold;
3156 lda=lda/100;
3157
3158//searching for peaks in deconvolved spectrum
3159 for(i = 1; i < size_ext - 1; i++){
3160 if(working_space[i] > working_space[i - 1] && working_space[i] > working_space[i + 1]){
3161 if(i >= shift && i < ssize + shift){
3162 if(working_space[i] > lda*maximum_decon && working_space[6 * size_ext + i] > threshold * maximum / 100.0){
3163 for(j = i - 1, a = 0, b = 0; j <= i + 1; j++){
3164 a += (double)(j - shift) * working_space[j];
3165 b += working_space[j];
3166 }
3167 a = a / b;
3168 if(a < 0)
3169 a = 0;
3170
3171 if(a >= ssize)
3172 a = ssize - 1;
3173 if(peak_index == 0){
3174 fPositionX[0] = a;
3175 peak_index = 1;
3176 }
3177
3178 else{
3179 for(j = 0, priz = 0; j < peak_index && priz == 0; j++){
3180 if(working_space[6 * size_ext + shift + (int)a] > working_space[6 * size_ext + shift + (int)fPositionX[j]])
3181 priz = 1;
3182 }
3183 if(priz == 0){
3184 if(j < fMaxPeaks){
3185 fPositionX[j] = a;
3186 }
3187 }
3188
3189 else{
3190 for(k = peak_index; k >= j; k--){
3191 if(k < fMaxPeaks){
3192 fPositionX[k] = fPositionX[k - 1];
3193 }
3194 }
3195 fPositionX[j - 1] = a;
3196 }
3197 if(peak_index < fMaxPeaks)
3198 peak_index += 1;
3199 }
3200 }
3201 }
3202 }
3203 }
3204
3205 for(i = 0; i < ssize; i++) destVector[i] = working_space[i + shift];
3206 delete [] working_space;
3207 fNPeaks = peak_index;
3208 if(peak_index == fMaxPeaks)
3209 Warning("SearchHighRes", "Peak buffer full");
3210 return fNPeaks;
3211}
3212
3213
3214//______________________________________________________________________________
3215Int_t NcSpectrum::Search1HighRes(float *source,float *destVector, int ssize,
3216 float sigma, double threshold,
3217 bool backgroundRemove,int deconIterations,
3218 bool markov, int averWindow)
3219{
3226
3227 return SearchHighRes(source,destVector,ssize,sigma,threshold,backgroundRemove,
3228 deconIterations,markov,averWindow);
3229}
3230
3231
3232//______________________________________________________________________________
3233Int_t NcSpectrum::StaticSearch(const TH1 *hist, Double_t sigma, Option_t *option, Double_t threshold)
3234{
3240
3241 NcSpectrum s;
3242 return s.Search(hist,sigma,option,threshold);
3243}
3244
3245
3246//______________________________________________________________________________
3247TH1* NcSpectrum::StaticBackground(const TH1 *hist,Int_t niter, Option_t *option)
3248{
3254
3255 NcSpectrum s;
3256 return s.Background(hist,niter,option);
3257}
PEAK_WINDOW
#define PEAK_WINDOW
Definition NcSpectrum.cxx:75
ClassImp
ClassImp(NcSpectrum)
NcSpectrum
Facilities for advanced spectral analysis.
Definition NcSpectrum.h:26
NcSpectrum::StaticBackground
static TH1 * StaticBackground(const TH1 *hist, Int_t niter=20, Option_t *option="")
Definition NcSpectrum.cxx:3247
NcSpectrum::StaticSearch
static Int_t StaticSearch(const TH1 *hist, Double_t sigma=2, Option_t *option="goff", Double_t threshold=0.05)
Definition NcSpectrum.cxx:3233
NcSpectrum::Search
virtual Int_t Search(const TH1 *hist, Double_t sigma=2, Option_t *option="", Double_t threshold=0.05)
Definition NcSpectrum.cxx:295
NcSpectrum::fgAverageWindow
static Int_t fgAverageWindow
Definition NcSpectrum.h:40
NcSpectrum::fResolution
Float_t fResolution
Definition NcSpectrum.h:38
NcSpectrum::fPositionX
Float_t * fPositionX
Definition NcSpectrum.h:36
NcSpectrum::SetDeconIterations
static void SetDeconIterations(Int_t n=3)
Definition NcSpectrum.cxx:160
NcSpectrum::kBackSmoothing11
@ kBackSmoothing11
Definition NcSpectrum.h:55
NcSpectrum::kBackSmoothing15
@ kBackSmoothing15
Definition NcSpectrum.h:57
NcSpectrum::kBackSmoothing7
@ kBackSmoothing7
Definition NcSpectrum.h:53
NcSpectrum::kBackSmoothing13
@ kBackSmoothing13
Definition NcSpectrum.h:56
NcSpectrum::kBackSmoothing9
@ kBackSmoothing9
Definition NcSpectrum.h:54
NcSpectrum::kBackIncreasingWindow
@ kBackIncreasingWindow
Definition NcSpectrum.h:49
NcSpectrum::kBackDecreasingWindow
@ kBackDecreasingWindow
Definition NcSpectrum.h:50
NcSpectrum::kBackSmoothing5
@ kBackSmoothing5
Definition NcSpectrum.h:52
NcSpectrum::kBackSmoothing3
@ kBackSmoothing3
Definition NcSpectrum.h:51
NcSpectrum::~NcSpectrum
virtual ~NcSpectrum()
Definition NcSpectrum.cxx:130
NcSpectrum::Deconvolution
const char * Deconvolution(float *source, const float *response, Int_t ssize, Int_t numberIterations, Int_t numberRepetitions, Double_t boost)
Definition NcSpectrum.cxx:1659
NcSpectrum::fgIterations
static Int_t fgIterations
Definition NcSpectrum.h:41
NcSpectrum::Search1HighRes
Int_t Search1HighRes(float *source, float *destVector, Int_t ssize, float sigma, Double_t threshold, bool backgroundRemove, Int_t deconIterations, bool markov, Int_t averWindow)
Definition NcSpectrum.cxx:3215
NcSpectrum::SearchHighRes
Int_t SearchHighRes(float *source, float *destVector, Int_t ssize, float sigma, Double_t threshold, bool backgroundRemove, Int_t deconIterations, bool markov, Int_t averWindow)
Definition NcSpectrum.cxx:2564
NcSpectrum::fPosition
Float_t * fPosition
Definition NcSpectrum.h:35
NcSpectrum::NcSpectrum
NcSpectrum(const NcSpectrum &)
NcSpectrum::Print
virtual void Print(Option_t *option="") const
Definition NcSpectrum.cxx:279
NcSpectrum::SetAverageWindow
static void SetAverageWindow(Int_t w=3)
Definition NcSpectrum.cxx:146
NcSpectrum::fMaxPeaks
Int_t fMaxPeaks
Definition NcSpectrum.h:33
NcSpectrum::fPositionY
Float_t * fPositionY
Definition NcSpectrum.h:37
NcSpectrum::SmoothMarkov
const char * SmoothMarkov(float *source, Int_t ssize, Int_t averWindow)
Definition NcSpectrum.cxx:1516
NcSpectrum::SetResolution
void SetResolution(Float_t resolution=1)
Definition NcSpectrum.cxx:414
NcSpectrum::Unfolding
const char * Unfolding(float *source, const float **respMatrix, Int_t ssizex, Int_t ssizey, Int_t numberIterations, Int_t numberRepetitions, Double_t boost)
Definition NcSpectrum.cxx:2312
NcSpectrum::DeconvolutionRL
const char * DeconvolutionRL(float *source, const float *response, Int_t ssize, Int_t numberIterations, Int_t numberRepetitions, Double_t boost)
Definition NcSpectrum.cxx:2063
NcSpectrum::fHistogram
TH1 * fHistogram
Definition NcSpectrum.h:39
NcSpectrum::fNPeaks
Int_t fNPeaks
Definition NcSpectrum.h:34
NcSpectrum::Background
virtual TH1 * Background(const TH1 *hist, Int_t niter=20, Option_t *option="")