Logo ROOT   6.16/01
Reference Guide
peaks.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_spectrum
3/// \notebook
4/// Illustrates how to find peaks in histograms.
5///
6/// This script generates a random number of gaussian peaks
7/// on top of a linear background.
8/// The position of the peaks is found via TSpectrum and injected
9/// as initial values of parameters to make a global fit.
10/// The background is computed and drawn on top of the original histogram.
11///
12/// This script can fit "peaks' heights" or "peaks' areas" (comment out
13/// or uncomment the line which defines `__PEAKS_C_FIT_AREAS__`).
14///
15/// To execute this example, do (in ROOT 5 or ROOT 6):
16///
17/// ~~~{.cpp}
18/// root > .x peaks.C (generate 10 peaks by default)
19/// root > .x peaks.C++ (use the compiler)
20/// root > .x peaks.C++(30) (generates 30 peaks)
21/// ~~~
22///
23/// To execute only the first part of the script (without fitting)
24/// specify a negative value for the number of peaks, eg
25///
26/// ~~~{.cpp}
27/// root > .x peaks.C(-20)
28/// ~~~
29///
30/// \macro_output
31/// \macro_image
32/// \macro_code
33///
34/// \author Rene Brun
35
36#include "TCanvas.h"
37#include "TMath.h"
38#include "TH1.h"
39#include "TF1.h"
40#include "TRandom.h"
41#include "TSpectrum.h"
42#include "TVirtualFitter.h"
43
44//
45// Comment out the line below, if you want "peaks' heights".
46// Uncomment the line below, if you want "peaks' areas".
47//
48// #define __PEAKS_C_FIT_AREAS__ 1 /* fit peaks' areas */
49
50Int_t npeaks = 30;
51Double_t fpeaks(Double_t *x, Double_t *par) {
52 Double_t result = par[0] + par[1]*x[0];
53 for (Int_t p=0;p<npeaks;p++) {
54 Double_t norm = par[3*p+2]; // "height" or "area"
55 Double_t mean = par[3*p+3];
56 Double_t sigma = par[3*p+4];
57#if defined(__PEAKS_C_FIT_AREAS__)
58 norm /= sigma * (TMath::Sqrt(TMath::TwoPi())); // "area"
59#endif /* defined(__PEAKS_C_FIT_AREAS__) */
60 result += norm*TMath::Gaus(x[0],mean,sigma);
61 }
62 return result;
63}
64void peaks(Int_t np=10) {
65 npeaks = TMath::Abs(np);
66 TH1F *h = new TH1F("h","test",500,0,1000);
67 // Generate n peaks at random
68 Double_t par[3000];
69 par[0] = 0.8;
70 par[1] = -0.6/1000;
71 Int_t p;
72 for (p=0;p<npeaks;p++) {
73 par[3*p+2] = 1; // "height"
74 par[3*p+3] = 10+gRandom->Rndm()*980; // "mean"
75 par[3*p+4] = 3+2*gRandom->Rndm(); // "sigma"
76#if defined(__PEAKS_C_FIT_AREAS__)
77 par[3*p+2] *= par[3*p+4] * (TMath::Sqrt(TMath::TwoPi())); // "area"
78#endif /* defined(__PEAKS_C_FIT_AREAS__) */
79 }
80 TF1 *f = new TF1("f",fpeaks,0,1000,2+3*npeaks);
81 f->SetNpx(1000);
82 f->SetParameters(par);
83 TCanvas *c1 = new TCanvas("c1","c1",10,10,1000,900);
84 c1->Divide(1,2);
85 c1->cd(1);
86 h->FillRandom("f",200000);
87 h->Draw();
88 TH1F *h2 = (TH1F*)h->Clone("h2");
89 // Use TSpectrum to find the peak candidates
90 TSpectrum *s = new TSpectrum(2*npeaks);
91 Int_t nfound = s->Search(h,2,"",0.10);
92 printf("Found %d candidate peaks to fit\n",nfound);
93 // Estimate background using TSpectrum::Background
94 TH1 *hb = s->Background(h,20,"same");
95 if (hb) c1->Update();
96 if (np <0) return;
97
98 //estimate linear background using a fitting method
99 c1->cd(2);
100 TF1 *fline = new TF1("fline","pol1",0,1000);
101 h->Fit("fline","qn");
102 // Loop on all found peaks. Eliminate peaks at the background level
103 par[0] = fline->GetParameter(0);
104 par[1] = fline->GetParameter(1);
105 npeaks = 0;
106 Double_t *xpeaks;
107 xpeaks = s->GetPositionX();
108 for (p=0;p<nfound;p++) {
109 Double_t xp = xpeaks[p];
110 Int_t bin = h->GetXaxis()->FindBin(xp);
111 Double_t yp = h->GetBinContent(bin);
112 if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
113 par[3*npeaks+2] = yp; // "height"
114 par[3*npeaks+3] = xp; // "mean"
115 par[3*npeaks+4] = 3; // "sigma"
116#if defined(__PEAKS_C_FIT_AREAS__)
117 par[3*npeaks+2] *= par[3*npeaks+4] * (TMath::Sqrt(TMath::TwoPi())); // "area"
118#endif /* defined(__PEAKS_C_FIT_AREAS__) */
119 npeaks++;
120 }
121 printf("Found %d useful peaks to fit\n",npeaks);
122 printf("Now fitting: Be patient\n");
123 TF1 *fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks);
124 // We may have more than the default 25 parameters
125 TVirtualFitter::Fitter(h2,10+3*npeaks);
126 fit->SetParameters(par);
127 fit->SetNpx(1000);
128 h2->Fit("fit");
129}
#define f(i)
Definition: RSha256.hxx:104
#define h(i)
Definition: RSha256.hxx:106
int Int_t
Definition: RtypesCore.h:41
double Double_t
Definition: RtypesCore.h:55
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
The Canvas class.
Definition: TCanvas.h:31
1-Dim function class
Definition: TF1.h:211
virtual void SetNpx(Int_t npx=100)
Set the number of points used to draw the function.
Definition: TF1.cxx:3426
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:628
virtual Double_t Eval(Double_t x, Double_t y=0, Double_t z=0, Double_t t=0) const
Evaluate this function.
Definition: TF1.cxx:1424
virtual Double_t GetParameter(Int_t ipar) const
Definition: TF1.h:496
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:571
The TH1 histogram class.
Definition: TH1.h:56
virtual TFitResultPtr Fit(const char *formula, Option_t *option="", Option_t *goption="", Double_t xmin=0, Double_t xmax=0)
Fit histogram with function fname.
Definition: TH1.cxx:3695
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom.cxx:533
Advanced Spectra Processing.
Definition: TSpectrum.h:18
static TVirtualFitter * Fitter(TObject *obj, Int_t maxpar=25)
Static function returning a pointer to the current fitter.
const Double_t sigma
return c1
Definition: legend1.C:41
Double_t x[n]
Definition: legend1.C:17
static constexpr double s
Double_t Gaus(Double_t x, Double_t mean=0, Double_t sigma=1, Bool_t norm=kFALSE)
Calculate a gaussian function with mean and sigma.
Definition: TMath.cxx:448
Double_t Sqrt(Double_t x)
Definition: TMath.h:679
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
constexpr Double_t TwoPi()
Definition: TMath.h:45