Version 17.0 example for multi-dimensional unfolding.
#include <iostream>
#include <cmath>
#include <map>
using namespace std;
void testUnfold5d()
{
TFile *outputFile=new TFile("testUnfold5_results.root","recreate");
TFile *inputFile=new TFile("testUnfold5_histograms.root");
outputFile->cd();
inputFile->GetObject("detector",detectorBinning);
inputFile->GetObject("generator",generatorBinning);
if((!detectorBinning)||(!generatorBinning)) {
cout<<"problem to read binning schemes\n";
}
detectorBinning->
Write();
generatorBinning->
Write();
TH1 *histDataReco,*histDataTruth;
inputFile->GetObject("histDataReco",histDataReco);
inputFile->GetObject("histDataTruth",histDataTruth);
inputFile->GetObject("histMCGenRec",histMCGenRec);
if((!histDataReco)||(!histDataTruth)||(!histMCGenRec)) {
cout<<"problem to read input histograms\n";
}
const char *REGULARISATION_DISTRIBUTION=0;
const char *REGULARISATION_AXISSTEERING="*[B]";
regMode,constraintMode,densityFlags,
generatorBinning,detectorBinning,
REGULARISATION_DISTRIBUTION,
REGULARISATION_AXISSTEERING);
unfold.SetInput(histDataReco );
#ifdef PRINT_MATRIX_L
TH2 *histL= unfold.GetL(
"L");
cout<<"L["<<unfold.GetLBinning()->GetBinName(j)<<"]";
if(c!=0.0) cout<<
" ["<<i<<
"]="<<
c;
}
cout<<"\n";
}
#endif
TGraph *lCurve=0;
const char *SCAN_DISTRIBUTION="signal";
const char *SCAN_AXISSTEERING=0;
Int_t iBest=unfold.ScanTau(nScan,0.,0.,&rhoLogTau,
SCAN_DISTRIBUTION,SCAN_AXISSTEERING,
&lCurve);
rhoLogTau->
GetKnot(iBest,t[0],rho[0]);
lCurve->GetPoint(iBest,x[0],y[0]);
TGraph *bestRhoLogTau=new TGraph(1,t,rho);
TGraph *bestLCurve=new TGraph(1,x,y);
for(
Int_t i=0;i<nScan;i++) {
rhoLogTau->
GetKnot(i,tAll[i],rhoAll[i]);
}
TGraph *knots=new TGraph(nScan,tAll,rhoAll);
cout<<"chi**2="<<unfold.GetChi2A()<<"+"<<unfold.GetChi2L()
<<" / "<<unfold.GetNdf()<<"\n";
TH1 *histDataUnfold=unfold.GetOutput(
"unfolded signal",0,0,0,
kFALSE);
histMCReco->
Scale(scaleFactor);
histMCTruth->
Scale(scaleFactor);
TH2 *histProbability=unfold.GetProbabilityMatrix(
"histProbability");
TH1 *histGlobalCorr=unfold.GetRhoItotal(
"histGlobalCorr",0,0,0,
kFALSE);
TH1 *histGlobalCorrScan=unfold.GetRhoItotal
(
"histGlobalCorrScan",0,SCAN_DISTRIBUTION,SCAN_AXISSTEERING,
kFALSE);
TH2 *histCorrCoeff=unfold.GetRhoIJtotal(
"histCorrCoeff",0,0,0,
kFALSE);
TCanvas *canvas = new TCanvas();
canvas->Print("testUnfold5.ps[");
canvas->Clear();
canvas->Divide(3,2);
canvas->cd(1);
histMCReco->
Draw(
"SAME HIST");
canvas->cd(2);
histProbability->
Draw(
"BOX");
canvas->cd(3);
histDataUnfold->
Draw(
"E");
histDataTruth->
Draw(
"SAME HIST");
histMCTruth->
Draw(
"SAME HIST");
canvas->cd(4);
knots->Draw("*");
bestRhoLogTau->SetMarkerColor(
kRed);
bestRhoLogTau->Draw("*");
canvas->cd(5);
histGlobalCorrScan->
Draw(
"HIST");
canvas->cd(6);
lCurve->Draw("AL");
bestLCurve->SetMarkerColor(
kRed);
bestLCurve->Draw("*");
canvas->Print("testUnfold5.ps");
canvas->Print("testUnfold5.ps]");
}