Gauss and Pol1 Fitting Results

[spyhunter01]

@TH1F variable bin width ROOT

{
  const Bool_t use_partial_covariance = kFALSE; // kFALSE or kTRUE
  
  const char *filename = "60Codata.txt";
  TGraph *g0 = new TGraph(filename, "%lg %lg %*lg"); // X axis = channels
  g0->SetTitle("g0;Channel number;Energy"); // energy calibration
  TGraph *g1 = new TGraph(filename, "%lg %*lg %lg"); // X axis = channels
  g1->SetTitle("g1;Channel number;Counts");
  TGraph *g2 = new TGraph(filename, "%*lg %lg %lg"); // X axis = energies
  g2->SetTitle("g2;Energy;Counts");
  if ((g1->GetN() < 2) || (g2->GetN() < 2)) return; // just a precaution
  
  TF1 *fg = new TF1("fg", "gaus");
  Double_t sb = 5.0; // sigma band for the fg integral (e.g. 3 or 5 or 7)
  
  Int_t ch_min = 360, ch_max = 460; // "channel numbers"
  TF1 *f1 = new TF1("f1", "gaus(0) + pol1(3)", ch_min, ch_max);
  f1->SetParameters(30000., 430., 10., 5000., 0.);
  TFitResultPtr r1 = g1->Fit(f1, "WEMRS");
  TMatrixD c1 = r1->GetCovarianceMatrix();
  TMatrixD c1g = c1.GetSub(0, 2, 0, 2);
  Double_t s1, s1e;
  fg->SetParameters(f1->GetParameters());
  // fg->SetParErrors(f1->GetParErrors());
  s1 = fg->Integral(fg->GetParameter(1) - sb * fg->GetParameter(2),
                    fg->GetParameter(1) + sb * fg->GetParameter(2));
  s1e = fg->IntegralError(fg->GetParameter(1) - sb * fg->GetParameter(2),
                          fg->GetParameter(1) + sb * fg->GetParameter(2),
                          fg->GetParameters(),
                          c1g.GetMatrixArray());
  std::cout << "gaus peak integral (mean +/- " << sb << " sigma) = "
            << s1 << " +/- " << s1e << std::endl;
  if (use_partial_covariance) {
    c1g -=
      c1.GetSub(0, 2, 3, 4) *
      c1.GetSub(3, 4, 3, 4).InvertFast() *
      c1.GetSub(3, 4, 0, 2);
    s1e = fg->IntegralError(fg->GetParameter(1) - sb * fg->GetParameter(2),
                            fg->GetParameter(1) + sb * fg->GetParameter(2),
                            fg->GetParameters(),
                            c1g.GetMatrixArray());
    std::cout << "note: using partial covariance" << std::endl;
    std::cout << "gaus peak integral (mean +/- " << sb << " sigma) = "
              << s1 << " +/- " << s1e << std::endl;
  }
  
  TF1 *f2 = new TF1("f2", "gaus(0) + pol1(3)",
                    g2->GetX()[ch_min], g2->GetX()[ch_max]);
  f2->SetParameters(30000., 1180., 28., 5000., 0.);
  TFitResultPtr r2 = g2->Fit(f2, "WEMRS");
  TMatrixD c2 = r2->GetCovarianceMatrix();
  TMatrixD c2g = c2.GetSub(0, 2, 0, 2);
  Double_t s2, s2e;
  fg->SetParameters(f2->GetParameters());
  // fg->SetParErrors(f2->GetParErrors());
  s2 = fg->Integral(fg->GetParameter(1) - sb * fg->GetParameter(2),
                    fg->GetParameter(1) + sb * fg->GetParameter(2));
  s2e = fg->IntegralError(fg->GetParameter(1) - sb * fg->GetParameter(2),
                          fg->GetParameter(1) + sb * fg->GetParameter(2),
                          fg->GetParameters(),
                          c2g.GetMatrixArray());
  std::cout << "gaus peak integral (mean +/- " << sb << " sigma) = "
            << s2 << " +/- " << s2e << std::endl;
  Double_t dx; // "average X bin width"
  dx = (g2->GetX()[ch_max] - g2->GetX()[ch_min]) / (ch_max - ch_min);
  std::cout << "gaus peak integral adjusted (mean +/- " << sb << " sigma) = "
            << s2 / dx << " +/- " << s2e / dx << std::endl;
  if (use_partial_covariance) {
    c2g -=
      c2.GetSub(0, 2, 3, 4) *
      c2.GetSub(3, 4, 3, 4).InvertFast() *
      c2.GetSub(3, 4, 0, 2);
    s2e = fg->IntegralError(fg->GetParameter(1) - sb * fg->GetParameter(2),
                            fg->GetParameter(1) + sb * fg->GetParameter(2),
                            fg->GetParameters(),
                            c2g.GetMatrixArray());
    std::cout << "note: using partial covariance" << std::endl;
    std::cout << "gaus peak integral (mean +/- " << sb << " sigma) = "
              << s2 << " +/- " << s2e << std::endl;
    std::cout << "gaus peak integral adjusted (mean +/- " << sb << " sigma) = "
              << s2 / dx << " +/- " << s2e / dx << std::endl;
  }
  
  TCanvas *c = new TCanvas("c", "c");
  c->Divide(1, 2);
  c->cd(1);
  g1->Draw("AP");
  c->cd(2);
  g2->Draw("AP");
  c->cd(0);
}

In your 2nd code, what exactly this part :" dx = (g2->GetX()[ch_max] - g2->GetX()[ch_min]) / (ch_max - ch_min);" adjusting? It looks like the slope of energy vs channel plot. First integral you find is divided by the dx, then the numbers are too different. You go 5sigma from the mean on both side to get the integral, do you take the portion in to account to retrieve the whole count. Let’s say if you used 1sigma, it would have covered 68.3% of the actual data.

In the first code, integral is taken from gauss+background, but in the second one it is taken from gauss only. My question is still the same: why we dont take the integral from gauss-background?

Fits look perfect, by the way.