Using the "60Codata.txt" file (attached in another post here), try this:
{
const Bool_t fit_with_fixed_background = 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
Int_t ch_min = 360, ch_max = 460; // "channel numbers"
TF1 *f1 = new TF1("f1", "gausn(0) + pol1(3)", ch_min, ch_max);
f1->SetParameters((30000. * 10. * 2.5), 430., 10., 5000., 0.);
g1->Fit(f1, "WEMR");
std::cout << "gaus peak integral = " << f1->GetParameter(0)
<< " +/- " << f1->GetParError(0)
<< std::endl;
if (fit_with_fixed_background) {
f1->FixParameter(3, f1->GetParameter(3));
f1->FixParameter(4, f1->GetParameter(4));
g1->Fit(f1, "WEMR");
std::cout << "note: background parameters are fixed" << std::endl;
std::cout << "gaus peak integral = " << f1->GetParameter(0)
<< " +/- " << f1->GetParError(0)
<< std::endl;
}
TF1 *f2 = new TF1("f2", "gausn(0) + pol1(3)",
g2->GetX()[ch_min], g2->GetX()[ch_max]);
f2->SetParameters((30000. * 28. * 2.5), 1180., 28., 5000., 0.);
g2->Fit(f2, "WEMR");
std::cout << "gaus peak integral = " << f2->GetParameter(0)
<< " +/- " << f2->GetParError(0)
<< 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 = " << f2->GetParameter(0) / dx
<< " +/- " << f2->GetParError(0) / dx
<< std::endl;
if (fit_with_fixed_background) {
f2->FixParameter(3, f2->GetParameter(3));
f2->FixParameter(4, f2->GetParameter(4));
g2->Fit(f2, "WEMR");
std::cout << "note: background parameters are fixed" << std::endl;
std::cout << "gaus peak integral = " << f2->GetParameter(0)
<< " +/- " << f2->GetParError(0)
<< std::endl;
std::cout << "gaus peak integral adjusted = " << f2->GetParameter(0) / dx
<< " +/- " << f2->GetParError(0) / 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);
}
and / or that:
{
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);
}
See also: Wikipedia → Covariance matrix → Partial covariance matrix