Facing Problem while fitting

arorastacy · March 7, 2020, 6:25am

Hi,
I drew a histogram by converting an Ascii file into root format.
Now i have to fit the peaks in the histogram.
I try to fit only two peaks from the histogram and only one of them is fitted quite nicely.
Why the other peak is not fitted properly?

Have a nice day.fithisto.C (776 Bytes)

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

arorastacy · March 7, 2020, 6:26am

This is the peak that i have represented by g2 in the above program.

arorastacy · March 7, 2020, 6:42am

What root is showing is given below:
I don’t know how it is taking all the parameters?

root [0] 
Processing fithisto.C...
Info in <TCanvas::MakeDefCanvas>:  created default TCanvas with name c1
x=128.000000, y=1.000000x=140.000000, y=1.000000x=141.000000, y=3.000000x=142.000000, y=2.000000x=143.000000, y=1.000000 FCN=1793.1 FROM MIGRAD    STATUS=CONVERGED     335 CALLS         336 TOTAL
                     EDM=3.43385e-09    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  Constant     2.41186e+01   6.74203e-01   9.38469e-03   1.23643e-04
   2  Mean         1.70820e+03   2.88060e+01   3.57311e-01  -1.83124e-07
   3  Sigma        1.16629e+03   1.85217e+01   2.05567e-01   3.67220e-06
 FCN=4603.23 FROM MIGRAD    STATUS=CONVERGED     118 CALLS         119 TOTAL
                     EDM=5.89992e-07    STRATEGY= 1  ERROR MATRIX UNCERTAINTY   3.0 per cent
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  Constant     2.53437e+03   1.34782e+03  -2.55656e+01  -1.42555e-06
   2  Mean         1.17827e+03   3.93463e-01  -1.01192e-03   1.21152e-03
   3  Sigma        2.65935e+00   3.16466e-01   7.23349e-03  -1.69006e-03
found 5212points

Wile_E_Coyote · March 7, 2020, 8:32am

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

arorastacy · March 7, 2020, 12:39pm

when i draw histograms i have a freedom of using bins which i require for the later part of my problem.
therefore i am not using graph. Also using the same code i didn’t get the proper fit.

Wile_E_Coyote · March 7, 2020, 1:53pm

What do you get when you run my example macro without any changes?

arorastacy · March 7, 2020, 2:09pm

Now, please tell me how to zoom it so that i could check the fitting.

Wile_E_Coyote · March 7, 2020, 2:12pm

So, it looks perfectly fine.
You can zoom using a mouse, for example (sometimes you first need to select a wider zoom range and then again zoom the newly displayed picture to get a more detailed view).

arorastacy · March 7, 2020, 2:39pm

Thanks for the help.
I have shown this to my supervisor but he said to draw it only on the histogram instead of the graph.
That’s why i was asking the same.

Wile_E_Coyote · March 7, 2020, 2:45pm

Change your supervisor.

{
  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", "gaus(0) + pol1(3)", 400., 435.); // linear background
  f1->SetParameters(2800., 415., 2., 0., 0.);
  h->Fit("f1", "BR+"); // e.g. "BR+" or "LBR+"
  TF1 *f2 = new TF1("f2", "gaus(0) + pol1(3)", 1160., 1200.);
  f2->SetParameters(2700., 1178., 2., 0., 0.);
  h->Fit("f2", "BR+"); // e.g. "BR+" or "LBR+"
  TF1 *f3 = new TF1("f3", "gaus", 4795., 4835.); // clean peak, no background
  f3->SetParameters(500., 4815., 3.);
  h->Fit("f3", "BR+"); // e.g. "BR+" or "LBR+"
}

arorastacy · March 7, 2020, 4:38pm

Thanks so much dear.
Your syntax is running quite nicely. As i am a newbie can you explain the syntax please.

Wile_E_Coyote · March 7, 2020, 4:50pm

I’ll let your supervisor shine.

arorastacy · March 7, 2020, 5:24pm

Please help me in understanding the syntax.

Wile_E_Coyote · March 7, 2020, 5:33pm

TGraph
TH1
TF1

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