TH1::IntegralError question

Dear root experts,

I have a gauss + landau TF1 fit to a TH1F histogram. I would like to extract the landau part of the function, draw it on the canvas over the histogram and perform a TF1::Integral over a given range. I have done this without a problem. A snippet of the code is shown below:

[code] TF1* combined = new TF1(“combined”,breitgauspluslandau,0.,1000.,8 );
combined->SetParameters(20,90,60,1,5,60,60);
h1->Fit(“combined”);

TF1* landau = new TF1(“landau”,"[0]TMath::Landau(x,[1],[2],[3])",0.,1000.);
Double_t
params = combined->GetParameters();
landau->SetParameters(params[4],params[5],params[6],params[7]);
Double_t errors[4] = {combined->GetParErrors()[4],combined->GetParErrors()[5],combined->GetParErrors()[6],combined->GetParErrors()[7]};
landau->SetParErrors(errors);[/code]

When i try to access the TF1::IntegralError of the landau, however, i get the following:

Can anyone describe what i need to do to access this information? Im in version 5.19.02.

Many thanks, Matt

Hi,

the function IntegralError works only for the fitted function, i.e in your case for your combined function. It needs to know the covariance matrix of all the parameters and this is not stored internally TF1 but in the TVirtualFitter.
You can alculate yourself as in the example

tutorials/fit/ErrorIntegral.C

but you would need to take the reduced covariance matrix for only the 4:7 (Landau) parameters.

Regards

Lorenzo

Thank you. It seems to work!

I have the same problem with a different sum of functions, in ROOT version 5.26. I was curious to know if it is now possible to set the covariance matrix from user input and an example on how to proceed could be made available. This would be so much simpler. I am now trying to look through the tutorial you suggested and find it quite cumbersome… It is a common issue in physics where one fits a sum of signal + background function, but only needs to access the integral and error on integral on the signal. It would therefore be useful to a lot of people to have this possibility.

Thanks in advance,

Andrée

I have some issue applying the ErrorIntegral.C tutorial to my specific situation (where I have 5 parameters). Basically, I don’t understand these lines:

   // calculate now the error (needs the derivatives of the function 
   // w..r.t the parameters)
   TF1 * deriv_par0 = new TF1("dfdp0",df_dPar,0,1,1);
   deriv_par0->SetParameter(0,0);

   TF1 * deriv_par1 = new TF1("dfdp1",df_dPar,0,1,1);
   deriv_par1->SetParameter(0,1.);

Your original function was a sin with 2 parameters. So I guess that you declare a TF1 for each derivative with n-1 params over the range of the integral you are interested in. But then I don’t get why you set the 0th parameter (always the same) to different values after declaring each derivative.

Maybe the other user who did set it up for a Landau function could post an example to this post. Thanks in advance!

Andrée

Hi,

in this particular example the parameter of the derivative function represents the index for the parameters of the original function. So The first partial derivative will have par0=0, the second one par0=1, etc…

I remind you that now (since ~ 2 years) we have implemented in TF1 the function TF1::IntegralError and the example has been updated and this code is not anymore used.

Lorenzo

Oh I see! That makes more sense, thanks.

If you look at my first post, you will see this was my original question. I am trying to use TF1::IntegralError on a signal part only, but I am fitting the signal and the background (exactly like the case of the original post). According to your answer, as of ROOT version 5.19, this would not work and one would have to adapt the tutorial, as you suggested yourself. Please let me know if it is now possible to use TF1::IntegralError on a subset of a fit, and if so, how to do it. I am using ROOT version 5.26.

Thanks,

Andrée

Hi,

I used a fit function on a histo histoSi :

histoSi->Fit(“gaus”);
TF1* fitSi = histoSi->GetFunction(“gaus”);

Then I want to do that :
fitSi->IntegralError(rangemin,rangemax)

I read this could be a covariant matrix probleme, that I have to re-adapt the code inside the fucntion IntegralError() …

Is there a simple way to do this, without dealing with covariant matrix myself (It is urgent, and while this is very interesting I don’t want to focus on this right now)

Thank you by advance

Hi,

could anyone answer please? I don’t understand since my function is fitted on my histogram, such that the covariant matrix should be assigned with some values …

Thanks

Hi,

I have similar problem, I am confused about what changed and how to use IntegralError now.

so I have a combined fit function:

TF1 *total = new TF1(“total”,“expo(0)+expo(2)+expo(4)”,0,5);
total->SetParameters(par);
hpT->Fit(total,“R”);

and it works great, the I want to get integral of the first exponent on some range, and its ERROR,
so I define new function expoGG:

