How to find the area error using peakfind.C

Dear all,

I am unable to understand how to find the error in area/integral using the code peakfind.C ([ROOT: tutorials/spectrum/peaks.C Source File]). Kindly suggest.

peaks.C (4.7 KB)

1 Like

Thanks a lot. I get area and the error associated. However, The area I got using this fitting is less than the value obtained from another software i.e. 1000 counts less. Could you suggest a reason for this. I have checked carefully the background subtraction for both look reasonable.

Try (see the TH1::Fit method description for available fit options):

   h2->Fit("fit", "L"); // "" or "L" or "WL" or "I" or "IL" or "IWL" or ...
1 Like

fit
The value didn’t improve after trying the fit options. I have attached the image showing the fitted peak with background line. The is for h2->Fit(“fit”) (seems the best one).

You could attach a root file with the “h_energy” histogram for inspection.

60Co_calibrated.root (21.9 KB)
60Co_uncalibrated.root (13.0 KB)
Area that I obtain after fitting in 60Co_calibrated.root is 8336 at 1330keV.

Using the "60Co_calibrated.root" file (attached in another post here), try this:

{
  const Bool_t fit_with_erfc_background = kTRUE; // kFALSE or kTRUE
  
  TFile *f = TFile::Open("60Co_calibrated.root");
  TCanvas *C2; f->GetObject("C2", C2);
  C2->Draw();
  
  TH1 *h = (TH1*)C2->FindObject("h_energy");
  h->SetDirectory(gROOT);
  h->GetListOfFunctions()->Delete(); // "remove" old stuff
  h->SetLineColor(kViolet);
  // h->Rebin(4); // try out the "average X bin width" != 1 case
  
  // Double_t x1 = 1125., x2 = 1215.; // the first peak
  // Double_t x1 = 1130., x2 = 1210.;
  Double_t x1 = 1285., x2 = 1375.; // the second peak
  // Double_t x1 = 1290., x2 = 1370.;
  
  Int_t i1 = h->FindFixBin(x1); i1 = TMath::Max(i1, 1);
  Int_t i2 = h->FindFixBin(x2); i2 = TMath::Min(i2, h->GetNbinsX());
  x1 = h->GetBinCenter(i1), x2 = h->GetBinCenter(i2);
  
  TF1 *ff; // the "fit function"
  if (fit_with_erfc_background) {
    ff = new TF1("ff",
                 "gausn(0) + ([3] + [4] * TMath::Erfc((x - [1]) / [2] / TMath::Sqrt2()))", x1, x2);
  } else {
    ff = new TF1("ff", "gausn(0) + pol1(3)", x1, x2);
  }
  
  // estimate "reasonable" initial values of parameters for the fits
  h->GetXaxis()->SetRangeUser(x1, x2);
  Double_t integral = h->Integral("width");
  Double_t mean = h->GetMean();
  Double_t sigma = h->GetStdDev();
  Double_t a = h->GetBinContent(i1), b = h->GetBinContent(i2);
  if (fit_with_erfc_background) {
    a = (a - b) / 2.; // y = a * erfc((x - mean) / sigma / sqrt(2)) + b
  } else {
    a = (b - a) / (x2 - x1); b = b - a * x2; // y = a * x + b
  }
  
  // note: the histogram is fine-grained so no need to try the "I" fit option
  
  if (gROOT->GetVersionInt() >= 60000) { // ROOT 6.00/00 or newer
    std::cout << std::endl << "... BLACK ... \"PR\" ..." << std::endl;
    ff->SetParameters(integral, mean, sigma, b, a);
    ff->SetLineColor(kBlack); ff->SetLineStyle(1);
    h->Fit(ff, "PR+"); // chi2 method (uses "expected errors" = Pearson chi2)
  }
  
  std::cout << std::endl << "... RED ... \"R\" ..." << std::endl;
  ff->SetParameters(integral, mean, sigma, b, a);
  ff->SetLineColor(kRed); ff->SetLineStyle(9);
  h->Fit(ff, "R+"); // chi2 method (uses "observed errors" = default chi2)
  
  std::cout << std::endl << "... GREEN ... \"WWR\" ..." << std::endl;
  ff->SetParameters(integral, mean, sigma, b, a);
  ff->SetLineColor(kGreen); ff->SetLineStyle(7);
  h->Fit(ff, "WWR+"); // chi2 method, all weights 1
  
  std::cout << std::endl << "... BLUE ... \"LR\" ..." << std::endl;
  ff->SetParameters(integral, mean, sigma, b, a);
  ff->SetLineColor(kBlue); ff->SetLineStyle(2);
  h->Fit(ff, "LR+"); // loglikelihood method
  
  // let's play with the last fit only ....
  
  // ... draw the background
  TF1 *fb; // the "background function"
  if (fit_with_erfc_background) {
    fb = new TF1("fb", "[2] + [3] * TMath::Erfc((x - [0]) / [1] / TMath::Sqrt2())", x1, x2);
    fb->SetParameters(ff->GetParameters() + 1);
  } else {
    fb = new TF1("fb", "pol1", x1, x2);
    fb->SetParameters(ff->GetParameters() + 3);
  }
  fb->SetLineColor(kYellow);
  fb->Draw("same");
  
  // ... if the "average X bin width" != 1, the "raw" and "adjusted" values differ
  std::cout << "gaus peak integral (\"raw\")  = "
            << ff->GetParameter(0) << " +- " << ff->GetParError(0)
            << std::endl;
  Double_t dx = (x2 - x1) / (i2 - i1); // "average X bin width"
  // std::cout << "average X bin width = " << dx << std::endl;
  std::cout << "gaus peak integral adjusted = "
            << ff->GetParameter(0) / dx << " +- " << ff->GetParError(0) / dx
            << std::endl;
}

