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)?

Newbie
1

Please read tips for efficient and successful posting and posting code

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");
  t->ReadFile("test.txt", "x/D");
  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");
  
}
2

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

3

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");
  t->ReadFile("a.txt", "x/D");
  // 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; 
}
4

Hi,
Can you show your histogram in a post please ?
You say it is not working. Do you have an error message ? What is not exactly working ?
Cheers,
R

5

It behaves abnormally, sometimes I get the error stating “par” and other undeclared identifiers and sometimes it runs giving the attached result.

c
c866×613 40.8 KB

This is what I want to do, as in this example (the only difference is that 1. I have .tex file, not the root file and 2. no. of peaks are not known, they need to be automatically detected. ) I’m attaching my sample histo

sample%20histo
sample%20histo1835×983 7.24 KB
whose Gaussian peak fitting is planned.

Maybe you can help, hopefully !!

7

Hi, Do you know how to perform this ? It’s quit urgent.