TF1 *expoGG = new TF1(“expoGG”,“expo”,0.0,5);

expoGG->SetParameter(0,total->GetParameter(0));
expoGG->SetParameter(1,total->GetParameter(1));
expoGG->SetParError(0,total->GetParError(0));
expoGG->SetParError(1,total->GetParError(1));

and I get the integral (works great):

float Ngg=expoGG->Integral(0.1,0.7)/((x2-x1)/nbins);

and the error:

float errNgg=expoGG->IntegralError(0.1,0.7);

and I get “Error in TF1Helper::IntegralError: Last used fitter is not compatible with the current TF1”

and it gives errNgg=0.

why is that?

Thanks for help!
Marek

Hi,
You need to pass the sub-correlation matrix for the function component you want to compute the errors.
The covarince matrix which is used by default, is the global one, and for this reason is not compatible with the function you want to compute the integral error

Lorenzo

Hi,

thanks, for response, could you please tell me how exactly do I do that? Or where can I find the information?

Thanks,
Marek

Hi,

An example is provided in this previous post

Basically you need to get the full covariance matrix from the total fit (using TFitResult::GetCovarianceMatrix) and then extract the sub matrix for the parameters 0 and 1, and use this matrix in TF1::IntegralError.

From your code do something like :

TF1 *total = new TF1("total","expo(0)+expo(2)+expo(4)",0,5);
total->SetParameters(par);
TFitResultPtr r = hpT->Fit(total,"R");
TMatrixDSym covTot = r->GetCovarianceMatrix();    
// make sub matrix with first two elements
TMatrixDSym covGG(0,1, covTot.GetMatrixArray() );

TF1 *expoGG = new TF1("expoGG","expo",0.0,5);

expoGG->SetParameter(0,total->GetParameter(0));
expoGG->SetParameter(1,total->GetParameter(1));
expoGG->SetParError(0,total->GetParError(0));
expoGG->SetParError(1,total->GetParError(1));

double errNgg=expoGG->IntegralError(0.1,0.7, expoGG->GetParameters(), covGG-GetMatrixArray());

Lorenzo

Dear Lorenzo,

Thank you for help! I did as you advised, but I changed

TFitResultPtr r = hpT->Fit(total,“R”);

to

TFitResultPtr r = hpT->Fit(total,“S”);

and changed “-” to “.” in:

double errNgg=expoGG->IntegralError(0.1,0.7, expoGG->GetParameters(), covGG.GetMatrixArray());

the code works now, the error message is gone.

and I get the error of integral != 0.

Are the two changes I made correct?

Below I is my code.


   Double_t par[6];
   TF1 *expo1    = new TF1("expo1","expo",0.0,0.1);
   TF1 *expo2    = new TF1("expo2","expo",0.1,0.7);
   TF1 *expo3    = new TF1("expo3","expo",0.7,5);

   TF1 *total = new TF1("total","expo(0)+expo(2)+expo(4)",0,5);

   hpT->Fit(expo1,"R");
   hpT->Fit(expo2,"R");
   hpT->Fit(expo3,"R");
   
   expo1->GetParameters(&par[0]);
   expo2->GetParameters(&par[2]);
   expo3->GetParameters(&par[4]);
   total->SetParameters(par);
   
   hpT->Fit(total,"R");
   
  TF1 *expoGG    = new TF1("expoGG","expo",0.0,5);

  expoGG->SetParameter(0,total->GetParameter(0));
  expoGG->SetParameter(1,total->GetParameter(1));
  expoGG->SetParError(0,total->GetParError(0));
  expoGG->SetParError(1,total->GetParError(1));
   
   TFitResultPtr r = hpT->Fit(total,"S"); // changed from "R" to "S"
   TMatrixDSym covTot = r->GetCovarianceMatrix();   
   // make sub matrix with first two elements
   TMatrixDSym covGG(0,1, covTot.GetMatrixArray() );

   
   double Ngg=expoGG->Integral(0.1,0.7)/((x2-x1)/nbins);
   double errNgg=expoGG->IntegralError(0.1,0.7, expoGG->GetParameters(), covGG.GetMatrixArray());

So I guess it is correct now. Thank you!

Marek

Hi,

Sorry, I could not test the code, since I did not have your full macro. Yes, it should be fine now.
Just if you need the fit option “R” you should keep it and use “R S” as fit option.

Lorenzo

Thank you,

it seems I still have some problem: depending on the hpT histogram I am using the code sometimes gives the integral error = -nan. I noticed that when this happens the covariance matrix has many very small values inside, roughly smaller (or greater) than 0.005 (-0.005).
In addiction the values are few times too small compared to what is expected.

