{
"cells": [
{
"cell_type": "markdown",
"id": "5c3a785b",
"metadata": {},
"source": [
"# ErrorIntegral\n",
"Estimate the error in the integral of a fitted function\n",
"taking into account the errors in the parameters resulting from the fit.\n",
"The error is estimated also using the correlations values obtained from\n",
"the fit\n",
"\n",
"run the macro doing:\n",
"\n",
"```cpp\n",
" .x ErrorIntegral.C\n",
"```\n",
"\n",
"After having computed the integral and its error using the integral and the integral\n",
"error using the generic functions TF1::Integral and TF1::IntegralError, we compute\n",
"the integrals and its error analytically using the fact that the fitting function is\n",
"\n",
"Therefore we have:\n",
" - integral in [0,1] : `ic = p[1]* (1-std::cos(p[0]) )/p[0]`\n",
" - derivative of integral with respect to p0:\n",
" `c0c = p[1] * (std::cos(p[0]) + p[0]*std::sin(p[0]) -1.)/p[0]/p[0]`\n",
" - derivative of integral with respect to p1:\n",
" `c1c = (1-std::cos(p[0]) )/p[0]`\n",
"\n",
"and then we can compute the integral error using error propagation and the covariance\n",
"matrix for the parameters p obtained from the fit.\n",
"\n",
"integral error : `sic = std::sqrt( c0c*c0c * covMatrix(0,0) + c1c*c1c * covMatrix(1,1) + 2.* c0c*c1c * covMatrix(0,1))`\n",
"\n",
"Note that, if possible, one should fit directly the function integral, which are the\n",
"number of events of the different components (e.g. signal and background).\n",
"In this way one obtains a better and more correct estimate of the integrals uncertainties,\n",
"since they are obtained directly from the fit without using the approximation of error propagation.\n",
"This is possible in ROOT. when using the TF1NormSum class, see the tutorial fitNormSum.C\n",
"\n",
"\n",
"\n",
"\n",
"**Author:** Lorenzo Moneta \n",
"This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Tuesday, March 19, 2024 at 07:08 PM."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "0183de2e",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:08:25.557400Z",
"iopub.status.busy": "2024-03-19T19:08:25.556983Z",
"iopub.status.idle": "2024-03-19T19:08:25.623611Z",
"shell.execute_reply": "2024-03-19T19:08:25.619109Z"
}
},
"outputs": [],
"source": [
"%%cpp -d\n",
"#include \"TF1.h\"\n",
"#include \"TH1D.h\"\n",
"#include \"TFitResult.h\"\n",
"#include \"TMath.h\"\n",
"#include \n",
"#include \n",
"#include \n",
"\n",
"TF1 * fitFunc; // fit function pointer\n",
"\n",
"const int NPAR = 2; // number of function parameters;"
]
},
{
"cell_type": "markdown",
"id": "4c44d19d",
"metadata": {},
"source": [
" ____________________________________________________________________\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "26d7967c",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:08:25.630772Z",
"iopub.status.busy": "2024-03-19T19:08:25.630142Z",
"iopub.status.idle": "2024-03-19T19:08:25.675298Z",
"shell.execute_reply": "2024-03-19T19:08:25.659735Z"
}
},
"outputs": [],
"source": [
"%%cpp -d\n",
"double f(double * x, double * p) {\n",
" // function used to fit the data\n",
" return p[1]*TMath::Sin( p[0] * x[0] );\n",
"}"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "a9bf92b7",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:08:25.690912Z",
"iopub.status.busy": "2024-03-19T19:08:25.690529Z",
"iopub.status.idle": "2024-03-19T19:08:27.475813Z",
"shell.execute_reply": "2024-03-19T19:08:27.474121Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"****************************************\n",
"Minimizer is Minuit2 / Migrad\n",
"Chi2 = 49.5952\n",
"NDf = 48\n",
"Edm = 1.61787e-06\n",
"NCalls = 58\n",
"p0 = 3.13205 +/- 0.0312726 \n",
"p1 = 29.7625 +/- 1.00876 \n",
"Covariance matrix from the fit \n",
"2x2 matrix is as follows\n",
"\n",
" | 0 | 1 |\n",
"-------------------------------\n",
" 0 | 0.000978 0.009147 \n",
" 1 | 0.009147 1.018 \n",
"\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Info in : created default TCanvas with name c1\n"
]
}
],
"source": [
"fitFunc = new TF1(\"f\",f,0,1,NPAR);\n",
"TH1D * h1 = new TH1D(\"h1\",\"h1\",50,0,1);\n",
"\n",
"double par[NPAR] = { 3.14, 1.};\n",
"fitFunc->SetParameters(par);\n",
"\n",
"h1->FillRandom(\"f\",1000); // fill histogram sampling fitFunc\n",
"fitFunc->SetParameter(0,3.); // vary a little the parameters\n",
"auto fitResult = h1->Fit(fitFunc,\"S\"); // fit the histogram and get fit result pointer\n",
"\n",
"h1->Draw();\n",
"\n",
"/* calculate the integral*/\n",
"double integral = fitFunc->Integral(0,1);\n",
"\n",
"auto covMatrix = fitResult->GetCovarianceMatrix();\n",
"std::cout << \"Covariance matrix from the fit \";\n",
"covMatrix.Print();"
]
},
{
"cell_type": "markdown",
"id": "c8459c8b",
"metadata": {},
"source": [
"need to pass covariance matrix to fit result.\n",
"Parameters values are are stored inside the function but we can also retrieve from TFitResult"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "40ddc5aa",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:08:27.482505Z",
"iopub.status.busy": "2024-03-19T19:08:27.482176Z",
"iopub.status.idle": "2024-03-19T19:08:27.734740Z",
"shell.execute_reply": "2024-03-19T19:08:27.688670Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Integral = 19.0047 +/- 0.616472\n"
]
}
],
"source": [
"double sigma_integral = fitFunc->IntegralError(0,1, fitResult->GetParams() , covMatrix.GetMatrixArray());\n",
"\n",
"std::cout << \"Integral = \" << integral << \" +/- \" << sigma_integral\n",
" << std::endl;"
]
},
{
"cell_type": "markdown",
"id": "f7151d92",
"metadata": {},
"source": [
"estimated integral and error analytically"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "fe0800d0",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:08:27.745114Z",
"iopub.status.busy": "2024-03-19T19:08:27.744727Z",
"iopub.status.idle": "2024-03-19T19:08:27.954634Z",
"shell.execute_reply": "2024-03-19T19:08:27.952074Z"
}
},
"outputs": [],
"source": [
"double * p = fitFunc->GetParameters();\n",
"double ic = p[1]* (1-std::cos(p[0]) )/p[0];\n",
"double c0c = p[1] * (std::cos(p[0]) + p[0]*std::sin(p[0]) -1.)/p[0]/p[0];\n",
"double c1c = (1-std::cos(p[0]) )/p[0];"
]
},
{
"cell_type": "markdown",
"id": "58d4599f",
"metadata": {},
"source": [
"estimated error with correlations"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "efa9ee9c",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:08:27.963745Z",
"iopub.status.busy": "2024-03-19T19:08:27.963323Z",
"iopub.status.idle": "2024-03-19T19:08:28.170616Z",
"shell.execute_reply": "2024-03-19T19:08:28.169412Z"
}
},
"outputs": [],
"source": [
"double sic = std::sqrt( c0c*c0c * covMatrix(0,0) + c1c*c1c * covMatrix(1,1)\n",
" + 2.* c0c*c1c * covMatrix(0,1));\n",
"\n",
"if ( std::fabs(sigma_integral-sic) > 1.E-6*sic )\n",
" std::cout << \" ERROR: test failed : different analytical integral : \"\n",
" << ic << \" +/- \" << sic << std::endl;"
]
},
{
"cell_type": "markdown",
"id": "92f13afe",
"metadata": {},
"source": [
"Draw all canvases "
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "5d64a611",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:08:28.182364Z",
"iopub.status.busy": "2024-03-19T19:08:28.182051Z",
"iopub.status.idle": "2024-03-19T19:08:28.890204Z",
"shell.execute_reply": "2024-03-19T19:08:28.889167Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"
\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%jsroot on\n",
"gROOT->GetListOfCanvases()->Draw()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "ROOT C++",
"language": "c++",
"name": "root"
},
"language_info": {
"codemirror_mode": "text/x-c++src",
"file_extension": ".C",
"mimetype": " text/x-c++src",
"name": "c++"
}
},
"nbformat": 4,
"nbformat_minor": 5
}