This is an example of unfolding a two-dimensional distribution also using an auxiliary measurement to constrain some background
#include <iostream>
#include <cmath>
#include <map>
using namespace std;
#define TEST_INPUT_COVARIANCE
void testUnfold5d()
{
TFile *outputFile=
new TFile(
"testUnfold5_results.root",
"recreate");
TFile *inputFile=
new TFile(
"testUnfold5_histograms.root");
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);
#ifdef TEST_ZERO_UNCORR_ERROR
for(
int i=0;i<=histMCGenRec->
GetNbinsX()+1;i++) {
for(
int j=0;j<=histMCGenRec->
GetNbinsY()+1;j++) {
}
}
#endif
if((!histDataReco)||(!histDataTruth)||(!histMCGenRec)) {
cout<<"problem to read input histograms\n";
}
const char *REGULARISATION_DISTRIBUTION=nullptr;
const char *REGULARISATION_AXISSTEERING="*[B]";
regMode,constraintMode,densityFlags,
generatorBinning,detectorBinning,
REGULARISATION_DISTRIBUTION,
REGULARISATION_AXISSTEERING);
#ifdef TEST_INPUT_COVARIANCE
for(
int i=1;i<=inputEmatrix->
GetNbinsX();i++) {
}
unfold.SetInput(histDataReco,0.0,0.0,inputEmatrix);
#else
unfold.SetInput(histDataReco );
#endif
#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
const char *SCAN_DISTRIBUTION="signal";
const char *SCAN_AXISSTEERING=nullptr;
Int_t iBest=unfold.ScanTau(nScan,0.,0.,&rhoLogTau,
SCAN_DISTRIBUTION,SCAN_AXISSTEERING,
&lCurve);
rhoLogTau->
GetKnot(iBest,t[0],rho[0]);
for(
Int_t i=0;i<nScan;i++) {
rhoLogTau->
GetKnot(i,tAll[i],rhoAll[i]);
}
cout<<"chi**2="<<unfold.GetChi2A()<<"+"<<unfold.GetChi2L()
<<" / "<<unfold.GetNdf()<<"\n";
TH1 *histDataUnfold=unfold.GetOutput(
"unfolded signal",
nullptr,
nullptr,
nullptr,
kFALSE);
histMCReco->
Scale(scaleFactor);
histMCTruth->
Scale(scaleFactor);
TH2 *histProbability=unfold.GetProbabilityMatrix(
"histProbability");
unfold.GetRhoItotal(
"histGlobalCorr",
nullptr,
nullptr,
nullptr,
kFALSE);
TH1 *histGlobalCorrScan=unfold.GetRhoItotal
(
"histGlobalCorrScan",
nullptr,SCAN_DISTRIBUTION,SCAN_AXISSTEERING,
kFALSE);
unfold.GetRhoIJtotal(
"histCorrCoeff",
nullptr,
nullptr,
nullptr,
kFALSE);
canvas.
Print(
"testUnfold5.ps[");
histMCReco->
Draw(
"SAME HIST");
histProbability->
Draw(
"BOX");
histDataUnfold->
Draw(
"E");
histDataTruth->
Draw(
"SAME HIST");
histMCTruth->
Draw(
"SAME HIST");
bestRhoLogTau->
Draw(
"*");
histGlobalCorrScan->
Draw(
"HIST");
canvas.
Print(
"testUnfold5.ps");
canvas.
Print(
"testUnfold5.ps]");
}
virtual void SetLineColor(Color_t lcolor)
Set the line color.
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
void Clear(Option_t *option="") override
Remove all primitives from the canvas.
TVirtualPad * cd(Int_t subpadnumber=0) override
Set current canvas & pad.
Bool_t cd() override
Change current directory to "this" directory.
void GetObject(const char *namecycle, T *&ptr)
Get an object with proper type checking.
A ROOT file is composed of a header, followed by consecutive data records (TKey instances) with a wel...
A TGraph is an object made of two arrays X and Y with npoints each.
void Draw(Option_t *chopt="") override
Draw this graph with its current attributes.
virtual Int_t GetPoint(Int_t i, Double_t &x, Double_t &y) const
Get x and y values for point number i.
TH1 is the base class of all histogram classes in ROOT.
virtual Int_t GetNbinsY() const
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
virtual Int_t GetNbinsX() const
virtual void SetBinError(Int_t bin, Double_t error)
Set the bin Error Note that this resets the bin eror option to be of Normal Type and for the non-empt...
void Draw(Option_t *option="") override
Draw this histogram with options.
virtual void SetMinimum(Double_t minimum=-1111)
static void SetDefaultSumw2(Bool_t sumw2=kTRUE)
When this static function is called with sumw2=kTRUE, all new histograms will automatically activate ...
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
virtual Double_t GetSumOfWeights() const
Return the sum of weights excluding under/overflows.
2-D histogram with a double per channel (see TH1 documentation)}
Service class for 2-D histogram classes.
TH1D * ProjectionY(const char *name="_py", Int_t firstxbin=0, Int_t lastxbin=-1, Option_t *option="") const
Project a 2-D histogram into a 1-D histogram along Y.
void SetBinContent(Int_t bin, Double_t content) override
Set bin content.
TH1D * ProjectionX(const char *name="_px", Int_t firstybin=0, Int_t lastybin=-1, Option_t *option="") const
Project a 2-D histogram into a 1-D histogram along X.
Double_t GetBinContent(Int_t binx, Int_t biny) const override
virtual Int_t Write(const char *name=nullptr, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0) override
Automatic pad generation by division.
void Print(const char *filename="") const override
This method is equivalent to SaveAs("filename"). See TPad::SaveAs for details.
Base class for spline implementation containing the Draw/Paint methods.
void Draw(Option_t *option="") override
Draw this function with its current attributes.
virtual void GetKnot(Int_t i, Double_t &x, Double_t &y) const =0
Binning schemes for use with the unfolding algorithm TUnfoldDensity.
TH2D * CreateErrorMatrixHistogram(const char *histogramName, Bool_t originalAxisBinning, Int_t **binMap=nullptr, const char *histogramTitle=nullptr, const char *axisSteering=nullptr) const
Create a TH2D histogram capable to hold a covariance matrix.
An algorithm to unfold distributions from detector to truth level.
@ kEScanTauRhoMax
maximum global correlation coefficient (from TUnfold::GetRhoI())
EDensityMode
choice of regularisation scale factors to cinstruct the matrix L
@ kDensityModeBinWidth
scale factors from multidimensional bin width
EConstraint
type of extra constraint
@ kEConstraintArea
enforce preservation of the area
ERegMode
choice of regularisation scheme
@ kRegModeCurvature
regularize the 2nd derivative of the output distribution
@ kHistMapOutputHoriz
truth level on x-axis of the response matrix
virtual void SetLogy(Int_t value=1)=0