On the plot I have 3 different slopes in 3 regions (0-0.1 , 0.1-0.7, 0.7-5), so each exponent describes well only one region. But I am interested on the influence of 1st and 3rd exponent on the second region (0.1-0.7). So I guess that it might be the problem for the calculation of the integral error, as I calculate it in the region where the function does not describe the data any more? When I calculate the error on the region, where the function fits the data well, I do not get error =-nan.

Might it be that in this case this method will not work properly? And I should estimate the error in some different way?

Hi,

I am not sure I have well understood exactly your problem. Assuming that the covariance matrix does not containing a nan, the nan could be caused by either a problem in calculating the integral or by a problem in computing the derivative of the function with respect to the parameter. Can these be computed in the region you are interested ?

The alternative and better method than use IntegralError, is to re-parametrize your fit by using as parameter directly the integral of the function (e.g. the number of events). This requires to normalise the input function components. It is possible now to do it in ROOT using the TF1NormSum class. See the tutorial,
$ROOTSYS/tutorials/fit/fitNormSum.C

root.cern.ch/root/htmldoc/tutor … Sum.C.html

Lorenzo

Thanks for reply,

This is my code

using namespace RooFit ;
//Adopt this style
gROOT->Reset();



// ==> Observables
RooRealVar m("m","",0,10) ;
RooRealVar pT("pT","",0,10);
RooRealVar y("y","",-5,5) ;


void blabla()
{
  // load all data
   RooDataSet* dataPtYM = RooDataSet::read("data/dimuVecPPbMarek_UPCpPbDiTrkTree.txt",RooArgList(m,pT,y)) ;

  // submit study for each different sample
  doAllCases(dataPtYM);

}

void doAllCases(RooDataSet *dataPtYM)
{


int nbins=100;
float x1=0;
float x2=5;

float r1=0.1;
float r2=0.7;
float r3=5.0;

  TH1F *hpT = new TH1F("histo","",nbins,x1,x2);



  for (Int_t i=0; i< dataPtYM->numEntries(); i++){

  
  hpT->Fill( dataPtYM->get(i)->getRealValue("pT") * dataPtYM->get(i)->getRealValue("pT")  );

}
	// expo functions used for initial parameters for total function
   Double_t par[6];
   TF1 *expo1    = new TF1("expo1","expo",0.0,r1);
   TF1 *expo2    = new TF1("expo2","expo",r1,r2);
   TF1 *expo3    = new TF1("expo3","expo",r2,r3);

   TF1 *total = new TF1("total","expo(0)+expo(2)+expo(4)",0,5);
   total->SetLineColor(2);

   hpT->Fit(expo1,"R");
   hpT->Fit(expo2,"R");
   hpT->Fit(expo3,"R");

   expo1->GetParameters(&par[0]);
   expo2->GetParameters(&par[2]);
   expo3->GetParameters(&par[4]);
   total->SetParameters(par);

   hpT->Fit(total,"R");
   
////////////// for signal extraction //////////////////
 
TFitResultPtr r = hpT->Fit(total,"R S");
TMatrixDSym covTot = r->GetCovarianceMatrix(); 
   

// make sub matrix for the second expo function
TMatrixDSym covEX(2,3, covTot.GetMatrixArray() );

         TF1 *expoEX    = new TF1("expoEX","expo",0.0,5);
        // exporange->SetLineStyle(5);

  expoEX->SetParameter(0,total->GetParameter(2));
  expoEX->SetParameter(1,total->GetParameter(3));
  expoEX->SetParError(0,total->GetParError(2));
  expoEX->SetParError(1,total->GetParError(3));
 // hpT->Fit(expoNEX,"RB+");

float Nex=expoEX->Integral(0,5)/((x2-x1)/nbins);
double errNex=expoEX->IntegralError(0,5, expoEX->GetParameters(), covEX.GetMatrixArray());



/////////////// canvas for plots /////////////////

TCanvas *c1 = new TCanvas("c1","c1",600,600);
  c1->SetLogy();

  hpT->Draw("E");

 // c1->Update();

     

/////////// TERMINAL OUTPUT /////////////

cout << endl << "integral and error : " << Nex << "+/-" << errNex <<endl;




}

I have 3 different slopes in 3 regions, and want to get the integral under the second one, but on whole rande (0,5) and the error. The integral is ok, the error is “nan”. Sometimes it has some value, generally one or two orders of mag to small (depending on txt file I use). What is wrong with my code? The problem seems quite simple.

Marek