# How to fit in a fast way multiple Gaussian distribution with TSpectrum

Hi ROOTers,

I’m trying to fit a multiphoton spectrum in a good and fast way. The code that I’m using in order to plot the data
(Run1_PHA_HG_0_50.txt) of the multiphoton spectrum is:

``````{
int    i=0;
float  x;
FILE   *fp;
fp = fopen("Run1_PHA_HG_0_50.txt","r");

auto janus = new TH1F("Janus Stats", "Histo from Janus", 4192,0,4192);

while (!feof(fp)) {
fscanf(fp,"%g\n",&x);
i++;
janus->SetBinContent(i,x);
}
janus->Draw();
}
``````

which produce this kind of plot:

Every gaussian fits was produced with the FitPanel. I like the first fit (pederestal) but I’am practically certain that the other fits are wrong as they do not cover the whole gaussian distribution. Why does this happen?

Also, is there a quick way to automatically perform this series of n Gaussian fits? Maybe the answer is releted to the TSpectrum.

It is possible to modify the file peaks.C in order to take in input my data and perform its job? Someone can help me?

Yes, in peaks.C you have:

``````   h->FillRandom("f",200000);
``````

replace that by the code you posted.

1 Like

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

You need to play with the “TSpectrum::Search” parameters (so that it finds all your peaks).

1 Like

I understand, thank you so much, i’m playing with it.

By the way, why for example the 3rd, 4th, 5th and 6th peaks, are not fitted very well? There are some reason? @Wile_E_Coyote

Then, the value printed in the terminal (mean, error and number of peaks) are not consistent with the physics:

``````Found 11 candidate peaks to fit
Found 11 useful peaks to fit
Now fitting: Be patient
FCN=11082.1 FROM MIGRAD    STATUS=CONVERGED    5017 CALLS        5018 TOTAL
EDM=5.44313e-07    STRATEGY= 1  ERROR MATRIX UNCERTAINTY   4.0 per cent
EXT PARAMETER                                   STEP         FIRST
NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE
1  p0          -2.31356e+01   1.03729e+00   5.46459e-04  -1.27494e-03
2  p1           1.99605e-02   8.26821e-04  -4.64592e-07  -9.47037e-01
3  p2           1.72452e+03   1.00114e+01  -1.60856e-03   3.75176e-05
4  p3           3.25092e+02   6.84930e-02  -1.60303e-07   3.22999e-03
5  p4           1.08561e+01   8.84586e-02   1.78936e-05   2.74314e-03
6  p5           1.63066e+03   1.14808e+01   1.66776e-03   2.44295e-06
7  p6           3.94844e+02   7.00740e-02  -1.13766e-05  -5.46475e-03
8  p7           1.01732e+01   8.54962e-02   1.74743e-05   5.38027e-03
9  p8           1.54816e+03   1.16074e+01  -5.84550e-04  -3.83515e-05
10  p9           2.55121e+02   7.19296e-02  -5.72782e-06  -4.55286e-03
11  p10          1.08632e+01   9.16013e-02  -3.27183e-06  -3.19309e-03
12  p11          1.28438e+03   9.34127e+00  -2.24053e-03   6.77367e-05
13  p12          4.65027e+02   7.95340e-02   5.93631e-06   2.37573e-04
14  p13          9.29726e+00   9.30303e-02   4.03660e-05   3.44510e-03
15  p14          1.11710e+03   9.05473e+00  -1.96995e-03   4.38044e-05
16  p15          1.86842e+02   8.55997e-02   9.30022e-06   9.84315e-04
17  p16          9.84343e+00   9.92434e-02   1.85257e-05   5.63803e-03
18  p17          8.88451e+02   8.42141e+00   1.96596e-03   1.40685e-05
19  p18          5.34955e+02   1.03965e-01   2.33208e-05   5.65390e-04
20  p19          8.46223e+00   9.90845e-02   2.96006e-05   2.66858e-04
21  p20          5.10919e+02   7.64177e+00   3.53640e-03  -1.36469e-04
22  p21          6.05296e+02   1.42041e-01   1.02773e-05   8.58372e-04
23  p22          7.55737e+00   1.26177e-01   1.80171e-05  -4.05531e-03
24  p23          5.75546e+02   8.08751e+00  -9.79202e-04   4.51454e-05
25  p24          1.20034e+02   1.03445e-01  -2.19113e-06   5.91139e-04
26  p25          7.73818e+00   1.07247e-01  -4.80358e-06  -4.97588e-03
27  p26          2.70839e+02   5.40511e+00   3.56555e-03  -6.64302e-05
28  p27          6.76077e+02   2.03411e-01  -7.87397e-05  -1.14759e-03
29  p28          7.15198e+00   1.72092e-01  -5.71297e-05  -2.61504e-03
30  p29          5.45974e+02   2.92205e+00  -1.68201e-03  -1.60541e-04
31  p30          4.43123e+02   5.96620e-01   2.28555e-04  -4.39655e-04
32  p31          1.98424e+02   6.71806e-01  -7.07737e-05  -1.11238e-03
33  p32          1.37485e+02   4.66569e+00   2.23185e-05   5.13067e-05
34  p33          5.61474e+01   1.69687e-01   6.00528e-05   1.31144e-03
35  p34          4.39090e+00   1.32999e-01   5.39130e-05   4.11720e-03
``````

