Simple Test program for the class TUnfoldDensity.
1-dimensional unfolding with background subtraction
the goal is to unfold the underlying "true" distribution of a variable Pt
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
- plus underflow bin from 0 to 6
- plus overflow bin above 26 there are two background sources bgr1 and bgr2 the signal has a finite trigger efficiency at a threshold of 8 GeV
one type of systematic error is studied, where the signal parameters are changed
Finally, the unfolding is compared to a "bin-by-bin" correction method
Processing /mnt/build/workspace/root-makedoc-v610/rootspi/rdoc/src/v6-10-00-patches/tutorials/unfold/testUnfold3.C...
TUnfold version is V17.6
chi**2=25.3912+4.28619 / 11
bin truth result (stat) (bgr) (sys)
===+=====+=========+========+========+=======
1 4002 4261.3 +/-137.0 +/- 8.9 +/-473.2 (unfolding)
4190.1 +/-128.7 +/-104.7 +/-266.3 (bin-by-bin)
2 1816 1775.1 +/- 82.0 +/- 7.3 +/- 25.2 (unfolding)
1714.1 +/- 61.0 +/- 17.0 +/- 70.8 (bin-by-bin)
3 1011 1020.9 +/- 60.1 +/- 6.7 +/- 24.8 (unfolding)
959.5 +/- 40.9 +/- 9.2 +/- 43.3 (bin-by-bin)
4 593 526.9 +/- 45.2 +/- 6.1 +/- 15.0 (unfolding)
506.2 +/- 27.9 +/- 6.6 +/- 45.2 (bin-by-bin)
5 379 371.6 +/- 38.6 +/- 6.5 +/- 19.8 (unfolding)
351.8 +/- 22.6 +/- 6.0 +/- 29.2 (bin-by-bin)
6 256 341.5 +/- 34.4 +/- 6.0 +/- 9.2 (unfolding)
282.5 +/- 19.2 +/- 5.2 +/- 46.6 (bin-by-bin)
7 193 182.0 +/- 30.1 +/- 5.7 +/- 7.5 (unfolding)
177.7 +/- 15.8 +/- 4.9 +/- 22.5 (bin-by-bin)
8 137 147.1 +/- 28.4 +/- 5.6 +/- 7.0 (unfolding)
133.6 +/- 13.9 +/- 4.5 +/- 20.5 (bin-by-bin)
9 123 109.6 +/- 26.3 +/- 5.4 +/- 10.4 (unfolding)
98.6 +/- 12.1 +/- 4.1 +/- 20.7 (bin-by-bin)
10 78 100.7 +/- 24.7 +/- 6.1 +/- 8.4 (unfolding)
89.7 +/- 13.2 +/- 5.0 +/- 17.1 (bin-by-bin)
{
do {
itype=0;
if(f<parm[0]) itype=1;
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";
}
}
Version 17.6, in parallel to changes in TUnfold
History:
- Version 17.5, in parallel to changes in TUnfold
- Version 17.4, in parallel to changes in TUnfold
- Version 17.3, in parallel to changes in TUnfold
- Version 17.2, in parallel to changes in TUnfold
- Version 17.1, in parallel to changes in TUnfold
- Version 17.0, change to use the class TUnfoldDensity
- Version 16.1, parallel to changes in TUnfold
- Version 16.0, parallel to changes in TUnfold
- Version 15, simple example including background subtraction
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/.
- Author
- Stefan Schmitt DESY, 14.10.2008
Definition in file testUnfold3.C.