GSI Object Oriented Online Offline (Go4)  GO4-6.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
TGo4FitPeakFinder.cxx
Go to the documentation of this file.
1 // $Id: TGo4FitPeakFinder.cxx 3494 2022-01-21 15:54:00Z linev $
2 //-----------------------------------------------------------------------
3 // The GSI Online Offline Object Oriented (Go4) Project
4 // Experiment Data Processing at EE department, GSI
5 //-----------------------------------------------------------------------
6 // Copyright (C) 2000- GSI Helmholtzzentrum fuer Schwerionenforschung GmbH
7 // Planckstr. 1, 64291 Darmstadt, Germany
8 // Contact: http://go4.gsi.de
9 //-----------------------------------------------------------------------
10 // This software can be used under the license agreements as stated
11 // in Go4License.txt file which is part of the distribution.
12 //-----------------------------------------------------------------------
13 
14 #include "TGo4FitPeakFinder.h"
15 
16 #include "TMath.h"
17 #include "TH1D.h"
18 #include "TMatrixD.h"
19 #include "TVectorD.h"
20 #include "TArrayD.h"
21 #include "TArrayI.h"
22 #include "TSpectrum.h"
23 
24 #include "TGo4Fitter.h"
25 #include "TGo4FitData.h"
26 
28 
29 }
30 
31 TGo4FitPeakFinder::TGo4FitPeakFinder(const char* Name, const char* DataName, Bool_t ClearModels, Int_t PolOrder) :
32  TGo4FitterAction(Name, "Peak finder action"),
33  fiPeakFinderType(0), fxDataName(DataName), fbClearModels(ClearModels),
34  fbUsePolynom(PolOrder>=0), fiPolynomOrder(PolOrder>=0 ? PolOrder : 0),
35  fd0MinWidth(5.), fd0MaxWidth(30.), fd0MaxAmplFactor(0.2),
36  fd1LineWidth(5.),
37  fd2NoiseFactor(2.), fd2NoiseMinimum(5.), fi2ChannelSum(2) {
38 }
39 
41 }
42 
44  if (PolynomOrder<0) SetUsePolynom(kFALSE);
45  else {
46  SetUsePolynom(kTRUE);
47  SetPolynomOrder(PolynomOrder);
48  }
49 }
50 
51 void TGo4FitPeakFinder::SetupForFirst(Double_t MaxAmplFactor, Double_t MinWidth, Double_t MaxWidth) {
53  Set0MaxAmplFactor(MaxAmplFactor);
54  Set0MinWidth(MinWidth);
55  Set0MaxWidth(MaxWidth);
56 }
57 
58 void TGo4FitPeakFinder::SetupForSecond(Double_t LineWidth) {
60  Set1LineWidth(LineWidth);
61 }
62 
63 void TGo4FitPeakFinder::SetupForThird(Double_t NoiseFactor, Double_t NoiseMinimum, Int_t ChannelSum) {
65  Set2NoiseFactor(NoiseFactor);
66  Set2NoiseMinimum(NoiseMinimum);
67  Set2ChannelSum(ChannelSum);
68 }
69 
71  TGo4Fitter* fitter = dynamic_cast<TGo4Fitter*> (Fitter);
72  if (fitter==0) return;
73 
74  TGo4FitData* data = fitter->FindData(GetDataName());
75  if (data==0) return;
76 
77  if (GetClearModels())
78  fitter->DeleteModelsAssosiatedTo(data->GetName());
79 
80  if (GetPeakFinderType()==0)
81  SergeyLinevPeakFinder(fitter,
82  data,
83  GetUsePolynom() ? GetPolynomOrder() : -1,
85  Get0MinWidth(),
86  Get0MaxWidth());
87 
88  if (GetPeakFinderType()==1)
89  ROOTPeakFinder(fitter,
90  data,
91  GetUsePolynom() ? GetPolynomOrder() : -1,
92  Get1LineWidth());
93 
94  if (GetPeakFinderType()==2)
95  HansEsselPeakFinder(fitter,
96  data,
97  500,
101  GetUsePolynom() ? GetPolynomOrder() : -1);
102 }
103 
104 void TGo4FitPeakFinder::Print(Option_t* option) const {
105  TGo4FitterAction::Print(option);
106 }
107 /*
108 void TGo4FitPeakFinder::ROOTPeakFinder(TGo4Fitter* fitter, TGo4FitData* data, Int_t PolynomOrder, Double_t Sigma) {
109  if ((fitter==0) || (data==0)) return;
110 
111  TGo4FitDataIter* iter = data->MakeIter();
112  if (iter==0) return;
113 
114  Int_t size = iter->CountPoints();
115 
116  if ((size<10) || (iter->ScalesSize()!=1)) { delete iter; return; }
117 
118  TArrayF Bins(size), Scales(size);
119 
120  Int_t nbin = 0;
121  if (iter->Reset()) do {
122  Bins[nbin] = data->GetAmplValue() * iter->Value();
123  Scales[nbin] = iter->x();
124  nbin++;
125  } while (iter->Next());
126 
127  delete iter;
128 
129  if (PolynomOrder>=0)
130  fitter->AddPolynomX(data->GetName(), "Pol", PolynomOrder);
131 
132  TSpectrum sp(100);
133  sp.Search1(Bins.GetArray(), size, Sigma);
134 
135  for(Int_t n=0;n<sp.GetNPeaks();n++) {
136  Double_t dindx = sp.GetPositionX()[n];
137  Int_t left = TMath::Nint(dindx-Sigma/2.);
138  Int_t right = TMath::Nint(dindx+Sigma/2.);
139 
140  Double_t pos = Scales[TMath::Nint(dindx)];
141  Double_t width = TMath::Abs(Scales[right]-Scales[left]);
142 
143  fitter->AddGauss1(data->GetName(), fitter->FindNextName("Gauss",n), pos, width);
144  }
145 }
146 */
147 
148 void TGo4FitPeakFinder::ROOTPeakFinder(TGo4Fitter* fitter, TGo4FitData* data, Int_t PolynomOrder, Double_t Sigma) {
149  if ((fitter==0) || (data==0)) return;
150 
151  TGo4FitDataIter* iter = data->MakeIter();
152  if (iter==0) return;
153 
154  Int_t size = iter->CountPoints();
155 
156  if ((size<10) || (iter->ScalesSize()!=1)) { delete iter; return; }
157 
158  TArrayD Bins(size), Scales(size), HScales(size+1);
159 
160  Int_t nbin = 0;
161  if (iter->Reset()) do {
162  Bins[nbin] = data->GetAmplValue() * iter->Value();
163  Double_t pos = iter->x();
164  Double_t width = iter->HasWidths() ? iter->Widths()[0] : 1.;
165  Scales[nbin] = pos;
166  HScales[nbin] = pos - width/2.;
167  HScales[nbin+1] = pos + width/2;
168  nbin++;
169  } while (iter->Next());
170 
171  delete iter;
172 
173  if (PolynomOrder>=0)
174  fitter->AddPolynomX(data->GetName(), "Pol", PolynomOrder);
175 
176  TH1D histo("ROOTPeakFinder", "ROOTPeakFinder spectrum", size, HScales.GetArray());
177  histo.SetDirectory(0);
178  for (Int_t i=0;i<size;i++)
179  histo.SetBinContent(i+1, Bins[i]);
180 
181  TSpectrum sp(100);
182  sp.Search(&histo, Sigma);
183 // sp.Search1(Bins.GetArray(), size, Sigma);
184 
185  for(Int_t n=0;n<sp.GetNPeaks();n++) {
186  Double_t pos = sp.GetPositionX()[n];
187  Double_t width = Sigma;
188 
189 
190 // this is old code for ROOT 3.*, may be usefull for somebody
191 /*
192  Double_t dindx = sp.GetPositionX()[n];
193  Int_t left = TMath::Nint(dindx-Sigma/2.);
194  Int_t right = TMath::Nint(dindx+Sigma/2.);
195  if (left<0) left = 0;
196  if (right>=size) right = size-1;
197 
198  Double_t pos = Scales[TMath::Nint(dindx)];
199  Double_t width = TMath::Abs(Scales[right]-Scales[left]);
200 */
201  fitter->AddGauss1(data->GetName(), fitter->FindNextName("Gauss",n), pos, width);
202  }
203 }
204 
205 
206 Double_t TGo4FitPeakFinder::CalcPolynom(const TArrayD& Coef, Double_t x) {
207  Double_t power = 1., sum = 0.;
208  for(Int_t n=0;n<Coef.GetSize();n++)
209  { sum+=Coef[n]*power; power*=x; }
210  return sum;
211 }
212 
213 void TGo4FitPeakFinder::DefinePolynom(Int_t size, Double_t* bins, Double_t* scales, TArrayD& Coef,
214  Double_t* weight, Double_t* backgr, Char_t* use) {
215  Int_t maxorder = Coef.GetSize()-1;
216 
217  TMatrixD matr(maxorder+1,maxorder+1);
218  TVectorD vector(maxorder+1);
219 
220  for(Int_t order1=0;order1<=maxorder;order1++)
221  for(Int_t order2=order1;order2<=maxorder;order2++) {
222  Double_t sum = 0;
223  for(Int_t i=0;i<size;i++)
224  if((use==0) || use[i]) {
225  Double_t zn = TMath::Power(scales[i],order1+order2);
226  if (weight) zn*=weight[i];
227  sum+=zn;
228  }
229  matr(order1,order2) = sum;
230  matr(order2,order1) = sum;
231  }
232 
233  for(Int_t order1=0;order1<=maxorder;order1++) {
234  Double_t sum = 0;
235  for(Int_t i=0;i<size;i++)
236  if((use==0) || use[i]) {
237  Double_t zn = bins[i];
238  if (backgr) zn-=backgr[i];
239  zn*=TMath::Power(scales[i],order1);
240  if (weight) zn*=weight[i];
241  sum+=zn;
242  }
243  vector(order1) = sum;
244  }
245 
246  matr.Invert();
247  vector*=matr;
248 
249  for(Int_t order=0;order<=maxorder;order++)
250  Coef[order]=vector(order);
251 }
252 
253 void TGo4FitPeakFinder::DefinePolynomEx(Int_t size, Double_t* bins, Double_t* scales, Double_t* weight, Double_t* backgr,
254  Int_t lbound, Int_t rbound, TArrayD& Coef) {
255  TArrayC use(size); use.Reset(0);
256  if (lbound<0) lbound = 0;
257  if (rbound>=size) rbound=size-1;
258  for(Int_t i=lbound;i<=rbound;i++) use[i]=1;
259  DefinePolynom(size, bins, scales, Coef, weight, backgr, use.GetArray());
260 }
261 
262 Bool_t FindValue(Double_t* bins, Double_t* backgr, Char_t* use, Int_t size, Int_t& p, Int_t dir, Double_t value) {
263  p += dir;
264  while ((p>=0) && (p<size) && use[p]) {
265  Double_t zn = bins[p]-backgr[p];
266  if (zn<value) return kTRUE;
267  p+=dir;
268  }
269  return kFALSE;
270 }
271 
273  TGo4FitData* data,
274  Int_t PolOrder,
275  Double_t AmplThreshold,
276  Double_t MinWidth,
277  Double_t MaxWidth) {
278  if ((fitter==0) || (data==0)) return;
279 
280  TGo4FitDataIter* iter = data->MakeIter();
281  if (iter==0) return;
282 
283  Int_t size = iter->CountPoints();
284 
285  if ((size<10) || (iter->ScalesSize()!=1)) { delete iter; return; }
286 
287  TArrayD Bins(size), Scales(size), Weights(size), Background(size);
288 
289  Weights.Reset(1.); Background.Reset(0.);
290 
291  Double_t* bins = Bins.GetArray();
292  Double_t* scales = Scales.GetArray();
293  Double_t* weights = Weights.GetArray();
294  Double_t* backgr = Background.GetArray();
295 
296  Int_t nbin = 0;
297  if (iter->Reset()) do {
298  bins[nbin] = data->GetAmplValue() * iter->Value();
299  if (iter->StandardDeviation()>0) weights[nbin] = 1./iter->StandardDeviation();
300  scales[nbin] = iter->x();
301  nbin++;
302  } while (iter->Next());
303 
304  delete iter;
305 
306  if ((PolOrder>=0) && (PolOrder<(size+1)*5)) {
307  Int_t nsegments = (PolOrder+1)*2;
308  TArrayC usage(size); usage.Reset(0);
309  for (Int_t nseg=0; nseg<nsegments; nseg++) {
310  Int_t lbound = Int_t((nseg+0.)/nsegments*size);
311  Int_t rbound = Int_t((nseg+1.)/nsegments*size-1.);
312  Int_t nselect = (rbound-lbound)/10;
313  if (nselect<1) nselect = 1;
314  Int_t nsel = 0;
315  while (nsel++<nselect) {
316  Int_t pmin=-1;
317  for(Int_t i=lbound;i<=rbound;i++)
318  if(!usage[i])
319  if ((pmin<0) || (bins[i]<bins[pmin])) pmin = i;
320  usage[pmin] = 1;
321  }
322  }
323 
324  TArrayD Coef(PolOrder+1);
325  DefinePolynom(size, bins, scales, Coef, weights, 0, usage.GetArray());
326 
327  for(Int_t i=0;i<size;i++)
328  backgr[i] = CalcPolynom(Coef, scales[i]);
329 
330  fitter->AddPolynomX(data->GetName(), "Pol", Coef);
331  }
332 
333  Double_t dmax = 0.;
334  for(Int_t i = 0; i < size; i++) {
335  Double_t zn = bins[i]-backgr[i];
336  if(zn > dmax) dmax = zn;
337  }
338  AmplThreshold = AmplThreshold*dmax;
339 
340  TArrayC usage(size); usage.Reset(1);
341  Char_t* use = usage.GetArray();
342 
343  Bool_t validline = kFALSE;
344 
345  do {
346 
347  validline = kFALSE;
348 
349  Int_t pmax = -1;
350  Double_t max = 0.;
351  for(Int_t i=0;i<size;i++)
352  if (use[i]) {
353  Double_t zn = bins[i]-backgr[i];
354  if((pmax<0) || (zn>max)) { pmax = i; max = zn; }
355  }
356 
357  Double_t sigma = TMath::Sqrt(1./weights[pmax]);
358 
359  if ((max<2.*sigma) || (max<AmplThreshold)) break;
360 
361  Int_t lbound = pmax, rbound = pmax;
362 
363  Bool_t lok = FindValue(bins, backgr, use, size, lbound, -1, max*0.5);
364  Bool_t rok = FindValue(bins, backgr, use, size, rbound, +1, max*0.5);
365 
366  if(lok || rok) {
367 
368  validline = kTRUE;
369 
370  Int_t w = (lok && rok) ? (rbound-lbound)/4 : ( lok ? (pmax-lbound)/2 : (rbound-pmax)/2 ) ;
371  if (w<2) w=2;
372 
373  TArrayD Coef1(3), Coef2(2);
374  DefinePolynomEx(size, bins, scales, weights, backgr, pmax-w, pmax+w, Coef1);
375  Double_t middle = -0.5*Coef1[1]/Coef1[2];
376  Double_t amplitude = CalcPolynom(Coef1, middle);
377 
378  Double_t left = 0, right = 0;
379  if (lok) {
380  DefinePolynomEx(size, bins, scales, weights, backgr, lbound-w, lbound+w, Coef2);
381  left = (amplitude*.5-Coef2[0])/Coef2[1];
382  }
383 
384  if (rok) {
385  DefinePolynomEx(size, bins, scales, weights, backgr, rbound-w, rbound+w, Coef2);
386  right = (amplitude*.5-Coef2[0])/Coef2[1];
387  }
388 
389  if (!lok) left = 2*middle-right; else
390  if (!rok) right = 2*middle-left;
391 
392  Double_t width = (right-left)*0.5;
393 
394  lbound = pmax-4*w; if (lbound<0) lbound=0;
395  rbound = pmax+4*w; if (rbound>=size) rbound=size-1;
396  for(Int_t i=lbound;i<=rbound;i++) use[i]=0;
397 
398  if(lok && rok && (width>=MinWidth) && (width<=MaxWidth) && (amplitude>AmplThreshold) && (amplitude>2*sigma))
399  fitter->AddGauss1(data->GetName(), fitter->FindNextName("Gauss",0), middle, width, amplitude);
400  }
401 
402  } while(validline);
403 }
404 
405 #include "f_find_peaks.c"
406 
408  TGo4FitData* data,
409  Int_t MaxNumPeaks,
410  Int_t ChannelSum,
411  Double_t NoiseFactor,
412  Double_t NoiseMinimum,
413  Int_t MinimasOrder) {
414  if ((fitter==0) || (data==0) || (MaxNumPeaks<1)) return;
415 
416  TGo4FitDataIter* iter = data->MakeIter();
417  if (iter==0) return;
418 
419  Int_t size = iter->CountPoints();
420 
421  if ((size<10) || (iter->ScalesSize()!=1)) { delete iter; return; }
422 
423  TArrayD Bins(size), Scales(size);
424  Double_t* bins = Bins.GetArray();
425  Double_t* scales = Scales.GetArray();
426 
427  Int_t nbin = 0;
428  if (iter->Reset()) do {
429  bins[nbin] = data->GetAmplValue() * iter->Value();
430  scales[nbin] = iter->x();
431  nbin++;
432  } while (iter->Next());
433 
434  delete iter;
435 
436 
437  int NumPeaks = 0;
438  TArrayI Minimas(MaxNumPeaks), Maximas(MaxNumPeaks), MaximaWidths(MaxNumPeaks);
439  TArrayI PeaksLeft(MaxNumPeaks), PeaksRight(MaxNumPeaks), PeaksMaximum(MaxNumPeaks);
440  TArrayD MinimaValues(MaxNumPeaks), MinimaScales(MaxNumPeaks), MaximaIntegrals(MaxNumPeaks);
441 
442  go4fit_find_peaks(bins,
443  TYPE__DOUBLE,
444  0,
445  size,
446  ChannelSum,
447  NoiseFactor,
448  NoiseMinimum,
449  0,
450  MaxNumPeaks,
451  &NumPeaks,
452  Minimas.GetArray(),
453  MinimaValues.GetArray(),
454  Maximas.GetArray(),
455  MaximaWidths.GetArray(),
456  MaximaIntegrals.GetArray(),
457  PeaksLeft.GetArray(),
458  PeaksMaximum.GetArray(),
459  PeaksRight.GetArray(),
460  0,
461  0);
462 
463  if (NumPeaks<=0) return;
464 
465  Int_t NumMinimas = (NumPeaks<MaxNumPeaks) ? NumPeaks+1 : MaxNumPeaks;
466  for(int n=0; n<NumMinimas; n++) {
467  MinimaScales[n] = scales[Minimas[n]];
468  }
469 
470  if ((MinimasOrder>=0) && (NumMinimas>1)) {
471  Int_t order = MinimasOrder;
472  if (order>=NumMinimas) order = NumMinimas-1;
473  TArrayD Coef(order+1);
474  DefinePolynom(NumMinimas, MinimaValues.GetArray(), MinimaScales.GetArray(), Coef);
475  fitter->AddPolynomX(data->GetName(), "Pol", Coef);
476  }
477 
478  for(int n=0;n<NumPeaks;n++) {
479  Double_t pos = scales[PeaksMaximum[n]];
480  Double_t width = 0.;
481  if (PeaksRight[n]-PeaksLeft[n]<1)
482  width = scales[PeaksMaximum[n]] - scales[PeaksMaximum[n]-1];
483  else
484  width = TMath::Abs(scales[PeaksRight[n]] - scales[PeaksLeft[n]])/2.;
485  Double_t ampl = bins[PeaksMaximum[n]];
486 
487  fitter->AddGauss1(data->GetName(), fitter->FindNextName("Gauss",0), pos, width, ampl);
488  }
489 }
490 
491 
492 
static void DefinePolynom(Int_t size, Double_t *bins, Double_t *scales, TArrayD &Coef, Double_t *weight=0, Double_t *backgr=0, Char_t *use=0)
void SetupForSecond(Double_t LineWidth)
Bool_t HasWidths() const
Definition: TGo4FitData.h:499
Int_t GetPeakFinderType() const
Double_t x() const
Definition: TGo4FitData.h:484
void Set2NoiseFactor(Double_t factor)
virtual TGo4FitDataIter * MakeIter()
Definition: TGo4FitData.h:165
Double_t Get0MinWidth() const
static void DefinePolynomEx(Int_t size, Double_t *bins, Double_t *scales, Double_t *weight, Double_t *backgr, Int_t lbound, Int_t rbound, TArrayD &Coef)
void AddPolynomX(const char *DataName, const char *NamePrefix, Int_t MaxOrder=1, Int_t GroupIndex=0, Double_t lrange=0., Double_t rrange=0.)
Definition: TGo4Fitter.cxx:223
void SergeyLinevPeakFinder(TGo4Fitter *fitter, TGo4FitData *data, Int_t PolOrder, Double_t AmplThreshold, Double_t MinWidth, Double_t MaxWidth)
void go4fit_find_peaks(void *pfData, int lType, int lFirstChan, int lPoints, int lAver, double dDeltaFact, double dDeltaMin, double *dNoise, int lPeaksMax, int *plPeaks, int *plMinima, double *pdMinima, int *plMean, int *plWidth, double *pdIntegral, int *plLeft, int *plMaximum, int *plRight, float *pfLowNoise, float *pfHighNoise)
Definition: f_find_peaks.c:28
void Set2NoiseMinimum(Double_t min)
void SetPeakFinderType(Int_t typ)
Bool_t GetUsePolynom() const
virtual Bool_t Reset(Bool_t UseRanges=kTRUE)
Double_t Get2NoiseMinimum() const
void DeleteModelsAssosiatedTo(const char *DataName)
Definition: TGo4Fitter.cxx:396
const Double_t * Widths() const
Definition: TGo4FitData.h:504
Int_t Get2ChannelSum() const
const char * GetDataName()
virtual void Print(Option_t *option) const
void Set0MaxWidth(Double_t max)
void Set0MaxAmplFactor(Double_t factor)
void SetupForFirst(Double_t MaxAmplFactor, Double_t MinWidth, Double_t MaxWidth)
TGo4FitModelGauss1 * AddGauss1(const char *DataName, const char *ModelName, Double_t iPosition, Double_t iWidth, Double_t iAmpl=1., Int_t Axis=0)
Definition: TGo4Fitter.cxx:345
void Set2ChannelSum(Int_t sum)
void Print(Option_t *option) const
Double_t Get2NoiseFactor() const
Bool_t FindValue(Double_t *bins, Double_t *backgr, Char_t *use, Int_t size, Int_t &p, Int_t dir, Double_t value)
#define TYPE__DOUBLE
Definition: f_find_peaks.c:19
Int_t CountPoints(Bool_t UseRanges=kTRUE)
virtual Bool_t Next(Bool_t UseRanges=kTRUE)
TString FindNextName(const char *Head, Int_t start, Bool_t isModel=kTRUE)
virtual void DoAction(TGo4FitterAbstract *Fitter)
void SetupPolynomialBackground(Int_t PolynomOrder)
Double_t Get1LineWidth() const
void SetUsePolynom(Bool_t use)
void SetupForThird(Double_t NoiseFactor, Double_t NoiseMinimum, Int_t ChannelSum)
void usage(const char *subtopic=0)
Double_t Get0MaxWidth() const
void Set1LineWidth(Double_t width)
Double_t Value() const
Definition: TGo4FitData.h:514
static Double_t CalcPolynom(const TArrayD &Coef, Double_t x)
Double_t Get0MaxAmplFactor() const
void ROOTPeakFinder(TGo4Fitter *fitter, TGo4FitData *data, Int_t PolynomOrder, Double_t Sigma)
TGo4FitData * FindData(const char *DataName)
Definition: TGo4Fitter.cxx:111
void SetPolynomOrder(Int_t order)
Double_t StandardDeviation() const
Definition: TGo4FitData.h:519
void HansEsselPeakFinder(TGo4Fitter *fitter, TGo4FitData *data, Int_t MaxNumPeaks=50, Int_t ChannelSum=1, Double_t NoiseFactor=2., Double_t NoiseMinimum=10., Int_t MinimasOrder=-1)
void Set0MinWidth(Double_t min)
Int_t ScalesSize() const
Definition: TGo4FitData.h:474
Int_t GetPolynomOrder() const