{
"cells": [
{
"cell_type": "markdown",
"id": "27894825",
"metadata": {},
"source": [
"# testUnfold1\n",
" Test program for the classes TUnfold and related.\n",
"\n",
" 1. Generate Monte Carlo and Data events\n",
" The events consist of\n",
" signal\n",
" background\n",
"\n",
" The signal is a resonance. It is generated with a Breit-Wigner,\n",
" smeared by a Gaussian\n",
"\n",
" 2. Unfold the data. The result is:\n",
" The background level\n",
" The shape of the resonance, corrected for detector effects\n",
"\n",
" Systematic errors from the MC shape variation are included\n",
" and propagated to the result\n",
"\n",
" 3. fit the unfolded distribution, including the correlation matrix\n",
"\n",
" 4. save six plots to a file `testUnfold1.ps`\n",
"```\n",
" 1 2 3\n",
" 4 5 6\n",
"```\n",
" 1. 2d-plot of the matrix describing the migrations\n",
" 2. generator-level distributions\n",
" - blue: unfolded data, total errors\n",
" - green: unfolded data, statistical errors\n",
" - red: generated data\n",
" - black: fit to green data points\n",
" 3. detector level distributions\n",
" - blue: unfolded data, folded back through the matrix\n",
" - black: Monte Carlo (with wrong peal position)\n",
" - blue: data\n",
" 4. global correlation coefficients\n",
" 5. $ \\chi^2 $ as a function of $ log(\\tau) $\n",
" the star indicates the final choice of $ \\tau $\n",
" 6. the L curve\n",
"\n",
"\n",
" **Version 17.6, in parallel to changes in TUnfold**\n",
"\n",
"#### History:\n",
" - Version 17.5, in parallel to changes in TUnfold\n",
" - Version 17.4, in parallel to changes in TUnfold\n",
" - Version 17.3, in parallel to changes in TUnfold\n",
" - Version 17.2, in parallel to changes in TUnfold\n",
" - Version 17.1, in parallel to changes in TUnfold\n",
" - Version 17.0, updated for using the classes TUnfoldDensity, TUnfoldBinning\n",
" - Version 16.1, parallel to changes in TUnfold\n",
" - Version 16.0, parallel to changes in TUnfold\n",
" - Version 15, with automated L-curve scan\n",
" - Version 14, with changes in TUnfoldSys.cxx\n",
" - Version 13, include test of systematic errors\n",
" - Version 12, catch error when defining the input\n",
" - Version 11, print chi**2 and number of degrees of freedom\n",
" - Version 10, with bug-fix in TUnfold.cxx\n",
" - Version 9, with bug-fix in TUnfold.cxx and TUnfold.h\n",
" - Version 8, with bug-fix in TUnfold.cxx and TUnfold.h\n",
" - Version 7, with bug-fix in TUnfold.cxx and TUnfold.h\n",
" - Version 6a, fix problem with dynamic array allocation under windows\n",
" - Version 6, bug-fixes in TUnfold.C\n",
" - Version 5, replace main() by testUnfold1()\n",
" - Version 4, with bug-fix in TUnfold.C\n",
" - Version 3, with bug-fix in TUnfold.C\n",
" - Version 2, with changed ScanLcurve() arguments\n",
" - Version 1, remove L curve analysis, use ScanLcurve() method instead\n",
" - Version 0, L curve analysis included here\n",
"\n",
" This file is part of TUnfold.\n",
"\n",
" TUnfold is free software: you can redistribute it and/or modify\n",
" it under the terms of the GNU General Public License as published by\n",
" the Free Software Foundation, either version 3 of the License, or\n",
" (at your option) any later version.\n",
"\n",
" TUnfold is distributed in the hope that it will be useful,\n",
" but WITHOUT ANY WARRANTY; without even the implied warranty of\n",
" MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n",
" GNU General Public License for more details.\n",
"\n",
" You should have received a copy of the GNU General Public License\n",
" along with TUnfold. If not, see .\n",
"\n",
"\n",
"\n",
"**Author:** Stefan Schmitt DESY, 14.10.2008 \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:21 PM."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d4e26161",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%%cpp -d\n",
"#include \n",
"#include \n",
"#include \n",
"#include \n",
"#include \n",
"#include \n",
"#include \n",
"#include \n",
"#include \n",
"\n",
"#include \"TUnfoldDensity.h\"\n",
"\n",
"\n",
"TRandom *rnd=nullptr;\n",
"\n",
"TH2 *gHistInvEMatrix;\n",
"\n",
"TVirtualFitter *gFitter=nullptr;"
]
},
{
"cell_type": "markdown",
"id": "55f181a7",
"metadata": {},
"source": [
" #define VERBOSE_LCURVE_SCAN\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4e621686",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%%cpp -d\n",
"void chisquare_corr(Int_t &npar, Double_t * /*gin */, Double_t &f, Double_t *u, Int_t /*flag */) {\n",
" // Minimization function for H1s using a Chisquare method\n",
" // only one-dimensional histograms are supported\n",
" // Correlated errors are taken from an external inverse covariance matrix\n",
" // stored in a 2-dimensional histogram\n",
" Double_t x;\n",
" TH1 *hfit = (TH1*)gFitter->GetObjectFit();\n",
" TF1 *f1 = (TF1*)gFitter->GetUserFunc();\n",
"\n",
" f1->InitArgs(&x,u);\n",
" npar = f1->GetNpar();\n",
" f = 0;\n",
"\n",
" Int_t npfit = 0;\n",
" Int_t nPoints=hfit->GetNbinsX();\n",
" Double_t *df=new Double_t[nPoints];\n",
" for (Int_t i=0;iGetBinCenter(i+1);\n",
" TF1::RejectPoint(kFALSE);\n",
" df[i] = f1->EvalPar(&x,u)-hfit->GetBinContent(i+1);\n",
" if (TF1::RejectedPoint()) df[i]=0.0;\n",
" else npfit++;\n",
" }\n",
" for (Int_t i=0;iGetBinContent(i+1,j+1);\n",
" }\n",
" }\n",
" delete[] df;\n",
" f1->SetNumberFitPoints(npfit);\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "dbc11be3",
"metadata": {},
"source": [
" Definition of a helper function: "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "630e2cb0",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%%cpp -d\n",
"Double_t bw_func(Double_t *x,Double_t *par) {\n",
" Double_t dm=x[0]-par[1];\n",
" return par[0]/(dm*dm+par[2]*par[2]);\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "4b169d60",
"metadata": {},
"source": [
" generate an event\n",
"output:\n",
"negative mass: background event\n",
"positive mass: signal event\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6a2b2903",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%%cpp -d\n",
"Double_t GenerateEvent(Double_t bgr, // relative fraction of background\n",
" Double_t mass, // peak position\n",
" Double_t gamma) // peak width\n",
"{\n",
" Double_t t;\n",
" if(rnd->Rndm()>bgr) {\n",
" // generate signal event\n",
" // with positive mass\n",
" do {\n",
" do {\n",
" t=rnd->Rndm();\n",
" } while(t>=1.0);\n",
" t=TMath::Tan((t-0.5)*TMath::Pi())*gamma+mass;\n",
" } while(t<=0.0);\n",
" return t;\n",
" } else {\n",
" // generate background event\n",
" // generate events following a power-law distribution\n",
" // f(E) = K * TMath::power((E0+E),N0)\n",
" static Double_t const E0=2.4;\n",
" static Double_t const N0=2.9;\n",
" do {\n",
" do {\n",
" t=rnd->Rndm();\n",
" } while(t>=1.0);\n",
" // the mass is returned negative\n",
" // In our example a convenient way to indicate it is a background event.\n",
" t= -(TMath::Power(1.-t,1./(1.-N0))-1.0)*E0;\n",
" } while(t>=0.0);\n",
" return t;\n",
" }\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "d42ede75",
"metadata": {},
"source": [
" smear the event to detector level\n",
"input:\n",
"mass on generator level (mTrue>0 !)\n",
"output:\n",
"mass on detector level\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "05cfd00e",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%%cpp -d\n",
"Double_t DetectorEvent(Double_t mTrue) {\n",
" // smear by double-gaussian\n",
" static Double_t frac=0.1;\n",
" static Double_t wideBias=0.03;\n",
" static Double_t wideSigma=0.5;\n",
" static Double_t smallBias=0.0;\n",
" static Double_t smallSigma=0.1;\n",
" if(rnd->Rndm()>frac) {\n",
" return rnd->Gaus(mTrue+smallBias,smallSigma);\n",
" } else {\n",
" return rnd->Gaus(mTrue+wideBias,wideSigma);\n",
" }\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "89d5ab56",
"metadata": {},
"source": [
"switch on histogram errors"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b85a7dd7",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"TH1::SetDefaultSumw2();"
]
},