Logo ROOT  
Reference Guide
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
66using namespace std;
67
68TRandom *rnd=0;
69
70Double_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
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[])
158void testUnfold3()
159{
160 // switch on histogram errors
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 =
402
403 // preserve the area
404 TUnfold::EConstraint constraintMode=
406
407 // bin content is divided by the bin width
408 TUnfoldDensity::EDensityMode densityFlags=
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",
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
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}
#define f(i)
Definition: RSha256.hxx:104
#define c(i)
Definition: RSha256.hxx:101
const Bool_t kFALSE
Definition: RtypesCore.h:90
@ kRed
Definition: Rtypes.h:64
@ kCyan
Definition: Rtypes.h:64
@ kBlue
Definition: Rtypes.h:64
double pow(double, double)
R__EXTERN TStyle * gStyle
Definition: TStyle.h:410
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition: TAttMarker.h:38
The Canvas class.
Definition: TCanvas.h:27
A TGraph is an object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
virtual void Draw(Option_t *chopt="")
Draw this graph with its current attributes.
Definition: TGraph.cxx:760
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:614
The TH1 histogram class.
Definition: TH1.h:56
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8519
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...
Definition: TH1.cxx:8006
virtual void SetMaximum(Double_t maximum=-1111)
Definition: TH1.h:394
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...
Definition: TH1.cxx:8662
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3275
virtual void SetMinimum(Double_t minimum=-1111)
Definition: TH1.h:395
static void SetDefaultSumw2(Bool_t sumw2=kTRUE)
When this static function is called with sumw2=kTRUE, all new histograms will automatically activate ...
Definition: TH1.cxx:6330
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...
Definition: TH1.cxx:8678
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2998
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4907
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:292
Service class for 2-Dim histogram classes.
Definition: TH2.h:30
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH2.cxx:294
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH2.h:88
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content.
Definition: TH2.cxx:2480
Random number generator class based on M.
Definition: TRandom3.h:27
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:27
virtual Double_t BreitWigner(Double_t mean=0, Double_t gamma=1)
Return a number distributed following a BreitWigner function with mean and gamma.
Definition: TRandom.cxx:214
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom.cxx:541
Base class for spline implementation containing the Draw/Paint methods.
Definition: TSpline.h:22
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition: TSpline.cxx:97
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.
Definition: TString.cxx:2311
void SetOptFit(Int_t fit=1)
The type of information about fit parameters printed in the histogram statistics box can be selected ...
Definition: TStyle.cxx:1542
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
Definition: TUnfoldSys.h:104
static const char * GetTUnfoldVersion(void)
Return a string describing the TUnfold version.
Definition: TUnfold.cxx:3680
EConstraint
type of extra constraint
Definition: TUnfold.h:109
@ kEConstraintArea
enforce preservation of the area
Definition: TUnfold.h:115
ERegMode
choice of regularisation scheme
Definition: TUnfold.h:119
@ kRegModeCurvature
regularize the 2nd derivative of the output distribution
Definition: TUnfold.h:131
@ kHistMapOutputHoriz
truth level on x-axis of the response matrix
Definition: TUnfold.h:142
const Double_t sigma
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
Double_t Exp(Double_t x)
Definition: TMath.h:717
Double_t Log(Double_t x)
Definition: TMath.h:750
Double_t Sqrt(Double_t x)
Definition: TMath.h:681
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:725
auto * a
Definition: textangle.C:12
static void output(int code)
Definition: gifencode.c:226