Go4Fit/TGo4FitPeakFinder.cxx
Go to the documentation of this file.
00001 // $Id: TGo4FitPeakFinder.cxx 478 2009-10-29 12:26:09Z linev $ 00002 //----------------------------------------------------------------------- 00003 // The GSI Online Offline Object Oriented (Go4) Project 00004 // Experiment Data Processing at EE department, GSI 00005 //----------------------------------------------------------------------- 00006 // Copyright (C) 2000- GSI Helmholtzzentrum für Schwerionenforschung GmbH 00007 // Planckstr. 1, 64291 Darmstadt, Germany 00008 // Contact: http://go4.gsi.de 00009 //----------------------------------------------------------------------- 00010 // This software can be used under the license agreements as stated 00011 // in Go4License.txt file which is part of the distribution. 00012 //----------------------------------------------------------------------- 00013 00014 #include "TGo4FitPeakFinder.h" 00015 00016 #include "TMath.h" 00017 #include "TH1D.h" 00018 #include "TMatrixD.h" 00019 #include "TVectorD.h" 00020 #include "TArrayD.h" 00021 #include "TArrayI.h" 00022 #include "TSpectrum.h" 00023 00024 #include "TGo4Fitter.h" 00025 #include "TGo4FitData.h" 00026 #include "TGo4FitModelPolynom.h" 00027 00028 TGo4FitPeakFinder::TGo4FitPeakFinder() : TGo4FitterAction() { 00029 00030 } 00031 00032 TGo4FitPeakFinder::TGo4FitPeakFinder(const char* Name, const char* DataName, Bool_t ClearModels, Int_t PolOrder) : 00033 TGo4FitterAction(Name, "Peak finder action"), 00034 fiPeakFinderType(0), fxDataName(DataName), fbClearModels(ClearModels), 00035 fbUsePolynom(PolOrder>=0), fiPolynomOrder(PolOrder>=0 ? PolOrder : 0), 00036 fd0MinWidth(5.), fd0MaxWidth(30.), fd0MaxAmplFactor(0.2), 00037 fd1LineWidth(5.), 00038 fd2NoiseFactor(2.), fd2NoiseMinimum(5.), fi2ChannelSum(2) { 00039 } 00040 00041 TGo4FitPeakFinder::~TGo4FitPeakFinder() { 00042 } 00043 00044 void TGo4FitPeakFinder::SetupPolynomialBackground(Int_t PolynomOrder) { 00045 if (PolynomOrder<0) SetUsePolynom(kFALSE); 00046 else { 00047 SetUsePolynom(kTRUE); 00048 SetPolynomOrder(PolynomOrder); 00049 } 00050 } 00051 00052 void TGo4FitPeakFinder::SetupForFirst(Double_t MaxAmplFactor, Double_t MinWidth, Double_t MaxWidth) { 00053 SetPeakFinderType(0); 00054 Set0MaxAmplFactor(MaxAmplFactor); 00055 Set0MinWidth(MinWidth); 00056 Set0MaxWidth(MaxWidth); 00057 } 00058 00059 void TGo4FitPeakFinder::SetupForSecond(Double_t LineWidth) { 00060 SetPeakFinderType(1); 00061 Set1LineWidth(LineWidth); 00062 } 00063 00064 void TGo4FitPeakFinder::SetupForThird(Double_t NoiseFactor, Double_t NoiseMinimum, Int_t ChannelSum) { 00065 SetPeakFinderType(2); 00066 Set2NoiseFactor(NoiseFactor); 00067 Set2NoiseMinimum(NoiseMinimum); 00068 Set2ChannelSum(ChannelSum); 00069 } 00070 00071 void TGo4FitPeakFinder::DoAction(TGo4FitterAbstract* Fitter) { 00072 TGo4Fitter* fitter = dynamic_cast<TGo4Fitter*> (Fitter); 00073 if (fitter==0) return; 00074 00075 TGo4FitData* data = fitter->FindData(GetDataName()); 00076 if (data==0) return; 00077 00078 if (GetClearModels()) 00079 fitter->DeleteModelsAssosiatedTo(data->GetName()); 00080 00081 if (GetPeakFinderType()==0) 00082 SergeyLinevPeakFinder(fitter, 00083 data, 00084 GetUsePolynom() ? GetPolynomOrder() : -1, 00085 Get0MaxAmplFactor(), 00086 Get0MinWidth(), 00087 Get0MaxWidth()); 00088 00089 if (GetPeakFinderType()==1) 00090 ROOTPeakFinder(fitter, 00091 data, 00092 GetUsePolynom() ? GetPolynomOrder() : -1, 00093 Get1LineWidth()); 00094 00095 if (GetPeakFinderType()==2) 00096 HansEsselPeakFinder(fitter, 00097 data, 00098 500, 00099 Get2ChannelSum(), 00100 Get2NoiseFactor(), 00101 Get2NoiseMinimum(), 00102 GetUsePolynom() ? GetPolynomOrder() : -1); 00103 } 00104 00105 void TGo4FitPeakFinder::Print(Option_t* option) const { 00106 TGo4FitterAction::Print(option); 00107 } 00108 /* 00109 void TGo4FitPeakFinder::ROOTPeakFinder(TGo4Fitter* fitter, TGo4FitData* data, Int_t PolynomOrder, Double_t Sigma) { 00110 if ((fitter==0) || (data==0)) return; 00111 00112 TGo4FitDataIter* iter = data->MakeIter(); 00113 if (iter==0) return; 00114 00115 Int_t size = iter->CountPoints(); 00116 00117 if ((size<10) || (iter->ScalesSize()!=1)) { delete iter; return; } 00118 00119 TArrayF Bins(size), Scales(size); 00120 00121 Int_t nbin = 0; 00122 if (iter->Reset()) do { 00123 Bins[nbin] = data->GetAmplValue() * iter->Value(); 00124 Scales[nbin] = iter->x(); 00125 nbin++; 00126 } while (iter->Next()); 00127 00128 delete iter; 00129 00130 if (PolynomOrder>=0) 00131 fitter->AddPolynomX(data->GetName(), "Pol", PolynomOrder); 00132 00133 TSpectrum sp(100); 00134 sp.Search1(Bins.GetArray(), size, Sigma); 00135 00136 for(Int_t n=0;n<sp.GetNPeaks();n++) { 00137 Double_t dindx = sp.GetPositionX()[n]; 00138 Int_t left = TMath::Nint(dindx-Sigma/2.); 00139 Int_t right = TMath::Nint(dindx+Sigma/2.); 00140 00141 Double_t pos = Scales[TMath::Nint(dindx)]; 00142 Double_t width = TMath::Abs(Scales[right]-Scales[left]); 00143 00144 fitter->AddGauss1(data->GetName(), fitter->FindNextName("Gauss",n), pos, width); 00145 } 00146 } 00147 */ 00148 00149 void TGo4FitPeakFinder::ROOTPeakFinder(TGo4Fitter* fitter, TGo4FitData* data, Int_t PolynomOrder, Double_t Sigma) { 00150 if ((fitter==0) || (data==0)) return; 00151 00152 TGo4FitDataIter* iter = data->MakeIter(); 00153 if (iter==0) return; 00154 00155 Int_t size = iter->CountPoints(); 00156 00157 if ((size<10) || (iter->ScalesSize()!=1)) { delete iter; return; } 00158 00159 TArrayD Bins(size), Scales(size), HScales(size+1); 00160 00161 Int_t nbin = 0; 00162 if (iter->Reset()) do { 00163 Bins[nbin] = data->GetAmplValue() * iter->Value(); 00164 Double_t pos = iter->x(); 00165 Double_t width = iter->HasWidths() ? iter->Widths()[0] : 1.; 00166 Scales[nbin] = pos; 00167 HScales[nbin] = pos - width/2.; 00168 HScales[nbin+1] = pos + width/2; 00169 nbin++; 00170 } while (iter->Next()); 00171 00172 delete iter; 00173 00174 if (PolynomOrder>=0) 00175 fitter->AddPolynomX(data->GetName(), "Pol", PolynomOrder); 00176 00177 TH1D histo("ROOTPeakFinder", "ROOTPeakFinder spectrum", size, HScales.GetArray()); 00178 histo.SetDirectory(0); 00179 for (Int_t i=0;i<size;i++) 00180 histo.SetBinContent(i+1, Bins[i]); 00181 00182 TSpectrum sp(100); 00183 sp.Search(&histo, Sigma); 00184 // sp.Search1(Bins.GetArray(), size, Sigma); 00185 00186 for(Int_t n=0;n<sp.GetNPeaks();n++) { 00187 Double_t pos = sp.GetPositionX()[n]; 00188 Double_t width = Sigma; 00189 00190 00191 // this is old code for ROOT 3.*, may be usefull for somebody 00192 /* 00193 Double_t dindx = sp.GetPositionX()[n]; 00194 Int_t left = TMath::Nint(dindx-Sigma/2.); 00195 Int_t right = TMath::Nint(dindx+Sigma/2.); 00196 if (left<0) left = 0; 00197 if (right>=size) right = size-1; 00198 00199 Double_t pos = Scales[TMath::Nint(dindx)]; 00200 Double_t width = TMath::Abs(Scales[right]-Scales[left]); 00201 */ 00202 fitter->AddGauss1(data->GetName(), fitter->FindNextName("Gauss",n), pos, width); 00203 } 00204 } 00205 00206 00207 Double_t TGo4FitPeakFinder::CalcPolynom(const TArrayD& Coef, Double_t x) { 00208 Double_t power = 1., sum = 0.; 00209 for(Int_t n=0;n<Coef.GetSize();n++) 00210 { sum+=Coef[n]*power; power*=x; } 00211 return sum; 00212 } 00213 00214 void TGo4FitPeakFinder::DefinePolynom(Int_t size, Double_t* bins, Double_t* scales, TArrayD& Coef, 00215 Double_t* weight, Double_t* backgr, Char_t* use) { 00216 Int_t maxorder = Coef.GetSize()-1; 00217 00218 TMatrixD matr(maxorder+1,maxorder+1); 00219 TVectorD vector(maxorder+1); 00220 00221 for(Int_t order1=0;order1<=maxorder;order1++) 00222 for(Int_t order2=order1;order2<=maxorder;order2++) { 00223 Double_t sum = 0; 00224 for(Int_t i=0;i<size;i++) 00225 if((use==0) || use[i]) { 00226 Double_t zn = TMath::Power(scales[i],order1+order2); 00227 if (weight) zn*=weight[i]; 00228 sum+=zn; 00229 } 00230 matr(order1,order2) = sum; 00231 matr(order2,order1) = sum; 00232 } 00233 00234 for(Int_t order1=0;order1<=maxorder;order1++) { 00235 Double_t sum = 0; 00236 for(Int_t i=0;i<size;i++) 00237 if((use==0) || use[i]) { 00238 Double_t zn = bins[i]; 00239 if (backgr) zn-=backgr[i]; 00240 zn*=TMath::Power(scales[i],order1); 00241 if (weight) zn*=weight[i]; 00242 sum+=zn; 00243 } 00244 vector(order1) = sum; 00245 } 00246 00247 matr.Invert(); 00248 vector*=matr; 00249 00250 for(Int_t order=0;order<=maxorder;order++) 00251 Coef[order]=vector(order); 00252 } 00253 00254 void TGo4FitPeakFinder::DefinePolynomEx(Int_t size, Double_t* bins, Double_t* scales, Double_t* weight, Double_t* backgr, 00255 Int_t lbound, Int_t rbound, TArrayD& Coef) { 00256 TArrayC use(size); use.Reset(0); 00257 if (lbound<0) lbound = 0; 00258 if (rbound>=size) rbound=size-1; 00259 for(Int_t i=lbound;i<=rbound;i++) use[i]=1; 00260 DefinePolynom(size, bins, scales, Coef, weight, backgr, use.GetArray()); 00261 } 00262 00263 Bool_t FindValue(Double_t* bins, Double_t* backgr, Char_t* use, Int_t size, Int_t& p, Int_t dir, Double_t value) { 00264 p += dir; 00265 while ((p>=0) && (p<size) && use[p]) { 00266 Double_t zn = bins[p]-backgr[p]; 00267 if (zn<value) return kTRUE; 00268 p+=dir; 00269 } 00270 return kFALSE; 00271 } 00272 00273 void TGo4FitPeakFinder::SergeyLinevPeakFinder(TGo4Fitter* fitter, 00274 TGo4FitData* data, 00275 Int_t PolOrder, 00276 Double_t AmplThreshold, 00277 Double_t MinWidth, 00278 Double_t MaxWidth) { 00279 if ((fitter==0) || (data==0)) return; 00280 00281 TGo4FitDataIter* iter = data->MakeIter(); 00282 if (iter==0) return; 00283 00284 Int_t size = iter->CountPoints(); 00285 00286 if ((size<10) || (iter->ScalesSize()!=1)) { delete iter; return; } 00287 00288 TArrayD Bins(size), Scales(size), Weights(size), Background(size); 00289 00290 Weights.Reset(1.); Background.Reset(0.); 00291 00292 Double_t* bins = Bins.GetArray(); 00293 Double_t* scales = Scales.GetArray(); 00294 Double_t* weights = Weights.GetArray(); 00295 Double_t* backgr = Background.GetArray(); 00296 00297 Int_t nbin = 0; 00298 if (iter->Reset()) do { 00299 bins[nbin] = data->GetAmplValue() * iter->Value(); 00300 if (iter->StandardDeviation()>0) weights[nbin] = 1./iter->StandardDeviation(); 00301 scales[nbin] = iter->x(); 00302 nbin++; 00303 } while (iter->Next()); 00304 00305 delete iter; 00306 00307 if ((PolOrder>=0) && (PolOrder<(size+1)*5)) { 00308 Int_t nsegments = (PolOrder+1)*2; 00309 TArrayC usage(size); usage.Reset(0); 00310 for (Int_t nseg=0; nseg<nsegments; nseg++) { 00311 Int_t lbound = Int_t((nseg+0.)/nsegments*size); 00312 Int_t rbound = Int_t((nseg+1.)/nsegments*size-1.); 00313 Int_t nselect = (rbound-lbound)/10; 00314 if (nselect<1) nselect = 1; 00315 Int_t nsel = 0; 00316 while (nsel++<nselect) { 00317 Int_t pmin=-1; 00318 for(Int_t i=lbound;i<=rbound;i++) 00319 if(!usage[i]) 00320 if ((pmin<0) || (bins[i]<bins[pmin])) pmin = i; 00321 usage[pmin] = 1; 00322 } 00323 } 00324 00325 TArrayD Coef(PolOrder+1); 00326 DefinePolynom(size, bins, scales, Coef, weights, 0, usage.GetArray()); 00327 00328 for(Int_t i=0;i<size;i++) 00329 backgr[i] = CalcPolynom(Coef, scales[i]); 00330 00331 fitter->AddPolynomX(data->GetName(), "Pol", Coef); 00332 } 00333 00334 Double_t max=0.; 00335 for(Int_t i=0;i<size;i++) { 00336 Double_t zn = bins[i]-backgr[i]; 00337 if(zn>max) max = zn; 00338 } 00339 AmplThreshold = AmplThreshold*max; 00340 00341 TArrayC usage(size); usage.Reset(1); 00342 Char_t* use = usage.GetArray(); 00343 00344 Bool_t validline = kFALSE; 00345 00346 do { 00347 00348 validline = kFALSE; 00349 00350 Int_t pmax=-1; Double_t max=0.; 00351 for(Int_t i=0;i<size;i++) 00352 if (use[i]) { 00353 Double_t zn = bins[i]-backgr[i]; 00354 if((pmax<0) || (zn>max)) { pmax = i; max = zn; } 00355 } 00356 00357 Double_t sigma = TMath::Sqrt(1./weights[pmax]); 00358 00359 if ((max<2.*sigma) || (max<AmplThreshold)) break; 00360 00361 Int_t lbound = pmax, rbound = pmax; 00362 00363 Bool_t lok = FindValue(bins, backgr, use, size, lbound, -1, max*0.5); 00364 Bool_t rok = FindValue(bins, backgr, use, size, rbound, +1, max*0.5); 00365 00366 if(lok || rok) { 00367 00368 validline = kTRUE; 00369 00370 Int_t w = (lok && rok) ? (rbound-lbound)/4 : ( lok ? (pmax-lbound)/2 : (rbound-pmax)/2 ) ; 00371 if (w<2) w=2; 00372 00373 TArrayD Coef1(3), Coef2(2); 00374 DefinePolynomEx(size, bins, scales, weights, backgr, pmax-w, pmax+w, Coef1); 00375 Double_t middle = -0.5*Coef1[1]/Coef1[2]; 00376 Double_t amplitude = CalcPolynom(Coef1, middle); 00377 00378 Double_t left = 0, right = 0; 00379 if (lok) { 00380 DefinePolynomEx(size, bins, scales, weights, backgr, lbound-w, lbound+w, Coef2); 00381 left = (amplitude*.5-Coef2[0])/Coef2[1]; 00382 } 00383 00384 if (rok) { 00385 DefinePolynomEx(size, bins, scales, weights, backgr, rbound-w, rbound+w, Coef2); 00386 right = (amplitude*.5-Coef2[0])/Coef2[1]; 00387 } 00388 00389 if (!lok) left = 2*middle-right; else 00390 if (!rok) right = 2*middle-left; 00391 00392 Double_t width = (right-left)*0.5; 00393 00394 lbound = pmax-4*w; if (lbound<0) lbound=0; 00395 rbound = pmax+4*w; if (rbound>=size) rbound=size-1; 00396 for(Int_t i=lbound;i<=rbound;i++) use[i]=0; 00397 00398 if(lok && rok && (width>=MinWidth) && (width<=MaxWidth) && (amplitude>AmplThreshold) && (amplitude>2*sigma)) 00399 fitter->AddGauss1(data->GetName(), fitter->FindNextName("Gauss",0), middle, width, amplitude); 00400 } 00401 00402 } while(validline); 00403 } 00404 00405 #include "f_find_peaks.c" 00406 00407 void TGo4FitPeakFinder::HansEsselPeakFinder(TGo4Fitter* fitter, 00408 TGo4FitData* data, 00409 Int_t MaxNumPeaks, 00410 Int_t ChannelSum, 00411 Double_t NoiseFactor, 00412 Double_t NoiseMinimum, 00413 Int_t MinimasOrder) { 00414 if ((fitter==0) || (data==0) || (MaxNumPeaks<1)) return; 00415 00416 TGo4FitDataIter* iter = data->MakeIter(); 00417 if (iter==0) return; 00418 00419 Int_t size = iter->CountPoints(); 00420 00421 if ((size<10) || (iter->ScalesSize()!=1)) { delete iter; return; } 00422 00423 TArrayD Bins(size), Scales(size); 00424 Double_t* bins = Bins.GetArray(); 00425 Double_t* scales = Scales.GetArray(); 00426 00427 Int_t nbin = 0; 00428 if (iter->Reset()) do { 00429 bins[nbin] = data->GetAmplValue() * iter->Value(); 00430 scales[nbin] = iter->x(); 00431 nbin++; 00432 } while (iter->Next()); 00433 00434 delete iter; 00435 00436 00437 int NumPeaks = 0; 00438 TArrayI Minimas(MaxNumPeaks), Maximas(MaxNumPeaks), MaximaWidths(MaxNumPeaks); 00439 TArrayI PeaksLeft(MaxNumPeaks), PeaksRight(MaxNumPeaks), PeaksMaximum(MaxNumPeaks); 00440 TArrayD MinimaValues(MaxNumPeaks), MinimaScales(MaxNumPeaks), MaximaIntegrals(MaxNumPeaks); 00441 00442 go4fit_find_peaks(bins, 00443 TYPE__DOUBLE, 00444 0, 00445 size, 00446 ChannelSum, 00447 NoiseFactor, 00448 NoiseMinimum, 00449 0, 00450 MaxNumPeaks, 00451 &NumPeaks, 00452 Minimas.GetArray(), 00453 MinimaValues.GetArray(), 00454 Maximas.GetArray(), 00455 MaximaWidths.GetArray(), 00456 MaximaIntegrals.GetArray(), 00457 PeaksLeft.GetArray(), 00458 PeaksMaximum.GetArray(), 00459 PeaksRight.GetArray(), 00460 0, 00461 0); 00462 00463 if (NumPeaks<=0) return; 00464 00465 Int_t NumMinimas = (NumPeaks<MaxNumPeaks) ? NumPeaks+1 : MaxNumPeaks; 00466 for(int n=0; n<NumMinimas; n++) { 00467 MinimaScales[n] = scales[Minimas[n]]; 00468 } 00469 00470 if ((MinimasOrder>=0) && (NumMinimas>1)) { 00471 Int_t order = MinimasOrder; 00472 if (order>=NumMinimas) order = NumMinimas-1; 00473 TArrayD Coef(order+1); 00474 DefinePolynom(NumMinimas, MinimaValues.GetArray(), MinimaScales.GetArray(), Coef); 00475 fitter->AddPolynomX(data->GetName(), "Pol", Coef); 00476 } 00477 00478 for(int n=0;n<NumPeaks;n++) { 00479 Double_t pos = scales[PeaksMaximum[n]]; 00480 Double_t width = 0.; 00481 if (PeaksRight[n]-PeaksLeft[n]<1) 00482 width = scales[PeaksMaximum[n]] - scales[PeaksMaximum[n]-1]; 00483 else 00484 width = TMath::Abs(scales[PeaksRight[n]] - scales[PeaksLeft[n]])/2.; 00485 Double_t ampl = bins[PeaksMaximum[n]]; 00486 00487 fitter->AddGauss1(data->GetName(), fitter->FindNextName("Gauss",0), pos, width, ampl); 00488 } 00489 } 00490 00491 00492
