TSpectrum and something else

Hi ROOTers,

We wrote this code in order to analize multiphoton spectrum:

FILE* file;
TH1F* hsipm;
float x;
int bins = 0;
TF1* full;
int sigma = 5;

const Int_t npeaks = 40;
Double_t par[3*npeaks];
Double_t fpeaks(Double_t *x, Double_t *par) {
Double_t result = 0.;

   for (Int_t p=0;p<npeaks;p++) {

      Double_t norm  = par[3*p+0];
      Double_t mean  = par[3*p+1];
      Double_t sigma = par[3*p+2];

      result += norm*TMath::Gaus(x[0],mean,sigma);
   }
   return result;
}

Double_t gain = 0.;
Double_t pedestal = 0.;

/*

//par[0]=mu -  par[1]=p  - par[2]=N (articolo vinogradov sipm1 pag. 6)
  Double_t Vino(Double_t* x, Double_t* par){
  Double_t mu = par[0];
  Double_t p = par[1];
//Double_t N = par[2];
  Int_t N = int((x[0]-pedestal)/gain);
  Double_t B;

  Double_t f0 = exp(-mu);
  Double_t term = mu*(1-p);
  Double_t sum = 0.;

  if(N==0)
    return f0;

  for(int i=1; i<N; ++i){

    B = TMath::Gamma(N)/(TMath::Gamma(i+1)*TMath::Gamma(i)*TMath::Gamma(N-i+1));
    sum += B*pow(term,i)*pow(p,N-i);

  }

  return sum*f0;

}

*/

int Multiphoton_spectrum_analysis() {

  file = fopen("Run0_PHA_HG_0_10.txt", "r");

  TCanvas* c1 = new TCanvas("c1","HG Multiphoton Spectrum",10,10,1920,1980);

  gStyle->SetOptFit(000);

  c1->SetGridx();
  c1->SetGridy();

  TH1F* htmp = new TH1F("htemp","htemp",4096,0,4096);

  while(!feof(file)){

    fscanf(file,"%g\n",&x);

    htmp->SetBinContent(bins,x);

    ++bins;

  }

  hsipm = new TH1F("HG Multiphoton Spectrum","",bins,0,bins);

  for(int a=0; a<bins; ++a)

  hsipm->SetBinContent(a,htmp->GetBinContent(a));

  sigma = bins/1024+1;

  hsipm->GetYaxis()->SetTitle("#Entries");
  hsipm->GetXaxis()->SetTitle("ADC channels");;
  hsipm->GetXaxis()->SetRange(0,1300);


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

  f->SetParameters(par);

  TSpectrum *s = new TSpectrum(npeaks);

  Int_t nfound;

  TString peak_name;
  TString full_func = "";

  TF1* g[npeaks];

  nfound = s->Search(hsipm,sigma,"",0.0007);

  printf("Found %d candidate peaks to fit\n", nfound);

  int peak_x_tmp[npeaks];
  double peak_x[npeaks];
  double peak_y[npeaks];

  TMath::Sort(nfound,s->GetPositionX(),peak_x_tmp,kFALSE);

  for(int a=0; a<nfound; ++a){

    peak_x[a] = s->GetPositionX()[peak_x_tmp[a]];
    peak_y[a] = s->GetPositionY()[peak_x_tmp[a]];

  }


  for(int a=0; a<nfound; ++a){

    peak_name.Form("g%d",a);
    g[a] = new TF1(peak_name,"gaus");

    TString func_tmp;

    func_tmp.Form("gaus(%d)+",3*a);
    full_func += func_tmp;

  }

  full_func.Remove(full_func.Length()-1);
  full = new TF1("full",full_func);

  for(int a=2; a<nfound; ++a){

    full->SetParLimits(3*a+0,peak_y[a]*0.5,peak_y[a]*1.1);
    full->SetParLimits(3*a+1,peak_x[a]-3*sigma, peak_x[a]+3*sigma);
    full->SetParLimits(3*a+2,0.5*sigma,6*sigma);

  }

  // pedestal

  full->SetParLimits(0,peak_y[0]*0.6,peak_y[0]*1.1);
  full->SetParLimits(1,peak_x[0]-1.2*sigma, peak_x[0]+1.2*sigma);
  full->SetParLimits(2,0.5*sigma,1*sigma);

  // 1 p.e.

  full->SetParLimits(3,peak_y[1]*0.5,peak_y[1]*1.1);
  full->SetParLimits(4,peak_x[1]-1.5*sigma, peak_x[1]+1.5*sigma);
  full->SetParLimits(5,0.5*sigma,4*sigma);
  full->SetNpx(4096*2);
  full->SetLineStyle(9);
  full->SetLineColor(2); //colore fit generale


  hsipm->Fit("full","","same",35.,peak_x[nfound-1]+4*sigma);
  hsipm->SetLineColor(1); //colore dei dati

  double fPar[npeaks] = {0.};
  double fErr[npeaks*3] = {0.};

  full->GetParameters(&fPar[0]);

  for(int a=0; a<nfound*3; ++a)

    fErr[a] = full->GetParErrors()[a];

  TF1* fg[npeaks];

  TString gname;

  double m_tmp = fPar[1];
  double s_tmp = fPar[2];

  fg[0] = new TF1("ped","gaus",m_tmp-3*s_tmp,m_tmp+3*s_tmp);
  fg[0]->SetParameters(&fPar[0]);
  fg[0]->SetLineColor(9); //colore piedistallo
  //fg[0]->Draw("same");

  for(int a=1; a<nfound; ++a){

    m_tmp = fPar[3*a+1];
    s_tmp = fPar[3*a+2];

    gname.Form("%d p.e.",a);

    fg[a] = new TF1(gname,"gaus",m_tmp-3*s_tmp,m_tmp+3*s_tmp);
    fg[a]->SetParameters(&fPar[3*a]);
    fg[a]->SetLineColor(9); //colore gaussiane
    //fg[a]->Draw("same");

  }

   double aNpeak[npeaks];
   double aConst[npeaks];
   double aMean[npeaks];
   double aSig[npeaks];
   double aConstErr[npeaks];
   double aMeanErr[npeaks];
   double aSigErr[npeaks];

   const double sqrt2pi = sqrt(2 * acos(-1.));

   for (int a = 1; a < nfound; ++a)

   {
     aNpeak[a - 1] = a;
     aMean[a - 1] = fPar[3 * a + 1];
     aMeanErr[a - 1] = fErr[3 * a + 1];
     aSig[a - 1] = fPar[3 * a + 2];
     aSigErr[a - 1] = fErr[3 * a + 2];
     aConst[a - 1] = fPar[3 * a] * sqrt2pi * aSig[a - 1];
     aConstErr[a - 1] = sqrt2pi * (fErr[3 * a] * aSig[a - 1] + aSigErr[a - 1] * fPar[3 * a]);
  }

  FILE* outfile;

  outfile = fopen("Gaussian_Parameters.txt", "w");

  for(int a=0; a<nfound; ++a)


  fprintf(outfile, "%d \t %f \t %f \t %f \t %f \t %f \t %f \n", a, fPar[3 * a], fErr[3*a], fPar[3 * a + 1], fErr[3*a+1], fPar[3 * a + 2], fErr[3*a+2]);
  fclose(outfile);

  TCanvas* c2 = new TCanvas("c2","HG Multiphoton Spectrum Analysis",10,10,1920,1980);

  c2->UseCurrentStyle();
  gStyle->SetOptFit(111);

  c2->Divide(2,2);
  c2->cd(1);
  c2->GetPad(1)->SetGridx();
  c2->GetPad(1)->SetGridy();
  c2->Update();

  Double_t xErr[npeaks] = {0.};

  TGraphErrors* mean = new TGraphErrors(nfound-1,aNpeak,aMean,xErr,aMeanErr);

  mean->SetMarkerStyle(22);
  mean->SetTitle("HG Calibration");
  mean->GetXaxis()->SetTitle("#p.e.");
  mean->GetYaxis()->SetTitle("ADC channels");
  mean->Draw("ap");

  TF1* fgain = new TF1("fgain","pol1");

  mean->Fit("fgain","","same",0.,nfound-2);
  mean->GetYaxis()->SetRangeUser(0.,aMean[nfound-2]+200.);

  c2->cd(1)->Update();

  pedestal = fgain->GetParameter(0);
  gain = fgain->GetParameter(1);

  c2->cd(2);
  c2->GetPad(2)->SetGridx();
  c2->GetPad(2)->SetGridy();

  Double_t aDpp[npeaks] = {0.};
  Double_t aDppErr[npeaks] = {0.};
  Double_t mean_dpp = 0.;

  for(int a=0; a<nfound-2; ++a){

    aDpp[a] = aMean[a+1]-aMean[a];
    aDppErr[a] = aMeanErr[a+1]+aMeanErr[a];
    mean_dpp += aDpp[a];
  }

  mean_dpp /= (nfound-2);

  TGraphErrors* gDpp = new TGraphErrors(nfound-2,aNpeak,aDpp,xErr,aDppErr);

  gDpp->SetMarkerStyle(22);
  gDpp->SetTitle("Delta peak to peak distribution");
  gDpp->GetXaxis()->SetTitle("#p.e.");
  gDpp->GetYaxis()->SetTitle("ADC channels");
  gDpp->Draw("ap");
  gDpp->Fit("pol1","","same",0.,nfound-3);
  gDpp->GetYaxis()->SetRangeUser(mean_dpp-50,mean_dpp+50.);


  c2->cd(3);
  c2->GetPad(3)->SetGridx();
  c2->GetPad(3)->SetGridy();

  TGraphErrors* gConst = new TGraphErrors(nfound-1,aNpeak,aConst,xErr,aConstErr);

  gConst->SetMarkerStyle(22);
  gConst->SetTitle("Gaussian Normalization par distribution");
  gConst->GetXaxis()->SetTitle("#p.e.");
  gConst->GetYaxis()->SetTitle("#Entries");
  gConst->Draw("ap");


  TF1 *poisson = new TF1("poisson","[0]*TMath::Power(([1]/[2]),(x/[2]))*(TMath::Exp(-([1]/[2])))/TMath::Gamma((x/[2])+1.)", 0, nfound);

  poisson->SetParameters(1.e6, 9., 2.);
  poisson->SetNpx(4096*2);

  gConst->Fit(poisson,"R","same",0.,nfound-2);


  c2->cd(4);
  c2->GetPad(4)->SetGridx();
  c2->GetPad(4)->SetGridy();

  TGraphErrors* gSig = new TGraphErrors(nfound-1,aNpeak,aSig,xErr,aSigErr);

  gSig->SetMarkerStyle(22);
  gSig->SetTitle("Gaussian Sigma par distribution");
  gSig->GetXaxis()->SetTitle("#p.e.");
  gSig->GetYaxis()->SetTitle("ADC channels");
  gSig->Draw("ap");
  gSig->GetYaxis()->SetRangeUser(0.,50.);

  TF1* pol1 = new TF1("pol1","pol1");

  pol1->SetParameters(20.,0.5);

  gSig->Fit("pol1","","same",0.,nfound-2);

  return 0;
}

