Perfect, your hint was very appreciate. Thanks.
I modified the file in this way:
/// \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;p<npeaks;p++) {
Double_t norm = par[3*p+2]; // "height" or "area"
Double_t mean = par[3*p+3];
Double_t sigma = par[3*p+4];
#if defined(__PEAKS_C_FIT_AREAS__)
norm /= sigma * (TMath::Sqrt(TMath::TwoPi())); // "area"
#endif /* defined(__PEAKS_C_FIT_AREAS__) */
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",4096,0,4096);
// 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; // "height"
par[3*p+3] = 10+gRandom->Rndm()*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,4096,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);
int i=0;
float x;
FILE *fp;
fp = fopen("Run1_PHA_HG_0_50.txt","r");
while (!feof(fp)) {
fscanf(fp,"%g\n",&x);
i++;
h->SetBinContent(i,x);
}
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,4096);
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;p<nfound;p++) {
Double_t xp = xpeaks[p];
Int_t bin = h->GetXaxis()->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");
}
The result that I obtained is the following:
Why the first peak was not fitted like the 11th? And why not all the fitted curve reach the peak? @couet
