/// \file /// \ingroup tutorial_spectrum /// \notebook /// Illustrates how to find peaks in histograms. /// /// This script generates a random number of gaussian peaks /// on top of a linear background. /// The position of the peaks is found via TSpectrum and injected /// as initial values of parameters to make a global fit. /// The background is computed and drawn on top of the original histogram. /// /// This script can fit "peaks' heights" or "peaks' areas" (comment out /// or uncomment the line which defines `__PEAKS_C_FIT_AREAS__`). /// /// To execute this example, do (in ROOT 5 or ROOT 6): /// /// ~~~{.cpp} /// root > .x peaks.C (generate 10 peaks by default) /// root > .x peaks.C++ (use the compiler) /// root > .x peaks.C++(30) (generates 30 peaks) /// ~~~ /// /// To execute only the first part of the script (without fitting) /// specify a negative value for the number of peaks, eg /// /// ~~~{.cpp} /// root > .x peaks.C(-20) /// ~~~ /// /// \macro_output /// \macro_image /// \macro_code /// /// \author Rene Brun #include "TCanvas.h" #include "TMath.h" #include "TH1.h" #include "TF1.h" #include "TRandom.h" #include "TSpectrum.h" #include "TVirtualFitter.h" // // Comment out the line below, if you want "peaks' heights". // Uncomment the line below, if you want "peaks' areas". // // #define __PEAKS_C_FIT_AREAS__ 1 /* fit peaks' areas */ Int_t npeaks = 30; Double_t fpeaks(Double_t *x, Double_t *par) { Double_t result = par[0] + par[1]*x[0]; for (Int_t p=0;pRndm()*980; // "mean" par[3*p+4] = 3+2*gRandom->Rndm(); // "sigma" #if defined(__PEAKS_C_FIT_AREAS__) par[3*p+2] *= par[3*p+4] * (TMath::Sqrt(TMath::TwoPi())); // "area" #endif /* defined(__PEAKS_C_FIT_AREAS__) */ } TF1 *f = new TF1("f",fpeaks,0,1000,2+3*npeaks); f->SetNpx(1000); f->SetParameters(par); TCanvas *c1 = new TCanvas("c1","c1",10,10,1000,900); c1->Divide(1,2); c1->cd(1); h->FillRandom("f",200000); h->Draw(); TH1F *h2 = (TH1F*)h->Clone("h2"); // Use TSpectrum to find the peak candidates TSpectrum *s = new TSpectrum(2*npeaks); Int_t nfound = s->Search(h,2,"",0.10); printf("Found %d candidate peaks to fit\n",nfound); // Estimate background using TSpectrum::Background TH1 *hb = s->Background(h,20,"same"); if (hb) c1->Update(); if (np <0) return; //estimate linear background using a fitting method c1->cd(2); TF1 *fline = new TF1("fline","pol1",0,1000); h->Fit("fline","qn"); // Loop on all found peaks. Eliminate peaks at the background level par[0] = fline->GetParameter(0); par[1] = fline->GetParameter(1); npeaks = 0; Double_t *xpeaks; xpeaks = s->GetPositionX(); for (p=0;pGetXaxis()->FindBin(xp); Double_t yp = h->GetBinContent(bin); if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue; par[3*npeaks+2] = yp; // "height" par[3*npeaks+3] = xp; // "mean" par[3*npeaks+4] = 3; // "sigma" #if defined(__PEAKS_C_FIT_AREAS__) par[3*npeaks+2] *= par[3*npeaks+4] * (TMath::Sqrt(TMath::TwoPi())); // "area" #endif /* defined(__PEAKS_C_FIT_AREAS__) */ npeaks++; } printf("Found %d useful peaks to fit\n",npeaks); printf("Now fitting: Be patient\n"); TF1 *fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks); // We may have more than the default 25 parameters TVirtualFitter::Fitter(h2,10+3*npeaks); fit->SetParameters(par); fit->SetNpx(1000); h2->Fit("fit"); }