And we get this results:

and

meanwhile in the terminal i get this warning:

Info in <ROOT::Math::ParameterSettings>: lower/upper bounds outside current parameter value. The value will be set to (low+up)/2

  1. How can I delete this warning? I suppose that i need to set parameterlimits but i did it.

  2. Why the chi squared are so big? The points lie on the lines…

  3. I saw that TSpectrum is deprecated, is there an alternative?

  4. There is a way to set two different style about statistic table? I used gStyle->SetOptFit(000) and gStyle->SetOptFit(111) in the two canvas but doesnt work.

  5. There is a way to do a similar fit (green) in ROOT? How?

I mean, do a multi gaussian fit like green line.

This setting is global.
@moneta surely knows about the other questions

1 Like

For the “lower/upper bounds outside current parameter value” warning, it seems to me that you call “SetParLimits” but there is no preceding “SetParameter” call.

For the “big” chi^2, notice the really “big” differences between the “data” (black) and the “fit” (red), especially near the highest peaks’ maxima and minima. You might try to play with fit options (e.g., try "L" or "P").

For the global style setting, see the TExec class reference and ROOT Forum → Search → TExec

1 Like

Ok, so there is no way to do locally, thanks @couet !

Yes as @Wile_E_Coyote said you can use TExec

  • I add the line: full->SetParameters(par); but nothing change.

  • I played with them, i obtain the best result with WW, but it is no acceptable (in this way I exclude the errors, right?). A “good” result is with S option.

  • I’m trying to use it, thanks.

Maybe the contents of “par” do not fully correspond to “full” parameters (note that you need to set their values before setting their limits).

The "S" option is not relevant to error handling.

I’m pretty sure that par fully correspond to full parameters. By the way, i did this:

full = new TF1("full",full_func);
full->SetParameters(par);

but nothing change.

I upload a data ile in order to run the code: Run0_PHA_HG_0_10.txt (15,6 KB)

It seems to me that your “par” contains only zeros (i.e., it’s uninitialized).

You mean that: Double_t par[3*npeaks]={0.}; ?

Why don’t you simply print the first several elements of your “par”.

par is a vector of zero, you are right.

A blast to the past for me and a nice fit for you. The magic word is hyper Gaussian

#include "TCanvas.h"
#include "TStyle.h"
#include "TGraphErrors.h"
#include "TArrayD.h"
#include "TH1F.h"
#include "TSpectrum.h"
#include "TF1.h"

const Int_t npeaks = 40;

Double_t HyperGaus(Double_t *x,Double_t *par)
{
   Double_t mean  = par[0];
   Double_t sigma = par[1];
   Double_t gamma = par[2];

   if (sigma == 0) return 1.e30;
   Double_t arg = TMath::Abs((x[0]-mean)/sigma);
   Double_t powarg = TMath::Power(arg,gamma);
   Double_t res = TMath::Exp(-0.5*powarg);
   return res;
}

Int_t nfound;
Double_t SumHyperGaus(Double_t *x,Double_t *par) {

   Double_t result = 0.;
   for (Int_t p=0;p<nfound;p++) {

      Double_t norm  = par[4*p+0];
      result += norm*HyperGaus(x,par+4*p+1);
   }
   return result;
}

int Multiphoton_spectrum_analysis() {

  FILE *file = fopen("./Run0_PHA_HG_0_10.txt", "r");

  TCanvas* c1 = new TCanvas("c1","HG Multiphoton Spectrum",10,10,1920,1980);

  gStyle->SetOptFit(000);

  c1->SetGridx();
  c1->SetGridy();

  TH1F* htmp = new TH1F("htemp","htemp",4096,0,4096);

  Int_t bins = 0;
  Float_t x;
  while(!feof(file)){
    fscanf(file,"%g\n",&x);
    htmp->SetBinContent(bins,x);
    ++bins;
  }

  TH1F *hsipm = new TH1F("HG Multiphoton Spectrum","",bins,0,bins);

  for(int a=0; a<bins; ++a)

  hsipm->SetBinContent(a,htmp->GetBinContent(a));

  Int_t sigma = bins/1024+1;

  hsipm->GetYaxis()->SetTitle("#Entries");
  hsipm->GetXaxis()->SetTitle("ADC channels");;
  hsipm->GetXaxis()->SetRange(0,1300);

  TSpectrum *s = new TSpectrum(npeaks);

  nfound = s->Search(hsipm,sigma,"",0.0007);

  printf("Found %d candidate peaks to fit\n", nfound);

  int peak_x_tmp[nfound];
  double peak_x[nfound];
  double peak_y[nfound];

  TMath::Sort(nfound,s->GetPositionX(),peak_x_tmp,kFALSE);

  for(int a=0; a<nfound; ++a){

    peak_x[a] = s->GetPositionX()[peak_x_tmp[a]];
    peak_y[a] = s->GetPositionY()[peak_x_tmp[a]];

  }

  TF1* f = new TF1("f",SumHyperGaus,0,4096,4*nfound);
  f->SetNpx(4096*2);
  f->SetLineStyle(9);
  f->SetLineColor(2); //colore fit generale

  for(int a=0; a<nfound; ++a){
    f->SetParameter(4*a+0,peak_y[a]);
    f->SetParameter(4*a+1,peak_x[a]);
    f->SetParameter(4*a+2,sigma);
    f->SetParameter(4*a+3,1.2);
  }

  for(int a=2; a<nfound; ++a){

    f->SetParLimits(4*a+0,peak_y[a]*0.5,peak_y[a]*1.1);
    f->SetParLimits(4*a+1,peak_x[a]-3*sigma, peak_x[a]+3*sigma);
    f->SetParLimits(4*a+2,0.5*sigma,6*sigma);
    f->SetParLimits(4*a+3,1.1,2.5);
    //f->SetParLimits(4*a+3,1.2,1.2);

  }

  // pedestal

  f->SetParLimits(0,peak_y[0]*0.6,peak_y[0]*1.1);
  f->SetParLimits(1,peak_x[0]-1.2*sigma, peak_x[0]+1.2*sigma);
  f->SetParLimits(2,0.5*sigma,1*sigma);
  f->SetParLimits(3,1.1,2.5);
  //f->SetParLimits(3,1.2,1.2);

  // 1 p.e.

  f->SetParLimits(4,peak_y[1]*0.5,peak_y[1]*1.1);
  f->SetParLimits(5,peak_x[1]-1.5*sigma, peak_x[1]+1.5*sigma);
  f->SetParLimits(6,0.5*sigma,4*sigma);
  f->SetParLimits(7,1.1,2.5);
  //f->SetParLimits(7,1.2,1.2);

  // option I makes sure the area under fitted curve is fitted, but this
  // takes a significant amount of CPU time
  //hsipm->Fit("f","I","same",35.,peak_x[nfound-1]+4*sigma);
  hsipm->Fit("f","","same",35.,peak_x[nfound-1]+4*sigma);
  hsipm->SetLineColor(1); //colore dei dati

  TArrayD fPar(4*nfound);
  TArrayD fErr(4*nfound);

  for(int a=0; a<nfound; ++a) {
    fPar[4*a+0] = f->GetParameter(4*a+0);
    fPar[4*a+1] = f->GetParameter(4*a+1);
    fPar[4*a+2] = f->GetParameter(4*a+2);
    fPar[4*a+3] = f->GetParameter(4*a+3);
    fErr[4*a+0] = f->GetParError(4*a+0);
    fErr[4*a+1] = f->GetParError(4*a+1);
    fErr[4*a+2] = f->GetParError(4*a+2);
    fErr[4*a+3] = f->GetParError(4*a+3);
  }

  TF1* fg[npeaks];

  double m_tmp = fPar[1];
  double s_tmp = fPar[2];

  fg[0] = new TF1("ped",HyperGaus,m_tmp-3*s_tmp,m_tmp+3*s_tmp);
  fg[0]->SetParameters(&fPar[0]);
  fg[0]->SetLineColor(9); //colore piedistallo
  //fg[0]->Draw("same");

  for(int a=1; a<nfound; ++a){

    m_tmp = fPar[4*a+1];
    s_tmp = fPar[4*a+2];

    TString gname;
    gname.Form("%d p.e.",a);

    fg[a] = new TF1(gname,HyperGaus,m_tmp-3*s_tmp,m_tmp+3*s_tmp);
    fg[a]->SetParameters(&fPar[4*a]);
    fg[a]->SetLineColor(9); //colore gaussiane
    //fg[a]->Draw("same");

  }

   TArrayD aNpeak(nfound-1);
   TArrayD aConst(nfound-1);
   TArrayD aMean (nfound-1);
   TArrayD aSig  (nfound-1);
   TArrayD aGamma(nfound-1);

   TArrayD aConstErr(nfound-1);
   TArrayD aMeanErr (nfound-1);
   TArrayD aSigErr  (nfound-1);
   TArrayD aGammaErr(nfound-1);

   const double sqrt2pi = sqrt(2 * acos(-1.));

   for (int a = 1; a < nfound; ++a)
   {
     aNpeak   [a-1] = a;
     aMean    [a-1] = fPar[4*a+1];
     aMeanErr [a-1] = fErr[4*a+1];
     aSig     [a-1] = fPar[4*a+2];
     aSigErr  [a-1] = fErr[4*a+2];
     aGamma   [a-1] = fPar[4*a+3];
     aGammaErr[a-1] = fErr[4*a+3];

     aConst   [a-1] = fPar[4*a]*sqrt2pi*aSig[a-1];
     aConstErr[a-1] = sqrt2pi*(fErr[4*a]*aSig[a-1]+aSigErr[a-1]*fPar[4*a]);
  }

  FILE* outfile;

  outfile = fopen("Gaussian_Parameters.txt", "w");

  for(int a=0; a<nfound; ++a)
    fprintf(outfile, "%d \t %f \t %f \t %f \t %f \t %f \t %f \t %f \t %f \n",
                       a,fPar[4*a],fErr[4*a],fPar[4*a+1],fErr[4*a+1],fPar[4*a+2],fErr[4*a+2],fPar[4*a+3],fErr[4*a+3]);
  fclose(outfile);

  TCanvas* c2 = new TCanvas("c2","HG Multiphoton Spectrum Analysis",10,10,1920,1980);

  c2->UseCurrentStyle();
  gStyle->SetOptFit(111);

  c2->Divide(2,2);
  c2->cd(1);
  c2->GetPad(1)->SetGridx();
  c2->GetPad(1)->SetGridy();
  c2->Update();

  TArrayD xErr(npeaks);

  TGraphErrors* mean = new TGraphErrors(nfound-1,aNpeak.GetArray(),aMean.GetArray(),xErr.GetArray(),aMeanErr.GetArray());

  mean->SetMarkerStyle(22);
  mean->SetTitle("HG Calibration");
  mean->GetXaxis()->SetTitle("#p.e.");
  mean->GetYaxis()->SetTitle("ADC channels");
  mean->Draw("ap");

  TF1* fgain = new TF1("fgain","pol1");

  mean->Fit("fgain","","same",0.,nfound-2);
  mean->GetYaxis()->SetRangeUser(0.,aMean[nfound-2]+200.);

  c2->cd(1)->Update();

  Double_t pedestal = fgain->GetParameter(0);
  Double_t gain = fgain->GetParameter(1);

  c2->cd(2);
  c2->GetPad(2)->SetGridx();
  c2->GetPad(2)->SetGridy();

  TArrayD aDpp(nfound-2);
  TArrayD aDppErr(nfound-2);
  Double_t mean_dpp = 0.;

  for(int a=0; a<nfound-2; ++a){

    aDpp[a] = aMean[a+1]-aMean[a];
    aDppErr[a] = aMeanErr[a+1]+aMeanErr[a];
    mean_dpp += aDpp[a];
  }

  mean_dpp /= (nfound-2);

  TGraphErrors* gDpp = new TGraphErrors(nfound-2,aNpeak.GetArray(),aDpp.GetArray(),xErr.GetArray(),aDppErr.GetArray());

  gDpp->SetMarkerStyle(22);
  gDpp->SetTitle("Delta peak to peak distribution");
  gDpp->GetXaxis()->SetTitle("#p.e.");
  gDpp->GetYaxis()->SetTitle("ADC channels");
  gDpp->Draw("ap");
  gDpp->Fit("pol1","","same",0.,nfound-3);
  gDpp->GetYaxis()->SetRangeUser(mean_dpp-50,mean_dpp+50.);

  c2->cd(3);
  c2->GetPad(3)->SetGridx();
  c2->GetPad(3)->SetGridy();

  TGraphErrors* gConst = new TGraphErrors(nfound-1,aNpeak.GetArray(),aConst.GetArray(),xErr.GetArray(),aConstErr.GetArray());

  gConst->SetMarkerStyle(22);
  gConst->SetTitle("Gaussian Normalization par distribution");
  gConst->GetXaxis()->SetTitle("#p.e.");
  gConst->GetYaxis()->SetTitle("#Entries");
  gConst->Draw("ap");


  TF1 *poisson = new TF1("poisson","[0]*TMath::Power(([1]/[2]),(x/[2]))*(TMath::Exp(-([1]/[2])))/TMath::Gamma((x/[2])+1.)", 0, nfound);

  poisson->SetParameters(1.e6, 9., 2.);
  poisson->SetNpx(4096*2);

  gConst->Fit(poisson,"R","same",0.,nfound-2);


  c2->cd(4);
  c2->GetPad(4)->SetGridx();
  c2->GetPad(4)->SetGridy();

  TGraphErrors* gSig = new TGraphErrors(nfound-1,aNpeak.GetArray(),aSig.GetArray(),xErr.GetArray(),aSigErr.GetArray());

  gSig->SetMarkerStyle(22);
  gSig->SetTitle("Gaussian Sigma par distribution");
  gSig->GetXaxis()->SetTitle("#p.e.");
  gSig->GetYaxis()->SetTitle("ADC channels");
  gSig->Draw("ap");
  gSig->GetYaxis()->SetRangeUser(0.,50.);

  TF1* pol1 = new TF1("pol1","pol1");

  pol1->SetParameters(20.,0.5);

  gSig->Fit("pol1","","same",0.,nfound-2);

  return 0;
}
1 Like

Your help is highly appreciated, thank you, really.

A few questions arise. It is true that graphically the fit is graphically perfect but looking at the fit parameters it seems that something has gone wrong:

Found 16 candidate peaks to fit
 FCN=105202 FROM MIGRAD    STATUS=CONVERGED   14374 CALLS       14375 TOTAL
                     EDM=2.93017e-07    STRATEGY= 1  ERROR MATRIX UNCERTAINTY   1.9 per cent
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  p0           9.47540e+03   2.19032e+01  -4.11213e-05** at limit **
   2  p1           1.02429e+02   1.45080e-02   1.40741e-06   2.55823e-01
   3  p2           2.50000e+00   3.37174e-04   1.61456e-06** at limit **
   4  p3           1.28264e+00   3.49022e-03  -2.88698e-06   1.62965e-03
   5  p4           1.90476e+04   1.90805e+00  -2.98346e-05** at limit **
   6  p5           1.60340e+02   1.66116e-02   1.01682e-07   3.11772e-03
   7  p6           5.91442e+00   2.05578e-02   4.09354e-07   1.54155e-02
   8  p7           1.18261e+00   3.63321e-03   1.41297e-06   1.38567e-02
   9  p8           4.77664e+04   1.40240e+00  -1.50920e-07** at limit **
  10  p9           2.31357e+02   1.01287e-02  -2.69824e-07  -4.78055e-02
  11  p10          5.06846e+00   1.08301e-02   9.70051e-08  -8.73756e-02
  12  p11          1.11403e+00   1.96513e-03   1.37306e-06  -3.88508e-02
  13  p12          7.11733e+04   2.03100e+00   9.25527e-06** at limit **
  14  p13          3.02345e+02   8.69980e-03  -7.51767e-08  -3.18896e-02
  15  p14          5.53092e+00   6.11178e-03   2.30727e-08  -1.07341e-01
  16  p15          1.10302e+00   8.35229e-04   2.79419e-06   5.29546e-02
  17  p16          8.29697e+04   1.75035e+00   2.61834e-06** at limit **
  18  p17          3.74988e+02   8.44828e-03  -8.71802e-08   4.74753e-02
  19  p18          6.02703e+00   7.11474e-03  -6.69980e-07  -4.64292e-01
  20  p19          1.10764e+00   1.01595e-03  -1.01879e-05   1.39439e-01
  21  p20          8.02835e+04   1.68127e+01   8.59465e-06** at limit **
  22  p21          4.47405e+02   9.27907e-03   9.67465e-08   9.76349e-02
  23  p22          6.61659e+00   1.05026e-02  -7.82051e-07  -1.52301e-01
  24  p23          1.13579e+00   1.63615e-03  -3.17363e-06  -3.16107e-02
  25  p24          6.67532e+04   1.21551e+02  -9.92028e-06   5.70403e-02
  26  p25          5.19250e+02   1.05850e-02   3.65760e-07   1.64454e-02
  27  p26          7.29914e+00   2.53526e-02   1.15507e-06   3.19413e-01
  28  p27          1.16073e+00   2.96816e-03   1.80163e-06  -7.47205e-02
  29  p28          4.87032e+04   9.23363e+01  -1.53043e-05  -2.68532e-04
  30  p29          5.90925e+02   1.33264e-02   7.95654e-08  -2.13725e-02
  31  p30          8.18717e+00   2.86863e-02   3.10540e-06   1.58467e-01
  32  p31          1.20280e+00   3.37619e-03   1.08631e-05  -5.71381e-02
  33  p32          3.21511e+04   6.50229e+01   1.22229e-06  -6.70512e-02
  34  p33          6.62495e+02   1.92062e-02  -4.15562e-07  -2.69225e-02
  35  p34          9.22962e+00   3.43608e-02  -1.73923e-07   1.56260e-01
  36  p35          1.25840e+00   4.20316e-03  -1.86229e-06   7.94666e-02
  37  p36          1.95294e+04   7.53299e+00   1.46239e-05** at limit **
  38  p37          7.33772e+02   2.36254e-02  -4.83796e-07  -9.10030e-02
  39  p38          1.02802e+01   2.67473e-02   2.97472e-07  -9.40134e-02
  40  p39          1.30847e+00   4.48459e-03   2.34653e-06   9.84635e-03
  41  p40          1.11882e+04   3.83210e+01   7.79975e-06   1.93249e-02
  42  p41          8.04638e+02   3.38542e-02  -1.34119e-06  -5.87357e-02
  43  p42          1.12209e+01   6.73357e-02  -1.10741e-06  -3.29220e-02
  44  p43          1.36827e+00   8.81434e-03  -5.39723e-06   9.79358e-04
  45  p44          6.16441e+03   2.17057e+01   5.29615e-06   1.01309e-02
  46  p45          8.75307e+02   4.82033e-02   9.93423e-09  -7.82219e-03
  47  p46          1.17251e+01   6.47382e-02  -1.44145e-06  -9.29764e-02
  48  p47          1.34955e+00   8.91423e-03  -2.52579e-06  -1.14385e-02
  49  p48          3.05249e+03   2.10481e+01  -6.04911e-06  -3.94238e-03
  50  p49          9.45655e+02   7.39211e-02   2.13036e-06   5.16220e-03
  51  p50          1.31712e+01   1.48457e-01   5.63265e-06  -3.37656e-02
  52  p51          1.47859e+00   2.26208e-02   1.39832e-05   1.00734e-02
  53  p52          1.52603e+03   1.51333e+01   3.76750e-05  -6.76617e-03
  54  p53          1.01593e+03   1.16509e-01   1.61320e-06   5.53068e-03
  55  p54          1.35772e+01   2.19277e-01   1.22562e-05  -1.70214e-02
  56  p55          1.40251e+00   3.65015e-02   3.01258e-05  -1.83848e-02
  57  p56          6.87811e+02   1.44462e+01   9.37282e-06  -2.39090e-03
  58  p57          1.08558e+03   2.08801e-01  -3.78648e-06   2.01303e-03
  59  p58          1.42689e+01   3.29219e-01   1.31926e-06   6.59303e-03
  60  p59          1.56284e+00   6.78403e-02   5.87985e-08  -2.02863e-03
  61  p60          3.34446e+02   9.68440e+00   1.00137e-04   1.03471e-04
  62  p61          1.15473e+03   3.31138e-01  -1.52573e-06  -2.87540e-03
  63  p62          1.47939e+01   7.14863e-01  -3.69000e-05   6.85185e-03
  64  p63          1.29355e+00   1.43411e-01   8.71889e-05  -3.43674e-03

Some parameters (e.g p9, p13, p17) have zero errors, which means that Root has given up because the parameter is at the limit of the range, and fixes it in the middle of the range, which unfortunately has no physical meaning and is reflected in the high chi-square/n.o.f (about 10^2). Other parameters are “at limits”.

Why this happen? if the fit is so perfect, I expect the chi-square to be just as perfect.

Thanks again for your big help.

Just remove the parameter limits and things are fine. In particular limits on height and position are not useful because you have from TSpectrum good starting values:

FCN=102635 FROM MIGRAD    STATUS=CONVERGED   20377 CALLS       20378 TOTAL
                     EDM=3.01172e-07    STRATEGY= 1  ERROR MATRIX UNCERTAINTY   1.2 per cent
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  p0           1.02541e+04   9.31222e+01   2.93468e-02  -4.00450e-06
   2  p1           1.02369e+02   1.56045e-02  -1.88708e-06  -9.60459e-04
   3  p2           2.08880e+00   4.05840e-02  -8.57820e-06  -8.86939e-03
   4  p3           1.13093e+00   1.47663e-02  -3.14096e-06  -1.90742e-02
   5  p4           1.98080e+04   6.90166e+01   4.96109e-02   4.37922e-06
   6  p5           1.60307e+02   1.62176e-02  -5.24457e-06  -9.25648e-03
   7  p6           5.62502e+00   3.74982e-02  -2.81694e-05   2.71375e-02
   8  p7           1.15677e+00   5.26796e-03  -3.22768e-06  -7.56588e-02
   9  p8           5.24787e+04   1.27003e+02  -3.03361e-02   9.13087e-08
  10  p9           2.31313e+02   9.33004e-03   5.11277e-07  -1.21506e-02
  11  p10          4.38126e+00   2.09624e-02   4.01907e-06   7.59486e-04
  12  p11          1.03412e+00   2.93985e-03   3.58288e-07  -2.21235e-02
  13  p12          7.42729e+04   1.41314e+02  -1.89999e-02   3.08412e-06
  14  p13          3.02322e+02   8.46518e-03  -8.79842e-07   6.41160e-03
  15  p14          5.19092e+00   1.91111e-02   7.90502e-07   4.33323e-02
  16  p15          1.07028e+00   2.44471e-03   3.15843e-07  -2.18970e-01
  17  p16          8.49393e+04   1.49020e+02  -5.99618e-02   1.72516e-06
  18  p17          3.74984e+02   8.33811e-03   2.29657e-06  -2.24234e-02
  19  p18          5.80169e+00   1.97632e-02   8.05751e-06  -3.28313e-03
  20  p19          1.08514e+00   2.37760e-03   1.01951e-06   2.63056e-02
  21  p20          8.05968e+04   1.37781e+02  -6.14525e-02  -8.81422e-07
  22  p21          4.47406e+02   9.37127e-03  -2.44565e-06   2.43563e-03
  23  p22          6.57894e+00   2.12416e-02   1.43575e-05  -7.23418e-02
  24  p23          1.13316e+00   2.53259e-03   1.16047e-06   5.32616e-01
  25  p24          6.67695e+04   1.21475e+02   1.83583e-02  -3.98217e-07
  26  p25          5.19249e+02   1.05880e-02   1.64783e-06   8.87442e-03
  27  p26          7.29428e+00   2.51076e-02   3.62715e-06  -1.97886e-02
  28  p27          1.15999e+00   2.90617e-03   1.37021e-07  -8.94262e-02
  29  p28          4.87018e+04   9.61724e+01  -3.38996e-02   2.92026e-06
  30  p29          5.90926e+02   1.31483e-02   7.29729e-06  -2.48014e-02
  31  p30          8.18759e+00   2.98714e-02   1.28891e-05   2.25773e-02
  32  p31          1.20286e+00   3.34548e-03   1.48882e-06   2.00579e-01
  33  p32          3.21473e+04   6.23603e+01  -3.39265e-03  -3.72801e-07
  34  p33          6.62491e+02   1.82036e-02  -5.95112e-06   7.59008e-03
  35  p34          9.23128e+00   2.92409e-02   1.96502e-05  -1.05662e-02
  36  p35          1.25903e+00   3.19849e-03   2.26718e-06  -6.72207e-01
  37  p36          1.96069e+04   5.75669e+01  -5.28892e-02  -2.48866e-08
  38  p37          7.33769e+02   2.25002e-02  -1.70649e-05   1.67897e-03
  39  p38          1.02151e+01   5.46677e-02   4.02298e-05   1.02870e-03
  40  p39          1.30150e+00   6.78401e-03   4.39096e-06   1.22155e-01
  41  p40          1.11834e+04   4.10681e+01  -1.29588e-02  -7.19223e-06
  42  p41          8.04644e+02   3.41669e-02   1.80084e-05  -4.23791e-03
  43  p42          1.12243e+01   6.84486e-02   3.45098e-05  -3.25288e-03
  44  p43          1.37001e+00   9.62359e-03   2.77547e-07   2.80403e-01
  45  p44          6.16709e+03   3.28097e+01  -2.74448e-02  -6.31651e-07
  46  p45          8.75303e+02   4.92328e-02  -2.02788e-05  -2.68062e-03
  47  p46          1.17185e+01   1.10730e-01   7.49818e-05  -4.07167e-03
  48  p47          1.34753e+00   1.51052e-02   1.45881e-05   3.29182e-02
  49  p48          3.05099e+03   2.22991e+01   3.83340e-04  -5.99980e-06
  50  p49          9.45659e+02   7.90957e-02   2.96576e-05  -1.55813e-03
  51  p50          1.31745e+01   1.53065e-01  -2.36602e-05  -4.78300e-03
  52  p51          1.48086e+00   2.43069e-02  -1.90073e-05   3.59031e-02
  53  p52          1.52680e+03   1.63800e+01  -1.08539e-02  -8.37003e-06
  54  p53          1.01593e+03   1.18765e-01  -7.44742e-05   1.84926e-03
  55  p54          1.35733e+01   2.43258e-01   8.14097e-05  -2.95641e-03
  56  p55          1.39966e+00   3.95722e-02   2.52129e-05   2.40325e-02
  57  p56          6.86886e+02   1.49710e+01   1.05628e-02  -2.97634e-05
  58  p57          1.08558e+03   2.09438e-01  -3.51877e-05   6.89297e-04
  59  p58          1.42634e+01   3.32508e-01   9.26532e-05  -1.27887e-03
  60  p59          1.56719e+00   6.97112e-02  -4.16264e-05   7.68160e-03
  61  p60          3.34848e+02   9.80470e+00  -4.87618e-03  -4.12954e-05
  62  p61          1.15472e+03   3.29875e-01   2.07314e-04  -6.31440e-04
  63  p62          1.47865e+01   7.32662e-01   2.84881e-04  -5.98310e-04
  64  p63          1.28515e+00   1.46152e-01   1.12130e-04   2.17236e-03

