Integral Error in Signal

[lamya86]

Hello,

I am trying to get the number of signal events with and without background events in 3 sigma range. I have 2D histogram for missing momentum vs missing mass which called (NDE_mmVSmp_exp) in my root histogram file. I plot missing mass and fit it with 2 functions (gaus+pol3). Then I got the number of tatal events within 3sigma with integral error(which contain signal +background) and also I get the number of events only in the signal. However, when I try to get the integral error in the signal, it is gaven me zero and it shown the following massage “Error in TF1Helper::IntegralError: Last used fitter is not compatible with the current TF1”. I know that I need to define conv matrix but I don’t know how I can add it to my code. Could please help me.

First off here is my code:
epip.py (1009 Bytes)

and here is my root histogram file:
histneutronefficincy.hipo.root (160.3 KB)

Here my plot with fitting

Thanks in advance for your help.

---

Please read tips for efficient and successful posting and posting code

_ROOT Version:v6-16-00
_Platform:pyroot
Compiler: Not Provided

---

[lamya86]

Here my code

#!/usr/bin/python

import sys
import ROOT

ROOT.gStyle.SetOptFit(11111)

ff = ROOT.TFile(sys.argv[1])
h2 = ff.Get("NDE_mmVSmp_exp")

c1 = ROOT.TCanvas('c1','c1',1100,800)
c1.Print('test.pdf[')
h2.Draw('colz')
c1.Print('test.pdf')
f1=ROOT.TF1("f1","gaus(0)+pol2(3)",0.81,1.02)
f1.SetParameters(1,1,0.2,1,1,1)
hy1 = h2.ProjectionY()
hy1.Draw()
hy1.Fit(f1,"R")
mu1,sig1 = f1.GetParameter(1), f1.GetParameter(2)
r11,r21 = mu1-3*sig1,mu1+3*sig1

fsig1 = ROOT.TF1("fsig","gaus",r11,r21)
fsig1.SetParameters(*[f1.GetParameter(i) for i in range(3)])
fsig1.SetLineColor(4)
fsig1.Draw('same')
ntot1 = abs(f1.Integral(r11,r21))/hy1.GetBinWidth(1)
nerrtot1 =abs(f1.IntegralError(r11,r21))/hy1.GetBinWidth(1)
nsig1 = abs(fsig1.Integral(r11,r21))/hy1.GetBinWidth(1)
nsigerr1 = abs(fsig1.IntegralError(r11,r21))/hy1.GetBinWidth(1)
hy1.SetTitle('ntot={:.2f}#pm {:.2f},nsig={:.2f},nsigerr={:.2f},ratio={:.2f}%;missing mass'.format(ntot1,nerrtot1,nsig1,nsigerr1,nsig1/ntot1*100))

c1.Print('test.pdf')

c1.Print('test.pdf]')

[Wile_E_Coyote]

A rough approach: epip.py (1.2 KB)
An exact approach (the second macro therein):

TH1F variable bin width

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");…

[lamya86]

Thank you for your reply. I have tried the way that you mentioned but I get confused to get the number of events in the background with error, the number of events in the signal with error and the total number of events in both signal and back ground. Could you please help me to get them correctly.

