How to find the area error using peakfind.C

Roofit and RooStats
15

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