Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
testUnfold3.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_unfold
3 /// \notebook
4 /// Simple Test program for the class TUnfoldDensity.
5 ///
6 /// 1-dimensional unfolding with background subtraction
7 ///
8 /// the goal is to unfold the underlying "true" distribution of a variable Pt
9 ///
10 /// the reconstructed Pt is measured in 24 bins from 4 to 28
11 /// the generator-level Pt is unfolded into 10 bins from 6 to 26
12 /// - plus underflow bin from 0 to 6
13 /// - plus overflow bin above 26
14 /// there are two background sources bgr1 and bgr2
15 /// the signal has a finite trigger efficiency at a threshold of 8 GeV
16 ///
17 /// one type of systematic error is studied, where the signal parameters are
18 /// changed
19 ///
20 /// Finally, the unfolding is compared to a "bin-by-bin" correction method
21 ///
22 /// \macro_output
23 /// \macro_code
24 ///
25 /// **Version 17.6, in parallel to changes in TUnfold**
26 ///
27 /// #### History:
28 /// - Version 17.5, in parallel to changes in TUnfold
29 /// - Version 17.4, in parallel to changes in TUnfold
30 /// - Version 17.3, in parallel to changes in TUnfold
31 /// - Version 17.2, in parallel to changes in TUnfold
32 /// - Version 17.1, in parallel to changes in TUnfold
33 /// - Version 17.0, change to use the class TUnfoldDensity
34 /// - Version 16.1, parallel to changes in TUnfold
35 /// - Version 16.0, parallel to changes in TUnfold
36 /// - Version 15, simple example including background subtraction
37 ///
38 /// This file is part of TUnfold.
39 ///
40 /// TUnfold is free software: you can redistribute it and/or modify
41 /// it under the terms of the GNU General Public License as published by
42 /// the Free Software Foundation, either version 3 of the License, or
43 /// (at your option) any later version.
44 ///
45 /// TUnfold is distributed in the hope that it will be useful,
46 /// but WITHOUT ANY WARRANTY; without even the implied warranty of
47 /// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
48 /// GNU General Public License for more details.
49 ///
50 /// You should have received a copy of the GNU General Public License
51 /// along with TUnfold. If not, see <http://www.gnu.org/licenses/>.
52 ///
53 /// \author Stefan Schmitt DESY, 14.10.2008
54 
55 #include <TMath.h>
56 #include <TCanvas.h>
57 #include <TRandom3.h>
58 #include <TFitter.h>
59 #include <TF1.h>
60 #include <TStyle.h>
61 #include <TVector.h>
62 #include <TGraph.h>
63 
64 #include "TUnfoldDensity.h"
65 
66 using namespace std;
67 
68 TRandom *rnd=0;
69 
70 Double_t GenerateEvent(const Double_t *parm,
71  const Double_t *triggerParm,
72  Double_t *intLumi,
73  Bool_t *triggerFlag,
74  Double_t *ptGen,Int_t *iType)
75 {
76  // generate an event
77  // input:
78  // parameters for the event generator
79  // return value:
80  // reconstructed Pt
81  // output to pointers:
82  // integrated luminosity
83  // several variables only accessible on generator level
84  //
85  // the parm array defines the physical parameters
86  // parm[0]: background source 1 fraction
87  // parm[1]: background source 2 fraction
88  // parm[2]: lower limit of generated Pt distribution
89  // parm[3]: upper limit of generated Pt distribution
90  // parm[4]: exponent for Pt distribution signal
91  // parm[5]: exponent for Pt distribution background 1
92  // parm[6]: exponent for Pt distribution background 2
93  // parm[7]: resolution parameter a goes with sqrt(Pt)
94  // parm[8]: resolution parameter b goes with Pt
95  // triggerParm[0]: trigger threshold turn-on position
96  // triggerParm[1]: trigger threshold turn-on width
97  // triggerParm[2]: trigger efficiency for high pt
98  //
99  // intLumi is advanced bu 1 for each *generated* event
100  // for data, several events may be generated, until one passes the trigger
101  //
102  // some generator-level quantities are also returned:
103  // triggerFlag: whether the event passed the trigger threshold
104  // ptGen: the generated pt
105  // iType: which type of process was simulated
106  //
107  // the "triggerFlag" also has another meaning:
108  // if(triggerFlag==0) only triggered events are returned
109  //
110  // Usage to generate data events
111  // ptObs=GenerateEvent(parm,triggerParm,0,0,0)
112  //
113  // Usage to generate MC events
114  // ptGen=GenerateEvent(parm,triggerParm,&triggerFlag,&ptGen,&iType);
115  //
116  Double_t ptObs;
117  Bool_t isTriggered=kFALSE;
118  do {
119  Int_t itype;
120  Double_t ptgen;
121  Double_t f=rnd->Rndm();
122  // decide whether this is background or signal
123  itype=0;
124  if(f<parm[0]) itype=1;
125  else if(f<parm[0]+parm[1]) itype=2;
126  // generate Pt according to distribution pow(ptgen,a)
127  // get exponent
128  Double_t a=parm[4+itype];
129  Double_t a1=a+1.0;
130  Double_t t=rnd->Rndm();
131  if(a1 == 0.0) {
132  Double_t x0=TMath::Log(parm[2]);
133  ptgen=TMath::Exp(t*(TMath::Log(parm[3])-x0)+x0);
134  } else {
135  Double_t x0=pow(parm[2],a1);
136  ptgen=pow(t*(pow(parm[3],a1)-x0)+x0,1./a1);
137  }
138  if(iType) *iType=itype;
139  if(ptGen) *ptGen=ptgen;
140 
141  // smearing in Pt with large asymmetric tail
142  Double_t sigma=
143  TMath::Sqrt(parm[7]*parm[7]*ptgen+parm[8]*parm[8]*ptgen*ptgen);
144  ptObs=rnd->BreitWigner(ptgen,sigma);
145 
146  // decide whether event was triggered
147  Double_t triggerProb =
148  triggerParm[2]/(1.+TMath::Exp((triggerParm[0]-ptObs)/triggerParm[1]));
149  isTriggered= rnd->Rndm()<triggerProb;
150  (*intLumi) ++;
151  } while((!triggerFlag) && (!isTriggered));
152  // return trigger decision
153  if(triggerFlag) *triggerFlag=isTriggered;
154  return ptObs;
155 }
156 
157 //int main(int argc, char *argv[])
158 void testUnfold3()
159 {
160  // switch on histogram errors
161  TH1::SetDefaultSumw2();
162 
163  // show fit result
164  gStyle->SetOptFit(1111);
165 
166  // random generator
167  rnd=new TRandom3();
168 
169  // data and MC luminosities
170  Double_t const lumiData= 30000;
171  Double_t const lumiMC =1000000;
172 
173  //========================
174  // Step 1: define binning, book histograms
175 
176  // reconstructed pt (fine binning)
177  Int_t const nDet=24;
178  Double_t const xminDet=4.0;
179  Double_t const xmaxDet=28.0;
180 
181  // generated pt (coarse binning)
182  Int_t const nGen=10;
183  Double_t const xminGen= 6.0;
184  Double_t const xmaxGen=26.0;
185 
186  //==================================================================
187  // book histograms
188  // (1) unfolding input: binning scheme is fine for detector,coarse for gen
189  // histUnfoldInput : reconstructed data, binning for unfolding
190  // histUnfoldMatrix : generated vs reconstructed distribution
191  // histUnfoldBgr1 : background source1, as predicted by MC
192  // histUnfoldBgr2 : background source2, as predicted by MC
193  // for systematic studies
194  // histUnfoldMatrixSys : generated vs reconstructed with different signal
195  //
196  // (2) histograms required for bin-by-bin method
197  // histDetDATAbbb : reconstructed data for bin-by-bin
198  // histDetMCbbb : reconstructed MC for bin-by-bin
199  // histDetMCBgr1bbb : reconstructed MC bgr1 for bin-by-bin
200  // histDetMCBgr2bbb : reconstructed MC bgr2 for bin-by-bin
201  // histDetMCBgrPtbbb : reconstructed MC bgr from low/high pt for bin-by-bin
202  // histGenMCbbb : generated MC truth
203  // for systematic studies
204  // histDetMCSysbbb : reconstructed with changed trigger
205  // histGenMCSysbbb : generated MC truth
206  //
207  // (3) monitoring and control
208  // histGenData : data truth for bias tests
209  // histDetMC : MC prediction
210 
211  // (1) create histograms required for unfolding
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);
224 
225  // (2) histograms required for bin-by-bin
226  TH1D *histBbbInput=
227  new TH1D("bbb input rec",";ptrec",nGen,xminGen,xmaxGen);
228  TH1D *histBbbSignalRec=
229  new TH1D("bbb signal rec",";ptrec",nGen,xminGen,xmaxGen);
230  TH1D *histBbbBgr1=
231  new TH1D("bbb bgr1 rec",";ptrec",nGen,xminGen,xmaxGen);
232  TH1D *histBbbBgr2=
233  new TH1D("bbb bgr2 rec",";ptrec",nGen,xminGen,xmaxGen);
234  TH1D *histBbbBgrPt=
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);
244 
245  // (3) control histograms
246  TH1D *histDataTruth=
247  new TH1D("DATA truth gen",";ptgen",nGen,xminGen,xmaxGen);
248  TH1D *histDetMC=
249  new TH1D("MC prediction rec",";ptrec",nDet,xminDet,xmaxDet);
250  // ==============================================================
251  // Step 2: generate data distributions
252  //
253  // data parameters: in real life these are unknown
254  static Double_t parm_DATA[]={
255  0.05, // fraction of background 1 (on generator level)
256  0.05, // fraction of background 2 (on generator level)
257  3.5, // lower Pt cut (generator level)
258  100.,// upper Pt cut (generator level)
259  -3.0,// signal exponent
260  0.1, // background 1 exponent
261  -0.5, // background 2 exponent
262  0.2, // energy resolution a term
263  0.01, // energy resolution b term
264  };
265  static Double_t triggerParm_DATA[]={8.0, // threshold is 8 GeV
266  4.0, // width is 4 GeV
267  0.95 // high Pt efficiency os 95%
268  };
269 
270  Double_t intLumi=0.0;
271  while(intLumi<lumiData) {
272  Int_t iTypeGen;
273  Bool_t isTriggered;
274  Double_t ptGen;
275  Double_t ptObs=GenerateEvent(parm_DATA,triggerParm_DATA,&intLumi,
276  &isTriggered,&ptGen,&iTypeGen);
277  if(isTriggered) {
278  // (1) histogram for unfolding
279  histUnfoldInput->Fill(ptObs);
280 
281  // (2) histogram for bin-by-bin
282  histBbbInput->Fill(ptObs);
283  }
284  // (3) monitoring
285  if(iTypeGen==0) histDataTruth->Fill(ptGen);
286  }
287 
288  // ==============================================================
289  // Step 3: generate default MC distributions
290  //
291  // MC parameters
292  // default settings
293  static Double_t parm_MC[]={
294  0.05, // fraction of background 1 (on generator level)
295  0.05, // fraction of background 2 (on generator level)
296  3.5, // lower Pt cut (generator level)
297  100.,// upper Pt cut (generator level)
298  -4.0,// signal exponent !!! steeper than in data
299  // to illustrate bin-by-bin bias
300  0.1, // background 1 exponent
301  -0.5, // background 2 exponent
302  0.2, // energy resolution a term
303  0.01, // energy resolution b term
304  };
305  static Double_t triggerParm_MC[]={8.0, // threshold is 8 GeV
306  4.0, // width is 4 GeV
307  0.95 // high Pt efficiency is 95%
308  };
309 
310  // weighting factor to accomodate for the different luminosity in data and MC
311  Double_t lumiWeight = lumiData/lumiMC;
312  intLumi=0.0;
313  while(intLumi<lumiMC) {
314  Int_t iTypeGen;
315  Bool_t isTriggered;
316  Double_t ptGen;
317  Double_t ptObs=GenerateEvent(parm_MC,triggerParm_MC,&intLumi,&isTriggered,
318  &ptGen,&iTypeGen);
319  if(!isTriggered) ptObs=0.0;
320 
321  // (1) distribution required for unfolding
322 
323  if(iTypeGen==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);
329  }
330 
331  // (2) distributions for bbb
332  if(iTypeGen==0) {
333  if((ptGen>=xminGen)&&(ptGen<xmaxGen)) {
334  histBbbSignalRec->Fill(ptObs,lumiWeight);
335  } else {
336  histBbbBgrPt->Fill(ptObs,lumiWeight);
337  }
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);
343  }
344 
345  // (3) control distribution
346  histDetMC ->Fill(ptObs,lumiWeight);
347  }
348 
349  // ==============================================================
350  // Step 4: generate MC distributions for systematic study
351  // example: study dependence on initial signal shape
352  // -> BGR distributions do not change
353  static Double_t parm_MC_SYS[]={
354  0.05, // fraction of background: unchanged
355  0.05, // fraction of background: unchanged
356  3.5, // lower Pt cut (generator level)
357  100.,// upper Pt cut (generator level)
358  -2.0, // signal exponent changed
359  0.1, // background 1 exponent
360  -0.5, // background 2 exponent
361  0.2, // energy resolution a term
362  0.01, // energy resolution b term
363  };
364 
365  intLumi=0.0;
366  while(intLumi<lumiMC) {
367  Int_t iTypeGen;
368  Bool_t isTriggered;
369  Double_t ptGen;
370  Double_t ptObs=GenerateEvent(parm_MC_SYS,triggerParm_MC,&intLumi,
371  &isTriggered,&ptGen,&iTypeGen);
372  if(!isTriggered) ptObs=0.0;
373 
374  // (1) for unfolding
375  if(iTypeGen==0) {
376  histUnfoldMatrixSys->Fill(ptGen,ptObs);
377  }
378 
379  // (2) for bin-by-bin
380  if(iTypeGen==0) {
381  if((ptGen>=xminGen)&&(ptGen<xmaxGen)) {
382  histBbbSignalRecSys->Fill(ptObs);
383  } else {
384  histBbbBgrPtSys->Fill(ptObs);
385  }
386  histBbbSignalGenSys->Fill(ptGen);
387  }
388  }
389 
390  // this method is new in version 16 of TUnfold
391  cout<<"TUnfold version is "<<TUnfold::GetTUnfoldVersion()<<"\n";
392 
393  //========================
394  // Step 5: unfolding
395  //
396 
397  // here one has to decide about the regularisation scheme
398 
399  // regularize curvature
400  TUnfold::ERegMode regMode =
401  TUnfold::kRegModeCurvature;
402 
403  // preserve the area
404  TUnfold::EConstraint constraintMode=
405  TUnfold::kEConstraintArea;
406 
407  // bin content is divided by the bin width
408  TUnfoldDensity::EDensityMode densityFlags=
409  TUnfoldDensity::kDensityModeBinWidth;
410 
411  // set up matrix of migrations
412  TUnfoldDensity unfold(histUnfoldMatrix,TUnfold::kHistMapOutputHoriz,
413  regMode,constraintMode,densityFlags);
414 
415  // define the input vector (the measured data distribution)
416  unfold.SetInput(histUnfoldInput);
417 
418  // subtract background, normalized to data luminosity
419  // with 10% scale error each
420  Double_t scale_bgr=1.0;
421  Double_t dscale_bgr=0.1;
422  unfold.SubtractBackground(histUnfoldBgr1,"background1",scale_bgr,dscale_bgr);
423  unfold.SubtractBackground(histUnfoldBgr2,"background2",scale_bgr,dscale_bgr);
424 
425  // add systematic error
426  unfold.AddSysError(histUnfoldMatrixSys,"signalshape_SYS",
427  TUnfold::kHistMapOutputHoriz,
428  TUnfoldSys::kSysErrModeMatrix);
429 
430  // run the unfolding
431  Int_t nScan=30;
432  TSpline *logTauX,*logTauY;
433  TGraph *lCurve;
434  // this method scans the parameter tau and finds the kink in the L curve
435  // finally, the unfolding is done for the best choice of tau
436  Int_t iBest=unfold.ScanLcurve(nScan,0.,0.,&lCurve,&logTauX,&logTauY);
437 
438  cout<<"chi**2="<<unfold.GetChi2A()<<"+"<<unfold.GetChi2L()
439  <<" / "<<unfold.GetNdf()<<"\n";
440 
441  // save graphs with one point to visualize best choice of tau
442  Double_t t[1],x[1],y[1];
443  logTauX->GetKnot(iBest,t[0],x[0]);
444  logTauY->GetKnot(iBest,t[0],y[0]);
445  TGraph *bestLcurve=new TGraph(1,x,y);
446  TGraph *bestLogTauLogChi2=new TGraph(1,t,x);
447 
448 
449  //===========================
450  // Step 6: retrieve unfolding results
451 
452  // get unfolding output
453  // includes the statistical and background errors
454  // but not the other systematic uncertainties
455  TH1 *histUnfoldOutput=unfold.GetOutput("PT(unfold,stat+bgrerr)");
456 
457  // retrieve error matrix of statistical errors
458  TH2 *histEmatStat=unfold.GetEmatrixInput("unfolding stat error matrix");
459 
460  // retrieve full error matrix
461  // This includes all systematic errors
462  TH2 *histEmatTotal=unfold.GetEmatrixTotal("unfolding total error matrix");
463 
464  // create two copies of the unfolded data, one with statistical errors
465  // the other with total errors
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++) {
471  Double_t c=histUnfoldOutput->GetBinContent(i);
472 
473  // histogram with unfolded data and stat errors
474  histUnfoldStat->SetBinContent(i,c);
475  histUnfoldStat->SetBinError
476  (i,TMath::Sqrt(histEmatStat->GetBinContent(i,i)));
477 
478  // histogram with unfolded data and total errors
479  histUnfoldTotal->SetBinContent(i,c);
480  histUnfoldTotal->SetBinError
481  (i,TMath::Sqrt(histEmatTotal->GetBinContent(i,i)));
482  }
483 
484  // create histogram with correlation matrix
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++) {
488  Double_t ei,ej;
489  ei=TMath::Sqrt(histEmatTotal->GetBinContent(i,i));
490  if(ei<=0.0) continue;
491  for(Int_t j=0;j<nGen+2;j++) {
492  ej=TMath::Sqrt(histEmatTotal->GetBinContent(j,j));
493  if(ej<=0.0) continue;
494  histCorr->SetBinContent(i,j,histEmatTotal->GetBinContent(i,j)/ei/ej);
495  }
496  }
497 
498  // retrieve bgr source 1
499  TH1 *histDetNormBgr1=unfold.GetBackground("bgr1 normalized",
500  "background1");
501  TH1 *histDetNormBgrTotal=unfold.GetBackground("bgr total normalized");
502 
503  //========================
504  // Step 7: plots
505 
506  TCanvas output;
507  output.Divide(3,2);
508  output.cd(1);
509  // data, MC prediction, background
510  histUnfoldInput->SetMinimum(0.0);
511  histUnfoldInput->Draw("E");
512  histDetMC->SetMinimum(0.0);
513  histDetMC->SetLineColor(kBlue);
514  histDetNormBgrTotal->SetLineColor(kRed);
515  histDetNormBgr1->SetLineColor(kCyan);
516  histDetMC->Draw("SAME HIST");
517  histDetNormBgr1->Draw("SAME HIST");
518  histDetNormBgrTotal->Draw("SAME HIST");
519 
520  output.cd(2);
521  // unfolded data, data truth, MC truth
522  histUnfoldTotal->SetMinimum(0.0);
523  histUnfoldTotal->SetMaximum(histUnfoldTotal->GetMaximum()*1.5);
524  // outer error: total error
525  histUnfoldTotal->Draw("E");
526  // middle error: stat+bgr
527  histUnfoldOutput->Draw("SAME E1");
528  // inner error: stat only
529  histUnfoldStat->Draw("SAME E1");
530  histDataTruth->Draw("SAME HIST");
531  histBbbSignalGen->SetLineColor(kBlue);
532  histBbbSignalGen->Draw("SAME HIST");
533 
534  output.cd(3);
535  // unfolding matrix
536  histUnfoldMatrix->SetLineColor(kBlue);
537  histUnfoldMatrix->Draw("BOX");
538 
539  // show tau as a function of chi**2
540  output.cd(4);
541  logTauX->Draw();
542  bestLogTauLogChi2->SetMarkerColor(kRed);
543  bestLogTauLogChi2->Draw("*");
544 
545  // show the L curve
546  output.cd(5);
547  lCurve->Draw("AL");
548  bestLcurve->SetMarkerColor(kRed);
549  bestLcurve->Draw("*");
550 
551  // show correlation matrix
552  output.cd(6);
553  histCorr->Draw("BOX");
554 
555  output.SaveAs("testUnfold3.ps");
556 
557  //============================================================
558  // step 8: compare results to the so-called bin-by-bin "correction"
559 
560  std::cout<<"bin truth result (stat) (bgr) (sys)\n";
561  std::cout<<"===+=====+=========+========+========+=======\n";
562  for(Int_t i=1;i<=nGen;i++) {
563  // data contribution in this bin
564  Double_t data=histBbbInput->GetBinContent(i);
565  Double_t errData=histBbbInput->GetBinError(i);
566 
567  // subtract background contributions
568  Double_t data_bgr=data
569  - scale_bgr*(histBbbBgr1->GetBinContent(i)
570  + histBbbBgr2->GetBinContent(i)
571  + histBbbBgrPt->GetBinContent(i));
572  Double_t errData_bgr=TMath::Sqrt
573  (errData*errData+
574  TMath::Power(dscale_bgr*histBbbBgr1->GetBinContent(i),2)+
575  TMath::Power(scale_bgr*histBbbBgr1->GetBinError(i),2)+
576  TMath::Power(dscale_bgr*histBbbBgr2->GetBinContent(i),2)+
577  TMath::Power(scale_bgr*histBbbBgr2->GetBinError(i),2)+
578  TMath::Power(dscale_bgr*histBbbBgrPt->GetBinContent(i),2)+
579  TMath::Power(scale_bgr*histBbbBgrPt->GetBinError(i),2));
580  // "correct" the data, using the Monte Carlo and neglecting off-diagonals
581  Double_t fCorr=(histBbbSignalGen->GetBinContent(i)/
582  histBbbSignalRec->GetBinContent(i));
583  Double_t data_bbb= data_bgr *fCorr;
584  // stat only error
585  Double_t errData_stat_bbb = errData*fCorr;
586  // stat plus background subtraction
587  Double_t errData_statbgr_bbb = errData_bgr*fCorr;
588 
589  // estimate systematic error by repeating the exercise
590  // using the MC with systematic shifts applied
591  Double_t fCorr_sys=(histBbbSignalGenSys->GetBinContent(i)/
592  histBbbSignalRecSys->GetBinContent(i));
593  Double_t shift_sys= data_bgr*fCorr_sys - data_bbb;
594 
595  // add systematic shift quadratically and get total error
596  Double_t errData_total_bbb=
597  TMath::Sqrt(errData_statbgr_bbb*errData_statbgr_bbb
598  +shift_sys*shift_sys);
599 
600  // get results from real unfolding
601  Double_t data_unfold= histUnfoldStat->GetBinContent(i);
602  Double_t errData_stat_unfold=histUnfoldStat->GetBinError(i);
603  Double_t errData_statbgr_unfold=histUnfoldOutput->GetBinError(i);
604  Double_t errData_total_unfold=histUnfoldTotal->GetBinError(i);
605 
606  // compare
607 
608  std::cout<<TString::Format
609  ("%3d %5.0f %8.1f +/-%5.1f +/-%5.1f +/-%5.1f (unfolding)",
610  i,histDataTruth->GetBinContent(i),data_unfold,
611  errData_stat_unfold,TMath::Sqrt(errData_statbgr_unfold*
612  errData_statbgr_unfold-
613  errData_stat_unfold*
614  errData_stat_unfold),
615  TMath::Sqrt(errData_total_unfold*
616  errData_total_unfold-
617  errData_statbgr_unfold*
618  errData_statbgr_unfold))<<"\n";
619  std::cout<<TString::Format
620  (" %8.1f +/-%5.1f +/-%5.1f +/-%5.1f (bin-by-bin)",
621  data_bbb,errData_stat_bbb,TMath::Sqrt(errData_statbgr_bbb*
622  errData_statbgr_bbb-
623  errData_stat_bbb*
624  errData_stat_bbb),
625  TMath::Sqrt(errData_total_bbb*
626  errData_total_bbb-
627  errData_statbgr_bbb*
628  errData_statbgr_bbb))
629  <<"\n\n";
630  }
631 }