105 if(f<parm[0]) itype=1;
106 else if(f<parm[0]+parm[1]) itype=2;
117 ptgen=
pow(t*(
pow(parm[3],a1)-x0)+x0,1./a1);
119 if(iType) *iType=itype;
120 if(ptGen) *ptGen=ptgen;
124 TMath::Sqrt(parm[7]*parm[7]*ptgen+parm[8]*parm[8]*ptgen*ptgen);
129 triggerParm[2]/(1.+
TMath::Exp((triggerParm[0]-ptObs)/triggerParm[1]));
130 isTriggered= rnd->
Rndm()<triggerProb;
132 }
while((!triggerFlag) && (!isTriggered));
134 if(triggerFlag) *triggerFlag=isTriggered;
193 TH1D *histUnfoldInput=
194 new TH1D(
"unfolding input rec",
";ptrec",nDet,xminDet,xmaxDet);
195 TH2D *histUnfoldMatrix=
196 new TH2D(
"unfolding matrix",
";ptgen;ptrec",
197 nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet);
198 TH1D *histUnfoldBgr1=
199 new TH1D(
"unfolding bgr1 rec",
";ptrec",nDet,xminDet,xmaxDet);
200 TH1D *histUnfoldBgr2=
201 new TH1D(
"unfolding bgr2 rec",
";ptrec",nDet,xminDet,xmaxDet);
202 TH2D *histUnfoldMatrixSys=
203 new TH2D(
"unfolding matrix sys",
";ptgen;ptrec",
204 nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet);
208 new TH1D(
"bbb input rec",
";ptrec",nGen,xminGen,xmaxGen);
209 TH1D *histBbbSignalRec=
210 new TH1D(
"bbb signal rec",
";ptrec",nGen,xminGen,xmaxGen);
212 new TH1D(
"bbb bgr1 rec",
";ptrec",nGen,xminGen,xmaxGen);
214 new TH1D(
"bbb bgr2 rec",
";ptrec",nGen,xminGen,xmaxGen);
216 new TH1D(
"bbb bgrPt rec",
";ptrec",nGen,xminGen,xmaxGen);
217 TH1D *histBbbSignalGen=
218 new TH1D(
"bbb signal gen",
";ptgen",nGen,xminGen,xmaxGen);
219 TH1D *histBbbSignalRecSys=
220 new TH1D(
"bbb signal sys rec",
";ptrec",nGen,xminGen,xmaxGen);
221 TH1D *histBbbBgrPtSys=
222 new TH1D(
"bbb bgrPt sys rec",
";ptrec",nGen,xminGen,xmaxGen);
223 TH1D *histBbbSignalGenSys=
224 new TH1D(
"bbb signal sys gen",
";ptgen",nGen,xminGen,xmaxGen);
228 new TH1D(
"DATA truth gen",
";ptgen",nGen,xminGen,xmaxGen);
230 new TH1D(
"MC prediction rec",
";ptrec",nDet,xminDet,xmaxDet);
246 static Double_t triggerParm_DATA[]={8.0,
252 while(intLumi<lumiData) {
256 Double_t ptObs=GenerateEvent(parm_DATA,triggerParm_DATA,&intLumi,
257 &isTriggered,&ptGen,&iTypeGen);
260 histUnfoldInput->
Fill(ptObs);
263 histBbbInput->
Fill(ptObs);
266 if(iTypeGen==0) histDataTruth->
Fill(ptGen);
286 static Double_t triggerParm_MC[]={8.0,
292 Double_t lumiWeight = lumiData/lumiMC;
294 while(intLumi<lumiMC) {
298 Double_t ptObs=GenerateEvent(parm_MC,triggerParm_MC,&intLumi,&isTriggered,
300 if(!isTriggered) ptObs=0.0;
305 histUnfoldMatrix->
Fill(ptGen,ptObs,lumiWeight);
306 }
else if(iTypeGen==1) {
307 histUnfoldBgr1->
Fill(ptObs,lumiWeight);
308 }
else if(iTypeGen==2) {
309 histUnfoldBgr2->
Fill(ptObs,lumiWeight);
314 if((ptGen>=xminGen)&&(ptGen<xmaxGen)) {
315 histBbbSignalRec->
Fill(ptObs,lumiWeight);
317 histBbbBgrPt->
Fill(ptObs,lumiWeight);
319 histBbbSignalGen->
Fill(ptGen,lumiWeight);
320 }
else if(iTypeGen==1) {
321 histBbbBgr1->
Fill(ptObs,lumiWeight);
322 }
else if(iTypeGen==2) {
323 histBbbBgr2->
Fill(ptObs,lumiWeight);
327 histDetMC ->
Fill(ptObs,lumiWeight);
347 while(intLumi<lumiMC) {
351 Double_t ptObs=GenerateEvent(parm_MC_SYS,triggerParm_MC,&intLumi,
352 &isTriggered,&ptGen,&iTypeGen);
353 if(!isTriggered) ptObs=0.0;
357 histUnfoldMatrixSys->
Fill(ptGen,ptObs);
362 if((ptGen>=xminGen)&&(ptGen<xmaxGen)) {
363 histBbbSignalRecSys->
Fill(ptObs);
365 histBbbBgrPtSys->
Fill(ptObs);
367 histBbbSignalGenSys->
Fill(ptGen);
394 regMode,constraintMode,densityFlags);
397 unfold.SetInput(histUnfoldInput);
403 unfold.SubtractBackground(histUnfoldBgr1,
"background1",scale_bgr,dscale_bgr);
404 unfold.SubtractBackground(histUnfoldBgr2,
"background2",scale_bgr,dscale_bgr);
407 unfold.AddSysError(histUnfoldMatrixSys,
"signalshape_SYS",
417 Int_t iBest=unfold.ScanLcurve(nScan,0.,0.,&lCurve,&logTauX,&logTauY);
419 cout<<
"chi**2="<<unfold.GetChi2A()<<
"+"<<unfold.GetChi2L()
420 <<
" / "<<unfold.GetNdf()<<
"\n";
424 logTauX->
GetKnot(iBest,t[0],x[0]);
425 logTauY->
GetKnot(iBest,t[0],y[0]);
436 TH1 *histUnfoldOutput=unfold.GetOutput(
"PT(unfold,stat+bgrerr)");
439 TH2 *histEmatStat=unfold.GetEmatrixInput(
"unfolding stat error matrix");
443 TH2 *histEmatTotal=unfold.GetEmatrixTotal(
"unfolding total error matrix");
447 TH1 *histUnfoldStat=
new TH1D(
"PT(unfold,staterr)",
";Pt(gen)",
448 nGen,xminGen,xmaxGen);
449 TH1 *histUnfoldTotal=
new TH1D(
"PT(unfold,totalerr)",
";Pt(gen)",
450 nGen,xminGen,xmaxGen);
451 for(
Int_t i=0;i<nGen+2;i++) {
466 TH2D *histCorr=
new TH2D(
"Corr(total)",
";Pt(gen);Pt(gen)",
467 nGen,xminGen,xmaxGen,nGen,xminGen,xmaxGen);
468 for(
Int_t i=0;i<nGen+2;i++) {
471 if(ei<=0.0)
continue;
472 for(
Int_t j=0;j<nGen+2;j++) {
474 if(ej<=0.0)
continue;
480 TH1 *histDetNormBgr1=unfold.GetBackground(
"bgr1 normalized",
482 TH1 *histDetNormBgrTotal=unfold.GetBackground(
"bgr total normalized");
492 histUnfoldInput->
Draw(
"E");
497 histDetMC->
Draw(
"SAME HIST");
498 histDetNormBgr1->
Draw(
"SAME HIST");
499 histDetNormBgrTotal->
Draw(
"SAME HIST");
506 histUnfoldTotal->
Draw(
"E");
508 histUnfoldOutput->
Draw(
"SAME E1");
510 histUnfoldStat->
Draw(
"SAME E1");
511 histDataTruth->
Draw(
"SAME HIST");
513 histBbbSignalGen->
Draw(
"SAME HIST");
518 histUnfoldMatrix->
Draw(
"BOX");
524 bestLogTauLogChi2->
Draw(
"*");
530 bestLcurve->
Draw(
"*");
534 histCorr->
Draw(
"BOX");
539 std::cout<<
"bin truth result (stat) (bgr) (sys)\n";
540 std::cout<<
"===+=====+=========+========+========+=======\n";
541 for(
Int_t i=1;i<=nGen;i++) {
564 Double_t errData_stat_bbb = errData*fCorr;
566 Double_t errData_statbgr_bbb = errData_bgr*fCorr;
572 Double_t shift_sys= data_bgr*fCorr_sys - data_bbb;
576 TMath::Sqrt(errData_statbgr_bbb*errData_statbgr_bbb
577 +shift_sys*shift_sys);
588 (
"%3d %5.0f %8.1f +/-%5.1f +/-%5.1f +/-%5.1f (unfolding)",
590 errData_stat_unfold,
TMath::Sqrt(errData_statbgr_unfold*
591 errData_statbgr_unfold-
593 errData_stat_unfold),
595 errData_total_unfold-
596 errData_statbgr_unfold*
597 errData_statbgr_unfold))<<
"\n";
599 (
" %8.1f +/-%5.1f +/-%5.1f +/-%5.1f (bin-by-bin)",
600 data_bbb,errData_stat_bbb,
TMath::Sqrt(errData_statbgr_bbb*
607 errData_statbgr_bbb))
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
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...
Random number generator class based on M.
virtual void SetMaximum(Double_t maximum=-1111)
R__EXTERN TStyle * gStyle
Base class for spline implementation containing the Draw/Paint methods //.
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
virtual void SetMinimum(Double_t minimum=-1111)
virtual void Draw(Option_t *chopt="")
Draw this graph with its current attributes.
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
static void SetDefaultSumw2(Bool_t sumw2=kTRUE)
When this static function is called with sumw2=kTRUE, all new histograms will automatically activate ...
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString...
This is the base class for the ROOT Random number generators.
double pow(double, double)
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
virtual void SetBinError(Int_t bin, Double_t error)
See convention for numbering bins in TH1::GetBin.
virtual Double_t Rndm()
Machine independent random number generator.
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Service class for 2-Dim histogram classes.
virtual void Draw(Option_t *option="")
Draw this histogram with options.
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...
void SetOptFit(Int_t fit=1)
The type of information about fit parameters printed in the histogram statistics box can be selected ...
static const char * GetTUnfoldVersion(void)
tomato 1-D histogram with a double per channel (see TH1 documentation)}
virtual void GetKnot(Int_t i, Double_t &x, Double_t &y) const =0
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
THist< 2, double, THistStatContent, THistStatUncertainty > TH2D
virtual Double_t BreitWigner(Double_t mean=0, Double_t gamma=1)
Return a number distributed following a BreitWigner function with mean and gamma. ...
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)
Automatic pad generation by division.
A Graph is a graphics object made of two arrays X and Y with npoints each.
THist< 1, double, THistStatContent, THistStatUncertainty > TH1D
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content.
Double_t Sqrt(Double_t x)
Int_t Fill(Double_t)
Invalid Fill method.
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
tomato 2-D histogram with a double per channel (see TH1 documentation)}