How to find the area error using peakfind.C

Roofit and RooStats
8

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