
{
"cell_type": "markdown",
"id": "9717c432",
"metadata": {},
"source": [
"show fit result"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3039a280",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"gStyle->SetOptFit(1111);"
]
},
{
"cell_type": "markdown",
"id": "5bd57735",
"metadata": {},
"source": [
"random generator"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8a6ba235",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"rnd=new TRandom3();"
]
},
{
"cell_type": "markdown",
"id": "38434fd0",
"metadata": {},
"source": [
"data and MC luminosity, cross-section"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "22f9d604",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"Double_t const luminosityData=100000;\n",
"Double_t const luminosityMC=1000000;\n",
"Double_t const crossSection=1.0;\n",
"\n",
"Int_t const nDet=250;\n",
"Int_t const nGen=100;\n",
"Double_t const xminDet=0.0;\n",
"Double_t const xmaxDet=10.0;\n",
"Double_t const xminGen=0.0;\n",
"Double_t const xmaxGen=10.0;"
]
},
{
"cell_type": "markdown",
"id": "893d0c1d",
"metadata": {},
"source": [
"============================================\n",
"generate MC distribution"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a26644f8",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"TH1D *histMgenMC=new TH1D(\"MgenMC\",\";mass(gen)\",nGen,xminGen,xmaxGen);\n",
"TH1D *histMdetMC=new TH1D(\"MdetMC\",\";mass(det)\",nDet,xminDet,xmaxDet);\n",
"TH2D *histMdetGenMC=new TH2D(\"MdetgenMC\",\";mass(det);mass(gen)\",\n",
" nDet,xminDet,xmaxDet,nGen,xminGen,xmaxGen);\n",
"Int_t neventMC=rnd->Poisson(luminosityMC*crossSection);\n",
"for(Int_t i=0;iFill(mGen,luminosityData/luminosityMC);\n",
" // reconstructed MC distribution (for comparison only)\n",
" histMdetMC->Fill(mDet,luminosityData/luminosityMC);\n",
"\n",
" // matrix describing how the generator input migrates to the\n",
" // reconstructed level. Unfolding input.\n",
" // NOTE on underflow/overflow bins:\n",
" // (1) the detector level under/overflow bins are used for\n",
" // normalisation (\"efficiency\" correction)\n",
" // in our toy example, these bins are populated from tails\n",
" // of the initial MC distribution.\n",
" // (2) the generator level underflow/overflow bins are\n",
" // unfolded. In this example:\n",
" // underflow bin: background events reconstructed in the detector\n",
" // overflow bin: signal events generated at masses > xmaxDet\n",
" // for the unfolded result these bins will be filled\n",
" // -> the background normalisation will be contained in the underflow bin\n",
" histMdetGenMC->Fill(mDet,mGen,luminosityData/luminosityMC);\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "88119557",
"metadata": {},
"source": [
"============================================\n",
"generate alternative MC\n",
"this will be used to derive a systematic error due to MC\n",
"parameter uncertainties"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "142eca70",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"TH2D *histMdetGenSysMC=new TH2D(\"MdetgenSysMC\",\";mass(det);mass(gen)\",\n",
" nDet,xminDet,xmaxDet,nGen,xminGen,xmaxGen);\n",
"neventMC=rnd->Poisson(luminosityMC*crossSection);\n",
"for(Int_t i=0;iFill(mDet,mGen,luminosityData/luminosityMC);\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "0f5ab3bf",
"metadata": {},
"source": [
"============================================\n",
"generate data distribution"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fec1fe8f",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"TH1D *histMgenData=new TH1D(\"MgenData\",\";mass(gen)\",nGen,xminGen,xmaxGen);\n",
"TH1D *histMdetData=new TH1D(\"MdetData\",\";mass(det)\",nDet,xminDet,xmaxDet);\n",
"Int_t neventData=rnd->Poisson(luminosityData*crossSection);\n",
"for(Int_t i=0;iFill(mGen);\n",
"\n",
" // reconstructed mass, unfolding input\n",
" histMdetData->Fill(mDet);\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "079036eb",
"metadata": {},
"source": [
"=========================================================================\n",
"divide by bin width to get density distributions"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c876191d",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"TH1D *histDensityGenData=new TH1D(\"DensityGenData\",\";mass(gen)\",\n",
" nGen,xminGen,xmaxGen);\n",
"TH1D *histDensityGenMC=new TH1D(\"DensityGenMC\",\";mass(gen)\",\n",
" nGen,xminGen,xmaxGen);\n",
"for(Int_t i=1;i<=nGen;i++) {\n",
" histDensityGenData->SetBinContent(i,histMgenData->GetBinContent(i)/\n",
" histMgenData->GetBinWidth(i));\n",
" histDensityGenMC->SetBinContent(i,histMgenMC->GetBinContent(i)/\n",
" histMgenMC->GetBinWidth(i));\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "c1034ff9",
"metadata": {},
"source": [
"=========================================================================\n",
"set up the unfolding\n",
"define migration matrix"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cdf5aa21",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"TUnfoldDensity unfold(histMdetGenMC,TUnfold::kHistMapOutputVert);"
]
},
{
"cell_type": "markdown",
"id": "7a26c8e9",
"metadata": {},
"source": [
"define input and bias scheme\n",
"do not use the bias, because MC peak may be at the wrong place\n",
"watch out for error codes returned by the SetInput method\n",
"errors larger or equal 10000 are fatal:\n",
"the data points specified as input are not sufficient to constrain the\n",
"unfolding process"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2fedd58d",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"if(unfold.SetInput(histMdetData)>=10000) {\n",
" std::cout<<\"Unfolding result may be wrong\\n\";\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "17a244af",
"metadata": {},
"source": [
"========================================================================\n",
"the unfolding is done here\n",
"\n",
"scan L curve and find best point"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "62172fa3",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"Int_t nScan=30;"
]
},
{
"cell_type": "markdown",
"id": "0c704aca",
"metadata": {},
"source": [
"use automatic L-curve scan: start with taumin=taumax=0.0"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "11c3e385",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"Double_t tauMin=0.0;\n",
"Double_t tauMax=0.0;\n",
"Int_t iBest;\n",
"TSpline *logTauX,*logTauY;\n",
"TGraph *lCurve;"
]
},
{
"cell_type": "markdown",
"id": "cd1a3364",
"metadata": {},
"source": [
"if required, report Info messages (for debugging the L-curve scan)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6e74b86a",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"#ifdef VERBOSE_LCURVE_SCAN\n",
"Int_t oldinfo=gErrorIgnoreLevel;\n",
"gErrorIgnoreLevel=kInfo;\n",
"#endif"
]
},
{
"cell_type": "markdown",
"id": "af45bc02",
"metadata": {},
"source": [
"this method scans the parameter tau and finds the kink in the L curve\n",
"finally, the unfolding is done for the best choice of tau"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "16f299a5",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"iBest=unfold.ScanLcurve(nScan,tauMin,tauMax,&lCurve,&logTauX,&logTauY);"
]
},
{
"cell_type": "markdown",
"id": "111df930",
"metadata": {},
"source": [
"if required, switch to previous log-level"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3ad7d251",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"#ifdef VERBOSE_LCURVE_SCAN\n",
"gErrorIgnoreLevel=oldinfo;\n",
"#endif"
]
},
{
"cell_type": "markdown",
"id": "bb2c6d53",
"metadata": {},
"source": [
"==========================================================================\n",
"define a correlated systematic error\n",
"for example, assume there is a 10% correlated error for all reconstructed\n",
"masses larger than 7"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "be26e98c",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"Double_t SYS_ERROR1_MSTART=6;\n",
"Double_t SYS_ERROR1_SIZE=0.1;\n",
"TH2D *histMdetGenSys1=new TH2D(\"Mdetgensys1\",\";mass(det);mass(gen)\",\n",
" nDet,xminDet,xmaxDet,nGen,xminGen,xmaxGen);\n",
"for(Int_t i=0;i<=nDet+1;i++) {\n",
" if(histMdetData->GetBinCenter(i)>=SYS_ERROR1_MSTART) {\n",
" for(Int_t j=0;j<=nGen+1;j++) {\n",
" histMdetGenSys1->SetBinContent(i,j,SYS_ERROR1_SIZE);\n",
" }\n",
" }\n",
"}\n",
"unfold.AddSysError(histMdetGenSysMC,\"SYSERROR_MC\",TUnfold::kHistMapOutputVert,\n",
" TUnfoldSys::kSysErrModeMatrix);\n",
"unfold.AddSysError(histMdetGenSys1,\"SYSERROR1\",TUnfold::kHistMapOutputVert,\n",
" TUnfoldSys::kSysErrModeRelative);"
]
},
{
"cell_type": "markdown",
"id": "f427d646",
"metadata": {},
"source": [
"==========================================================================\n",
"print some results"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f6a511a5",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"std::cout<<\"tau=\"<GetKnot(iBest,t[0],x[0]);\n",
"logTauY->GetKnot(iBest,t[0],y[0]);\n",
"TGraph *bestLcurve=new TGraph(1,x,y);\n",
"TGraph *bestLogTauLogChi2=new TGraph(1,t,x);"
]
},
{
"cell_type": "markdown",
"id": "000ca521",
"metadata": {},
"source": [
"==========================================================================\n",
"retrieve results into histograms"
]
},
{
"cell_type": "markdown",
"id": "7e24b6f9",
"metadata": {},
"source": [
"get unfolded distribution"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ad40c3e6",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"TH1 *histMunfold=unfold.GetOutput(\"Unfolded\");"
]
},
{
"cell_type": "markdown",
"id": "2575454d",
"metadata": {},
"source": [
"get unfolding result, folded back"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bc34b778",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"TH1 *histMdetFold=unfold.GetFoldedOutput(\"FoldedBack\");"
]
},
{
"cell_type": "markdown",
"id": "88fa651a",
"metadata": {},
"source": [
"get error matrix (input distribution [stat] errors only)\n",
"TH2D *histEmatData=unfold.GetEmatrix(\"EmatData\");"
]
},
{
"cell_type": "markdown",
"id": "b8eaae4b",
"metadata": {},
"source": [
"get total error matrix:\n",
"migration matrix uncorrelated and correlated systematic errors\n",
"added in quadrature to the data statistical errors"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "04d6bcbc",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"TH2 *histEmatTotal=unfold.GetEmatrixTotal(\"EmatTotal\");"
]
},
{
"cell_type": "markdown",
"id": "dbefade3",
"metadata": {},
"source": [
"create data histogram with the total errors"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9b114dde",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"TH1D *histTotalError=\n",
" new TH1D(\"TotalError\",\";mass(gen)\",nGen,xminGen,xmaxGen);\n",
"for(Int_t bin=1;bin<=nGen;bin++) {\n",
" histTotalError->SetBinContent(bin,histMunfold->GetBinContent(bin));\n",
" histTotalError->SetBinError\n",
" (bin,TMath::Sqrt(histEmatTotal->GetBinContent(bin,bin)));\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "10c13b58",
"metadata": {},
"source": [
"get global correlation coefficients\n",
"for this calculation one has to specify whether the\n",
"underflow/overflow bins are included or not\n",
"default: include all bins\n",
"here: exclude underflow and overflow bins"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "10254ac0",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"TH1 *histRhoi=unfold.GetRhoItotal(\"rho_I\",\n",
" nullptr, // use default title\n",
" nullptr, // all distributions\n",
" \"*[UO]\", // discard underflow and overflow bins on all axes\n",
" kTRUE, // use original binning\n",
" &gHistInvEMatrix // store inverse of error matrix\n",
" );"
]
},