It seems to me that all cases return reasonable “area” values (e.g., for the second peak, between 8200 and 8330, but note that their errors are between 100 and 140).

While thinking about a possible source of differences in values returned by “another software”, it came to my mind that it could be the “average X bin width”, if it is not 1 (see the end of the macro above how to “adjust” the integral).

BTW. You may also be interested in this old post (it also shows how to “adjust” the integral when using a graph instead of a histogram).

1 Like

The binwidth that I have chosen is 1. But I guess the other software has may be 0.3. But in my code, when I play with the bin width for calibrated histogram:
TH1* h_energy = new TH1F(“h_energy”, “Calibrated Energy Spectra”,4096,0,4096);
the calibration always shifts. For reference, I am attaching the code that I wrote with the help of Peakfind.C and I am including the input file also.
root.zip (6.7 KB)

I’m afraid the way you create histograms in the “call_spectra_1.C” routine may create problems with getting proper “area” values later.

Using the 05cm_60Co.Spe.txt file (extracted and renamed from the "root.zip" file, attached in another post here), try to play with the macro given below:

Spe_60Co.cxx

1 Like

Thanks a lot for this optimised and beautiful code, This certainly works. I have few doubts:
Do we need to change dE value always according to the peak?
The result that I obtained from “i1 = h->FindFixBin(x1);” and “i1 = TMath::Max(i1, 1);” is same, are there some exceptions to this when the value may change?
Another issue is that I wanted to have generalized peak fitting because I have many source measurements and some spectra have two peaks (not completely resolved), those peaks can not be fitted with this method, right as it fits one peak at a time? Example file is attached.
5cm_57Co.zip (5.1 KB)

I now added some small comments in the “Spe_60Co.cxx” macro regarding the “i1” and “i2” (just download it again).

The “dEn” value is your choice (it’s up to you how “wide” / “narrow” you want to fit the peak, just try different values). Note that different lines may have different peak widths, so it makes sense to have “dEn” as an additional function parameter.

If you have two overlapping peaks, then you need a function that has two Gaussians, e.g., "gausn(0) + gausn(3) + pol1(6)" (note that you need to set “reasonable” initial values for all parameters, before you try to fit it). That’s actually what the “peaksC” tutorial does.

