104 if(f<parm[0]) itype=1;
105 else if(f<parm[0]+parm[1]) itype=2;
116 ptgen=
pow(t*(
pow(parm[3],a1)-x0)+x0,1./a1);
118 if(iType) *iType=itype;
119 if(ptGen) *ptGen=ptgen;
123 TMath::Sqrt(parm[7]*parm[7]*ptgen+parm[8]*parm[8]*ptgen*ptgen);
124 ptObs=rnd->BreitWigner(ptgen,sigma);
128 triggerParm[2]/(1.+
TMath::Exp((triggerParm[0]-ptObs)/triggerParm[1]));
129 isTriggered= rnd->Rndm()<triggerProb;
131 }
while((!triggerFlag) && (!isTriggered));
133 if(triggerFlag) *triggerFlag=isTriggered;
192 TH1D *histUnfoldInput=
193 new TH1D(
"unfolding input rec",
";ptrec",nDet,xminDet,xmaxDet);
194 TH2D *histUnfoldMatrix=
195 new TH2D(
"unfolding matrix",
";ptgen;ptrec",
196 nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet);
197 TH1D *histUnfoldBgr1=
198 new TH1D(
"unfolding bgr1 rec",
";ptrec",nDet,xminDet,xmaxDet);
199 TH1D *histUnfoldBgr2=
200 new TH1D(
"unfolding bgr2 rec",
";ptrec",nDet,xminDet,xmaxDet);
201 TH2D *histUnfoldMatrixSys=
202 new TH2D(
"unfolding matrix sys",
";ptgen;ptrec",
203 nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet);
207 new TH1D(
"bbb input rec",
";ptrec",nGen,xminGen,xmaxGen);
208 TH1D *histBbbSignalRec=
209 new TH1D(
"bbb signal rec",
";ptrec",nGen,xminGen,xmaxGen);
211 new TH1D(
"bbb bgr1 rec",
";ptrec",nGen,xminGen,xmaxGen);
213 new TH1D(
"bbb bgr2 rec",
";ptrec",nGen,xminGen,xmaxGen);
215 new TH1D(
"bbb bgrPt rec",
";ptrec",nGen,xminGen,xmaxGen);
216 TH1D *histBbbSignalGen=
217 new TH1D(
"bbb signal gen",
";ptgen",nGen,xminGen,xmaxGen);
218 TH1D *histBbbSignalRecSys=
219 new TH1D(
"bbb signal sys rec",
";ptrec",nGen,xminGen,xmaxGen);
220 TH1D *histBbbBgrPtSys=
221 new TH1D(
"bbb bgrPt sys rec",
";ptrec",nGen,xminGen,xmaxGen);
222 TH1D *histBbbSignalGenSys=
223 new TH1D(
"bbb signal sys gen",
";ptgen",nGen,xminGen,xmaxGen);
227 new TH1D(
"DATA truth gen",
";ptgen",nGen,xminGen,xmaxGen);
229 new TH1D(
"MC prediction rec",
";ptrec",nDet,xminDet,xmaxDet);
245 static Double_t triggerParm_DATA[]={8.0,
251 while(intLumi<lumiData) {
255 Double_t ptObs=GenerateEvent(parm_DATA,triggerParm_DATA,&intLumi,
256 &isTriggered,&ptGen,&iTypeGen);
259 histUnfoldInput->
Fill(ptObs);
262 histBbbInput->
Fill(ptObs);
265 if(iTypeGen==0) histDataTruth->
Fill(ptGen);
285 static Double_t triggerParm_MC[]={8.0,
291 Double_t lumiWeight = lumiData/lumiMC;
293 while(intLumi<lumiMC) {
297 Double_t ptObs=GenerateEvent(parm_MC,triggerParm_MC,&intLumi,&isTriggered,
299 if(!isTriggered) ptObs=0.0;
304 histUnfoldMatrix->
Fill(ptGen,ptObs,lumiWeight);
305 }
else if(iTypeGen==1) {
306 histUnfoldBgr1->
Fill(ptObs,lumiWeight);
307 }
else if(iTypeGen==2) {
308 histUnfoldBgr2->
Fill(ptObs,lumiWeight);
313 if((ptGen>=xminGen)&&(ptGen<xmaxGen)) {
314 histBbbSignalRec->
Fill(ptObs,lumiWeight);
316 histBbbBgrPt->
Fill(ptObs,lumiWeight);
318 histBbbSignalGen->
Fill(ptGen,lumiWeight);
319 }
else if(iTypeGen==1) {
320 histBbbBgr1->
Fill(ptObs,lumiWeight);
321 }
else if(iTypeGen==2) {
322 histBbbBgr2->
Fill(ptObs,lumiWeight);
326 histDetMC ->
Fill(ptObs,lumiWeight);
346 while(intLumi<lumiMC) {
350 Double_t ptObs=GenerateEvent(parm_MC_SYS,triggerParm_MC,&intLumi,
351 &isTriggered,&ptGen,&iTypeGen);
352 if(!isTriggered) ptObs=0.0;
356 histUnfoldMatrixSys->
Fill(ptGen,ptObs);
361 if((ptGen>=xminGen)&&(ptGen<xmaxGen)) {
362 histBbbSignalRecSys->
Fill(ptObs);
364 histBbbBgrPtSys->
Fill(ptObs);
366 histBbbSignalGenSys->
Fill(ptGen);
393 regMode,constraintMode,densityFlags);
396 unfold.SetInput(histUnfoldInput);
402 unfold.SubtractBackground(histUnfoldBgr1,
"background1",scale_bgr,dscale_bgr);
403 unfold.SubtractBackground(histUnfoldBgr2,
"background2",scale_bgr,dscale_bgr);
406 unfold.AddSysError(histUnfoldMatrixSys,
"signalshape_SYS",
416 Int_t iBest=unfold.ScanLcurve(nScan,0.,0.,&lCurve,&logTauX,&logTauY);
418 cout<<
"chi**2="<<unfold.GetChi2A()<<
"+"<<unfold.GetChi2L()
419 <<
" / "<<unfold.GetNdf()<<
"\n";
423 logTauX->
GetKnot(iBest,t[0],x[0]);
424 logTauY->
GetKnot(iBest,t[0],y[0]);
425 TGraph *bestLcurve=
new TGraph(1,x,y);
426 TGraph *bestLogTauLogChi2=
new TGraph(1,t,x);
435 TH1 *histUnfoldOutput=unfold.GetOutput(
"PT(unfold,stat+bgrerr)");
438 TH2 *histEmatStat=unfold.GetEmatrixInput(
"unfolding stat error matrix");
442 TH2 *histEmatTotal=unfold.GetEmatrixTotal(
"unfolding total error matrix");
446 TH1 *histUnfoldStat=
new TH1D(
"PT(unfold,staterr)",
";Pt(gen)",
447 nGen,xminGen,xmaxGen);
448 TH1 *histUnfoldTotal=
new TH1D(
"PT(unfold,totalerr)",
";Pt(gen)",
449 nGen,xminGen,xmaxGen);
450 for(
Int_t i=0;i<nGen+2;i++) {
465 TH2D *histCorr=
new TH2D(
"Corr(total)",
";Pt(gen);Pt(gen)",
466 nGen,xminGen,xmaxGen,nGen,xminGen,xmaxGen);
467 for(
Int_t i=0;i<nGen+2;i++) {
470 if(ei<=0.0)
continue;
471 for(
Int_t j=0;j<nGen+2;j++) {
473 if(ej<=0.0)
continue;
479 TH1 *histDetNormBgr1=unfold.GetBackground(
"bgr1 normalized",
481 TH1 *histDetNormBgrTotal=unfold.GetBackground(
"bgr total normalized");
486 TCanvas *output =
new TCanvas();
491 histUnfoldInput->
Draw(
"E");
496 histDetMC->
Draw(
"SAME HIST");
497 histDetNormBgr1->
Draw(
"SAME HIST");
498 histDetNormBgrTotal->
Draw(
"SAME HIST");
505 histUnfoldTotal->
Draw(
"E");
507 histUnfoldOutput->
Draw(
"SAME E1");
509 histUnfoldStat->
Draw(
"SAME E1");
510 histDataTruth->
Draw(
"SAME HIST");
512 histBbbSignalGen->
Draw(
"SAME HIST");
517 histUnfoldMatrix->
Draw(
"BOX");
522 bestLogTauLogChi2->SetMarkerColor(
kRed);
523 bestLogTauLogChi2->Draw(
"*");
528 bestLcurve->SetMarkerColor(
kRed);
529 bestLcurve->Draw(
"*");
533 histCorr->
Draw(
"BOX");
538 std::cout<<
"bin truth result (stat) (bgr) (sys)\n";
539 std::cout<<
"===+=====+=========+========+========+=======\n";
540 for(
Int_t i=1;i<=nGen;i++) {
563 Double_t errData_stat_bbb = errData*fCorr;
565 Double_t errData_statbgr_bbb = errData_bgr*fCorr;
571 Double_t shift_sys= data_bgr*fCorr_sys - data_bbb;
575 TMath::Sqrt(errData_statbgr_bbb*errData_statbgr_bbb
576 +shift_sys*shift_sys);
587 (
"%3d %5.0f %8.1f +/-%5.1f +/-%5.1f +/-%5.1f (unfolding)",
589 errData_stat_unfold,
TMath::Sqrt(errData_statbgr_unfold*
590 errData_statbgr_unfold-
592 errData_stat_unfold),
594 errData_total_unfold-
595 errData_statbgr_unfold*
596 errData_statbgr_unfold))<<
"\n";
598 (
" %8.1f +/-%5.1f +/-%5.1f +/-%5.1f (bin-by-bin)",
599 data_bbb,errData_stat_bbb,
TMath::Sqrt(errData_statbgr_bbb*
606 errData_statbgr_bbb))
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
static const char * GetTUnfoldVersion(void)
Random number generator class based on M.
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
virtual void SetMaximum(Double_t maximum=-1111)
R__EXTERN TStyle * gStyle
Base class for spline implementation containing the Draw/Paint methods //.
virtual void SetMinimum(Double_t minimum=-1111)
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
static void SetDefaultSumw2(Bool_t sumw2=kTRUE)
static function.
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 SetBinError(Int_t bin, Double_t error)
see convention for numbering bins in TH1::GetBin
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
virtual void SetLineColor(Color_t lcolor)
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 ...
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.
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content.
Double_t Sqrt(Double_t x)
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...
Int_t Fill(Double_t)
Invalid Fill method.
2-D histogram with a double per channel (see TH1 documentation)}