the reconstructed Pt is measured in 24 bins from 4 to 28 the generator-level Pt is unfolded into 10 bins from 6 to 26
{
do {
itype=0;
else if(
f<parm[0]+parm[1]) itype=2;
if(a1 == 0.0) {
} else {
ptgen=pow(t*(pow(parm[3],a1)-x0)+x0,1./a1);
}
if(iType) *iType=itype;
if(ptGen) *ptGen=ptgen;
TMath::Sqrt(parm[7]*parm[7]*ptgen+parm[8]*parm[8]*ptgen*ptgen);
triggerParm[2]/(1.+
TMath::Exp((triggerParm[0]-ptObs)/triggerParm[1]));
isTriggered= rnd->
Rndm()<triggerProb;
(*intLumi) ++;
} while((!triggerFlag) && (!isTriggered));
if(triggerFlag) *triggerFlag=isTriggered;
return ptObs;
}
void testUnfold3()
{
new TH1D(
"unfolding input rec",
";ptrec",nDet,xminDet,xmaxDet);
new TH2D(
"unfolding matrix",
";ptgen;ptrec",
nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet);
new TH1D(
"unfolding bgr1 rec",
";ptrec",nDet,xminDet,xmaxDet);
new TH1D(
"unfolding bgr2 rec",
";ptrec",nDet,xminDet,xmaxDet);
TH2D *histUnfoldMatrixSys=
new TH2D(
"unfolding matrix sys",
";ptgen;ptrec",
nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet);
new TH1D(
"bbb input rec",
";ptrec",nGen,xminGen,xmaxGen);
new TH1D(
"bbb signal rec",
";ptrec",nGen,xminGen,xmaxGen);
new TH1D(
"bbb bgr1 rec",
";ptrec",nGen,xminGen,xmaxGen);
new TH1D(
"bbb bgr2 rec",
";ptrec",nGen,xminGen,xmaxGen);
new TH1D(
"bbb bgrPt rec",
";ptrec",nGen,xminGen,xmaxGen);
new TH1D(
"bbb signal gen",
";ptgen",nGen,xminGen,xmaxGen);
TH1D *histBbbSignalRecSys=
new TH1D(
"bbb signal sys rec",
";ptrec",nGen,xminGen,xmaxGen);
new TH1D(
"bbb bgrPt sys rec",
";ptrec",nGen,xminGen,xmaxGen);
TH1D *histBbbSignalGenSys=
new TH1D(
"bbb signal sys gen",
";ptgen",nGen,xminGen,xmaxGen);
new TH1D(
"DATA truth gen",
";ptgen",nGen,xminGen,xmaxGen);
new TH1D(
"MC prediction rec",
";ptrec",nDet,xminDet,xmaxDet);
0.05,
0.05,
3.5,
100.,
-3.0,
0.1,
-0.5,
0.2,
0.01,
};
static Double_t triggerParm_DATA[]={8.0,
4.0,
0.95
};
while(intLumi<lumiData) {
Double_t ptObs=GenerateEvent(parm_DATA,triggerParm_DATA,&intLumi,
&isTriggered,&ptGen,&iTypeGen);
if(isTriggered) {
histUnfoldInput->
Fill(ptObs);
histBbbInput->
Fill(ptObs);
}
if(iTypeGen==0) histDataTruth->
Fill(ptGen);
}
0.05,
0.05,
3.5,
100.,
-4.0,
0.1,
-0.5,
0.2,
0.01,
};
4.0,
0.95
};
intLumi=0.0;
while(intLumi<lumiMC) {
Double_t ptObs=GenerateEvent(parm_MC,triggerParm_MC,&intLumi,&isTriggered,
&ptGen,&iTypeGen);
if(!isTriggered) ptObs=0.0;
if(iTypeGen==0) {
histUnfoldMatrix->
Fill(ptGen,ptObs,lumiWeight);
} else if(iTypeGen==1) {
histUnfoldBgr1->
Fill(ptObs,lumiWeight);
} else if(iTypeGen==2) {
histUnfoldBgr2->
Fill(ptObs,lumiWeight);
}
if(iTypeGen==0) {
if((ptGen>=xminGen)&&(ptGen<xmaxGen)) {
histBbbSignalRec->
Fill(ptObs,lumiWeight);
} else {
histBbbBgrPt->
Fill(ptObs,lumiWeight);
}
histBbbSignalGen->
Fill(ptGen,lumiWeight);
} else if(iTypeGen==1) {
histBbbBgr1->
Fill(ptObs,lumiWeight);
} else if(iTypeGen==2) {
histBbbBgr2->
Fill(ptObs,lumiWeight);
}
histDetMC ->
Fill(ptObs,lumiWeight);
}
0.05,
0.05,
3.5,
100.,
-2.0,
0.1,
-0.5,
0.2,
0.01,
};
intLumi=0.0;
while(intLumi<lumiMC) {
Double_t ptObs=GenerateEvent(parm_MC_SYS,triggerParm_MC,&intLumi,
&isTriggered,&ptGen,&iTypeGen);
if(!isTriggered) ptObs=0.0;
if(iTypeGen==0) {
histUnfoldMatrixSys->
Fill(ptGen,ptObs);
}
if(iTypeGen==0) {
if((ptGen>=xminGen)&&(ptGen<xmaxGen)) {
histBbbSignalRecSys->
Fill(ptObs);
} else {
histBbbBgrPtSys->
Fill(ptObs);
}
histBbbSignalGenSys->
Fill(ptGen);
}
}
regMode,constraintMode,densityFlags);
unfold.SetInput(histUnfoldInput);
unfold.SubtractBackground(histUnfoldBgr1,"background1",scale_bgr,dscale_bgr);
unfold.SubtractBackground(histUnfoldBgr2,"background2",scale_bgr,dscale_bgr);
unfold.AddSysError(histUnfoldMatrixSys,"signalshape_SYS",
Int_t iBest=unfold.ScanLcurve(nScan,0.,0.,&lCurve,&logTauX,&logTauY);
cout<<"chi**2="<<unfold.GetChi2A()<<"+"<<unfold.GetChi2L()
<<" / "<<unfold.GetNdf()<<"\n";
TH1 *histUnfoldOutput=unfold.GetOutput(
"PT(unfold,stat+bgrerr)");
TH2 *histEmatStat=unfold.GetEmatrixInput(
"unfolding stat error matrix");
TH2 *histEmatTotal=unfold.GetEmatrixTotal(
"unfolding total error matrix");
TH1 *histUnfoldStat=
new TH1D(
"PT(unfold,staterr)",
";Pt(gen)",
nGen,xminGen,xmaxGen);
TH1 *histUnfoldTotal=
new TH1D(
"PT(unfold,totalerr)",
";Pt(gen)",
nGen,xminGen,xmaxGen);
for(
Int_t i=0;i<nGen+2;i++) {
}
TH2D *histCorr=
new TH2D(
"Corr(total)",
";Pt(gen);Pt(gen)",
nGen,xminGen,xmaxGen,nGen,xminGen,xmaxGen);
for(
Int_t i=0;i<nGen+2;i++) {
if(ei<=0.0) continue;
for(
Int_t j=0;j<nGen+2;j++) {
if(ej<=0.0) continue;
}
}
TH1 *histDetNormBgr1=unfold.GetBackground(
"bgr1 normalized",
"background1");
TH1 *histDetNormBgrTotal=unfold.GetBackground(
"bgr total normalized");
histUnfoldInput->
Draw(
"E");
histDetMC->
Draw(
"SAME HIST");
histDetNormBgr1->
Draw(
"SAME HIST");
histDetNormBgrTotal->
Draw(
"SAME HIST");
histUnfoldTotal->
Draw(
"E");
histUnfoldOutput->
Draw(
"SAME E1");
histUnfoldStat->
Draw(
"SAME E1");
histDataTruth->
Draw(
"SAME HIST");
histBbbSignalGen->
Draw(
"SAME HIST");
histUnfoldMatrix->
Draw(
"BOX");
bestLogTauLogChi2->
Draw(
"*");
output.
SaveAs(
"testUnfold3.ps");
std::cout<<"bin truth result (stat) (bgr) (sys)\n";
std::cout<<"===+=====+=========+========+========+=======\n";
for(
Int_t i=1;i<=nGen;i++) {
(errData*errData+
Double_t errData_stat_bbb = errData*fCorr;
Double_t errData_statbgr_bbb = errData_bgr*fCorr;
Double_t shift_sys= data_bgr*fCorr_sys - data_bbb;
+shift_sys*shift_sys);
("%3d %5.0f %8.1f +/-%5.1f +/-%5.1f +/-%5.1f (unfolding)",
errData_stat_unfold,
TMath::Sqrt(errData_statbgr_unfold*
errData_statbgr_unfold-
errData_stat_unfold*
errData_stat_unfold),
errData_total_unfold-
errData_statbgr_unfold*
errData_statbgr_unfold))<<"\n";
(" %8.1f +/-%5.1f +/-%5.1f +/-%5.1f (bin-by-bin)",
data_bbb,errData_stat_bbb,
TMath::Sqrt(errData_statbgr_bbb*
errData_statbgr_bbb-
errData_stat_bbb*
errData_stat_bbb),
errData_total_bbb-
errData_statbgr_bbb*
errData_statbgr_bbb))
<<"\n\n";
}
}
int Int_t
Signed integer 4 bytes (int).
bool Bool_t
Boolean (0=false, 1=true) (bool).
double Double_t
Double 8 bytes.
virtual void SetLineColor(Color_t lcolor)
Set the line color.
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
TVirtualPad * cd(Int_t subpadnumber=0) override
Set current canvas & pad.
void Draw(Option_t *chopt="") override
Default Draw method for all objects.
1-D histogram with a double per channel (see TH1 documentation)
TH1 is the base class of all histogram classes in ROOT.
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
virtual Double_t GetMaximum(Double_t maxval=FLT_MAX) const
Return maximum value smaller than maxval of bins in the range, unless the value has been overridden b...
virtual void SetMaximum(Double_t maximum=-1111)
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...
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
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 SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
2-D histogram with a double per channel (see TH1 documentation)
Service class for 2-D histogram classes.
void SetBinContent(Int_t bin, Double_t content) override
Set bin content.
Double_t GetBinContent(Int_t binx, Int_t biny) const override
Int_t Fill(Double_t) override
Invalid Fill method.
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 SaveAs(const char *filename="", Option_t *option="") const override
Save the pad content in a file.
Random number generator class based on M.
This is the base class for the ROOT Random number generators.
Double_t Rndm() override
Machine independent random number generator.
virtual Double_t BreitWigner(Double_t mean=0, Double_t gamma=1)
Return a number distributed following a BreitWigner function with mean and gamma.
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
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString.
An algorithm to unfold distributions from detector to truth level.
EDensityMode
choice of regularisation scale factors to cinstruct the matrix L
@ kDensityModeBinWidth
scale factors from multidimensional bin width
@ kSysErrModeMatrix
matrix is an alternative to the default matrix, the errors are the difference to the original matrix
static const char * GetTUnfoldVersion(void)
return a string describing the TUnfold version
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
Double_t Exp(Double_t x)
Returns the base-e exponential function of x, which is e raised to the power x.
Double_t Log(Double_t x)
Returns the natural logarithm of x.
Double_t Sqrt(Double_t x)
Returns the square root of x.
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
This file is part of TUnfold.
TUnfold is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
TUnfold is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with TUnfold. If not, see http://www.gnu.org/licenses/.