Facing Problem while fitting

Wile_E_Coyote · March 7, 2020, 5:33pm

arorastacy · March 7, 2020, 7:10pm

  Int_t n = g->GetN();
  Double_t *x = new Double_t[(n + 1)];
  x[0] = (3.0 * g->GetX()[0] - g->GetX()[1]) / 2.0;
  x[n] = (3.0 * g->GetX()[(n - 1)] - g->GetX()[(n - 2)]) / 2.0;
  for (Int_t i = 1; i < n; i++)
    { x[i] = (g->GetX()[(i - 1)] + g->GetX()[i]) / 2.0; }
  TH1D *h = new TH1D("h", "spectrum;x-value;#counts", n, x);
  h->FillN(n, g->GetX(), g->GetY());
  h->Sumw2(kFALSE); // restore proper errors
  delete [] x; // no longer needed
  delete g; // no longer needed

Please tell me about this portion.

arorastacy · March 7, 2020, 7:49pm

Also tell me

how to find area under all the peaks in the histogram?
how to locate centroid in each peak by writing a code?

Wile_E_Coyote · March 7, 2020, 8:37pm

The first parameter of “gaus” is the “peak’s height” while the first parameter of “gausn” is the “peak’s area”.
In both cases, the second parameter is the “peak’s mean” (i.e. the “peak’s position”).

{
  TGraph *g = new TGraph("eu_gu_1.dat");
  // g->Sort(); // just a precaution
  TF1 *f1 = new TF1("f1", "gausn(0) + pol1(3)", 400., 435.); // linear background
  f1->SetParameters(14000., 415., 2., 0., 0.);
  g->Fit("f1", "BR+");
  TF1 *f2 = new TF1("f2", "gausn(0) + pol1(3)", 1160., 1200.);
  f2->SetParameters(15000., 1178., 2., 0., 0.);
  g->Fit("f2", "BR+");
  TF1 *f3 = new TF1("f3", "gausn", 4795., 4835.); // clean peak, no background
  f3->SetParameters(4400., 4815., 3.);
  g->Fit("f3", "BR+");
  g->Draw("AL");
}

{
  TGraph *g = new TGraph("eu_gu_1.dat");
  // g->Sort(); // just a precaution
  Int_t n = g->GetN();
  Double_t *x = new Double_t[(n + 1)];
  x[0] = (3.0 * g->GetX()[0] - g->GetX()[1]) / 2.0;
  x[n] = (3.0 * g->GetX()[(n - 1)] - g->GetX()[(n - 2)]) / 2.0;
  for (Int_t i = 1; i < n; i++)
    { x[i] = (g->GetX()[(i - 1)] + g->GetX()[i]) / 2.0; }
  TH1D *h = new TH1D("h", "spectrum;x-value;#counts", n, x);
  h->FillN(n, g->GetX(), g->GetY());
  h->Sumw2(kFALSE); // restore proper errors
  delete [] x; // no longer needed
  delete g; // no longer needed
  TF1 *f1 = new TF1("f1", "gausn(0) + pol1(3)", 400., 435.); // linear background
  f1->SetParameters(14000., 415., 2., 0., 0.);
  h->Fit("f1", "BR+"); // e.g. "BR+" or "LBR+"
  TF1 *f2 = new TF1("f2", "gausn(0) + pol1(3)", 1160., 1200.);
  f2->SetParameters(15000., 1178., 2., 0., 0.);
  h->Fit("f2", "BR+"); // e.g. "BR+" or "LBR+"
  TF1 *f3 = new TF1("f3", "gausn", 4795., 4835.); // clean peak, no background
  f3->SetParameters(4400., 4815., 3.);
  h->Fit("f3", "BR+"); // e.g. "BR+" or "LBR+"
}

arorastacy · March 8, 2020, 6:56am

How do you get this area?
I have to calculate it for a number of peaks from the code in the program. How to do that?

Wile_E_Coyote · March 8, 2020, 7:09am

For the definitions of the “gaus”, “gausn” and “polN” functions see, for example, the underlaying: TFormula

You can retrieve the current values of function’s parameters and their errors (e.g. after fitting) using: TF1::GetParameter and TF1::GetParError

arorastacy · March 8, 2020, 7:46am

As you have mentioned earlier that for first peak area = 14000 how do you get that?
How can i use below code in my program.
https://root.cern/doc/master/classTH1.html#aaf053a4b487c46b9e20c0bf534b34012

Wile_E_Coyote · March 8, 2020, 8:17am

You need to estimate initial function’s parameters looking at the spectrum, e.g. (estimated_sigma = estimated_fwhm / 2.35): estimated_area = 2.51 * estimated_height * estimated_sigma = 1.06 * estimated_height * estimated_fwhm

arorastacy · March 8, 2020, 9:57am

Analytically i understood how am i getting values of the area under the peaks.

But when i try to retrieve its value from the original code i got an error:

 Double_t TH1::Integral(Int_t binx1 = 400., Int_t binx2 = 435., Option_t *option) const
 {
    double err = 0;
    return DoIntegral(binx1,binx2,0,-1,0,-1,err,option);
 }


  TGraph *g = new TGraph("eu_gu_1.dat");
  // g->Sort(); // just a precaution
  Int_t n = g->GetN();
  Double_t *x = new Double_t[(n + 1)];
  x[0] = (3.0 * g->GetX()[0] - g->GetX()[1]) / 2.0;
  x[n] = (3.0 * g->GetX()[(n - 1)] - g->GetX()[(n - 2)]) / 2.0;
  for (Int_t i = 1; i < n; i++)
    { x[i] = (g->GetX()[(i - 1)] + g->GetX()[i]) / 2.0; }
  TH1D *h = new TH1D("h", "spectrum;x-value;#counts", n, x);
  h->FillN(n, g->GetX(), g->GetY());
  h->Sumw2(kFALSE); // restore proper errors
  delete [] x; // no longer needed
  delete g; // no longer needed
  TF1 *f1 = new TF1("f1", "gausn(0) + pol1(3)", 400., 435.); // linear background
  f1->SetParameters(14000., 415., 2., 0., 0.);
  h->Fit("f1", "BR+"); // e.g. "BR+" or "LBR+"
  TF1 *f2 = new TF1("f2", "gausn(0) + pol1(3)", 1160., 1200.);
  f2->SetParameters(15000., 1178., 2., 0., 0.);
  h->Fit("f2", "BR+"); // e.g. "BR+" or "LBR+"
  TF1 *f3 = new TF1("f3", "gausn", 4795., 4835.); // clean peak, no background
  f3->SetParameters(4400., 4815., 3.);
  h->Fit("f3", "BR+"); // e.g. "BR+" or "LBR+"

 error: function definition is not allowed here
 {
 ^
<<< cling interactive line includer >>>:1:1: error: expected '}'

Wile_E_Coyote · March 8, 2020, 10:01am

std::cout << h->Integral(h->FindFixBin(400.), h->FindFixBin(435.)) << std::endl;

arorastacy · March 8, 2020, 10:07am

So if i have to find it out for a number of peaks then each time i have to change values of the range. Is it possible to get the area of all the peaks simultaneously.

Wile_E_Coyote · March 8, 2020, 10:14am

${ROOTSYS}/tutorials/spectrum

arorastacy · March 8, 2020, 10:39am

Let’s say i choose the code:

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

In line 2, I am unable to understand what does this signifies? par[0] = ? par[1] = ?
How are the values of p[0] & p[1] taken?

Wile_E_Coyote · March 8, 2020, 11:06am

Well, it’s a good question for your supervisor.

arorastacy · March 8, 2020, 11:10am

Do I need to change par[0] + par[1] x[0] into gaus function? Am i thinking in right direction?

arorastacy · March 8, 2020, 3:27pm

I calculated the area under all the peaks by both the methods

by sigma2.5H
2.integral fix bin method
The two values are not exactly similar and are given below:
------------------------ Integration ------sigma2.5H
1st peak (400-435) – 17412 ------- 14346
2nd peak (1160-1200) 17388 ------- 15724
3rd peak (4795-4835) - 4598 --------- 4423 (comparably close)

Wile_E_Coyote · March 8, 2020, 4:06pm

Try with “LBR+” fit options (you’ll get better “agreement” for the last peak).

Then show it to your supervisor, who should easily be able to explain the nonnegligible differences between the fitted peak area and the corresponding histogram integral.

arorastacy · March 8, 2020, 4:48pm

Thanks for helping my throughout my problem but please tell me about how should I combine my code with calculating the no. of peaks codes.
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;
}

vaidgourav2995 · March 8, 2020, 7:31pm

Sir, Please explain this to me also. I am also unable to solve this.

system · March 22, 2020, 7:41pm

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