The HyperGaussians reduced the chi-squared roughly an order of magnitude but your high-statistics spectrum shows that the shape might need further improvements. Looking at the higher ADC channels, one can see that the peak shape is asymmetric, it is wider on the right side. One could incorporate this by introducing a fifth parameter in the peak shape, the asymmetry: x < peak position: sigma*(1-asym), x > peak position: sigma*(1+asym) .

You have parametrized each peak with its own set of parameters. In order to reduce the
number of fit parameters and stabilize/constrain the fit you might want to describe sigma, gamma and asym as polynomials as a function of the ADC channel.

1 Like

Here the code with the asymmetry, another 30% improvement in ch-sqaured

#include "TCanvas.h"
#include "TStyle.h"
#include "TGraphErrors.h"
#include "TArrayD.h"
#include "TH1F.h"
#include "TSpectrum.h"
#include "TF1.h"

const Int_t npeaks = 40;

Double_t HyperGaus(Double_t *x,Double_t *par)
{  
   Double_t mean  = par[0];
   Double_t sigma = par[1];
   Double_t gamma = par[2];
   Double_t asym  = par[3];
   
   if (sigma == 0) return 1.e30;
   Double_t arg = (x[0]-mean > 0) ? TMath::Abs((x[0]-mean)/sigma/(1-asym))
                                  : TMath::Abs((x[0]-mean)/sigma/(1+asym));
   Double_t powarg = TMath::Power(arg,gamma);
   Double_t res = TMath::Exp(-0.5*powarg);
   return res;
}

Int_t nfound;
Double_t SumHyperGaus(Double_t *x,Double_t *par) {

   Double_t result = 0.;
   for (Int_t p=0;p<nfound;p++) {

      Double_t norm  = par[5*p+0];
      result += norm*HyperGaus(x,par+5*p+1);
   }
   return result;
}

