/Go4Fit/TGo4FitPeakFinder.cxx Source File

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