Get Covariance Matrix without printing

Hey, all! I need some help. I am fitting histograms with a 2nd-order polynomial plus a gaussian, and also integrating the polynomial. In order to extract the error in the integral, I am using GetCovarianceMatrix, and getting the submatrix corresponding to the specific parameters of the polynomial, to get the error in the integral by using the function IntegralError. My problem is that every time I get the covariance matrix, I am also printing the entire matrix, which may become memory-costly on my simple laptop.

How can I avoid printing anything when getting the covariance matrix?

Alternatively, is there another way to get the integral error? IntegralError can’t seem to extract the error in the integral unless I integrate the entire fitted function, as opposed to just the polynomial function.

Below is my function definition that fits a histogram with a (Background) Polynomial + (Signal) Gaussian (pre-defined functions), and integrates (if I tell it to):

Double_t fitGAUS(TH1D *histo, Double_t low, Double_t high, Double_t term0, Double_t term1,Double_t term2, Double_t mean, Double_t width, Double_t dleft, Double_t dright, TString intgrl)
{
  /// create a TF1 with the range from tlo to thi and 6 parameters
  TF1 *fitFcn = new TF1("fitFcn",fitFunction, low, high, 6);

  /// Set Parameters
  fitFcn->SetParameter(0, term0); // 					pol c0
  fitFcn->SetParameter(1, term1); // pol c1
  fitFcn->SetParameter(2, term2); // pol c2
  fitFcn->SetParLimits(3, 0,100000); // peak
  fitFcn->SetParameter(4, mean); // mean
  fitFcn->SetParLimits(4,mean-width, mean+width); // mean limits
  fitFcn->SetParameter(5, width);   // uncertainty
  fitFcn->SetParLimits(5, 0.003, 6*width);   // uncertainty limits

  // Perform Fitting
  histo->Fit("fitFcn","NQ","",low,high);

  // Get fit results
  TFitResultPtr fitResult = histo->Fit("fitFcn","NRS");

  // Get total covariance matrix (will be used for integral error)
  TMatrixDSym covFit = fitResult->GetCovarianceMatrix();

  // Get covariance sub-matrices for signal and background
  TMatrixDSym covBackground(0, 2, covFit.GetMatrixArray() );
  TMatrixDSym covSignal(3, 5, covFit.GetMatrixArray() );

  // Define Gaussian and Background functions, and extract parameters and par-errors
  TF1 *backFcn = new TF1(Form("backFcn"), background,low,high,3);
  TF1 *signalFcn = new TF1(Form("signalFcn"), GaussianPeak,low,high,3);
  Double_t par[6];

  fitFcn->GetParameters(par);
  backFcn->SetParameter(0, par[0]);
  backFcn->SetParameter(1, par[1]);
  backFcn->SetParameter(2, par[2]);
  signalFcn->SetParameter(0, par[3]);
  signalFcn->SetParameter(1, par[4]);
  signalFcn->SetParameter(2, par[5]);

  backFcn->SetParError(0, fitFcn->GetParError(0) );
  backFcn->SetParError(1, fitFcn->GetParError(1) );
  backFcn->SetParError(2, fitFcn->GetParError(2) );
  signalFcn->SetParError(0, fitFcn->GetParError(3) );
  signalFcn->SetParError(1, fitFcn->GetParError(4) );
  signalFcn->SetParError(2, fitFcn->GetParError(5) );

  /////// Integrals /////////

  Double_t binw        = histo->GetBinWidth(1);
  Double_t IntGauss    = Int_t (signalFcn->Integral(dleft,dright)/binw);
  Double_t IntGaussErr = Int_t (signalFcn->IntegralError(dleft, dright, signalFcn->GetParameters(), covSignal.GetMatrixArray() )/binw);
  Double_t IntPoly     = Int_t (backFcn->Integral(dleft,dright)/binw);
  Double_t IntPolyErr  = Int_t (backFcn->IntegralError(dleft, dright, backFcn->GetParameters(), covBackground.GetMatrixArray() )/binw);
  Double_t IntFit      = Int_t (fitFcn->Integral(dleft,dright)/binw);
  Double_t IntFitErr   = Int_t (fitFcn->IntegralError(dleft, dright)/binw);

  if (intgrl == "G") {return IntGaus;}
  else if (intgrl == "Gerror" ) {return IntGaussErr;}

  else if (intgrl == "P") {return IntPoly;}
  else if (intgrl == "Perror" ) {return IntPolyErr;}

  else if (intgrl == "Fit") {return IntFit;}
  else if (intgrl == "FitError" ) {return IntFitErr;}

  else {return 0;}
}

Hi,

I don’t understand which function print the covariance matrix, FitResult::GetCovarianceMatrix() ?

The best way to get the error on the background events is to fit directly them instead of computing the integral.
This is possible by normalising the functions. This can be done automatically in ROOT now using the NSUM operator of TF1.
Example is that he tutorial tutorials/fit/fitNormSum.C or simpler by just doing (e.g. for a gaussian + exponential)

auto f1 = new TF1("f1","NSUM([sg] * gaus, [bg] * expo", xmin, xmax);

Best Regards

Lorenzo

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.