Here is my code using exact approach

 void electronPionID(){
    gStyle->SetOptFit(11111);
    const Bool_t use_partial_covariances = kFALSE; // kFALSE or kTRUE
    TFile *f= TFile::Open("histneutronefficincy.hipo.root");
    TH2D *h1 = new TH2D("h1","",150 , 0, 10,150,0,2);
    h1= (TH2D*)f->Get("NDE_mmVSmp_exp");

    TH1D *g1 = h1->ProjectionY();
    TF1 *fg = new TF1("fg", "gaus(0)");
    //Int_t ch_min = 0.8, ch_max = 1.1; // "channel numbers"
    TF1 *f1 = new TF1("f1", "gaus(0) + pol2(3)", 0.83,1.13);
    
    f1->SetParameters(1, 1, 0.2, 1, 1,1);
    TFitResultPtr r1 = g1->Fit(f1, "WEMRS");
    std::cout << "r1) = " << r1  << std::endl;
    TMatrixD c1 = r1->GetCovarianceMatrix();
    c1.Print();
    TMatrixD c1g = c1.GetSub(0, 2, 0, 2);
    c1g.Print();
    if (use_partial_covariances) c1g -=
                                 c1.GetSub(0, 2, 3, 4) *
                                 c1.GetSub(3, 4, 3, 4).InvertFast() *
                                 c1.GetSub(3, 4, 0, 2);
    c1g.Print();
    Double_t s1, s1e, mu,sigma, N_s, N_s_Err;
    fg->SetParameters(f1->GetParameters());
    // fg->SetParErrors(f1->GetParErrors());
    mu = fg->GetParameter(1);
    sigma = fg->GetParameter(2);
    std::cout << "mu = " << mu << std::endl;
    std::cout << "sigma = " << sigma << std::endl;
    Double_t dx; // "average X bin width"
    dx = g1->GetBinWidth(1);
    s1 = fg->Integral(mu -(3.0*sigma), mu+(3.0 * sigma));
    N_s=s1/dx;
    s1e = fg->IntegralError(mu -(3.0*sigma), mu+(3.0 * sigma),fg->GetParameters(),c1g.GetMatrixArray());
    N_s_Err = s1e/dx;
    std::cout << "gaus peak integral (mean +/- 3 sigma) = " << s1
              << " +/- " << s1e
              << std::endl;
    
    Double_t N_sb, N_sb_Err, mu_total,sigma_total;
    mu_total = f1->GetParameter(1);
    sigma_total = f1->GetParameter(2);
    N_sb = f1->Integral(mu_total-(3*sigma_total), mu_total+(3*sigma_total))/dx;
    N_sb_Err = f1-> IntegralError(mu_total-(3*sigma_total), mu_total+(3*sigma_total))/dx;
    std::cout <<  "   N_sb = " << N_sb
              << " +/- " << N_sb_Err
              << std::endl;
    g1->SetTitle(Form("N_sb = %4.3f, N_sb_Err = %4.3f , N_s = %4.3f , N_s_Err = %4.3f ", N_sb,N_sb_Err,N_s,N_s_Err)); 
    
    
    

    
    TCanvas *c = new TCanvas("c", "c",1100,800);
    c->Divide(1, 2);
    c->Print("test1.pdf[");   
    c->cd(1);
    g1->Draw();
    fsig->Draw("same");
    c->Print("test1.pdf");
    c->cd(2);
    h1->Draw("colz");
    //c->Print("test1.pdf");
    c->cd(0);
    c->Print("test1.pdf]");
}

[Wile_E_Coyote]

The original example uses “pol1(3)” so there are 2 parameters, but you use “pol2(3)” so you have 3 parameters. That means that you need to change “3, 4” into “3, 5” in all relevant places.
In the first macro therein, you also need to add two lines (needed when “fit_with_fixed_background” = “kTRUE”):

    // ...
    f1->FixParameter(4, f1->GetParameter(4)); // this line already exists
    f1->FixParameter(5, f1->GetParameter(5));
    // ...
    f2->FixParameter(4, f2->GetParameter(4)); // this line already exists
    f2->FixParameter(5, f2->GetParameter(5));
    // ...

You can get the integral and error of your “signal + background” directly from your fitted “f1”. For the “signal” alone, you then have “fg”. If you want the “background” alone, then you will need to create another function (using “pol2”) and play the same game as for your “fg”.

[lamya86]

Thank you Wile_E_Coyote for helping me. I still have issue with background integral error which is given me zero and the number of the events in the background is larger than the total number. Also I notice that chi2/ndf is too large. That mean there is some thing wrong I made it which I don’t know what exactly is. Could you please help me

void electronPionID(){
    gStyle->SetOptFit(11111);
    const Bool_t use_partial_covariances = kFALSE; // kFALSE or kTRUE
    TFile *f= TFile::Open("histneutronefficincy.hipo.root");
    TH2D *h1 = new TH2D("h1","",150 , 0, 10,150,0,2);
    h1= (TH2D*)f->Get("NDE_mmVSmp_exp");

    TH1D *g1 = h1->ProjectionY();
    TH1D *g2 = h1->ProjectionX();
    TF1 *fg = new TF1("fg", "gaus(0)");
    TF1 *fbc = new TF1("fbc", "pol2(3)");
    //Int_t ch_min = 0.8, ch_max = 1.1; // "channel numbers"
    TF1 *f1 = new TF1("f1", "gaus(0) + pol2(3)", 0.70,1.13);
    
    f1->SetParameters(1, 1, 0.2, 1, 1,1);
    TFitResultPtr r1 = g1->Fit(f1, "WEMRS");
    std::cout << "r1) = " << r1  << std::endl;
    TMatrixD c1 = r1->GetCovarianceMatrix();
    c1.Print();
 
    TMatrixD c1g = c1.GetSub(0, 2, 0, 2);
    c1g.Print();

    if (use_partial_covariances) c1g -=
                                 c1.GetSub(0, 2, 3, 5) *
                                 c1.GetSub(3, 5, 3, 5).InvertFast() *
                                 c1.GetSub(3, 5, 0, 2);

    Double_t N_sig, N_sigErr, N_bg, N_bgErr, N_tot, N_totErr;
    // Getting Number of events in the signal peak
    fg->SetParameters(f1->GetParameters());
    N_sig = fg->Integral(fg->GetParameter(1) -(3.0*fg->GetParameter(2)), fg->GetParameter(1)+(3.0 * fg->GetParameter(2)))/g1->GetBinWidth(1);
    N_sigErr= fg->IntegralError(fg->GetParameter(1) -(3.0*fg->GetParameter(2)), fg->GetParameter(1)+(3.0 * fg->GetParameter(2)),fg->GetParameters(),c1g.GetMatrixArray())/g1->GetBinWidth(1);
    std::cout << "Number of events in signal peak (mean +/- 3 sigma) = " << N_sig
              << " +/- " << N_sigErr
              << std::endl;
    
    // Getting Number of events in the back ground 
    fbc->SetParameters(f1->GetParameters());
    N_bg = fbc->Integral(fg->GetParameter(1)-(3.0*fg->GetParameter(2)), fg->GetParameter(1)+(3.0 * fg->GetParameter(2)))/g1->GetBinWidth(1);
    N_bgErr = fbc->IntegralError(fg->GetParameter(1) -(3.0*fg->GetParameter(2)), fg->GetParameter(1)+(3.0 * fg->GetParameter(2)),fbc->GetParameters(),c1g.GetMatrixArray())/g1->GetBinWidth(1);
    std::cout << "Number of events in background (mean +/- 3 sigma) = " << N_bg
              << " +/- " << N_bgErr
              << std::endl;

    // Getting Number of events for Both 
    N_tot= f1->Integral(fg->GetParameter(1)-(3.0*fg->GetParameter(2)), fg->GetParameter(1)+(3.0 * fg->GetParameter(2)))/g1->GetBinWidth(1);
    N_totErr = f1->IntegralError(fg->GetParameter(1)-(3.0*fg->GetParameter(2)), fg->GetParameter(1)+(3.0 * fg->GetParameter(2)),f1->GetParameters(),c1g.GetMatrixArray())/g1->GetBinWidth(1);
    std::cout <<  "  Number of the total events = " << N_tot
              << " +/- " << N_totErr
              << std::endl;
   
    double integral = f1->GetParameter(0) * f1->GetParameter(2) * TMath::Sqrt(TMath::TwoPi()) / g1->GetBinWidth(0);
    double integralErr = integral * sqrt(pow(f1->GetParError(0) / f1->GetParameter(0), 2) + pow(f1->GetParError(2) / f1->GetParameter(2), 2));
    std::cout <<  "  Number of sig = " << integral
              << " +/- " << integralErr
              << std::endl;
    
    double x1,x2;
    x1= f1->GetParameter(1)-(3.0*f1->GetParameter(2));
    x2= f1->GetParameter(1)+(3.0 * f1->GetParameter(2));
    TF1 * fSig = new TF1("fSig","gaus",x1, x2);
     // Set the values of the gaussian and parabola
    for (int ipar=0;ipar<3;ipar++){
        fSig->SetParameter(ipar,f1->GetParameter(ipar));
        //fbc->SetParameter(ipar,f1->GetParameter(ipar+3));
    }
   
  
    TCanvas *c = new TCanvas("c", "c",1100,800);
    c->Divide(1, 2);
    TLatex *latex = new TLatex();
    TLatex *latex1 = new TLatex();
    TLatex *latex2 = new TLatex();
    c->Print("test1.pdf[");   
    c->cd(1);
    g1->Draw();

    latex->DrawLatexNDC(0.13,0.80,Form("total: %.2f +/- %.2f", N_tot, N_totErr));
    latex1->DrawLatexNDC(0.13,0.74,Form(" sig: %.2f +/- %.2f",N_sig, N_sigErr));
    latex2->DrawLatexNDC(0.13,0.69,Form("bacground : %.2f +/- %.2f",N_bg, N_bgErr));
    fSig->SetLineColor(kBlue);
    fSig->DrawClone("Same"); 
    //fbc->DrawClone("Same");
    
    c->Print("test1.pdf");
    c->cd(2);
    h1->Draw("colz");
    //c->Print("test1.pdf");
    c->cd(0);
    c->Print("test1.pdf]");
  
  
}

[Wile_E_Coyote]

At least:

// ...
TF1 *fbc = new TF1("fbc", "pol2");
// ...
TMatrixD c1bc = c1.GetSub(3, 5, 3, 5);
if (use_partial_covariances) c1bc -=
                               c1.GetSub(3, 5, 0, 2) *
                               c1.GetSub(0, 2, 0, 2).InvertFast() *
                               c1.GetSub(0, 2, 3, 5);
// ...
fbc->SetParameters(f1->GetParameters() + 3);
// ...