int Multiphoton_spectrum_analysis() {

  FILE *file = fopen("./Run0_PHA_HG_0_10.txt", "r");

  TCanvas* c1 = new TCanvas("c1","HG Multiphoton Spectrum",10,10,1920,1980);

  gStyle->SetOptFit(000);

  c1->SetGridx();
  c1->SetGridy();

  TH1F* htmp = new TH1F("htemp","htemp",4096,0,4096);

  Int_t bins = 0;
  Float_t x;
  while(!feof(file)){
    fscanf(file,"%g\n",&x);
    htmp->SetBinContent(bins,x);
    ++bins;
  }

  TH1F *hsipm = new TH1F("HG Multiphoton Spectrum","",bins,0,bins);

  for(int a=0; a<bins; ++a)

  hsipm->SetBinContent(a,htmp->GetBinContent(a));

  Int_t sigma = bins/1024+1;

  hsipm->GetYaxis()->SetTitle("#Entries");
  hsipm->GetXaxis()->SetTitle("ADC channels");;
  hsipm->GetXaxis()->SetRange(0,1300);

  TSpectrum *s = new TSpectrum(npeaks);

  nfound = s->Search(hsipm,sigma,"",0.0007);

  printf("Found %d candidate peaks to fit\n", nfound);

  int peak_x_tmp[nfound];
  double peak_x[nfound];
  double peak_y[nfound];

  TMath::Sort(nfound,s->GetPositionX(),peak_x_tmp,kFALSE);

  for(int a=0; a<nfound; ++a){

    peak_x[a] = s->GetPositionX()[peak_x_tmp[a]];
    peak_y[a] = s->GetPositionY()[peak_x_tmp[a]];

  }

  TF1* f = new TF1("f",SumHyperGaus,0,4096,5*nfound);
  f->SetNpx(4096*2);
  f->SetLineStyle(9);
  f->SetLineColor(2); //colore fit generale

  for(int a=0; a<nfound; ++a){
    f->SetParameter(5*a+0,peak_y[a]);
    f->SetParameter(5*a+1,peak_x[a]);
    f->SetParameter(5*a+2,sigma);
    f->SetParameter(5*a+3,1.2);
    f->SetParameter(5*a+4,0.0);
  }

  for(int a=2; a<nfound; ++a){

    //f->SetParLimits(5*a+0,peak_y[a]*0.5,peak_y[a]*1.1);
    //f->SetParLimits(5*a+1,peak_x[a]-3*sigma, peak_x[a]+3*sigma);
    //f->SetParLimits(5*a+2,0.5*sigma,6*sigma);
    //f->SetParLimits(5*a+3,1.1,2.5);

  }

  // pedestal

  //f->SetParLimits(0,peak_y[0]*0.6,peak_y[0]*1.1);
  //f->SetParLimits(1,peak_x[0]-1.2*sigma, peak_x[0]+1.2*sigma);
  //f->SetParLimits(2,0.5*sigma,1*sigma);
  //f->SetParLimits(3,1.1,2.5);

  // 1 p.e.

  //f->SetParLimits(4,peak_y[1]*0.5,peak_y[1]*1.1);
  //f->SetParLimits(5,peak_x[1]-1.5*sigma, peak_x[1]+1.5*sigma);
  //f->SetParLimits(6,0.5*sigma,4*sigma);
  //f->SetParLimits(7,1.1,2.5);

  // option I makes sure the area under fitted curve is fitted, but this
  // takes a significant amount of CPU time
  //hsipm->Fit("f","I","same",35.,peak_x[nfound-1]+4*sigma);
  hsipm->Fit("f","","same",35.,peak_x[nfound-1]+4*sigma);
  hsipm->SetLineColor(1); //colore dei dati

  TArrayD fPar(5*nfound);
  TArrayD fErr(5*nfound);

  for(int a=0; a<nfound; ++a) {
    fPar[5*a+0] = f->GetParameter(5*a+0);
    fPar[5*a+1] = f->GetParameter(5*a+1);
    fPar[5*a+2] = f->GetParameter(5*a+2);
    fPar[5*a+3] = f->GetParameter(5*a+3);
    fPar[5*a+4] = f->GetParameter(5*a+4);
    fErr[5*a+0] = f->GetParError(5*a+0);
    fErr[5*a+1] = f->GetParError(5*a+1);
    fErr[5*a+2] = f->GetParError(5*a+2);
    fErr[5*a+3] = f->GetParError(5*a+3);
    fErr[5*a+4] = f->GetParError(5*a+4);
  }

  TF1* fg[npeaks];

  double m_tmp = fPar[1];
  double s_tmp = fPar[2];

  fg[0] = new TF1("ped",HyperGaus,m_tmp-3*s_tmp,m_tmp+3*s_tmp);
  fg[0]->SetParameters(&fPar[0]);
  fg[0]->SetLineColor(9); //colore piedistallo
  //fg[0]->Draw("same");

  for(int a=1; a<nfound; ++a){

    m_tmp = fPar[5*a+1];
    s_tmp = fPar[5*a+2];

    TString gname;
    gname.Form("%d p.e.",a);

    fg[a] = new TF1(gname,HyperGaus,m_tmp-3*s_tmp,m_tmp+3*s_tmp);
    fg[a]->SetParameters(&fPar[5*a]);
    fg[a]->SetLineColor(9); //colore gaussiane
    //fg[a]->Draw("same");

  }

   TArrayD aNpeak(nfound-1);
   TArrayD aConst(nfound-1);
   TArrayD aMean (nfound-1);
   TArrayD aSig  (nfound-1);
   TArrayD aGamma(nfound-1);
   TArrayD aASymm(nfound-1);

   TArrayD aConstErr(nfound-1);
   TArrayD aMeanErr (nfound-1);
   TArrayD aSigErr  (nfound-1);
   TArrayD aGammaErr(nfound-1);
   TArrayD aASymmErr(nfound-1);

   const double sqrt2pi = sqrt(2 * acos(-1.));

   for (int a = 1; a < nfound; ++a)
   {
     aNpeak   [a-1] = a;
     aMean    [a-1] = fPar[5*a+1];
     aMeanErr [a-1] = fErr[5*a+1];
     aSig     [a-1] = fPar[5*a+2];
     aSigErr  [a-1] = fErr[5*a+2];
     aGamma   [a-1] = fPar[5*a+3];
     aGammaErr[a-1] = fErr[5*a+3];
     aASymm   [a-1] = fPar[5*a+4];
     aASymmErr[a-1] = fErr[5*a+4];

     aConst   [a-1] = fPar[5*a]*sqrt2pi*aSig[a-1];
     aConstErr[a-1] = sqrt2pi*(fErr[5*a]*aSig[a-1]+aSigErr[a-1]*fPar[5*a]);
  }

  FILE* outfile;

  outfile = fopen("Gaussian_Parameters.txt", "w");

  for(int a=0; a<nfound; ++a)
    fprintf(outfile, "%d \t %f \t %f \t %f \t %f \t %f \t %f \t %f \t %f \n",
                       a,fPar[5*a],fErr[5*a],fPar[5*a+1],fErr[5*a+1],fPar[5*a+2],fErr[5*a+2],fPar[5*a+3],fErr[5*a+3]);
  fclose(outfile);

  TCanvas* c2 = new TCanvas("c2","HG Multiphoton Spectrum Analysis",10,10,1920,1980);

  c2->UseCurrentStyle();
  gStyle->SetOptFit(111);

  c2->Divide(2,2);
  c2->cd(1);
  c2->GetPad(1)->SetGridx();
  c2->GetPad(1)->SetGridy();
  c2->Update();

  TArrayD xErr(npeaks);

  TGraphErrors* mean = new TGraphErrors(nfound-1,aNpeak.GetArray(),aMean.GetArray(),xErr.GetArray(),aMeanErr.GetArray());

  mean->SetMarkerStyle(22);
  mean->SetTitle("HG Calibration");
  mean->GetXaxis()->SetTitle("#p.e.");
  mean->GetYaxis()->SetTitle("ADC channels");
  mean->Draw("ap");

  TF1* fgain = new TF1("fgain","pol1");

  mean->Fit("fgain","","same",0.,nfound-2);
  mean->GetYaxis()->SetRangeUser(0.,aMean[nfound-2]+200.);

  c2->cd(1)->Update();

  Double_t pedestal = fgain->GetParameter(0);
  Double_t gain = fgain->GetParameter(1);

  c2->cd(2);
  c2->GetPad(2)->SetGridx();
  c2->GetPad(2)->SetGridy();

  TArrayD aDpp(nfound-2);
  TArrayD aDppErr(nfound-2);
  Double_t mean_dpp = 0.;

  for(int a=0; a<nfound-2; ++a){

    aDpp[a] = aMean[a+1]-aMean[a];
    aDppErr[a] = aMeanErr[a+1]+aMeanErr[a];
    mean_dpp += aDpp[a];
  }

  mean_dpp /= (nfound-2);

  TGraphErrors* gDpp = new TGraphErrors(nfound-2,aNpeak.GetArray(),aDpp.GetArray(),xErr.GetArray(),aDppErr.GetArray());

  gDpp->SetMarkerStyle(22);
  gDpp->SetTitle("Delta peak to peak distribution");
  gDpp->GetXaxis()->SetTitle("#p.e.");
  gDpp->GetYaxis()->SetTitle("ADC channels");
  gDpp->Draw("ap");
  gDpp->Fit("pol1","","same",0.,nfound-3);
  gDpp->GetYaxis()->SetRangeUser(mean_dpp-50,mean_dpp+50.);


  c2->cd(3);
  c2->GetPad(3)->SetGridx();
  c2->GetPad(3)->SetGridy();

  TGraphErrors* gConst = new TGraphErrors(nfound-1,aNpeak.GetArray(),aConst.GetArray(),xErr.GetArray(),aConstErr.GetArray());

  gConst->SetMarkerStyle(22);
  gConst->SetTitle("Gaussian Normalization par distribution");
  gConst->GetXaxis()->SetTitle("#p.e.");
  gConst->GetYaxis()->SetTitle("#Entries");
  gConst->Draw("ap");


  TF1 *poisson = new TF1("poisson","[0]*TMath::Power(([1]/[2]),(x/[2]))*(TMath::Exp(-([1]/[2])))/TMath::Gamma((x/[2])+1.)", 0, nfound);

  poisson->SetParameters(1.e6, 9., 2.);
  poisson->SetNpx(4096*2);

  gConst->Fit(poisson,"R","same",0.,nfound-2);


  c2->cd(4);
  c2->GetPad(4)->SetGridx();
  c2->GetPad(4)->SetGridy();

  TGraphErrors* gSig = new TGraphErrors(nfound-1,aNpeak.GetArray(),aSig.GetArray(),xErr.GetArray(),aSigErr.GetArray());

  gSig->SetMarkerStyle(22);
  gSig->SetTitle("Gaussian Sigma par distribution");
  gSig->GetXaxis()->SetTitle("#p.e.");
  gSig->GetYaxis()->SetTitle("ADC channels");
  gSig->Draw("ap");
  gSig->GetYaxis()->SetRangeUser(0.,50.);

  TF1* pol1 = new TF1("pol1","pol1");

  pol1->SetParameters(20.,0.5);

  gSig->Fit("pol1","","same",0.,nfound-2);

  c2->cd(5);
  c2->GetPad(5)->SetGridx();
  c2->GetPad(5)->SetGridy();
  TGraphErrors* gGamma = new TGraphErrors(nfound-1,aNpeak.GetArray(),aGamma.GetArray(),xErr.GetArray(),aGammaErr.GetArray());
  gGamma->SetMarkerStyle(22);
  gGamma->SetTitle("Gaussian Gamma par distribution");
  gGamma->GetXaxis()->SetTitle("#p.e.");
  gGamma->GetYaxis()->SetTitle("ADC channels");
  gGamma->Draw("ap");

  c2->cd(6);
  c2->GetPad(6)->SetGridx();
  c2->GetPad(6)->SetGridy();
  TGraphErrors* gASymm = new TGraphErrors(nfound-1,aNpeak.GetArray(),aASymm.GetArray(),xErr.GetArray(),aASymmErr.GetArray());
  gASymm->SetMarkerStyle(22);
  gASymm->SetTitle("Gaussian ASymm par distribution");
  gASymm->GetXaxis()->SetTitle("#p.e.");
  gASymm->GetYaxis()->SetTitle("ADC channels");
  gASymm->Draw("ap");

  return 0;
}
1 Like

What you wrote is absolutely true. There is a problem of asymmetry in spectra placed at a high number of ADC channels. This behaviour seems to be caused by the experimental setup and we are trying very hard to minimize and even abate it.

We can now say that the curve fits the data “well”, although it might be interesting to apply this asymmetry parameter to a multi-Gaussian fit, what do you think?

The use of Hypergaussians is interesting but it distorts the shape of the peak, as I show here in the figure:

Just put parameter limits on the parameter gamma:

f->SetParameter(5*a+3,2.0);
f->SetParLimits(5*a+3,2.0,2.0);

and you fit a normal Gaussian. The peak you show is the one with gamma = 1.009 !!

Realize with HyperGaussians the following:

  • gamma > 2: the peak shape gets squarer
  • gamma < 2: peak gets sharper and the tails “fatter”

Your peaks seem to have fatter tail on the left side than a Gaussian peak, the HyperGaussian is compensating for that with gamma < 2 but therefore makes it too sharp in the center. So one of the things you could try is to add to the peaks description some exp(-x)
tail on the left side starting at a few sigma from the center. Then the fitted gamma will be
much closer to 2.

In the mean time I see a reasonable fit with gamma=1.2 (fixed)

