I’ll let your supervisor shine.
Please help me in understanding the syntax.
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.
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?
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+"
}
- 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?
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
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
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
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 '}'
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.
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?
Well, it’s a good question for your supervisor.
Do I need to change par[0] + par[1] x[0] into gaus function? Am i thinking in right direction?
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)
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.
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;
}