Integral for a fit function

ROOT
1

Hi,

  I have a fucntion which i used to fit  my signal hisogram (The fitted histgram is attached).

Now I want to calculate the number of events using this fit funciton by its integration. The relevant part of code which does this is:

=====================
…
…
…
arglist[0]=0.;
gMinuit->mnexcm(“SIMPLEX”, arglist ,0,ierflg);
gMinuit->mnexcm(“IMPROVE”, arglist ,0,ierflg);

//Get the parameters
Double_t p[5],errp[5];

gMinuit->GetParameter(0, p[0], errp[0]);
gMinuit->GetParameter(1, p[1], errp[1]);
gMinuit->GetParameter(2, p[2], errp[2]);
gMinuit->GetParameter(3, p[3], errp[3]);
gMinuit->GetParameter(4, p[4], errp[4]);

TF1 *f1 = new TF1(“plotfunc”,plotfunc,xmin,xmax,5);

f1->SetParameters§;
f1->Print();

double Tot_SB=0.0;
Tot_SB= f1->Integral(2000.0,3000.0); //Get the interal in the range 2000-3000 GeV
cout<<Tot_SB<<endl;

The “Tot_SB” comes out to be 1817.54 which does not seems to be correct if one see the fitted histogram.

Could someone help me how to get correct integral value from the fitted function???

Thanks in advance.
with best,
sushil

1 Like
2

sorry forgot to attached the plot. Please find it here.

3

sorry not able to attach the plot… but one can see it here:

schauhan.web.cern.ch/schauhan/B … Region.eps

with best,
sushil

4

Could you post the shortest possible RUNNING script reproducing your problem, together with a file containing your histogram?
Why do you use Minuit directly instead of TH1::Fit?

Rene

5

Hello Rene,

I have put the root file which has the histogram and the code (SigFit.C) at the following link

schauhan.web.cern.ch/schauhan/B … /Integral/

Thanks and with best,
sushil

6

I assume that you want to divide the result of TF1::Integral by the x range (ie 1000).

You can also simplify your code as shown below

Rene

Double_t Sigfit(Double_t *x,Double_t *p) {
            Double_t sig1  = exp(p[0]+(p[1]*x[0]));
            Double_t arg1 = (x[0]-p[3])/p[4];
            Double_t sig2= p[2]*exp(-0.5*pow(arg1,2.0));          
            return (sig1+sig2);
}
void sigfit() { 
  TFile *f=new TFile("Data_1000.root","r");
  TH1F *histo1=(TH1F*)f->Get("hsb_data_0");
   Double_t xmin=1000.,xmax=6000.;  //range to be plotted 
   TF1 *f1 = new TF1("Sigfit",Sigfit,xmin,xmax,5);
   f1->SetParameters(0.001,0.0001,0.001,0.1,0.01);
   histo1->Fit(f1,"L","",1000,3000);
   printf("Integral(2000->3000)=%g\n",f1->Integral(2000,3000)/1000.);
}
7

Hello Rene,

Thanks for the code.  

From my version of code I was getting 1818 (  I have not divided by 1000 ) and from  your example I get 1.9 (after dividing by 1000).

   I am confused why we need to divide it by the integration range here. I thnk I am missing something very basic here???

Thanks for your prompt help.
with best,
sushil