ROOT logo

From $ROOTSYS/tutorials/spectrum/peaks.C

// 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.
//
// To execute this example, do
//  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
//  root > .x peaks.C(-20)
//
//Author: Rene Brun

#include "TCanvas.h"
#include "TMath.h"
#include "TH1.h"
#include "TF1.h"
#include "TRandom.h"
#include "TSpectrum.h"
#include "TVirtualFitter.h"
   
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;p<npeaks;p++) {
      Double_t norm  = par[3*p+2];
      Double_t mean  = par[3*p+3];
      Double_t sigma = par[3*p+4];
      result += norm*TMath::Gaus(x[0],mean,sigma);
   }
   return result;
}
void peaks(Int_t np=10) {
   npeaks = TMath::Abs(np);
   TH1F *h = new TH1F("h","test",500,0,1000);
   //generate n peaks at random
   Double_t par[3000];
   par[0] = 0.8;
   par[1] = -0.6/1000;
   Int_t p;
   for (p=0;p<npeaks;p++) {
      par[3*p+2] = 1;
      par[3*p+3] = 10+gRandom->Rndm()*980;
      par[3*p+4] = 3+2*gRandom->Rndm();
   }
   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;
   Float_t *xpeaks = s->GetPositionX();
   for (p=0;p<nfound;p++) {
      Float_t xp = xpeaks[p];
      Int_t bin = h->GetXaxis()->FindBin(xp);
      Float_t yp = h->GetBinContent(bin);
      if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
      par[3*npeaks+2] = yp;
      par[3*npeaks+3] = xp;
      par[3*npeaks+4] = 3;
      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");             
}
 peaks.C:1
 peaks.C:2
 peaks.C:3
 peaks.C:4
 peaks.C:5
 peaks.C:6
 peaks.C:7
 peaks.C:8
 peaks.C:9
 peaks.C:10
 peaks.C:11
 peaks.C:12
 peaks.C:13
 peaks.C:14
 peaks.C:15
 peaks.C:16
 peaks.C:17
 peaks.C:18
 peaks.C:19
 peaks.C:20
 peaks.C:21
 peaks.C:22
 peaks.C:23
 peaks.C:24
 peaks.C:25
 peaks.C:26
 peaks.C:27
 peaks.C:28
 peaks.C:29
 peaks.C:30
 peaks.C:31
 peaks.C:32
 peaks.C:33
 peaks.C:34
 peaks.C:35
 peaks.C:36
 peaks.C:37
 peaks.C:38
 peaks.C:39
 peaks.C:40
 peaks.C:41
 peaks.C:42
 peaks.C:43
 peaks.C:44
 peaks.C:45
 peaks.C:46
 peaks.C:47
 peaks.C:48
 peaks.C:49
 peaks.C:50
 peaks.C:51
 peaks.C:52
 peaks.C:53
 peaks.C:54
 peaks.C:55
 peaks.C:56
 peaks.C:57
 peaks.C:58
 peaks.C:59
 peaks.C:60
 peaks.C:61
 peaks.C:62
 peaks.C:63
 peaks.C:64
 peaks.C:65
 peaks.C:66
 peaks.C:67
 peaks.C:68
 peaks.C:69
 peaks.C:70
 peaks.C:71
 peaks.C:72
 peaks.C:73
 peaks.C:74
 peaks.C:75
 peaks.C:76
 peaks.C:77
 peaks.C:78
 peaks.C:79
 peaks.C:80
 peaks.C:81
 peaks.C:82
 peaks.C:83
 peaks.C:84
 peaks.C:85
 peaks.C:86
 peaks.C:87
 peaks.C:88
 peaks.C:89
 peaks.C:90
 peaks.C:91
 peaks.C:92
 peaks.C:93
 peaks.C:94
 peaks.C:95
 peaks.C:96
 peaks.C:97
thumb