# Gaussian fitting of mutiple peaks in the case of unknown number of peaks in histogram generated from .tex file (for example: test.tex). What are the change required in following code (this code is for 5 peaks)?

ROOT Version: 6.14
Platform: Not Provided
Compiler: Not Provided

``````TF1 *fpol;
Int_t npol;

const Int_t npeaks = 5;
Double_t fpeaks(Double_t *x, Double_t *par) {
Double_t result = fpol->EvalPar(x,par);
for (Int_t p=0;p<npeaks;p++) {
Double_t norm  = par[3*p+npol];
Double_t mean  = par[3*p+npol+1];
Double_t sigma = par[3*p+npol+2];
result += norm*TMath::Gaus(x[0],mean,sigma);
}
return result;
}
void fitg() {
TTree *t = new TTree("t", "my tree");
TH1D *h = new TH1D("h", "Histo", , ., .);
//fit background with a polynomial
fpol = new TF1("fpol","pol3",60,550);
npol = fpol->GetNpar();
h->Fit("fpol","0");
Double_t xpeak[npeaks] = {200, , , , };
Double_t sigma[npeaks] = {2, , , , };
Double_t amp[npeaks]   = {100,100,100,100,100};
Double_t par[100];
for (Int_t ip=0;ip<npol;ip++) par[ip] = fpol->GetParameter(ip);
TF1 *fp = new TF1("fp",fpeaks,60,550,npol+3*npeaks);
for (Int_t i=0;i<npeaks;i++) {
fp->SetParameter(3*i+npol,1000*amp[i]);
fp->SetParameter(3*i+npol+1, xpeak[i]);
fp->SetParLimits(3*i+npol+1, xpeak[i]-5,xpeak[i]+5);
fp->SetParameter(3*i+npol+2, sigma[i]);
fp->SetParLimits(3*i+npol+2, 0.5,3*sigma[i]);
}
fp->SetNpx(5000);
TVirtualFitter::SetMaxIterations(10000);
h->Fit(fp,"rw");

}
``````

Hi,
you can give a look to TSpectrum class.

After declaring

``````TSpectrum *spectrum= new TSpectrum();
spectrum->Search(h);
``````

Then you can use the method TSpectrum::GetNPeaks() and TSpectrum::GetPositionX() to have the number of peaks and vector with the peak position.

Then your code should work fine using the informations get from TSpectrum.

Cheers,
S

Thanks, but i already tried with the following code (modified from > tutorial>Tspectrum>peaks.C) which didn’t work.

``````#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=5)
{
TTree *t = new TTree("t", "my tree");
// creating a histogram ...
TH1D *h = new TH1D("h", "my histo", 200, 100., 350.);
// fill histo ...
t->Project("h", "x");

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->cd(1);
h->Draw();
//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;
}
``````

Hi,