
{
"cell_type": "markdown",
"id": "def39e1a",
"metadata": {},
"source": [
"======================================================================\n",
"fit Breit-Wigner shape to unfolded data, using the full error matrix\n",
"here we use a \"user\" chi**2 function to take into account\n",
"the full covariance matrix"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "15963135",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"gFitter=TVirtualFitter::Fitter(histMunfold);\n",
"gFitter->SetFCN(chisquare_corr);\n",
"\n",
"TF1 *bw=new TF1(\"bw\",bw_func,xminGen,xmaxGen,3);\n",
"bw->SetParameter(0,1000.);\n",
"bw->SetParameter(1,3.8);\n",
"bw->SetParameter(2,0.2);"
]
},
{
"cell_type": "markdown",
"id": "a4b4b5c5",
"metadata": {},
"source": [
"for (wrong!) fitting without correlations, drop the option \"U\"\n",
"here."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "57996f21",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"histMunfold->Fit(bw,\"UE\");"
]
},
{
"cell_type": "markdown",
"id": "8af2d969",
"metadata": {},
"source": [
"=====================================================================\n",
"plot some histograms"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ff0f8c46",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"TCanvas output;\n",
"output.Divide(3,2);"
]
},
{
"cell_type": "markdown",
"id": "df02c30f",
"metadata": {},
"source": [
"Show the matrix which connects input and output\n",
"There are overflow bins at the bottom, not shown in the plot\n",
"These contain the background shape.\n",
"The overflow bins to the left and right contain\n",
"events which are not reconstructed. These are necessary for proper MC\n",
"normalisation"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c8264925",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"output.cd(1);\n",
"histMdetGenMC->Draw(\"BOX\");"
]
},
{
"cell_type": "markdown",
"id": "ad7574d1",
"metadata": {},
"source": [
"draw generator-level distribution:\n",
"data (red) [for real data this is not available]\n",
"MC input (black) [with completely wrong peak position and shape]\n",
"unfolded data (blue)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0b88575f",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"output.cd(2);\n",
"histTotalError->SetLineColor(kBlue);\n",
"histTotalError->Draw(\"E\");\n",
"histMunfold->SetLineColor(kGreen);\n",
"histMunfold->Draw(\"SAME E1\");\n",
"histDensityGenData->SetLineColor(kRed);\n",
"histDensityGenData->Draw(\"SAME\");\n",
"histDensityGenMC->Draw(\"SAME HIST\");"
]
},
{
"cell_type": "markdown",
"id": "98c2842e",
"metadata": {},
"source": [
"show detector level distributions\n",
"data (red)\n",
"MC (black) [with completely wrong peak position and shape]\n",
"unfolded data (blue)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "06baef6d",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"output.cd(3);\n",
"histMdetFold->SetLineColor(kBlue);\n",
"histMdetFold->Draw();\n",
"histMdetMC->Draw(\"SAME HIST\");\n",
"\n",
"TH1 *histInput=unfold.GetInput(\"Minput\",\";mass(det)\");\n",
"\n",
"histInput->SetLineColor(kRed);\n",
"histInput->Draw(\"SAME\");"
]
},
{
"cell_type": "markdown",
"id": "75e4a260",
"metadata": {},
"source": [
"show correlation coefficients"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ced0d3d0",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"output.cd(4);\n",
"histRhoi->Draw();"
]
},
{
"cell_type": "markdown",
"id": "98b98d8f",
"metadata": {},
"source": [
"show tau as a function of chi**2"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "dab69ff1",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"output.cd(5);\n",
"logTauX->Draw();\n",
"bestLogTauLogChi2->SetMarkerColor(kRed);\n",
"bestLogTauLogChi2->Draw(\"*\");"
]
},
{
"cell_type": "markdown",
"id": "7fd434aa",
"metadata": {},
"source": [
"show the L curve"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9c3170c9",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"output.cd(6);\n",
"lCurve->Draw(\"AL\");\n",
"bestLcurve->SetMarkerColor(kRed);\n",
"bestLcurve->Draw(\"*\");\n",
"\n",
"output.SaveAs(\"testUnfold1.ps\");\n",
"\n",
"return 0;"
]
},
{
"cell_type": "markdown",
"id": "59c1f1d4",
"metadata": {},
"source": [
"Draw all canvases "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "914813f3",
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"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
}