124 if(
f<parm[0]) itype=1;
125 else if(
f<parm[0]+parm[1]) itype=2;
136 ptgen=
pow(t*(
pow(parm[3],a1)-x0)+x0,1./a1);
138 if(iType) *iType=itype;
139 if(ptGen) *ptGen=ptgen;
143 TMath::Sqrt(parm[7]*parm[7]*ptgen+parm[8]*parm[8]*ptgen*ptgen);
148 triggerParm[2]/(1.+
TMath::Exp((triggerParm[0]-ptObs)/triggerParm[1]));
149 isTriggered= rnd->
Rndm()<triggerProb;
151 }
while((!triggerFlag) && (!isTriggered));
153 if(triggerFlag) *triggerFlag=isTriggered;
212 TH1D *histUnfoldInput=
213 new TH1D(
"unfolding input rec",
";ptrec",nDet,xminDet,xmaxDet);
214 TH2D *histUnfoldMatrix=
215 new TH2D(
"unfolding matrix",
";ptgen;ptrec",
216 nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet);
217 TH1D *histUnfoldBgr1=
218 new TH1D(
"unfolding bgr1 rec",
";ptrec",nDet,xminDet,xmaxDet);
219 TH1D *histUnfoldBgr2=
220 new TH1D(
"unfolding bgr2 rec",
";ptrec",nDet,xminDet,xmaxDet);
221 TH2D *histUnfoldMatrixSys=
222 new TH2D(
"unfolding matrix sys",
";ptgen;ptrec",
223 nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet);
227 new TH1D(
"bbb input rec",
";ptrec",nGen,xminGen,xmaxGen);
228 TH1D *histBbbSignalRec=
229 new TH1D(
"bbb signal rec",
";ptrec",nGen,xminGen,xmaxGen);
231 new TH1D(
"bbb bgr1 rec",
";ptrec",nGen,xminGen,xmaxGen);
233 new TH1D(
"bbb bgr2 rec",
";ptrec",nGen,xminGen,xmaxGen);
235 new TH1D(
"bbb bgrPt rec",
";ptrec",nGen,xminGen,xmaxGen);
236 TH1D *histBbbSignalGen=
237 new TH1D(
"bbb signal gen",
";ptgen",nGen,xminGen,xmaxGen);
238 TH1D *histBbbSignalRecSys=
239 new TH1D(
"bbb signal sys rec",
";ptrec",nGen,xminGen,xmaxGen);
240 TH1D *histBbbBgrPtSys=
241 new TH1D(
"bbb bgrPt sys rec",
";ptrec",nGen,xminGen,xmaxGen);
242 TH1D *histBbbSignalGenSys=
243 new TH1D(
"bbb signal sys gen",
";ptgen",nGen,xminGen,xmaxGen);
247 new TH1D(
"DATA truth gen",
";ptgen",nGen,xminGen,xmaxGen);
249 new TH1D(
"MC prediction rec",
";ptrec",nDet,xminDet,xmaxDet);
265 static Double_t triggerParm_DATA[]={8.0,
271 while(intLumi<lumiData) {
275 Double_t ptObs=GenerateEvent(parm_DATA,triggerParm_DATA,&intLumi,
276 &isTriggered,&ptGen,&iTypeGen);
279 histUnfoldInput->
Fill(ptObs);
282 histBbbInput->
Fill(ptObs);
285 if(iTypeGen==0) histDataTruth->
Fill(ptGen);
305 static Double_t triggerParm_MC[]={8.0,
311 Double_t lumiWeight = lumiData/lumiMC;
313 while(intLumi<lumiMC) {
317 Double_t ptObs=GenerateEvent(parm_MC,triggerParm_MC,&intLumi,&isTriggered,
319 if(!isTriggered) ptObs=0.0;
324 histUnfoldMatrix->
Fill(ptGen,ptObs,lumiWeight);
325 }
else if(iTypeGen==1) {
326 histUnfoldBgr1->
Fill(ptObs,lumiWeight);
327 }
else if(iTypeGen==2) {
328 histUnfoldBgr2->
Fill(ptObs,lumiWeight);
333 if((ptGen>=xminGen)&&(ptGen<xmaxGen)) {
334 histBbbSignalRec->
Fill(ptObs,lumiWeight);
336 histBbbBgrPt->
Fill(ptObs,lumiWeight);
338 histBbbSignalGen->
Fill(ptGen,lumiWeight);
339 }
else if(iTypeGen==1) {
340 histBbbBgr1->
Fill(ptObs,lumiWeight);
341 }
else if(iTypeGen==2) {
342 histBbbBgr2->
Fill(ptObs,lumiWeight);
346 histDetMC ->
Fill(ptObs,lumiWeight);
366 while(intLumi<lumiMC) {
370 Double_t ptObs=GenerateEvent(parm_MC_SYS,triggerParm_MC,&intLumi,
371 &isTriggered,&ptGen,&iTypeGen);
372 if(!isTriggered) ptObs=0.0;
376 histUnfoldMatrixSys->
Fill(ptGen,ptObs);
381 if((ptGen>=xminGen)&&(ptGen<xmaxGen)) {
382 histBbbSignalRecSys->
Fill(ptObs);
384 histBbbBgrPtSys->
Fill(ptObs);
386 histBbbSignalGenSys->
Fill(ptGen);
413 regMode,constraintMode,densityFlags);
416 unfold.SetInput(histUnfoldInput);
422 unfold.SubtractBackground(histUnfoldBgr1,
"background1",scale_bgr,dscale_bgr);
423 unfold.SubtractBackground(histUnfoldBgr2,
"background2",scale_bgr,dscale_bgr);
426 unfold.AddSysError(histUnfoldMatrixSys,
"signalshape_SYS",
436 Int_t iBest=unfold.ScanLcurve(nScan,0.,0.,&lCurve,&logTauX,&logTauY);
438 cout<<
"chi**2="<<unfold.GetChi2A()<<
"+"<<unfold.GetChi2L()
439 <<
" / "<<unfold.GetNdf()<<
"\n";
455 TH1 *histUnfoldOutput=unfold.GetOutput(
"PT(unfold,stat+bgrerr)");
458 TH2 *histEmatStat=unfold.GetEmatrixInput(
"unfolding stat error matrix");
462 TH2 *histEmatTotal=unfold.GetEmatrixTotal(
"unfolding total error matrix");
466 TH1 *histUnfoldStat=
new TH1D(
"PT(unfold,staterr)",
";Pt(gen)",
467 nGen,xminGen,xmaxGen);
468 TH1 *histUnfoldTotal=
new TH1D(
"PT(unfold,totalerr)",
";Pt(gen)",
469 nGen,xminGen,xmaxGen);
470 for(
Int_t i=0;i<nGen+2;i++) {
485 TH2D *histCorr=
new TH2D(
"Corr(total)",
";Pt(gen);Pt(gen)",
486 nGen,xminGen,xmaxGen,nGen,xminGen,xmaxGen);
487 for(
Int_t i=0;i<nGen+2;i++) {
490 if(ei<=0.0)
continue;
491 for(
Int_t j=0;j<nGen+2;j++) {
493 if(ej<=0.0)
continue;
499 TH1 *histDetNormBgr1=unfold.GetBackground(
"bgr1 normalized",
501 TH1 *histDetNormBgrTotal=unfold.GetBackground(
"bgr total normalized");
511 histUnfoldInput->
Draw(
"E");
516 histDetMC->
Draw(
"SAME HIST");
517 histDetNormBgr1->
Draw(
"SAME HIST");
518 histDetNormBgrTotal->
Draw(
"SAME HIST");
525 histUnfoldTotal->
Draw(
"E");
527 histUnfoldOutput->
Draw(
"SAME E1");
529 histUnfoldStat->
Draw(
"SAME E1");
530 histDataTruth->
Draw(
"SAME HIST");
532 histBbbSignalGen->
Draw(
"SAME HIST");
537 histUnfoldMatrix->
Draw(
"BOX");
543 bestLogTauLogChi2->
Draw(
"*");
549 bestLcurve->
Draw(
"*");
553 histCorr->
Draw(
"BOX");
555 output.SaveAs(
"testUnfold3.ps");
560 std::cout<<
"bin truth result (stat) (bgr) (sys)\n";
561 std::cout<<
"===+=====+=========+========+========+=======\n";
562 for(
Int_t i=1;i<=nGen;i++) {
585 Double_t errData_stat_bbb = errData*fCorr;
587 Double_t errData_statbgr_bbb = errData_bgr*fCorr;
593 Double_t shift_sys= data_bgr*fCorr_sys - data_bbb;
597 TMath::Sqrt(errData_statbgr_bbb*errData_statbgr_bbb
598 +shift_sys*shift_sys);
609 (
"%3d %5.0f %8.1f +/-%5.1f +/-%5.1f +/-%5.1f (unfolding)",
611 errData_stat_unfold,
TMath::Sqrt(errData_statbgr_unfold*
612 errData_statbgr_unfold-
614 errData_stat_unfold),
616 errData_total_unfold-
617 errData_statbgr_unfold*
618 errData_statbgr_unfold))<<
"\n";
620 (
" %8.1f +/-%5.1f +/-%5.1f +/-%5.1f (bin-by-bin)",
621 data_bbb,errData_stat_bbb,
TMath::Sqrt(errData_statbgr_bbb*
628 errData_statbgr_bbb))
double pow(double, double)
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.
virtual void Draw(Option_t *chopt="")
Draw this graph with its current attributes.
1-D histogram with a double per channel (see TH1 documentation)}
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.
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 void Draw(Option_t *option="")
Draw this histogram with options.
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-Dim histogram classes.
Int_t Fill(Double_t)
Invalid Fill method.
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.
Random number generator class based on M.
This is the base class for the ROOT Random number generators.
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 Double_t Rndm()
Machine independent random number generator.
Base class for spline implementation containing the Draw/Paint methods.
virtual void Draw(Option_t *option="")
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 Sqrt(Double_t x)
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
static void output(int code)