1 Like

Now I understande perfectly what you mean. Fantastic!

Just the last question, probably is my fault but:

Maybe i put something wronge in the code and the other fit get non sense and:

#0  0x00007f8c0d5f4dba in __GI___wait4 (pid=5949, stat_loc=stat_loc
entry=0x7ffc0ab25168, options=options
entry=0, usage=usage
entry=0x0) at ../sysdeps/unix/sysv/linux/wait4.c:27
#1  0x00007f8c0d5f4d7b in __GI___waitpid (pid=<optimized out>, stat_loc=stat_loc
entry=0x7ffc0ab25168, options=options
entry=0) at waitpid.c:38
#2  0x00007f8c0d5640e7 in do_system (line=<optimized out>) at ../sysdeps/posix/system.c:172
#3  0x00007f8c0dc31a3e in TUnixSystem::StackTrace() () from /home/mllsml92e13c351e/ROOT6_24_02/lib/libCore.so
#4  0x00007f8c09223480 in cling::MultiplexInterpreterCallbacks::PrintStackTrace() () from /home/mllsml92e13c351e/ROOT6_24_02/lib/libCling.so
#5  0x00007f8c09216bc2 in cling_runtime_internal_throwIfInvalidPointer () from /home/mllsml92e13c351e/ROOT6_24_02/lib/libCling.so
#6  0x00007f8bf7cfba72 in ?? ()
#7  0x40c0d30000000000 in ?? ()
#8  0x40d0e90000000000 in ?? ()
#9  0x40e5340000000000 in ?? ()
#10 0x40ef97e000000000 in ?? ()
#11 0x40f26a3000000000 in ?? ()
#12 0x40f1d19000000000 in ?? ()
#13 0x40ee07c000000000 in ?? ()
#14 0x40e6896000000000 in ?? ()
#15 0x40dd728000000000 in ?? ()
#16 0x40d1568000000000 in ?? ()
#17 0x40c4c08000000000 in ?? ()
#18 0x40b7e50000000000 in ?? ()
#19 0x40a8080000000000 in ?? ()
#20 0x4096e00000000000 in ?? ()
#21 0x4085580000000000 in ?? ()
#22 0x4074700000000000 in ?? ()
#23 0x4059a00000000000 in ?? ()
#24 0x4064100000000000 in ?? ()
#25 0x406cf00000000000 in ?? ()
#26 0x4072e80000000000 in ?? ()
#27 0x4077680000000000 in ?? ()
#28 0x407be80000000000 in ?? ()
#29 0x4080340000000000 in ?? ()
#30 0x4082740000000000 in ?? ()
#31 0x4084ac0000000000 in ?? ()
#32 0x4086e40000000000 in ?? ()
#33 0x40891c0000000000 in ?? ()
#34 0x408b540000000000 in ?? ()
#35 0x408d840000000000 in ?? ()
#36 0x408fb40000000000 in ?? ()
#37 0x4090f20000000000 in ?? ()
#38 0x4092060000000000 in ?? ()
#39 0x0000000700000009 in ?? ()
#40 0x0000000200000004 in ?? ()
#41 0x0000000100000000 in ?? ()
#42 0x0000000500000003 in ?? ()
#43 0x0000000800000006 in ?? ()
#44 0x0000000b0000000a in ?? ()
#45 0x0000000d0000000c in ?? ()
#46 0x0000000f0000000e in ?? ()
#47 0x00007f8bf97254d0 in ?? () from /home/mllsml92e13c351e/ROOT6_24_02/lib/libHist.so
#48 0x00007ffc0ab27aa0 in ?? ()
#49 0x0000563b5ab33af8 in ?? ()
#50 0x0000000000000059 in ?? ()
#51 0x0000000000000000 in ?? ()
Error in <HandleInterpreterException>: Trying to dereference null pointer or trying to call routine taking non-null arguments.
Execution of your code was aborted.
In file included from input_line_104:1:
/home/mllsml92e13c351e/MEGA/SSD/Università/Dottorato di Ricerca in Fisica/CAEN Front End Readout System 5200 /Test FERS 5200/Test SiPM Hamamatsu S12, S13, S14 e Ketek/Test sul Hamamatsu S13/Run 23_02_2022/Multiphoton_spectrum_analysis.C:314:3: warning: null passed to a callee that requires a non-null argument [-Wnonnull]
  c2->GetPad(5)->SetGridx();

something of bad is written in my terminal.

Can you share with us the final code? Do you see this strange line in your terminal?

The code where I fixed gamma on 1.2:

#include "TCanvas.h"
#include "TStyle.h"
#include "TGraphErrors.h"
#include "TArrayD.h"
#include "TH1F.h"
#include "TSpectrum.h"
#include "TF1.h"

const Int_t npeaks = 40;

Double_t HyperGaus(Double_t *x,Double_t *par)
{  
   Double_t mean  = par[0];
   Double_t sigma = par[1];
   Double_t gamma = par[2];
   Double_t asym  = par[3];
   
   if (sigma == 0) return 1.e30;
   Double_t arg = (x[0]-mean > 0) ? TMath::Abs((x[0]-mean)/sigma/(1-asym))
                                  : TMath::Abs((x[0]-mean)/sigma/(1+asym));
   Double_t powarg = TMath::Power(arg,gamma);
   Double_t res = TMath::Exp(-0.5*powarg);
   return res;
}

Int_t nfound;
Double_t SumHyperGaus(Double_t *x,Double_t *par) {

   Double_t result = 0.;
   for (Int_t p=0;p<nfound;p++) {

      Double_t norm  = par[5*p+0];
      result += norm*HyperGaus(x,par+5*p+1);
   }
   return result;
}