For example, I have only 11 peaks, meanwhile in the terminal are displayed 35 peaks and if we see the 1st line, the value are not coerent with the reality…I don’t understand why…

I’m using the file of Rene Brum that we can find here: ROOT: tutorials/spectrum/peaks.C File Reference

I edited it in order to fit my multiphon peaks, as you can read in the preavius post here. I change the npeaks but the results are not coherent…

You either need a better background model in your fit function, or you need to subtract the background in advance:

``````{
TTree *t = new TTree("t", "t"); // a temporary tree
TH1D *h = new TH1D("h", "Histo;ADC [ch];Counts",
t->GetEntries(), 0., (Double_t)t->GetEntries());
t->Project("h", "Entry\$", "v");
delete t; // no longer needed
h->Sumw2(kFALSE); h->ResetStats();
TSpectrum *s = new TSpectrum();
TH1 *bg = s->Background(h, 40, "nosmoothing"); // estimate the background
bg->SetLineColor(kRed);
TCanvas *c = new TCanvas("c", "c");
c->Divide(1, 2);
c->cd(1);
h->DrawCopy("");
bg->Draw("SAME");
c->cd(2);
h->Add(bg, -1.); // subtract the estimated background
Int_t nfound = s->Search(h, 7., "", 0.001); // search for peaks
// h->Draw();
c->cd(0);
}
``````
1 Like

Thanks for the time that you are spendind. Now I did manually the fit and the results seems to be much better than before:

So, the trick for the “bst fit” is subtrack the correct background?

@Wile_E_Coyote, I’m trying to make the procedecure of fit automatic, the code that I wrote is the following:

``````{

TTree *t = new TTree("t", "t"); // A temporary tree

TH1D *h = new TH1D("h", "Multiphoton Spectrum; ADC [ch]; Counts", t->GetEntries(), 0., (Double_t)t->GetEntries()); // A new histogram

t->Project("h", "Entry\$", "v");

delete t; // No longer needed

h->Sumw2(kFALSE);
h->ResetStats();

TSpectrum *s = new TSpectrum(); // Use TSpectrum to find the peak candidates

TH1 *bg = s->Background(h, 40, "nosmoothing"); // Estimate the background

bg->SetLineColor(kRed); // Set the color

TCanvas *c = new TCanvas("c", "c");

c->Divide(1, 2); // The pad is 1 coloumn and two lines
c->cd(1); // Move to line one

h->DrawCopy("");

bg->Draw("SAME"); // Draw the bg in the first line

c->cd(2); // Move to line two

h->Add(bg, -1.); // Subtract the estimated background

Int_t nfound = s->Search(h, 7., "", 0.05); // Search for peaks

printf("We found %d candidate peaks to fit.\n",nfound);

npeaks = 0;
Double_t *xpeaks;
Double_t par[3000];
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);

par[3*npeaks] = yp; // "height"
par[3*npeaks+1] = xp; // "mean"
par[3*npeaks+2] = 3; // "sigma"

npeaks++;
}

printf("Found %d useful peaks to fit\n",npeaks);
printf("Now fitting: Be patient\n");

TF1 *fit = new TF1("fit",fpeaks,0,4096,3*npeaks);

TVirtualFitter::Fitter(h,3*npeaks);
fit->SetParameters(par);
fit->SetNpx(1000);

h->Fit("fit");
h->Draw(); // Draw the plot without bg

c->cd(0);

}
``````

and the result is:

`````` *** Break *** segmentation violation

``````

I tried to incorporate part of the previus code (that works for some aspect) but probably i made some errors…

For the (original) “`fpeaks`” function to work, the “`npeaks`” must be a global variable.
Also, the (original) “`fpeaks`” function expects the first two parameters to be “responsible” for a (background) linear function.

1 Like

If I understand well… Now `npeaks` is a global variable but nothing seems to be changed…

``````{

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;
}

TTree *t = new TTree("t", "t"); // A temporary tree

TH1D *h = new TH1D("h", "Multiphoton Spectrum; ADC [ch]; Counts", t->GetEntries(), 0., (Double_t)t->GetEntries()); // A new histogram

t->Project("h", "Entry\$", "v");

delete t; // No longer needed

h->Sumw2(kFALSE);
h->ResetStats();

TSpectrum *s = new TSpectrum(); // Use TSpectrum to find the peak candidates

TH1 *bg = s->Background(h, 40, "nosmoothing"); // Estimate the background

bg->SetLineColor(kRed); // Set the color

TCanvas *c = new TCanvas("c", "c");

c->Divide(1, 2); // The pad is 1 coloumn and two lines
c->cd(1);        // Move to line one

h->DrawCopy("");

bg->Draw("SAME"); // Draw the bg in the first line

c->cd(2); // Move to line two

h->Add(bg, -1.); // Subtract the estimated background

Int_t nfound = s->Search(h, 7., "", 0.05); // Search for peaks
Int_t npeaks = 0;

printf("We found %d candidate peaks to fit.\n", nfound);

Double_t *xpeaks;
Double_t par[3000];
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);

par[3 * npeaks] = yp;     // "height"
par[3 * npeaks + 1] = xp; // "mean"
par[3 * npeaks + 2] = 3;  // "sigma"

npeaks++;
}

printf("Found %d useful peaks to fit.\n", npeaks);
printf("Now fitting: Be patient.\n");

TF1 *fit = new TF1("fit", fpeaks, 0, 4096, 3 * npeaks);

TVirtualFitter::Fitter(h, 3 * npeaks);
fit->SetParameters(par);
fit->SetNpx(1000);

h->Fit("fit");
h->Draw(); // Draw the plot without bg

c->cd(0);

}

``````

@Wile_E_Coyote , What you mean for:

``````Also, the (original) “fpeaks” function expects the first two parameters to be
“responsible” for a (background) linear function.
``````

?

I’m not understanding…

@dastudillo I defined `npeaks` outside all the functions:

``````Int_t npeaks = 0;

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 Prova(){

TTree *t = new TTree("t", "t"); // A temporary tree

TH1D *h = new TH1D("h", "Multiphoton Spectrum; ADC [ch]; Counts", t->GetEntries(), 0., (Double_t)t->GetEntries()); // A new histogram

t->Project("h", "Entry\$", "v");

delete t; // No longer needed

h->Sumw2(kFALSE);
h->ResetStats();

TSpectrum *s = new TSpectrum(); // Use TSpectrum to find the peak candidates

TH1 *bg = s->Background(h, 40, "nosmoothing"); // Estimate the background

bg->SetLineColor(kRed); // Set the color

TCanvas *c = new TCanvas("c", "c");

c->Divide(1, 2); // The pad is 1 coloumn and two lines
c->cd(1);        // Move to line one

h->DrawCopy("");

bg->Draw("SAME"); // Draw the bg in the first line

c->cd(2); // Move to line two

h->Add(bg, -1.); // Subtract the estimated background

Int_t nfound = s->Search(h, 7., "", 0.05); // Search for peaks

printf("We found %d candidate peaks to fit.\n", nfound);

Double_t *xpeaks;
Double_t par[3000];
xpeaks = s->GetPositionX();

for (Int_t p = 0; p < nfound; p++)
{
Double_t xp = xpeaks[p];
Int_t bin = h->GetXaxis()->FindBin(xp);
Double_t yp = h->GetBinContent(bin);

par[3 * npeaks] = yp;     // "height"
par[3 * npeaks + 1] = xp; // "mean"
par[3 * npeaks + 2] = 3;  // "sigma"

npeaks++;
}

printf("Found %d useful peaks to fit.\n", npeaks);
printf("Now fitting: Be patient.\n");

TF1 *fit = new TF1("fit", fpeaks, 0, 4096, 3 * npeaks);

TVirtualFitter::Fitter(h, 3 * npeaks);
fit->SetParameters(par);
fit->SetNpx(1000);

h->Fit("fit");
h->Draw(); // Draw the plot without bg

c->cd(0);
}
``````

The result is the following:

What’s happening?

In the “`fpeaks`”, you have (note the “`+ 2`”):
`Double_t norm = par[3 * p + 2]; // "height" or "area"`
but in your “`Prova`”, you set:
`par[3 * npeaks] = yp; // "height"`

Right, was an error. Now i restored the original code of Brune:

``````Int_t npeaks = 0;

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 Multiphoton(){

TTree *t = new TTree("t", "t"); // A temporary tree

TH1D *h = new TH1D("h", "Multiphoton Spectrum; ADC [ch]; Counts", t->GetEntries(), 0., (Double_t)t->GetEntries()); // A new histogram

t->Project("h", "Entry\$", "v");

delete t; // No longer needed

h->Sumw2(kFALSE);
h->ResetStats();

TSpectrum *s = new TSpectrum(); // Use TSpectrum to find the peak candidates

TH1 *bg = s->Background(h, 40, "nosmoothing"); // Estimate the background

bg->SetLineColor(kRed); // Set the color

TCanvas *c = new TCanvas("c", "c");

c->Divide(1, 2); // The pad is 1 coloumn and two lines
c->cd(1);        // Move to line one

h->DrawCopy("");

bg->Draw("SAME"); // Draw the bg in the first line

c->cd(2); // Move to line two

h->Add(bg, -1.); // Subtract the estimated background

Int_t nfound = s->Search(h, 7., "", 0.05); // Search for peaks

printf("We found %d candidate peaks to fit.\n", nfound);

Double_t *xpeaks;
Double_t par[3000];
xpeaks = s->GetPositionX();

for (Int_t p = 0; p < nfound; p++)
{
Double_t xp = xpeaks[p];
Int_t bin = h->GetXaxis()->FindBin(xp);
Double_t yp = h->GetBinContent(bin);

par[3 * npeaks + 2] = yp;     // "height"
par[3 * npeaks + 3] = xp; // "mean"
par[3 * npeaks + 4] = 3;  // "sigma"

npeaks++;
}

printf("Now fitting: Be patient.\n");

TF1 *fit = new TF1("fit", fpeaks, 0, 4096, 3 * npeaks);

TVirtualFitter::Fitter(h, 3 * npeaks);
fit->SetParameters(par);
fit->SetNpx(1000);

h->Fit("fit");
h->Draw(); // Draw the plot without bg

c->cd(0);
}
``````

The parameters are now egual but…

Ps, if I re-run the code the output change…it is so strange:

Where did this strange “`+ 2`” come from
Can it be that the total number of parameters (in “`fpeaks`”) is not equal to “`3 * npeaks`

I suppose that the `+2,+3 and +4` comes out from the fact that the parameter `p[0], p[1]` are already used in `result`.

In the original file, the number of parameters are `2+3*npeaks`, in my case I suppose that are `3*npeaks+2`.

But I’m not sure…

I also tried to play with the parameters but I get worse results…where I wrong? @Wile_E_Coyote

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.