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
using namespace std;
{
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";
}
}
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
R__EXTERN TStyle * gStyle
virtual void SetLineColor(Color_t lcolor)
Set the line color.
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
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.
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.
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.
void SetOptFit(Int_t fit=1)
The type of information about fit parameters printed in the histogram statistics box can be selected ...
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.