{ "cells": [ { "cell_type": "markdown", "id": "9136481c", "metadata": {}, "source": [ "# testUnfold7c\n", "Test program for the classes TUnfoldDensity and TUnfoldBinning.\n", "\n", "A toy test of the TUnfold package\n", "\n", "\n", "This example is documented in conference proceedings:\n", "\n", " arXiv:1611.01927\n", " 12th Conference on Quark Confinement and the Hadron Spectrum (Confinement XII)\n", "\n", "This is an example of unfolding a one-dimensional distribution. It compares\n", "various unfolding methods:\n", "\n", " matrix inversion, template fit, L-curve scan,\n", " iterative unfolding, etc\n", "\n", "Further details can be found in talk by S.Schmitt at:\n", "\n", " XII Quark Confinement and the Hadron Spectrum\n", " 29.8. - 3.9.2016 Thessaloniki, Greece\n", " statictics session (+proceedings)\n", "\n", "The example comprises several macros\n", " - testUnfold7a.C create root files with TTree objects for\n", " signal, background and data\n", " - write files testUnfold7_signal.root\n", " testUnfold7_background.root\n", " testUnfold7_data.root\n", "\n", " - testUnfold7b.C loop over trees and fill histograms based on the\n", " TUnfoldBinning objects\n", " - read testUnfold7binning.xml\n", " testUnfold7_signal.root\n", " testUnfold7_background.root\n", " testUnfold7_data.root\n", "\n", " - write testUnfold7_histograms.root\n", "\n", " - testUnfold7c.C run the unfolding\n", " - read testUnfold7_histograms.root\n", " - write testUnfold7_result.root\n", " - write many histograms, to compare various unfolding methods\n", "\n", "\n", " **Version 17.6, in parallel to changes in TUnfold**\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, May 19, 2026 at 08:11 PM." ] }, { "cell_type": "code", "execution_count": null, "id": "3b82e323", "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", "#include \n", "#include \n", "#include \n", "#include \n", "#include \n", "#include \n", "#include \n", "#include \n", "#include \n", "#include \n", "#include \n", "#include \n", "#include \n", "#include \"TUnfoldDensity.h\"\n", "\n", "using std::vector, std::pair, std::cout;\n", "\n", "\n", "#define TEST_INPUT_COVARIANCE\n", "\n", "void CreateHistogramCopies(TH1 *h[3],TUnfoldBinning const *binning);\n", "void CreateHistogramCopies(TH2 *h[3],TUnfoldBinning const *binningX);\n", "\n", "TH2 *AddOverflowXY(TH2 *h,double widthX,double widthY);\n", "TH1 *AddOverflowX(TH1 *h,double width);\n", "\n", "void DrawOverflowX(TH1 *h,double posy);\n", "void DrawOverflowY(TH1 *h,double posx);\n", "\n", "\n", "double const kLegendFontSize=0.05;\n", "int kNbinC=0;\n", "\n", "void DrawPadProbability(TH2 *h);\n", "void DrawPadEfficiency(TH1 *h);\n", "void DrawPadReco(TH1 *histMcRec,TH1 *histMcbgrRec,TH1 *histDataRec,\n", " TH1 *histDataUnfold,TH2 *histProbability,TH2 *histRhoij);\n", "void DrawPadTruth(TH1 *histMcsigGen,TH1 *histDataGen,TH1 *histDataUnfold,\n", " char const *text=nullptr,double tau=0.0,vector const *r=nullptr,\n", " TF1 *f=nullptr);\n", "void DrawPadCorrelations(TH2 *h,\n", " vector > > const *table);\n", "\n", "TFitResultPtr DoFit(TH1 *h,TH2 *rho,TH1 *truth,char const *text,\n", " vector > > &table,int niter=0);\n", "\n", "void GetNiterGraphs(int iFirst,int iLast,vector > > const &table,int color,\n", " TGraph *graph[4],int style);\n", "void GetNiterHist(int ifit,vector > > const &table,\n", " TH1 *hist[4],int color,int style,int fillStyle);\n", "\n", "#ifdef WITH_IDS\n", "void IDSfirst(TVectorD *data, TVectorD *dataErr, TMatrixD *A_, Double_t lambdaL_, TVectorD* &unfres1IDS_,TVectorD *&soustr);\n", "\n", "void IDSiterate(TVectorD *data, TVectorD *dataErr, TMatrixD *A_,TMatrixD *Am_,\n", " Double_t lambdaU_, Double_t lambdaM_, Double_t lambdaS_,\n", " TVectorD* &unfres2IDS_ ,TVectorD *&soustr);\n", "#endif\n", "\n", "TRandom3 *g_rnd;\n", "\n", "TH1 *g_fcnHist=nullptr;\n", "TMatrixD *g_fcnMatrix=nullptr;\n", "\n", "#ifdef WITH_IDS" ] }, { "cell_type": "markdown", "id": "e82bd003", "metadata": {}, "source": [ " Definition of a helper function: " ] }, { "cell_type": "code", "execution_count": null, "id": "8fe7d736", "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%cpp -d\n", "void GetNiterGraphs(int iFirst,int iLast,vector > > const &table,int color,\n", " TGraph *graph[4],int style) {\n", " TVectorD niter(iLast-iFirst);\n", " TVectorD eniter(iLast-iFirst);\n", " TVectorD chi2(iLast-iFirst);\n", " TVectorD gcor(iLast-iFirst);\n", " TVectorD mean(iLast-iFirst);\n", " TVectorD emean(iLast-iFirst);\n", " TVectorD sigma(iLast-iFirst);\n", " TVectorD esigma(iLast-iFirst);\n", " for(int ifit=iFirst;ifit const &r=table[ifit].second;\n", " niter(ifit-iFirst)=r[4];\n", " chi2(ifit-iFirst)=r[2]/r[3];\n", " gcor(ifit-iFirst)=r[5];\n", " TF1 const *f=table[ifit].first;\n", " mean(ifit-iFirst)=f->GetParameter(1);\n", " emean(ifit-iFirst)=f->GetParError(1);\n", " sigma(ifit-iFirst)=f->GetParameter(2);\n", " esigma(ifit-iFirst)=f->GetParError(2);\n", " }\n", " graph[0]=new TGraph(niter,chi2);\n", " graph[1]=new TGraph(niter,gcor);\n", " graph[2]=new TGraphErrors(niter,mean,eniter,emean);\n", " graph[3]=new TGraphErrors(niter,sigma,eniter,esigma);\n", " for(int g=0;g<4;g++) {\n", " if(graph[g]) {\n", " graph[g]->SetLineColor(color);\n", " graph[g]->SetMarkerColor(color);\n", " graph[g]->SetMarkerStyle(style);\n", " }\n", " }\n", "}" ] }, { "cell_type": "markdown", "id": "377d37ab", "metadata": {}, "source": [ " Definition of a helper function: " ] }, { "cell_type": "code", "execution_count": null, "id": "666ea7cf", "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%cpp -d\n", "void GetNiterHist(int ifit,vector > > const &table,\n", " TH1 *hist[4],int color,int style,int fillStyle) {\n", " vector const &r=table[ifit].second;\n", " TF1 const *f=table[ifit].first;\n", " hist[0]=new TH1D(table[ifit].first->GetName()+TString(\"_chi2\"),\n", " \";iteration;unfold-truth #chi^{2}/N_{D.F.}\",1,0.2,1500.);\n", " hist[0]->SetBinContent(1,r[2]/r[3]);\n", "\n", " hist[1]=new TH1D(table[ifit].first->GetName()+TString(\"_gcor\"),\n", " \";iteration;avg(#rho_{i})\",1,0.2,1500.);\n", " hist[1]->SetBinContent(1,r[5]);\n", " hist[2]=new TH1D(table[ifit].first->GetName()+TString(\"_mu\"),\n", " \";iteration;parameter #mu\",1,0.2,1500.);\n", " hist[2]->SetBinContent(1,f->GetParameter(1));\n", " hist[2]->SetBinError(1,f->GetParError(1));\n", " hist[3]=new TH1D(table[ifit].first->GetName()+TString(\"_sigma\"),\n", " \";iteration;parameter #sigma\",1,0.2,1500.);\n", " hist[3]->SetBinContent(1,f->GetParameter(2));\n", " hist[3]->SetBinError(1,f->GetParError(2));\n", " for(int h=0;h<4;h++) {\n", " if(hist[h]) {\n", " hist[h]->SetLineColor(color);\n", " hist[h]->SetLineStyle(style);\n", " if( hist[h]->GetBinError(1)>0.0) {\n", " hist[h]->SetFillColor(color-10);\n", " hist[h]->SetFillStyle(fillStyle);\n", " }\n", " hist[h]->SetMarkerStyle(0);\n", " }\n", " }\n", "}" ] }, { "cell_type": "markdown", "id": "86b756e6", "metadata": {}, "source": [ " Definition of a helper function: " ] }, { "cell_type": "code", "execution_count": null, "id": "7146cb46", "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%cpp -d\n", "void CreateHistogramCopies(TH1 *h[3],TUnfoldBinning const *binning) {\n", " TString baseName(h[0]->GetName());\n", " Int_t *binMap;\n", " h[1]=binning->CreateHistogram(baseName+\"_axis\",kTRUE,&binMap);\n", " h[2]=(TH1 *)h[1]->Clone(baseName+\"_binw\");\n", " Int_t nMax=binning->GetEndBin()+1;\n", " for(Int_t iSrc=0;iSrcGetBinContent(iSrc)+h[1]->GetBinContent(iDest);\n", " double e=TMath::Hypot(h[0]->GetBinError(iSrc),h[1]->GetBinError(iDest));\n", " h[1]->SetBinContent(iDest,c);\n", " h[1]->SetBinError(iDest,e);\n", " h[2]->SetBinContent(iDest,c);\n", " h[2]->SetBinError(iDest,e);\n", " }\n", " for(int iDest=0;iDest<=h[2]->GetNbinsX()+1;iDest++) {\n", " double c=h[2]->GetBinContent(iDest);\n", " double e=h[2]->GetBinError(iDest);\n", " double bw=binning->GetBinSize(iDest);\n", " /* if(bw!=h[2]->GetBinWidth(iDest)) {\n", " cout<<\"bin \"<GetBinWidth(iDest)\n", " <<\"\\n\";\n", " } */\n", " if(bw>0.0) {\n", " h[2]->SetBinContent(iDest,c/bw);\n", " h[2]->SetBinError(iDest,e/bw);\n", " } else {\n", " }\n", " }\n", "}" ] }, { "cell_type": "markdown", "id": "81cbac41", "metadata": {}, "source": [ " Definition of a helper function: " ] }, { "cell_type": "code", "execution_count": null, "id": "42e4d1ad", "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%cpp -d\n", "void CreateHistogramCopies(TH2 *h[3],TUnfoldBinning const *binningX) {\n", " h[1]=nullptr;\n", " h[2]=nullptr;\n", "}" ] }, { "cell_type": "markdown", "id": "b4ed9874", "metadata": {}, "source": [ " Definition of a helper function: " ] }, { "cell_type": "code", "execution_count": null, "id": "5f397570", "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%cpp -d\n", "TH2 *AddOverflowXY(TH2 *h,double widthX,double widthY) {\n", " // add overflow bin to X-axis\n", " int nx=h->GetNbinsX();\n", " int ny=h->GetNbinsY();\n", " double *xBins=new double[nx+2];\n", " double *yBins=new double[ny+2];\n", " for(int i=1;i<=nx;i++) {\n", " xBins[i-1]=h->GetXaxis()->GetBinLowEdge(i);\n", " }\n", " xBins[nx]=h->GetXaxis()->GetBinUpEdge(nx);\n", " xBins[nx+1]=xBins[nx]+widthX;\n", " for(int i=1;i<=ny;i++) {\n", " yBins[i-1]=h->GetYaxis()->GetBinLowEdge(i);\n", " }\n", " yBins[ny]=h->GetYaxis()->GetBinUpEdge(ny);\n", " yBins[ny+1]=yBins[ny]+widthY;\n", " TString name(h->GetName());\n", " name+=\"U\";\n", " TH2 *r=new TH2D(name,h->GetTitle(),nx+1,xBins,ny+1,yBins);\n", " for(int ix=0;ix<=nx+1;ix++) {\n", " for(int iy=0;iy<=ny+1;iy++) {\n", " r->SetBinContent(ix,iy,h->GetBinContent(ix,iy));\n", " r->SetBinError(ix,iy,h->GetBinError(ix,iy));\n", " }\n", " }\n", " delete [] yBins;\n", " delete [] xBins;\n", " return r;\n", "}" ] }, { "cell_type": "markdown", "id": "37e56426", "metadata": {}, "source": [ " Definition of a helper function: " ] }, { "cell_type": "code", "execution_count": null, "id": "31583209", "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%cpp -d\n", "TH1 *AddOverflowX(TH1 *h,double widthX) {\n", " // add overflow bin to X-axis\n", " int nx=h->GetNbinsX();\n", " double *xBins=new double[nx+2];\n", " for(int i=1;i<=nx;i++) {\n", " xBins[i-1]=h->GetXaxis()->GetBinLowEdge(i);\n", " }\n", " xBins[nx]=h->GetXaxis()->GetBinUpEdge(nx);\n", " xBins[nx+1]=xBins[nx]+widthX;\n", " TString name(h->GetName());\n", " name+=\"U\";\n", " TH1 *r=new TH1D(name,h->GetTitle(),nx+1,xBins);\n", " for(int ix=0;ix<=nx+1;ix++) {\n", " r->SetBinContent(ix,h->GetBinContent(ix));\n", " r->SetBinError(ix,h->GetBinError(ix));\n", " }\n", " delete [] xBins;\n", " return r;\n", "}" ] }, { "cell_type": "markdown", "id": "10c8567a", "metadata": {}, "source": [ " Definition of a helper function: " ] }, { "cell_type": "code", "execution_count": null, "id": "ecd0704d", "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%cpp -d\n", "void DrawOverflowX(TH1 *h,double posy) {\n", " double x1=h->GetXaxis()->GetBinLowEdge(h->GetNbinsX());\n", " double x2=h->GetXaxis()->GetBinUpEdge(h->GetNbinsX());\n", " double y0=h->GetYaxis()->GetBinLowEdge(1);\n", " double y2=h->GetYaxis()->GetBinUpEdge(h->GetNbinsY());;\n", " if(h->GetDimension()==1) {\n", " y0=h->GetMinimum();\n", " y2=h->GetMaximum();\n", " }\n", " double w1=-0.3;\n", " TText *textX=new TText((1.+w1)*x2-w1*x1,(1.-posy)*y0+posy*y2,\"Overflow bin\");\n", " textX->SetNDC(kFALSE);\n", " textX->SetTextSize(0.05);\n", " textX->SetTextAngle(90.);\n", " textX->Draw();\n", " TLine *lineX=new TLine(x1,y0,x1,y2);\n", " lineX->Draw();\n", "}" ] }, { "cell_type": "markdown", "id": "f7f8b8d4", "metadata": {}, "source": [ " Definition of a helper function: " ] }, { "cell_type": "code", "execution_count": null, "id": "c3829f44", "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%cpp -d\n", "void DrawOverflowY(TH1 *h,double posx) {\n", " double x0=h->GetXaxis()->GetBinLowEdge(1);\n", " double x2=h->GetXaxis()->GetBinUpEdge(h->GetNbinsX());\n", " double y1=h->GetYaxis()->GetBinLowEdge(h->GetNbinsY());;\n", " double y2=h->GetYaxis()->GetBinUpEdge(h->GetNbinsY());;\n", " double w1=-0.3;\n", " TText *textY=new TText((1.-posx)*x0+posx*x2,(1.+w1)*y1-w1*y2,\"Overflow bin\");\n", " textY->SetNDC(kFALSE);\n", " textY->SetTextSize(0.05);\n", " textY->Draw();\n", " TLine *lineY=new TLine(x0,y1,x2,y1);\n", " lineY->Draw();\n", "}" ] }, { "cell_type": "markdown", "id": "8969ad69", "metadata": {}, "source": [ " Definition of a helper function: " ] }, { "cell_type": "code", "execution_count": null, "id": "e4f0aa1f", "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%cpp -d\n", "void DrawPadProbability(TH2 *h) {\n", " h->Draw(\"COLZ\");\n", " h->SetTitle(\"migration probabilities;P_{T}(gen) [GeV];P_{T}(rec) [GeV]\");\n", " DrawOverflowX(h,0.05);\n", " DrawOverflowY(h,0.35);\n", "}" ] }, { "cell_type": "markdown", "id": "1e19587e", "metadata": {}, "source": [ " Definition of a helper function: " ] }, { "cell_type": "code", "execution_count": null, "id": "bacda729", "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%cpp -d\n", "void DrawPadEfficiency(TH1 *h) {\n", " h->SetTitle(\"efficiency;P_{T}(gen) [GeV];#epsilon\");\n", " h->SetLineColor(kBlue);\n", " h->SetMinimum(0.75);\n", " h->SetMaximum(1.0);\n", " h->Draw();\n", " DrawOverflowX(h,0.05);\n", " TLegend *legEfficiency=new TLegend(0.3,0.58,0.6,0.75);\n", " legEfficiency->SetBorderSize(0);\n", " legEfficiency->SetFillStyle(0);\n", " legEfficiency->SetTextSize(kLegendFontSize);\n", " legEfficiency->AddEntry(h,\"reconstruction\",\"l\");\n", " legEfficiency->AddEntry((TObject*)nullptr,\" efficiency\",\"\");\n", " legEfficiency->Draw();\n", "}" ] }, { "cell_type": "markdown", "id": "b059cc56", "metadata": {}, "source": [ " Definition of a helper function: " ] }, { "cell_type": "code", "execution_count": null, "id": "bb87917e", "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%cpp -d\n", "void DrawPadReco(TH1 *histMcRec,TH1 *histMcbgrRec,TH1 *histDataRec,\n", " TH1 *histDataUnfold,TH2 *histProbability,TH2 *histRhoij) {\n", " //gPad->SetLogy(kTRUE);\n", " double amax=0.0;\n", " for(int i=1;i<=histMcRec->GetNbinsX();i++) {\n", " amax=TMath::Max(amax,histMcRec->GetBinContent(i)\n", " +5.0*histMcRec->GetBinError(i));\n", " amax=TMath::Max(amax,histDataRec->GetBinContent(i)\n", " +2.0*histDataRec->GetBinError(i));\n", " }\n", " histMcRec->SetTitle(\"Reconstructed;P_{T}(rec);Nevent / GeV\");\n", " histMcRec->Draw(\"HIST\");\n", " histMcRec->SetLineColor(kBlue);\n", " histMcRec->SetMinimum(1.0);\n", " histMcRec->SetMaximum(amax);\n", " //histMcbgrRec->SetFillMode(1);\n", " histMcbgrRec->SetLineColor(kBlue-6);\n", " histMcbgrRec->SetFillColor(kBlue-10);\n", " histMcbgrRec->Draw(\"SAME HIST\");\n", "\n", " TH1 * histFoldBack=nullptr;\n", " if(histDataUnfold && histProbability && histRhoij) {\n", " histFoldBack=(TH1 *)\n", " histMcRec->Clone(histDataUnfold->GetName()+TString(\"_folded\"));\n", " int nrec=histFoldBack->GetNbinsX();\n", " if((nrec==histProbability->GetNbinsY())&&\n", " (nrec==histMcbgrRec->GetNbinsX())&&\n", " (nrec==histDataRec->GetNbinsX())\n", " ) {\n", " for(int ix=1;ix<=nrec;ix++) {\n", " double sum=0.0;\n", " double sume2=0.0;\n", " for(int iy=0;iy<=histProbability->GetNbinsX()+1;iy++) {\n", " sum += histDataUnfold->GetBinContent(iy)*\n", " histDataUnfold->GetBinWidth(iy)*\n", " histProbability->GetBinContent(iy,ix);\n", " for(int iy2=0;iy2<=histProbability->GetNbinsX()+1;iy2++) {\n", " sume2 += histDataUnfold->GetBinError(iy)*\n", " histDataUnfold->GetBinWidth(iy)*\n", " histProbability->GetBinContent(iy,ix)*\n", " histDataUnfold->GetBinError(iy2)*\n", " histDataUnfold->GetBinWidth(iy2)*\n", " histProbability->GetBinContent(iy2,ix)*\n", " histRhoij->GetBinContent(iy,iy2);\n", " }\n", " }\n", " sum /= histFoldBack->GetBinWidth(ix);\n", " sum += histMcbgrRec->GetBinContent(ix);\n", " histFoldBack->SetBinContent(ix,sum);\n", " histFoldBack->SetBinError(ix,TMath::Sqrt(sume2)\n", " /histFoldBack->GetBinWidth(ix));\n", " }\n", " } else {\n", " cout<<\"can not fold back: \"<GetNbinsY()\n", " <<\" \"<GetNbinsX()\n", " <<\" \"<GetNbinsX()\n", " <<\"\\n\";\n", " exit(0);\n", " }\n", "\n", " histFoldBack->SetLineColor(kBlack);\n", " histFoldBack->SetMarkerStyle(0);\n", " histFoldBack->Draw(\"SAME HIST\");\n", " }\n", "\n", " histDataRec->SetLineColor(kRed);\n", " histDataRec->SetMarkerColor(kRed);\n", " histDataRec->Draw(\"SAME\");\n", " DrawOverflowX(histMcRec,0.5);\n", "\n", " TLegend *legRec=new TLegend(0.4,0.5,0.68,0.85);\n", " legRec->SetBorderSize(0);\n", " legRec->SetFillStyle(0);\n", " legRec->SetTextSize(kLegendFontSize);\n", " legRec->AddEntry(histMcRec,\"MC total\",\"l\");\n", " legRec->AddEntry(histMcbgrRec,\"background\",\"f\");\n", " if(histFoldBack) {\n", " int ndf=-kNbinC;\n", " double sumD=0.,sumF=0.,chi2=0.;\n", " for(int i=1;i<=histDataRec->GetNbinsX();i++) {\n", " //cout<GetBinContent(i)<<\" \"<GetBinContent(i)<<\" \"<<\" w=\"<GetBinWidth(i)<<\"\\n\";\n", " sumD+=histDataRec->GetBinContent(i)*histDataRec->GetBinWidth(i);\n", " sumF+=histFoldBack->GetBinContent(i)*histFoldBack->GetBinWidth(i);\n", " double pull=(histFoldBack->GetBinContent(i)-histDataRec->GetBinContent(i))/histDataRec->GetBinError(i);\n", " chi2+= pull*pull;\n", " ndf+=1;\n", " }\n", " legRec->AddEntry(histDataRec,TString::Format(\"data N_{evt}=%.0f\",sumD),\"lp\");\n", " legRec->AddEntry(histFoldBack,TString::Format(\"folded N_{evt}=%.0f\",sumF),\"l\");\n", " legRec->AddEntry((TObject*)nullptr,TString::Format(\"#chi^{2}=%.1f ndf=%d\",chi2,ndf),\"\");\n", " //exit(0);\n", " } else {\n", " legRec->AddEntry(histDataRec,\"data\",\"lp\");\n", " }\n", " legRec->Draw();\n", "}" ] }, { "cell_type": "markdown", "id": "55716f01", "metadata": {}, "source": [ " Definition of a helper function: " ] }, { "cell_type": "code", "execution_count": null, "id": "faba8c9f", "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%cpp -d\n", "void DrawPadTruth(TH1 *histMcsigGen,TH1 *histDataGen,TH1 *histDataUnfold,\n", " char const *text,double tau,vector const *r,\n", " TF1 *f) {\n", " //gPad->SetLogy(kTRUE);\n", " double amin=0.;\n", " double amax=0.;\n", " for(int i=1;i<=histMcsigGen->GetNbinsX();i++) {\n", " if(histDataUnfold) {\n", " amin=TMath::Min(amin,histDataUnfold->GetBinContent(i)\n", " -1.1*histDataUnfold->GetBinError(i));\n", " amax=TMath::Max(amax,histDataUnfold->GetBinContent(i)\n", " +1.1*histDataUnfold->GetBinError(i));\n", " }\n", " amin=TMath::Min(amin,histMcsigGen->GetBinContent(i)\n", " -histMcsigGen->GetBinError(i));\n", " amin=TMath::Min(amin,histDataGen->GetBinContent(i)\n", " -histDataGen->GetBinError(i));\n", " amax=TMath::Max(amax,histMcsigGen->GetBinContent(i)\n", " +10.*histMcsigGen->GetBinError(i));\n", " amax=TMath::Max(amax,histDataGen->GetBinContent(i)\n", " +2.*histDataGen->GetBinError(i));\n", " }\n", " histMcsigGen->SetMinimum(amin);\n", " histMcsigGen->SetMaximum(amax);\n", "\n", " histMcsigGen->SetTitle(\"Truth;P_{T};Nevent / GeV\");\n", " histMcsigGen->SetLineColor(kBlue);\n", " histMcsigGen->Draw(\"HIST\");\n", " histDataGen->SetLineColor(kRed);\n", " histDataGen->SetMarkerColor(kRed);\n", " histDataGen->SetMarkerSize(1.0);\n", " histDataGen->Draw(\"SAME HIST\");\n", " if(histDataUnfold) {\n", " histDataUnfold->SetMarkerStyle(21);\n", " histDataUnfold->SetMarkerSize(0.7);\n", " histDataUnfold->Draw(\"SAME\");\n", " }\n", " DrawOverflowX(histMcsigGen,0.5);\n", "\n", " if(f) {\n", " f->SetLineStyle(kSolid);\n", " f->Draw(\"SAME\");\n", " }\n", "\n", " TLegend *legTruth=new TLegend(0.32,0.65,0.6,0.9);\n", " legTruth->SetBorderSize(0);\n", " legTruth->SetFillStyle(0);\n", " legTruth->SetTextSize(kLegendFontSize);\n", " legTruth->AddEntry(histMcsigGen,\"MC\",\"l\");\n", " if(!histDataUnfold) legTruth->AddEntry((TObject *)nullptr,\" Landau(5,2)\",\"\");\n", " legTruth->AddEntry(histDataGen,\"data\",\"l\");\n", " if(!histDataUnfold) legTruth->AddEntry((TObject *)nullptr,\" Landau(6,1.8)\",\"\");\n", " if(histDataUnfold) {\n", " TString t;\n", " if(text) t=text;\n", " else t=histDataUnfold->GetName();\n", " if(tau>0) {\n", " t+=TString::Format(\" #tau=%.2g\",tau);\n", " }\n", " legTruth->AddEntry(histDataUnfold,t,\"lp\");\n", " if(r) {\n", " legTruth->AddEntry((TObject *)nullptr,\"test wrt data:\",\"\");\n", " legTruth->AddEntry((TObject *)nullptr,TString::Format\n", " (\"#chi^{2}/%d=%.1f prob=%.3f\",\n", " (int)(*r)[3],(*r)[2]/(*r)[3],\n", " TMath::Prob((*r)[2],(*r)[3])),\"\");\n", " }\n", " }\n", " if(f) {\n", " legTruth->AddEntry(f,\"fit\",\"l\");\n", " }\n", " legTruth->Draw();\n", " if(histDataUnfold ) {\n", " TPad *subpad = new TPad(\"subpad\",\"\",0.35,0.29,0.88,0.68);\n", " subpad->SetFillStyle(0);\n", " subpad->Draw();\n", " subpad->cd();\n", " amin=0.;\n", " amax=0.;\n", " int istart=11;\n", " for(int i=istart;i<=histMcsigGen->GetNbinsX();i++) {\n", " amin=TMath::Min(amin,histMcsigGen->GetBinContent(i)\n", " -histMcsigGen->GetBinError(i));\n", " amin=TMath::Min(amin,histDataGen->GetBinContent(i)\n", " -histDataGen->GetBinError(i));\n", " amin=TMath::Min(amin,histDataUnfold->GetBinContent(i)\n", " -histDataUnfold->GetBinError(i));\n", " amax=TMath::Max(amax,histMcsigGen->GetBinContent(i)\n", " +histMcsigGen->GetBinError(i));\n", " amax=TMath::Max(amax,histDataGen->GetBinContent(i)\n", " +histDataGen->GetBinError(i));\n", " amax=TMath::Max(amax,histDataUnfold->GetBinContent(i)\n", " +histDataUnfold->GetBinError(i));\n", " }\n", " TH1 *copyMcsigGen=(TH1*)histMcsigGen->Clone();\n", " TH1 *copyDataGen=(TH1*)histDataGen->Clone();\n", " TH1 *copyDataUnfold=(TH1*)histDataUnfold->Clone();\n", " copyMcsigGen->GetXaxis()->SetRangeUser\n", " (copyMcsigGen->GetXaxis()->GetBinLowEdge(istart),\n", " copyMcsigGen->GetXaxis()->GetBinUpEdge(copyMcsigGen->GetNbinsX()-1));\n", " copyMcsigGen->SetTitle(\";;\");\n", " copyMcsigGen->GetYaxis()->SetRangeUser(amin,amax);\n", " copyMcsigGen->Draw(\"HIST\");\n", " copyDataGen->Draw(\"SAME HIST\");\n", " copyDataUnfold->Draw(\"SAME\");\n", " if(f) {\n", " ((TF1 *)f->Clone())->Draw(\"SAME\");\n", " }\n", " }\n", "}" ] }, { "cell_type": "markdown", "id": "3501ca94", "metadata": {}, "source": [ " Definition of a helper function: " ] }, { "cell_type": "code", "execution_count": null, "id": "0b1d0300", "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%cpp -d\n", "void DrawPadCorrelations(TH2 *h,\n", " vector > > const *table) {\n", " h->SetMinimum(-1.);\n", " h->SetMaximum(1.);\n", " h->SetTitle(\"correlation coefficients;P_{T}(gen) [GeV];P_{T}(gen) [GeV]\");\n", " h->Draw(\"COLZ\");\n", " DrawOverflowX(h,0.05);\n", " DrawOverflowY(h,0.05);\n", " if(table) {\n", " TLegend *legGCor=new TLegend(0.13,0.6,0.5,0.8);\n", " legGCor->SetBorderSize(0);\n", " legGCor->SetFillStyle(0);\n", " legGCor->SetTextSize(kLegendFontSize);\n", " vector const &r=(*table)[table->size()-1].second;\n", " legGCor->AddEntry((TObject *)nullptr,TString::Format(\"min(#rho_{ij})=%5.2f\",r[6]),\"\");\n", " legGCor->AddEntry((TObject *)nullptr,TString::Format(\"max(#rho_{ij})=%5.2f\",r[7]),\"\");\n", " legGCor->AddEntry((TObject *)nullptr,TString::Format(\"avg(#rho_i)=%5.2f\",r[5]),\"\");\n", " legGCor->Draw();\n", " }\n", "}" ] }, { "cell_type": "markdown", "id": "d590aad9", "metadata": {}, "source": [ " Definition of a helper function: " ] }, { "cell_type": "code", "execution_count": null, "id": "b88f84be", "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%cpp -d\n", "void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag) {\n", " if(flag==0) {\n", " cout<<\"fcn flag=0: npar=\"<GetNrows();\n", " TVectorD dy(n);\n", " double x0=0,y0=0.;\n", " for(int i=0;i<=n;i++) {\n", " double x1;\n", " if(i<1) x1=g_fcnHist->GetXaxis()->GetBinLowEdge(i+1);\n", " else x1=g_fcnHist->GetXaxis()->GetBinUpEdge(i);\n", " double y1=TMath::LandauI((x1-u[1])/u[2]);\n", " if(i>0) {\n", " double iy=u[0]*u[2]*(y1-y0)/(x1-x0);\n", " dy(i-1)=iy-g_fcnHist->GetBinContent(i);\n", " //cout<<\"i=\"< > > &table,int niter) {\n", " TString option=\"IESN\";\n", " cout<GetName()<<\"\\n\";\n", " double gcorAvg=0.;\n", " double rhoMin=0.;\n", " double rhoMax=0.;\n", " if(rho) {\n", " g_fcnHist=h;\n", " int n=h->GetNbinsX()-1; // overflow is included as extra bin, exclude in fit\n", " TMatrixDSym v(n);\n", " //g_fcnMatrix=new TMatrixD(n,n);\n", " for(int i=0;iGetBinContent(i+1,j+1)*\n", " (h->GetBinError(i+1)*h->GetBinError(j+1));\n", " }\n", " }\n", " TMatrixDSymEigen ev(v);\n", " TMatrixD d(n,n);\n", " TVectorD di(ev.GetEigenValues());\n", " for(int i=0;i0.0) {\n", " d(i,i)=1./di(i);\n", " } else {\n", " cout<<\"bad eigenvalue i=\"<1.E-7) {\n", " error++;\n", " }\n", " for(int j=0;j1.E-7)) error++;\n", " }\n", " }\n", " // calculate global correlation coefficient (all bins)\n", " TMatrixDSym v1(n+1);\n", " rhoMin=1.;\n", " rhoMax=-1.;\n", " for(int i=0;i<=n;i++) {\n", " for(int j=0;j<=n;j++) {\n", " double rho_ij=rho->GetBinContent(i+1,j+1);\n", " v1(i,j)=rho_ij*\n", " (h->GetBinError(i+1)*h->GetBinError(j+1));\n", " if(i!=j) {\n", " if(rho_ijrhoMax) rhoMax=rho_ij;\n", " }\n", " }\n", " }\n", " TMatrixDSymEigen ev1(v1);\n", " TMatrixD d1(n+1,n+1);\n", " TVectorD di1(ev1.GetEigenValues());\n", " for(int i=0;i<=n;i++) {\n", " if(di1(i)>0.0) {\n", " d1(i,i)=1./di1(i);\n", " } else {\n", " cout<<\"bad eigenvalue i=\"<=0.0) {\n", " double gcor=TMath::Sqrt(gcor2);\n", " gcorAvg += gcor;\n", " } else {\n", " cout<<\"bad global correlation \"<Print();\n", " exit(0);\n", " } */\n", " //g_fcnMatrix->Invert();\n", " //from: HFitImpl.cxx\n", " // TVirtualFitter::FCNFunc_t userFcn = 0;\n", " // typedef void (* FCNFunc_t )(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);\n", " // userFcn = (TVirtualFitter::GetFitter())->GetFCN();\n", " // (TVirtualFitter::GetFitter())->SetUserFunc(f1);\n", " //...\n", " //fitok = fitter->FitFCN( userFcn );\n", "\n", " TVirtualFitter::Fitter(h)->SetFCN(fcn);\n", " option += \"U\";\n", "\n", " }\n", " double xmax=h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()-1);\n", " TF1 *landau=new TF1(text,\"[0]*TMath::Landau(x,[1],[2],0)\",\n", " 0.,xmax);\n", " landau->SetParameter(0,6000.);\n", " landau->SetParameter(1,5.);\n", " landau->SetParameter(2,2.);\n", " landau->SetParError(0,10.);\n", " landau->SetParError(1,0.5);\n", " landau->SetParError(2,0.1);\n", " TFitResultPtr s=h->Fit(landau,option,nullptr,0.,xmax);\n", " vector r(8);\n", " int np=landau->GetNpar();\n", " fcn(np,nullptr,r[0],landau->GetParameters(),0);\n", " r[1]=h->GetNbinsX()-1-landau->GetNpar();\n", " for(int i=0;iGetNbinsX()-1;i++) {\n", " double di=h->GetBinContent(i+1)-truth->GetBinContent(i+1);\n", " if(g_fcnMatrix) {\n", " for(int j=0;jGetNbinsX()-1;j++) {\n", " double dj=h->GetBinContent(j+1)-truth->GetBinContent(j+1);\n", " r[2]+=di*dj*(*g_fcnMatrix)(i,j);\n", " }\n", " } else {\n", " double pull=di/h->GetBinError(i+1);\n", " r[2]+=pull*pull;\n", " }\n", " r[3]+=1.0;\n", " }\n", " r[4]=niter;\n", " if(!niter) r[4]=0.25;\n", " r[5]=gcorAvg;\n", " r[6]=rhoMin;\n", " r[7]=rhoMax;\n", " if(rho) {\n", " g_fcnHist=nullptr;\n", " delete g_fcnMatrix;\n", " g_fcnMatrix=nullptr;\n", " }\n", " table.push_back(make_pair(landau,r));\n", " return s;\n", "}" ] }, { "cell_type": "markdown", "id": "ad00a142", "metadata": {}, "source": [ " ===================== interface to IDS unfolding code follows here\n", "contact Bogdan Malescu to find it\n", " " ] }, { "cell_type": "code", "execution_count": null, "id": "85b14ff7", "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%cpp -d\n", "#include \"ids_code.cc\"\n", "\n", "void IDSfirst(TVectorD *data, TVectorD *dataErr, TMatrixD *A_, Double_t lambdaL_, TVectorD* &unfres1IDS_,TVectorD *&soustr){\n", "\n", " int N_=data->GetNrows();\n", " soustr = new TVectorD(N_);\n", " for( Int_t i=0; iGetNrows();\n", " ModifyMatrix( Am_, A_, unfres2IDS_, dataErr, N_, lambdaM_, soustr, lambdaS_ );\n", " delete unfres2IDS_;\n", " unfres2IDS_ = Unfold( data, dataErr, Am_, N_, lambdaU_, soustr );\n", "}" ] }, { "cell_type": "code", "execution_count": null, "id": "d5faec10", "metadata": { "collapsed": false }, "outputs": [], "source": [ " gErrorIgnoreLevel=kInfo;" ] }, { "cell_type": "markdown", "id": "039420b0", "metadata": {}, "source": [ "switch on histogram errors" ] }, { "cell_type": "code", "execution_count": null, "id": "b335826c", "metadata": { "collapsed": false }, "outputs": [], "source": [ "TH1::SetDefaultSumw2();\n", "\n", "gStyle->SetOptStat(0);\n", "\n", "g_rnd=new TRandom3(4711);" ] }, { "cell_type": "markdown", "id": "e11e4b13", "metadata": {}, "source": [ "==============================================\n", "step 1 : open output file" ] }, { "cell_type": "code", "execution_count": null, "id": "f111e951", "metadata": { "collapsed": false }, "outputs": [], "source": [ "TFile *outputFile=new TFile(\"testUnfold7_results.root\",\"recreate\");" ] }, { "cell_type": "markdown", "id": "85da5a50", "metadata": {}, "source": [ "==============================================\n", "step 2 : read binning schemes and input histograms" ] }, { "cell_type": "code", "execution_count": null, "id": "477639ec", "metadata": { "collapsed": false }, "outputs": [], "source": [ "TFile *inputFile=new TFile(\"testUnfold7_histograms.root\");\n", "\n", "outputFile->cd();\n", "\n", "TUnfoldBinning *fineBinning,*coarseBinning;\n", "\n", "inputFile->GetObject(\"fine\",fineBinning);\n", "inputFile->GetObject(\"coarse\",coarseBinning);\n", "\n", "if((!fineBinning)||(!coarseBinning)) {\n", " cout<<\"problem to read binning schemes\\n\";\n", "}" ] }, { "cell_type": "markdown", "id": "91bb06fc", "metadata": {}, "source": [ "save binning schemes to output file" ] }, { "cell_type": "code", "execution_count": null, "id": "9f4332ef", "metadata": { "collapsed": false }, "outputs": [], "source": [ "fineBinning->Write();\n", "coarseBinning->Write();" ] }, { "cell_type": "markdown", "id": "e0c21503", "metadata": {}, "source": [ "read histograms" ] }, { "cell_type": "code", "execution_count": null, "id": "62720700", "metadata": { "collapsed": false }, "outputs": [], "source": [ "#define READ(TYPE,binning,name) \\\n", "TYPE *name[3]; inputFile->GetObject(#name,name[0]); \\\n", "name[0]->Write(); \\\n", "if(!name[0]) cout<<\"Error reading \" #name \"\\n\"; \\\n", "CreateHistogramCopies(name,binning);\n", "\n", "outputFile->cd();\n", "\n", "READ(TH1,fineBinning,histDataRecF);\n", "READ(TH1,coarseBinning,histDataRecC);\n", "READ(TH1,fineBinning,histDataBgrF);\n", "READ(TH1,coarseBinning,histDataBgrC);\n", "READ(TH1,coarseBinning,histDataGen);\n", "\n", "READ(TH2,fineBinning,histMcsigGenRecF);\n", "READ(TH2,coarseBinning,histMcsigGenRecC);\n", "READ(TH1,fineBinning,histMcsigRecF);\n", "READ(TH1,coarseBinning,histMcsigRecC);\n", "READ(TH1,coarseBinning,histMcsigGen);\n", "\n", "READ(TH1,fineBinning,histMcbgrRecF);\n", "READ(TH1,coarseBinning,histMcbgrRecC);\n", "\n", "TH1 *histOutputCtau0[3];\n", "TH2 *histRhoCtau0;\n", "TH1 *histOutputCLCurve[3];" ] }, { "cell_type": "markdown", "id": "f9543a66", "metadata": {}, "source": [ "TH2 *histRhoCLCurve;" ] }, { "cell_type": "code", "execution_count": null, "id": "37bfbf1f", "metadata": { "collapsed": false }, "outputs": [], "source": [ "TH2 *histProbC;\n", "double tauMin=1.e-4;\n", "double tauMax=1.e-1;\n", "double fBgr=1.0; // 0.2/0.25;\n", "double biasScale=1.0;\n", "TUnfold::ERegMode mode=TUnfold::kRegModeSize; //Derivative;" ] }, { "cell_type": "markdown", "id": "a0433d9b", "metadata": {}, "source": [ "double tauC;" ] }, { "cell_type": "code", "execution_count": null, "id": "fb7a72dd", "metadata": { "collapsed": false }, "outputs": [], "source": [ "{\n", "TUnfoldDensity *tunfoldC=\n", " new TUnfoldDensity(histMcsigGenRecC[0],\n", " TUnfold::kHistMapOutputHoriz,\n", " mode,\n", " TUnfold::kEConstraintNone,//Area,\n", " TUnfoldDensity::kDensityModeNone,\n", " coarseBinning,\n", " coarseBinning);\n", "tunfoldC->SetInput(histDataRecC[0],biasScale);\n", "tunfoldC->SubtractBackground(histMcbgrRecC[0],\"BGR\",fBgr,0.0);\n", "tunfoldC->DoUnfold(0.);\n", "histOutputCtau0[0]=tunfoldC->GetOutput(\"histOutputCtau0\");\n", "histRhoCtau0=tunfoldC->GetRhoIJtotal(\"histRhoCtau0\");\n", "CreateHistogramCopies(histOutputCtau0,coarseBinning);\n", "tunfoldC->ScanLcurve(50,tauMin,tauMax,nullptr);\n", "/* tauC= */tunfoldC->GetTau();\n", "//tunfoldC->ScanTau(50,1.E-7,1.E-1,0,TUnfoldDensity::kEScanTauRhoAvg);\n", "histOutputCLCurve[0]=tunfoldC->GetOutput(\"histOutputCLCurve\");\n", "/* histRhoCLCurve= */tunfoldC->GetRhoIJtotal(\"histRhoCLCurve\");\n", "CreateHistogramCopies(histOutputCLCurve,coarseBinning);\n", "histProbC=tunfoldC->GetProbabilityMatrix(\"histProbC\",\";P_T(gen);P_T(rec)\");\n", "}\n", "TH1 *histOutputFtau0[3];\n", "TH2 *histRhoFtau0;\n", "TH1 *histOutputFLCurve[3];" ] }, { "cell_type": "markdown", "id": "0886ce50", "metadata": {}, "source": [ "TH2 *histRhoFLCurve;" ] }, { "cell_type": "code", "execution_count": null, "id": "aa207c1d", "metadata": { "collapsed": false }, "outputs": [], "source": [ "TH2 *histProbF;\n", "TGraph *lCurve;\n", "TSpline *logTauX,*logTauY;\n", "tauMin=3.E-4;\n", "tauMax=3.E-2;" ] }, { "cell_type": "markdown", "id": "eba04c6b", "metadata": {}, "source": [ "double tauF;" ] }, { "cell_type": "code", "execution_count": null, "id": "8f31e0b5", "metadata": { "collapsed": false }, "outputs": [], "source": [ "{\n", "TUnfoldDensity *tunfoldF=\n", " new TUnfoldDensity(histMcsigGenRecF[0],\n", " TUnfold::kHistMapOutputHoriz,\n", " mode,\n", " TUnfold::kEConstraintNone,//Area,\n", " TUnfoldDensity::kDensityModeNone,\n", " coarseBinning,\n", " fineBinning);\n", "tunfoldF->SetInput(histDataRecF[0],biasScale);\n", "tunfoldF->SubtractBackground(histMcbgrRecF[0],\"BGR\",fBgr,0.0);\n", "tunfoldF->DoUnfold(0.);\n", "histOutputFtau0[0]=tunfoldF->GetOutput(\"histOutputFtau0\");\n", "histRhoFtau0=tunfoldF->GetRhoIJtotal(\"histRhoFtau0\");\n", "CreateHistogramCopies(histOutputFtau0,coarseBinning);\n", "tunfoldF->ScanLcurve(50,tauMin,tauMax,nullptr);\n", "//tunfoldF->DoUnfold(tauC);\n", "/* tauF= */tunfoldF->GetTau();\n", "//tunfoldF->ScanTau(50,1.E-7,1.E-1,0,TUnfoldDensity::kEScanTauRhoAvg);\n", "histOutputFLCurve[0]=tunfoldF->GetOutput(\"histOutputFLCurve\");\n", "/* histRhoFLCurve= */tunfoldF->GetRhoIJtotal(\"histRhoFLCurve\");\n", "CreateHistogramCopies(histOutputFLCurve,coarseBinning);\n", "histProbF=tunfoldF->GetProbabilityMatrix(\"histProbF\",\";P_T(gen);P_T(rec)\");\n", "}\n", "TH1 *histOutputFAtau0[3];\n", "TH2 *histRhoFAtau0;\n", "TH1 *histOutputFALCurve[3];\n", "TH2 *histRhoFALCurve;\n", "TH1 *histOutputFArho[3];\n", "TH2 *histRhoFArho;\n", "TSpline *rhoScan=nullptr;\n", "TSpline *logTauCurvature=nullptr;\n", "\n", "double tauFA,tauFArho;\n", "{\n", "TUnfoldDensity *tunfoldFA=\n", " new TUnfoldDensity(histMcsigGenRecF[0],\n", " TUnfold::kHistMapOutputHoriz,\n", " mode,\n", " TUnfold::kEConstraintArea,\n", " TUnfoldDensity::kDensityModeNone,\n", " coarseBinning,\n", " fineBinning);\n", "tunfoldFA->SetInput(histDataRecF[0],biasScale);\n", "tunfoldFA->SubtractBackground(histMcbgrRecF[0],\"BGR\",fBgr,0.0);\n", "tunfoldFA->DoUnfold(0.);\n", "histOutputFAtau0[0]=tunfoldFA->GetOutput(\"histOutputFAtau0\");\n", "histRhoFAtau0=tunfoldFA->GetRhoIJtotal(\"histRhoFAtau0\");\n", "CreateHistogramCopies(histOutputFAtau0,coarseBinning);\n", "tunfoldFA->ScanTau(50,tauMin,tauMax,&rhoScan,TUnfoldDensity::kEScanTauRhoAvg);\n", "tauFArho=tunfoldFA->GetTau();\n", "histOutputFArho[0]=tunfoldFA->GetOutput(\"histOutputFArho\");\n", "histRhoFArho=tunfoldFA->GetRhoIJtotal(\"histRhoFArho\");\n", "CreateHistogramCopies(histOutputFArho,coarseBinning);\n", "\n", "tunfoldFA->ScanLcurve(50,tauMin,tauMax,&lCurve,&logTauX,&logTauY,&logTauCurvature);\n", "tauFA=tunfoldFA->GetTau();\n", "histOutputFALCurve[0]=tunfoldFA->GetOutput(\"histOutputFALCurve\");\n", "histRhoFALCurve=tunfoldFA->GetRhoIJtotal(\"histRhoFALCurve\");\n", "CreateHistogramCopies(histOutputFALCurve,coarseBinning);\n", "}\n", "lCurve->Write();\n", "logTauX->Write();\n", "logTauY->Write();\n", "\n", "\n", "double widthC=coarseBinning->GetBinSize(histProbC->GetNbinsY()+1);\n", "double widthF=fineBinning->GetBinSize(histProbF->GetNbinsY()+1);\n", "\n", "TH2 *histProbCO=AddOverflowXY(histProbC,widthC,widthC);\n", "TH2 *histProbFO=AddOverflowXY(histProbF,widthC,widthF);" ] }, { "cell_type": "markdown", "id": "2efb414e", "metadata": {}, "source": [ "efficiency" ] }, { "cell_type": "code", "execution_count": null, "id": "e8b17ec4", "metadata": { "collapsed": false }, "outputs": [], "source": [ "TH1 *histEfficiencyC=histProbCO->ProjectionX(\"histEfficiencyC\");\n", "kNbinC=histProbCO->GetNbinsX();" ] }, { "cell_type": "markdown", "id": "57bd9cce", "metadata": {}, "source": [ "reconstructed quantities with overflow (coarse binning)\n", "MC: add signal and bgr" ] }, { "cell_type": "code", "execution_count": null, "id": "15ffc955", "metadata": { "collapsed": false }, "outputs": [], "source": [ "TH1 *histMcsigRecCO=AddOverflowX(histMcsigRecC[2],widthC);\n", "TH1 *histMcbgrRecCO=AddOverflowX(histMcbgrRecC[2],widthC);\n", "histMcbgrRecCO->Scale(fBgr);\n", "TH1 *histMcRecCO=(TH1 *)histMcsigRecCO->Clone(\"histMcRecC0\");\n", "histMcRecCO->Add(histMcsigRecCO,histMcbgrRecCO);\n", "TH1 *histDataRecCO=AddOverflowX(histDataRecC[2],widthC);" ] }, { "cell_type": "markdown", "id": "ef91ac84", "metadata": {}, "source": [ "TH1 *histDataRecCNO=AddOverflowX(histDataRecC[1],widthC);" ] }, { "cell_type": "code", "execution_count": null, "id": "da4fb125", "metadata": { "collapsed": false }, "outputs": [], "source": [ "TH1 *histMcsigRecFO=AddOverflowX(histMcsigRecF[2],widthF);\n", "TH1 *histMcbgrRecFO=AddOverflowX(histMcbgrRecF[2],widthF);\n", "histMcbgrRecFO->Scale(fBgr);\n", "TH1 *histMcRecFO=(TH1 *)histMcsigRecFO->Clone(\"histMcRecF0\");\n", "histMcRecFO->Add(histMcsigRecFO,histMcbgrRecFO);\n", "TH1 *histDataRecFO=AddOverflowX(histDataRecF[2],widthF);" ] }, { "cell_type": "markdown", "id": "510f71b7", "metadata": {}, "source": [ "truth level with overflow" ] }, { "cell_type": "code", "execution_count": null, "id": "2b245bf4", "metadata": { "collapsed": false }, "outputs": [], "source": [ "TH1 *histMcsigGenO=AddOverflowX(histMcsigGen[2],widthC);\n", "TH1 *histDataGenO=AddOverflowX(histDataGen[2],widthC);" ] }, { "cell_type": "markdown", "id": "1b262d3a", "metadata": {}, "source": [ "unfolding result with overflow" ] }, { "cell_type": "code", "execution_count": null, "id": "4f367a8f", "metadata": { "collapsed": false }, "outputs": [], "source": [ "TH1 *histOutputCtau0O=AddOverflowX(histOutputCtau0[2],widthC);\n", "TH2 *histRhoCtau0O=AddOverflowXY(histRhoCtau0,widthC,widthC);" ] }, { "cell_type": "markdown", "id": "ae3bddc7", "metadata": {}, "source": [ "TH1 *histOutputCLCurveO=AddOverflowX(histOutputCLCurve[2],widthC);\n", "TH2 *histRhoCLCurveO=AddOverflowXY(histRhoCLCurve,widthC,widthC);" ] }, { "cell_type": "code", "execution_count": null, "id": "81597a65", "metadata": { "collapsed": false }, "outputs": [], "source": [ "TH1 *histOutputFtau0O=AddOverflowX(histOutputFtau0[2],widthC);\n", "TH2 *histRhoFtau0O=AddOverflowXY(histRhoFtau0,widthC,widthC);\n", "TH1 *histOutputFAtau0O=AddOverflowX(histOutputFAtau0[2],widthC);\n", "TH2 *histRhoFAtau0O=AddOverflowXY(histRhoFAtau0,widthC,widthC);" ] }, { "cell_type": "markdown", "id": "1289b3b7", "metadata": {}, "source": [ "TH1 *histOutputFLCurveO=AddOverflowX(histOutputFLCurve[2],widthC);\n", "TH2 *histRhoFLCurveO=AddOverflowXY(histRhoFLCurve,widthC,widthC);" ] }, { "cell_type": "code", "execution_count": null, "id": "ae4f0741", "metadata": { "collapsed": false }, "outputs": [], "source": [ "TH1 *histOutputFALCurveO=AddOverflowX(histOutputFALCurve[2],widthC);\n", "TH2 *histRhoFALCurveO=AddOverflowXY(histRhoFALCurve,widthC,widthC);\n", "TH1 *histOutputFArhoO=AddOverflowX(histOutputFArho[2],widthC);\n", "TH2 *histRhoFArhoO=AddOverflowXY(histRhoFArho,widthC,widthC);" ] }, { "cell_type": "markdown", "id": "3b1ba5e3", "metadata": {}, "source": [ "bin-by-bin" ] }, { "cell_type": "code", "execution_count": null, "id": "7613c131", "metadata": { "collapsed": false }, "outputs": [], "source": [ "TH2 *histRhoBBBO=(TH2 *)histRhoCtau0O->Clone(\"histRhoBBBO\");\n", "for(int i=1;i<=histRhoBBBO->GetNbinsX();i++) {\n", " for(int j=1;j<=histRhoBBBO->GetNbinsX();j++) {\n", " histRhoBBBO->SetBinContent(i,j,(i==j)?1.:0.);\n", " }\n", "}\n", "TH1 *histDataBgrsub=(TH1 *)histDataRecCO->Clone(\"histDataBgrsub\");\n", "histDataBgrsub->Add(histMcbgrRecCO,-fBgr);\n", "TH1 *histOutputBBBO=(TH1 *)histDataBgrsub->Clone(\"histOutputBBBO\");\n", "histOutputBBBO->Divide(histMcsigRecCO);\n", "histOutputBBBO->Multiply(histMcsigGenO);" ] }, { "cell_type": "markdown", "id": "a889507b", "metadata": {}, "source": [ "iterative" ] }, { "cell_type": "code", "execution_count": null, "id": "f08760cc", "metadata": { "collapsed": false }, "outputs": [], "source": [ "int niter=1000;\n", "cout<<\"maximum number of iterations: \"<histOutputAgo,histOutputAgorep;\n", "vector histRhoAgo,histRhoAgorep;\n", "vector nIter;\n", "histOutputAgo.push_back((TH1*)histMcsigGenO->Clone(\"histOutputAgo-1\"));\n", "histOutputAgorep.push_back((TH1*)histMcsigGenO->Clone(\"histOutputAgorep-1\"));\n", "histRhoAgo.push_back((TH2*)histRhoBBBO->Clone(\"histRhoAgo-1\"));\n", "histRhoAgorep.push_back((TH2*)histRhoBBBO->Clone(\"histRhoAgorep-1\"));\n", "nIter.push_back(-1);\n", "\n", "int nx=histProbCO->GetNbinsX();\n", "int ny=histProbCO->GetNbinsY();\n", "TMatrixD covAgo(nx+ny,nx+ny);\n", "TMatrixD A(ny,nx);\n", "TMatrixD AToverEps(nx,ny);\n", "for(int i=0;iGetBinContent(i+1,j+1);\n", " }\n", " for(int j=0;jGetBinContent(i+1,j+1);\n", " A(j,i)=aji;\n", " AToverEps(i,j)=aji/epsilonI;\n", " }\n", "}\n", "for(int i=0;iGetBinError(i+1)\n", " *histOutputAgo[0]->GetXaxis()->GetBinWidth(i+1),2.);\n", "}\n", "for(int i=0;iGetBinError(i+1)\n", " *histDataRecCO->GetXaxis()->GetBinWidth(i+1),2.);\n", "}\n", "#define NREPLICA 300\n", "vector y(NREPLICA);\n", "vector yMb(NREPLICA);\n", "vector yErr(NREPLICA);\n", "vector x(NREPLICA);\n", "TVectorD b(nx);\n", "for(int nr=0;nrGetBinContent(i+1)\n", " *histOutputAgo[0]->GetXaxis()->GetBinWidth(i+1);\n", " for(int nr=1;nrGetBinContent(i+1)\n", " *histDataRecCO->GetXaxis()->GetBinWidth(i+1);\n", " for(int nr=1;nrPoisson((*y[0])(i));\n", " (*yErr[nr])(i)=TMath::Sqrt((*y[nr])(i));\n", " }\n", " b(i)=histMcbgrRecCO->GetBinContent(i+1)*\n", " histMcbgrRecCO->GetXaxis()->GetBinWidth(i+1);\n", " for(int nr=0;nrClone\n", " (TString::Format(\"histOutputAgo%d\",iter));\n", " histOutputAgo.push_back(h);\n", " for(int i=0;iGetXaxis()->GetBinWidth(i+1);\n", " h->SetBinContent(i+1,(*x[0])(i)/bw);\n", " h->SetBinError(i+1,TMath::Sqrt(covAgo(i,i))/bw);\n", " }\n", " TH2 *h2=(TH2*)histRhoAgo[0]->Clone\n", " (TString::Format(\"histRhoAgo%d\",iter));\n", " histRhoAgo.push_back(h2);\n", " for(int i=0;i=1.0)) {\n", " cout<<\"bad error matrix: iter=\"<SetBinContent(i+1,j+1,rho);\n", " }\n", " }\n", " // error and correlations from replica analysis\n", " h=(TH1*)histOutputAgo[0]->Clone\n", " (TString::Format(\"histOutputAgorep%d\",iter));\n", " h2=(TH2*)histRhoAgo[0]->Clone\n", " (TString::Format(\"histRhoAgorep%d\",iter));\n", " histOutputAgorep.push_back(h);\n", " histRhoAgorep.push_back(h2);\n", "\n", " TVectorD mean(nx);\n", " double w=1./(NREPLICA-1.);\n", " for(int nr=1;nrGetXaxis()->GetBinWidth(i+1);\n", " h->SetBinContent(i+1,(*x[0])(i)/bw);\n", " h->SetBinError(i+1,TMath::Sqrt(covAgorep(i,i))/bw);\n", " // cout<=1.0)) {\n", " cout<<\"bad error matrix: iter=\"<SetBinContent(i+1,j+1,rho);\n", " }\n", " }\n", " }\n", "}\n", "\n", "#ifdef WITH_IDS" ] }, { "cell_type": "markdown", "id": "e4f62ce8", "metadata": {}, "source": [ "IDS Malaescu" ] }, { "cell_type": "code", "execution_count": null, "id": "af4bec03", "metadata": { "collapsed": false }, "outputs": [], "source": [ "int niterIDS=100;\n", "vector unfresIDS(NREPLICA),soustr(NREPLICA);\n", "cout<<\"IDS number of iterations: \"<GetBinContent(ix+1,iy+1);\n", " }\n", "}\n", "double lambdaL=0.;\n", "Double_t lambdaUmin = 1.0000002;\n", "Double_t lambdaMmin = 0.0000001;\n", "Double_t lambdaS = 0.000001;\n", "double lambdaU=lambdaUmin;\n", "double lambdaM=lambdaMmin;\n", "vector histOutputIDS;\n", "vector histRhoIDS;\n", "histOutputIDS.push_back((TH1*)histOutputAgo[0]->Clone(\"histOutputIDS-1\"));\n", "histRhoIDS.push_back((TH2*)histRhoAgo[0]->Clone(\"histRhoIDS-1\"));\n", "histOutputIDS.push_back((TH1*)histOutputAgo[0]->Clone(\"histOutputIDS0\"));\n", "histRhoIDS.push_back((TH2*)histRhoAgo[0]->Clone(\"histRhoIDS0\"));\n", "for(int iter=1;iter<=niterIDS;iter++) {\n", " if(!(iter %10)) cout<Clone\n", " (TString::Format(\"histOutputIDS%d\",iter));\n", " TH2 *h2=(TH2*)histRhoIDS[0]->Clone\n", " (TString::Format(\"histRhoIDS%d\",iter));\n", " histOutputIDS.push_back(h);\n", " histRhoIDS.push_back(h2);\n", " TVectorD mean(nx);\n", " double w=1./(NREPLICA-1.);\n", " for(int nr=1;nrGetXaxis()->GetBinWidth(i+1);\n", " h->SetBinContent(i+1,(*unfresIDS[0])(i)/bw/\n", " histEfficiencyC->GetBinContent(i+1));\n", " h->SetBinError(i+1,TMath::Sqrt(covIDSrep(i,i))/bw/\n", " histEfficiencyC->GetBinContent(i+1));\n", " // cout<=1.0)) {\n", " cout<<\"bad error matrix: iter=\"<SetBinContent(i+1,j+1,rho);\n", " }\n", " }\n", " }\n", "}\n", "#endif" ] }, { "cell_type": "markdown", "id": "a50f6ec9", "metadata": {}, "source": [ "double NEdSmc=histDataBgrsub->GetSumOfWeights();" ] }, { "cell_type": "code", "execution_count": null, "id": "b86d7b10", "metadata": { "collapsed": false }, "outputs": [], "source": [ "vector > > table;\n", "\n", "TCanvas *c1=new TCanvas(\"c1\",\"\",600,600);\n", "TCanvas *c2sq=new TCanvas(\"c2sq\",\"\",600,600);\n", "c2sq->Divide(1,2);\n", "TCanvas *c2w=new TCanvas(\"c2w\",\"\",600,300);\n", "c2w->Divide(2,1);\n", "TCanvas *c4=new TCanvas(\"c4\",\"\",600,600);\n", "c4->Divide(2,2);" ] }, { "cell_type": "markdown", "id": "0d636076", "metadata": {}, "source": [ "TCanvas *c3n=new TCanvas(\"c3n\",\"\",600,600);" ] }, { "cell_type": "code", "execution_count": null, "id": "60decd7e", "metadata": { "collapsed": false }, "outputs": [], "source": [ "TPad *subn[3];" ] }, { "cell_type": "markdown", "id": "ed447382", "metadata": {}, "source": [ "gROOT->SetStyle(\"xTimes2\");" ] }, { "cell_type": "code", "execution_count": null, "id": "90889a7b", "metadata": { "collapsed": false }, "outputs": [], "source": [ "subn[0]= new TPad(\"subn0\",\"\",0.,0.5,1.,1.);" ] }, { "cell_type": "markdown", "id": "300b07bc", "metadata": {}, "source": [ "gROOT->SetStyle(\"square\");" ] }, { "cell_type": "code", "execution_count": null, "id": "6bed214b", "metadata": { "collapsed": false }, "outputs": [], "source": [ "subn[1]= new TPad(\"subn1\",\"\",0.,0.,0.5,0.5);\n", "subn[2]= new TPad(\"subn2\",\"\",0.5,0.0,1.,0.5);\n", "for(int i=0;i<3;i++) {\n", " subn[i]->SetFillStyle(0);\n", " subn[i]->Draw();\n", "}\n", "TCanvas *c3c=new TCanvas(\"c3c\",\"\",600,600);\n", "TPad *subc[3];" ] }, { "cell_type": "markdown", "id": "8ea92262", "metadata": {}, "source": [ "gROOT->SetStyle(\"xTimes2\");" ] }, { "cell_type": "code", "execution_count": null, "id": "08dd182c", "metadata": { "collapsed": false }, "outputs": [], "source": [ "subc[0]= new TPad(\"sub0\",\"\",0.,0.5,1.,1.);" ] }, { "cell_type": "markdown", "id": "f359599d", "metadata": {}, "source": [ "gROOT->SetStyle(\"squareCOLZ\");" ] }, { "cell_type": "code", "execution_count": null, "id": "313eefc3", "metadata": { "collapsed": false }, "outputs": [], "source": [ "subc[1]= new TPad(\"sub1\",\"\",0.,0.,0.5,0.5);" ] }, { "cell_type": "markdown", "id": "3fe7ba59", "metadata": {}, "source": [ "gROOT->SetStyle(\"square\");" ] }, { "cell_type": "code", "execution_count": null, "id": "535ab4da", "metadata": { "collapsed": false }, "outputs": [], "source": [ "subc[2]= new TPad(\"sub2\",\"\",0.5,0.0,1.,0.5);\n", "for(int i=0;i<3;i++) {\n", " subc[i]->SetFillStyle(0);\n", " subc[i]->Draw();\n", "}" ] }, { "cell_type": "markdown", "id": "22c8f83e", "metadata": {}, "source": [ "=========================== example ==================================" ] }, { "cell_type": "code", "execution_count": null, "id": "8f2caecd", "metadata": { "collapsed": false }, "outputs": [], "source": [ "c2w->cd(1);\n", "DrawPadTruth(histMcsigGenO,histDataGenO,nullptr);\n", "c2w->cd(2);\n", "DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,nullptr,nullptr,nullptr);\n", "c2w->SaveAs(\"exampleTR.eps\");" ] }, { "cell_type": "markdown", "id": "02b05fb3", "metadata": {}, "source": [ "=========================== example ==================================" ] }, { "cell_type": "code", "execution_count": null, "id": "6e4770f3", "metadata": { "collapsed": false }, "outputs": [], "source": [ "c2w->cd(1);\n", "DrawPadProbability(histProbCO);\n", "c2w->cd(2);\n", "DrawPadEfficiency(histEfficiencyC);\n", "c2w->SaveAs(\"exampleAE.eps\");\n", "\n", "int iFitInversion=table.size();\n", "DoFit(histOutputCtau0O,histRhoCtau0O,histDataGenO,\"inversion\",table);" ] }, { "cell_type": "markdown", "id": "0776d2c5", "metadata": {}, "source": [ "=========================== inversion ==================================" ] }, { "cell_type": "code", "execution_count": null, "id": "d9ed14f3", "metadata": { "collapsed": false }, "outputs": [], "source": [ "subc[0]->cd();\n", "DrawPadTruth(histMcsigGenO,histDataGenO,histOutputCtau0O,\"inversion\",0.,\n", " &table[table.size()-1].second);\n", "subc[1]->cd();\n", "DrawPadCorrelations(histRhoCtau0O,&table);\n", "subc[2]->cd();\n", "DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,\n", " histOutputCtau0O,histProbCO,histRhoCtau0O);\n", "c3c->SaveAs(\"inversion.eps\");\n", "\n", "\n", "DoFit(histOutputFtau0O,histRhoFtau0O,histDataGenO,\"template\",table);" ] }, { "cell_type": "markdown", "id": "a523d0d4", "metadata": {}, "source": [ "=========================== template ==================================" ] }, { "cell_type": "code", "execution_count": null, "id": "4908762f", "metadata": { "collapsed": false }, "outputs": [], "source": [ "subc[0]->cd();\n", "DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFtau0O,\"fit\",0.,\n", " &table[table.size()-1].second);\n", "subc[1]->cd();\n", "DrawPadCorrelations(histRhoFtau0O,&table);\n", "subc[2]->cd();\n", "DrawPadReco(histMcRecFO,histMcbgrRecFO,histDataRecFO,\n", " histOutputFtau0O,histProbFO,histRhoFtau0O);\n", "c3c->SaveAs(\"template.eps\");\n", "\n", "DoFit(histOutputFAtau0O,histRhoFAtau0O,histDataGenO,\"template+area\",table);" ] }, { "cell_type": "markdown", "id": "5dd0368e", "metadata": {}, "source": [ "=========================== template+area ==================================" ] }, { "cell_type": "code", "execution_count": null, "id": "3249c7e7", "metadata": { "collapsed": false }, "outputs": [], "source": [ "subc[0]->cd();\n", "DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFAtau0O,\"fit\",0.,\n", " &table[table.size()-1].second);\n", "subc[1]->cd();\n", "DrawPadCorrelations(histRhoFAtau0O,&table);\n", "subc[2]->cd();\n", "DrawPadReco(histMcRecFO,histMcbgrRecFO,histDataRecFO,\n", " histOutputFAtau0O,histProbFO,histRhoFAtau0O);\n", "c3c->SaveAs(\"templateA.eps\");\n", "\n", "int iFitFALCurve=table.size();\n", "DoFit(histOutputFALCurveO,histRhoFALCurveO,histDataGenO,\"Tikhonov+area\",table);" ] }, { "cell_type": "markdown", "id": "ca0d138c", "metadata": {}, "source": [ "=========================== template+area+tikhonov =====================" ] }, { "cell_type": "code", "execution_count": null, "id": "264b4426", "metadata": { "collapsed": false }, "outputs": [], "source": [ "subc[0]->cd();\n", "DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFALCurveO,\"Tikhonov\",tauFA,\n", " &table[table.size()-1].second);\n", "subc[1]->cd();\n", "DrawPadCorrelations(histRhoFALCurveO,&table);\n", "subc[2]->cd();\n", "DrawPadReco(histMcRecFO,histMcbgrRecFO,histDataRecFO,\n", " histOutputFALCurveO,histProbFO,histRhoFALCurveO);\n", "c3c->SaveAs(\"lcurveFA.eps\");\n", "\n", "int iFitFArho=table.size();\n", "DoFit(histOutputFArhoO,histRhoFArhoO,histDataGenO,\"min(rhomax)\",table);" ] }, { "cell_type": "markdown", "id": "d6cf3591", "metadata": {}, "source": [ "=========================== template+area+tikhonov =====================" ] }, { "cell_type": "code", "execution_count": null, "id": "c5d014b7", "metadata": { "collapsed": false }, "outputs": [], "source": [ "subc[0]->cd();\n", "DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFArhoO,\"Tikhonov\",tauFArho,\n", " &table[table.size()-1].second);\n", "subc[1]->cd();\n", "DrawPadCorrelations(histRhoFArho,&table);\n", "subc[2]->cd();\n", "DrawPadReco(histMcRecFO,histMcbgrRecFO,histDataRecFO,\n", " histOutputFArhoO,histProbFO,histRhoFArhoO);\n", "c3c->SaveAs(\"rhoscanFA.eps\");\n", "\n", "int iFitBinByBin=table.size();\n", "DoFit(histOutputBBBO,histRhoBBBO,histDataGenO,\"bin-by-bin\",table);" ] }, { "cell_type": "markdown", "id": "826a4765", "metadata": {}, "source": [ "=========================== bin-by-bin =================================\n", "c->cd(1);\n", "DrawPadProbability(histProbCO);\n", "c->cd(2);\n", "DrawPadCorrelations(histRhoBBBO,&table);" ] }, { "cell_type": "code", "execution_count": null, "id": "b253c2a4", "metadata": { "collapsed": false }, "outputs": [], "source": [ "c2sq->cd(1);\n", "DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,\n", " histOutputBBBO,histProbCO,histRhoBBBO);\n", "c2sq->cd(2);\n", "DrawPadTruth(histMcsigGenO,histDataGenO,histOutputBBBO,\"bin-by-bin\",0.,\n", " &table[table.size()-1].second);\n", "c2sq->SaveAs(\"binbybin.eps\");" ] }, { "cell_type": "markdown", "id": "7ffa8d68", "metadata": {}, "source": [ "=========================== iterative ===================================" ] }, { "cell_type": "code", "execution_count": null, "id": "d62efafc", "metadata": { "collapsed": false }, "outputs": [], "source": [ "int iAgoFirstFit=table.size();\n", "for(size_t i=1;icd();\n", " DrawPadTruth(histMcsigGenO,histDataGenO,histOutputAgorep[i],\n", " TString::Format(\"iterative N=%d\",nIter[i]),0.,\n", " isFitted ? &table[table.size()-1].second : nullptr);\n", " subc[1]->cd();\n", " DrawPadCorrelations(histRhoAgorep[i],&table);\n", " subc[2]->cd();\n", " DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,\n", " histOutputAgorep[i],histProbCO,histRhoAgorep[i]);\n", " c3c->SaveAs(TString::Format(\"iterative%d.eps\",nIter[i]));\n", "}\n", "int iAgoLastFit=table.size();\n", "\n", "#ifdef WITH_IDS\n", "int iIDSFirstFit=table.size();" ] }, { "cell_type": "markdown", "id": "2d1204ad", "metadata": {}, "source": [ "=========================== IDS ===================================" ] }, { "cell_type": "code", "execution_count": null, "id": "4d5db608", "metadata": { "collapsed": false }, "outputs": [], "source": [ "for(size_t i=2;icd();\n", " DrawPadTruth(histMcsigGenO,histDataGenO,histOutputIDS[i],\n", " TString::Format(\"IDS N=%d\",nIter[i]),0.,\n", " isFitted ? &table[table.size()-1].second : 0);\n", " subc[1]->cd();\n", " DrawPadCorrelations(histRhoIDS[i],&table);\n", " subc[2]->cd();\n", " DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,\n", " histOutputIDS[i],histProbCO,histRhoIDS[i]);\n", " c3c->SaveAs(TString::Format(\"ids%d.eps\",nIter[i]));\n", "}\n", "int iIDSLastFit=table.size();\n", "#endif\n", "\n", "int nfit=table.size();\n", "TH1D *fitChindf=new TH1D(\"fitChindf\",\";algorithm;#chi^{2}/NDF\",nfit,0,nfit);\n", "TH1D *fitNorm=new TH1D(\"fitNorm\",\";algorithm;Landau amplitude [1/GeV]\",nfit,0,nfit);\n", "TH1D *fitMu=new TH1D(\"fitMu\",\";algorithm;Landau #mu [GeV]\",nfit,0,nfit);\n", "TH1D *fitSigma=new TH1D(\"fitSigma\",\";algorithm;Landau #sigma [GeV]\",nfit,0,nfit);\n", "for(int fit=0;fit const &r=table[fit].second;\n", " fitChindf->GetXaxis()->SetBinLabel(fit+1,f->GetName());\n", " fitNorm->GetXaxis()->SetBinLabel(fit+1,f->GetName());\n", " fitMu->GetXaxis()->SetBinLabel(fit+1,f->GetName());\n", " fitSigma->GetXaxis()->SetBinLabel(fit+1,f->GetName());\n", " double chi2=r[0];\n", " double ndf=r[1];\n", " fitChindf->SetBinContent(fit+1,chi2/ndf);\n", " fitChindf->SetBinError(fit+1,TMath::Sqrt(2./ndf));\n", " fitNorm->SetBinContent(fit+1,f->GetParameter(0));\n", " fitNorm->SetBinError(fit+1,f->GetParError(0));\n", " fitMu->SetBinContent(fit+1,f->GetParameter(1));\n", " fitMu->SetBinError(fit+1,f->GetParError(1));\n", " fitSigma->SetBinContent(fit+1,f->GetParameter(2));\n", " fitSigma->SetBinError(fit+1,f->GetParError(2));\n", " cout<<\"\\\"\"<GetName()<<\"\\\",\"<GetParameter(i)<<\",\\\"\\302\\261\\\",\"<GetParError(i);\n", " }\n", " cout<<\"\\n\";\n", "}" ] }, { "cell_type": "markdown", "id": "f1c0a8e7", "metadata": {}, "source": [ "=========================== L-curve ==========================" ] }, { "cell_type": "code", "execution_count": null, "id": "b8f150f0", "metadata": { "collapsed": false }, "outputs": [], "source": [ "c4->cd(1);\n", "lCurve->SetTitle(\"L curve;log_{10} L_{x};log_{10} L_{y}\");\n", "lCurve->SetLineColor(kRed);\n", "lCurve->Draw(\"AL\");\n", "c4->cd(2);\n", "gPad->Clear();\n", "c4->cd(3);\n", "logTauX->SetTitle(\";log_{10} #tau;log_{10} L_{x}\");\n", "logTauX->SetLineColor(kBlue);\n", "logTauX->Draw();\n", "c4->cd(4);\n", "logTauY->SetTitle(\";log_{10} #tau;log_{10} L_{y}\");\n", "logTauY->SetLineColor(kBlue);\n", "logTauY->Draw();\n", "c4->SaveAs(\"lcurveL.eps\");" ] }, { "cell_type": "markdown", "id": "39f95312", "metadata": {}, "source": [ "========================= rho and L-curve scan ===============" ] }, { "cell_type": "code", "execution_count": null, "id": "56ec58c8", "metadata": { "collapsed": false }, "outputs": [], "source": [ "c4->cd(1);\n", "logTauCurvature->SetTitle(\";log_{10}(#tau);L curve curvature\");\n", "logTauCurvature->SetLineColor(kRed);\n", "logTauCurvature->Draw();\n", "c4->cd(2);\n", "rhoScan->SetTitle(\";log_{10}(#tau);average(#rho_{i})\");\n", "rhoScan->SetLineColor(kRed);\n", "rhoScan->Draw();\n", "c4->cd(3);\n", "DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFALCurveO,\"Tikhonov\",tauFA,\n", " &table[iFitFALCurve].second);\n", "c4->cd(4);\n", "DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFArhoO,\"Tikhonov\",tauFArho,\n", " &table[iFitFArho].second);\n", "\n", "c4->SaveAs(\"scanTau.eps\");\n", "\n", "\n", "TGraph *graphNiterAgoBay[4];\n", "GetNiterGraphs(iAgoFirstFit,iAgoFirstFit+1,table,kRed-2,graphNiterAgoBay,20);\n", "TGraph *graphNiterAgo[4];\n", "GetNiterGraphs(iAgoFirstFit,iAgoLastFit,table,kRed,graphNiterAgo,24);\n", "#ifdef WITH_IDS\n", "TGraph *graphNiterIDS[4];\n", "GetNiterGraphs(iIDSFirstFit,iIDSLastFit,table,kMagenta,graphNiterIDS,21);\n", "#endif\n", "TH1 *histNiterInversion[4];\n", "GetNiterHist(iFitInversion,table,histNiterInversion,kCyan,1,1001);\n", "TH1 *histNiterFALCurve[4];\n", "GetNiterHist(iFitFALCurve,table,histNiterFALCurve,kBlue,2,3353);\n", "TH1 *histNiterFArho[4];\n", "GetNiterHist(iFitFArho,table,histNiterFArho,kAzure-4,3,3353);\n", "TH1 *histNiterBinByBin[4];\n", "GetNiterHist(iFitBinByBin,table,histNiterBinByBin,kOrange+1,3,3335);\n", "\n", "histNiterInversion[0]->GetYaxis()->SetRangeUser(0.3,500.);\n", "histNiterInversion[1]->GetYaxis()->SetRangeUser(-0.1,0.9);\n", "histNiterInversion[2]->GetYaxis()->SetRangeUser(5.6,6.3);\n", "histNiterInversion[3]->GetYaxis()->SetRangeUser(1.6,2.4);\n", "\n", "TLine *line=nullptr;\n", "c1->cd();\n", "for(int i=0;i<2;i++) {\n", " gPad->Clear();\n", " gPad->SetLogx();\n", " gPad->SetLogy((i<1));\n", " if(! histNiterInversion[i]) continue;\n", " histNiterInversion[i]->Draw(\"][\");\n", " histNiterFALCurve[i]->Draw(\"SAME ][\");\n", " histNiterFArho[i]->Draw(\"SAME ][\");\n", " histNiterBinByBin[i]->Draw(\"SAME ][\");\n", " graphNiterAgo[i]->Draw(\"LP\");\n", " graphNiterAgoBay[i]->SetMarkerStyle(20);\n", " graphNiterAgoBay[i]->Draw(\"P\");\n", "#ifdef WITH_IDS\n", " graphNiterIDS[i]->Draw(\"LP\");\n", "#endif\n", " TLegend *legend;\n", " if(i==1) {\n", " legend=new TLegend(0.48,0.28,0.87,0.63);\n", " } else {\n", " legend=new TLegend(0.45,0.5,0.88,0.88);\n", " }\n", " legend->SetBorderSize(0);\n", " legend->SetFillStyle(1001);\n", " legend->SetFillColor(kWhite);\n", " legend->SetTextSize(kLegendFontSize*0.7);\n", " legend->AddEntry( histNiterInversion[0],\"inversion\",\"l\");\n", " legend->AddEntry( histNiterFALCurve[0],\"Tikhonov L-curve\",\"l\");\n", " legend->AddEntry( histNiterFArho[0],\"Tikhonov global cor.\",\"l\");\n", " legend->AddEntry( histNiterBinByBin[0],\"bin-by-bin\",\"l\");\n", " legend->AddEntry( graphNiterAgoBay[0],\"\\\"Bayesian\\\"\",\"p\");\n", " legend->AddEntry( graphNiterAgo[0],\"iterative\",\"p\");\n", "#ifdef WITH_IDS\n", " legend->AddEntry( graphNiterIDS[0],\"IDS\",\"p\");\n", "#endif\n", " legend->Draw();\n", "\n", " c1->SaveAs(TString::Format(\"niter%d.eps\",i));\n", "}" ] }, { "cell_type": "markdown", "id": "20f51b9e", "metadata": {}, "source": [ "c4->cd(1);\n", "DrawPadCorrelations(histRhoFALCurveO,&table);" ] }, { "cell_type": "code", "execution_count": null, "id": "d4a17ee8", "metadata": { "collapsed": false }, "outputs": [], "source": [ "c2sq->cd(1);\n", "DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFALCurveO,\"Tikhonov\",tauFA,\n", " &table[iFitFALCurve].second,table[iFitFALCurve].first);\n", "\n", "c2sq->cd(2);\n", "gPad->SetLogx();\n", "gPad->SetLogy(false);\n", "histNiterInversion[3]->DrawClone(\"E2\");\n", "histNiterInversion[3]->SetFillStyle(0);\n", "histNiterInversion[3]->Draw(\"SAME HIST ][\");\n", "histNiterFALCurve[3]->DrawClone(\"SAME E2\");\n", "histNiterFALCurve[3]->SetFillStyle(0);\n", "histNiterFALCurve[3]->Draw(\"SAME HIST ][\");\n", "histNiterFArho[3]->DrawClone(\"SAME E2\");\n", "histNiterFArho[3]->SetFillStyle(0);\n", "histNiterFArho[3]->Draw(\"SAME HIST ][\");\n", "histNiterBinByBin[3]->DrawClone(\"SAME E2\");\n", "histNiterBinByBin[3]->SetFillStyle(0);\n", "histNiterBinByBin[3]->Draw(\"SAME HIST ][\");\n", "double yTrue=1.8;\n", "line=new TLine(histNiterInversion[3]->GetXaxis()->GetXmin(),\n", " yTrue,\n", " histNiterInversion[3]->GetXaxis()->GetXmax(),\n", " yTrue);\n", "line->SetLineWidth(3);\n", "line->Draw();\n", "\n", "graphNiterAgo[3]->Draw(\"LP\");\n", "graphNiterAgoBay[3]->SetMarkerStyle(20);\n", "graphNiterAgoBay[3]->Draw(\"P\");\n", "#ifdef WITH_IDS\n", "graphNiterIDS[3]->Draw(\"LP\");\n", "#endif\n", "\n", "TLegend *legend;\n", "legend=new TLegend(0.55,0.53,0.95,0.97);\n", "legend->SetBorderSize(0);\n", "legend->SetFillStyle(1001);\n", "legend->SetFillColor(kWhite);\n", "legend->SetTextSize(kLegendFontSize);\n", "legend->AddEntry( line,\"truth\",\"l\");\n", "legend->AddEntry( histNiterInversion[3],\"inversion\",\"l\");\n", "legend->AddEntry( histNiterFALCurve[3],\"Tikhonov L-curve\",\"l\");\n", "legend->AddEntry( histNiterFArho[3],\"Tikhonov global cor.\",\"l\");\n", "legend->AddEntry( histNiterBinByBin[3],\"bin-by-bin\",\"l\");\n", "legend->AddEntry( graphNiterAgoBay[3],\"\\\"Bayesian\\\"\",\"p\");\n", "legend->AddEntry( graphNiterAgo[3],\"iterative\",\"p\");\n", "#ifdef WITH_IDS\n", "legend->AddEntry( graphNiterIDS[3],\"IDS\",\"p\");\n", "#endif\n", "legend->Draw();\n", "\n", "c2sq->SaveAs(\"fitSigma.eps\");\n", "\n", "outputFile->Write();\n", "\n", "delete outputFile;" ] }, { "cell_type": "markdown", "id": "2099af6e", "metadata": {}, "source": [ "Draw all canvases " ] }, { "cell_type": "code", "execution_count": null, "id": "ad62a71a", "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 }