BTW. If you want to have a “generalized peak fitting”, you will need to make sure that the spectrum’s “background” is subtracted first. There exist some “Spectrum tutorials” and some examples on the ROOT Forum which deal with it.

1 Like

Thank you @Wile_E_Coyote . I have written a small macro (by taking some online help) to fit double gaussian peak using parameters obtained from single fit data.
double_peak_fit_57Co.root (41.0 KB)
Here are few output statements:
gaus peak integral1 = 590800 ± 774.567
gaus peak integral2 = 63938.9 ± 250.478
ChiSquare= 570.468 NDf= 39 ChiSquare/NDf=14.6274

Does the values look reasonable for chisquare/NDF and errors?

Your results seem reasonable (the drawn “yellow” curve seems to be the “background” as it was before some “background subtraction procedure” was applied).

The “background subtraction procedure” created some bins with content < 0. Unfortunately, ROOT has problems with such cases so, when fitting, you should set the “xmax” so that the bin content is still >= 0 (e.g., use “x2 = 146.” in the macro below as, if you try “x2 >= 148.” and “restore_background = kFALSE”, the chi2 will immediately “explode”). Maybe @moneta can comment on this.

Another problem is that the “background subtraction procedure” did not handle bin errors, and then the “chi2 / ndf” is not really very meaningful (well, one can also question the errors of integrals, of course).

Using the "double_peak_fit_57Co.root" file (attached in another post here), try this:

{
  const Bool_t restore_background = kTRUE; // kFALSE or kTRUE
  
  gStyle->SetOptFit();
  TFile *f = TFile::Open("double_peak_fit_57Co.root");
  TCanvas *c; f->GetObject("c1", c);
  c->Draw();
  gROOT->cd(); // newly created histograms should go here
  // TH1F *h = (TH1F*)c->FindObject("h_energy")->Clone();
  TH1F *h = (TH1F*)c->GetPrimitive("h_energy")->Clone();
  if (restore_background) {
    h->Add((TH1F*)c->GetPrimitive("h_energy_background"));
    h->Sumw2(kFALSE);
  }
  h->GetListOfFunctions()->Delete(); // "remove" old stuff
  h->SetLineColor(kViolet);
  h->Draw();
  Double_t x1 = 109., x2 = 146.;
  Int_t i1 = h->FindFixBin(x1); i1 = TMath::Max(i1, 1);
  Int_t i2 = h->FindFixBin(x2); i2 = TMath::Min(i2, h->GetNbinsX());
  x1 = h->GetBinCenter(i1), x2 = h->GetBinCenter(i2);
  Double_t dx = (x2 - x1) / (i2 - i1); // "average X bin width"
  TF1 *ff = new TF1("ff", "gausn(0) + gausn(3) + pol1(6)", x1, x2);
  ff->SetParameters(450000., 122.1, 3.5, 45000., 136.5, 3.5, 500., 0.);
  h->GetXaxis()->SetRangeUser(x1, x2);
  h->Fit(ff, "L", "", x1, x2); // e.g. "L" or "P" or ""
  std::cout << "gaus peak 1 integral (\"raw\")  = "
            << ff->GetParameter(0) << " +- " << ff->GetParError(0)
            << std::endl;
  std::cout << "gaus peak 1 integral adjusted = "
            << ff->GetParameter(0) / dx << " +- " << ff->GetParError(0) / dx
            << std::endl;
  std::cout << "gaus peak 2 integral (\"raw\")  = "
            << ff->GetParameter(3) << " +- " << ff->GetParError(3)
            << std::endl;
  std::cout << "gaus peak 2 integral adjusted = "
            << ff->GetParameter(3) / dx << " +- " << ff->GetParError(3) / dx
            << std::endl;
}

BTW. When using the “Maximum Likelihood Fit”, the displayed “chi2” should be equal to “2 * FCN” (if it is not, you are facing a known problem).

1 Like

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