int Multiphoton_spectrum_analysis() {

  FILE *file = fopen("./Run0_PHA_HG_0_10.txt", "r");

  TCanvas* c1 = new TCanvas("c1","HG Multiphoton Spectrum",10,10,1920,1980);

  gStyle->SetOptFit(000);

  c1->SetGridx();
  c1->SetGridy();

  TH1F* htmp = new TH1F("htemp","htemp",4096,0,4096);

  Int_t bins = 0;
  Float_t x;
  while(!feof(file)){
    fscanf(file,"%g\n",&x);
    htmp->SetBinContent(bins,x);
    ++bins;
  }

  TH1F *hsipm = new TH1F("HG Multiphoton Spectrum","",bins,0,bins);

  for(int a=0; a<bins; ++a)

  hsipm->SetBinContent(a,htmp->GetBinContent(a));

  Int_t sigma = bins/1024+1;

  hsipm->GetYaxis()->SetTitle("#Entries");
  hsipm->GetXaxis()->SetTitle("ADC channels");;
  hsipm->GetXaxis()->SetRange(0,1300);

  TSpectrum *s = new TSpectrum(npeaks);

  nfound = s->Search(hsipm,sigma,"",0.0007);

  printf("Found %d candidate peaks to fit\n", nfound);

  int peak_x_tmp[nfound];
  double peak_x[nfound];
  double peak_y[nfound];

  TMath::Sort(nfound,s->GetPositionX(),peak_x_tmp,kFALSE);

  for(int a=0; a<nfound; ++a){

    peak_x[a] = s->GetPositionX()[peak_x_tmp[a]];
    peak_y[a] = s->GetPositionY()[peak_x_tmp[a]];

  }

  TF1* f = new TF1("f",SumHyperGaus,0,4096,5*nfound);
  f->SetNpx(4096*2);
  f->SetLineStyle(9);
  f->SetLineColor(2); //colore fit generale

  for(int a=0; a<nfound; ++a){
    f->SetParameter(5*a+0,peak_y[a]);
    f->SetParameter(5*a+1,peak_x[a]);
    f->SetParameter(5*a+2,sigma);
    f->SetParameter(5*a+3,1.2);
    f->SetParameter(5*a+4,0.0);
  }

  for(int a=2; a<nfound; ++a){

    //f->SetParLimits(5*a+0,peak_y[a]*0.5,peak_y[a]*1.1);
    //f->SetParLimits(5*a+1,peak_x[a]-3*sigma, peak_x[a]+3*sigma);
    //f->SetParLimits(5*a+2,0.5*sigma,6*sigma);
    f->SetParLimits(5*a+3,1.2,1.2);

  }

  // pedestal

  //f->SetParLimits(0,peak_y[0]*0.6,peak_y[0]*1.1);
  //f->SetParLimits(1,peak_x[0]-1.2*sigma, peak_x[0]+1.2*sigma);
  //f->SetParLimits(2,0.5*sigma,1*sigma);
  f->SetParLimits(3,1.2,1.2);

  // 1 p.e.

  //f->SetParLimits(5,peak_y[1]*0.5,peak_y[1]*1.1);
  //f->SetParLimits(6,peak_x[1]-1.5*sigma, peak_x[1]+1.5*sigma);
  //f->SetParLimits(7,0.5*sigma,4*sigma);
  f->SetParLimits(8,1.2,1.2);

  // option I makes sure the area under fitted curve is fitted, but this
  // takes a significant amount of CPU time
  //hsipm->Fit("f","I","same",35.,peak_x[nfound-1]+4*sigma);
  hsipm->Fit("f","","same",35.,peak_x[nfound-1]+4*sigma);
  hsipm->SetLineColor(1); //colore dei dati

  TArrayD fPar(5*nfound);
  TArrayD fErr(5*nfound);

  for(int a=0; a<nfound; ++a) {
    fPar[5*a+0] = f->GetParameter(5*a+0);
    fPar[5*a+1] = f->GetParameter(5*a+1);
    fPar[5*a+2] = f->GetParameter(5*a+2);
    fPar[5*a+3] = f->GetParameter(5*a+3);
    fPar[5*a+4] = f->GetParameter(5*a+4);
    fErr[5*a+0] = f->GetParError(5*a+0);
    fErr[5*a+1] = f->GetParError(5*a+1);
    fErr[5*a+2] = f->GetParError(5*a+2);
    fErr[5*a+3] = f->GetParError(5*a+3);
    fErr[5*a+4] = f->GetParError(5*a+4);
  }

  TF1* fg[npeaks];

  double m_tmp = fPar[1];
  double s_tmp = fPar[2];

  fg[0] = new TF1("ped",HyperGaus,m_tmp-3*s_tmp,m_tmp+3*s_tmp);
  fg[0]->SetParameters(&fPar[0]);
  fg[0]->SetLineColor(9); //colore piedistallo
  //fg[0]->Draw("same");

  for(int a=1; a<nfound; ++a){

    m_tmp = fPar[5*a+1];
    s_tmp = fPar[5*a+2];

    TString gname;
    gname.Form("%d p.e.",a);

    fg[a] = new TF1(gname,HyperGaus,m_tmp-3*s_tmp,m_tmp+3*s_tmp);
    fg[a]->SetParameters(&fPar[5*a]);
    fg[a]->SetLineColor(9); //colore gaussiane
    //fg[a]->Draw("same");

  }

   TArrayD aNpeak(nfound-1);
   TArrayD aConst(nfound-1);
   TArrayD aMean (nfound-1);
   TArrayD aSig  (nfound-1);
   TArrayD aGamma(nfound-1);
   TArrayD aASymm(nfound-1);

   TArrayD aConstErr(nfound-1);
   TArrayD aMeanErr (nfound-1);
   TArrayD aSigErr  (nfound-1);
   TArrayD aGammaErr(nfound-1);
   TArrayD aASymmErr(nfound-1);

   const double sqrt2pi = sqrt(2 * acos(-1.));

   for (int a = 1; a < nfound; ++a)
   {
     aNpeak   [a-1] = a;
     aMean    [a-1] = fPar[5*a+1];
     aMeanErr [a-1] = fErr[5*a+1];
     aSig     [a-1] = fPar[5*a+2];
     aSigErr  [a-1] = fErr[5*a+2];
     aGamma   [a-1] = fPar[5*a+3];
     aGammaErr[a-1] = fErr[5*a+3];
     aASymm   [a-1] = fPar[5*a+4];
     aASymmErr[a-1] = fErr[5*a+4];

     aConst   [a-1] = fPar[5*a]*sqrt2pi*aSig[a-1];
     aConstErr[a-1] = sqrt2pi*(fErr[5*a]*aSig[a-1]+aSigErr[a-1]*fPar[5*a]);

  }

  FILE* outfile;

  outfile = fopen("Gaussian_Parameters.txt", "w");

  for(int a=0; a<nfound; ++a)
    fprintf(outfile, "%d \t %f \t %f \t %f \t %f \t %f \t %f \t %f \t %f \n",
                       a,fPar[5*a],fErr[5*a],fPar[5*a+1],fErr[5*a+1],fPar[5*a+2],fErr[5*a+2],fPar[5*a+3],fErr[5*a+3]);
  fclose(outfile);

  TCanvas* c2 = new TCanvas("c2","HG Multiphoton Spectrum Analysis",10,10,1920,1980);

  c2->UseCurrentStyle();
  gStyle->SetOptFit(111);

  c2->Divide(2,3);
  c2->cd(1);
  c2->GetPad(1)->SetGridx();
  c2->GetPad(1)->SetGridy();
  c2->Update();

  TArrayD xErr(npeaks);

  TGraphErrors* mean = new TGraphErrors(nfound-1,aNpeak.GetArray(),aMean.GetArray(),xErr.GetArray(),aMeanErr.GetArray());

  mean->SetMarkerStyle(22);
  mean->SetTitle("HG Calibration");
  mean->GetXaxis()->SetTitle("#p.e.");
  mean->GetYaxis()->SetTitle("ADC channels");
  mean->Draw("ap");

  TF1* fgain = new TF1("fgain","pol1");

  mean->Fit("fgain","","same",0.,nfound-2);
  mean->GetYaxis()->SetRangeUser(0.,aMean[nfound-2]+200.);

  c2->cd(1)->Update();

  Double_t pedestal = fgain->GetParameter(0);
  Double_t gain = fgain->GetParameter(1);

  c2->cd(2);
  c2->GetPad(2)->SetGridx();
  c2->GetPad(2)->SetGridy();

  TArrayD aDpp(nfound-2);
  TArrayD aDppErr(nfound-2);
  Double_t mean_dpp = 0.;

  for(int a=0; a<nfound-2; ++a){

    aDpp[a] = aMean[a+1]-aMean[a];
    aDppErr[a] = aMeanErr[a+1]+aMeanErr[a];
    mean_dpp += aDpp[a];
  }

  mean_dpp /= (nfound-2);

  TGraphErrors* gDpp = new TGraphErrors(nfound-2,aNpeak.GetArray(),aDpp.GetArray(),xErr.GetArray(),aDppErr.GetArray());

  gDpp->SetMarkerStyle(22);
  gDpp->SetTitle("Delta peak to peak distribution");
  gDpp->GetXaxis()->SetTitle("#p.e.");
  gDpp->GetYaxis()->SetTitle("ADC channels");
  gDpp->Draw("ap");
  gDpp->Fit("pol1","","same",0.,nfound-3);
  gDpp->GetYaxis()->SetRangeUser(mean_dpp-50,mean_dpp+50.);


  c2->cd(3);
  c2->GetPad(3)->SetGridx();
  c2->GetPad(3)->SetGridy();

  TGraphErrors* gConst = new TGraphErrors(nfound-1,aNpeak.GetArray(),aConst.GetArray(),xErr.GetArray(),aConstErr.GetArray());

  gConst->SetMarkerStyle(22);
  gConst->SetTitle("Gaussian Normalization par distribution");
  gConst->GetXaxis()->SetTitle("#p.e.");
  gConst->GetYaxis()->SetTitle("#Entries");
  gConst->Draw("ap");


  TF1 *poisson = new TF1("poisson","[0]*TMath::Power(([1]/[2]),(x/[2]))*(TMath::Exp(-([1]/[2])))/TMath::Gamma((x/[2])+1.)", 0, nfound);

  poisson->SetParameters(1.e6, 9., 2.);
  poisson->SetNpx(4096*2);

  gConst->Fit(poisson,"R","same",0.,nfound-2);


  c2->cd(4);
  c2->GetPad(4)->SetGridx();
  c2->GetPad(4)->SetGridy();

  TGraphErrors* gSig = new TGraphErrors(nfound-1,aNpeak.GetArray(),aSig.GetArray(),xErr.GetArray(),aSigErr.GetArray());

  gSig->SetMarkerStyle(22);
  gSig->SetTitle("Gaussian Sigma par distribution");
  gSig->GetXaxis()->SetTitle("#p.e.");
  gSig->GetYaxis()->SetTitle("ADC channels");
  gSig->Draw("ap");
  gSig->GetYaxis()->SetRangeUser(0.,50.);

  TF1* pol1 = new TF1("pol1","pol1");

  pol1->SetParameters(20.,0.5);

  gSig->Fit("pol1","","same",0.,nfound-2);

  c2->cd(5);
  c2->GetPad(5)->SetGridx();
  c2->GetPad(5)->SetGridy();
  TGraphErrors* gGamma = new TGraphErrors(nfound-1,aNpeak.GetArray(),aGamma.GetArray(),xErr.GetArray(),aGammaErr.GetArray());
  gGamma->SetMarkerStyle(22);
  gGamma->SetTitle("Gaussian Gamma par distribution");
  gGamma->GetXaxis()->SetTitle("#p.e.");
  gGamma->GetYaxis()->SetTitle("ADC channels");
  gGamma->Draw("ap");

  c2->cd(6);
  c2->GetPad(6)->SetGridx();
  c2->GetPad(6)->SetGridy();
  TGraphErrors* gASymm = new TGraphErrors(nfound-1,aNpeak.GetArray(),aASymm.GetArray(),xErr.GetArray(),aASymmErr.GetArray());
  gASymm->SetMarkerStyle(22);
  gASymm->SetTitle("Gaussian ASymm par distribution");
  gASymm->GetXaxis()->SetTitle("#p.e.");
  gASymm->GetYaxis()->SetTitle("ADC channels");
  gASymm->Draw("ap");
}
1 Like