Logo ROOT  

Reference Guide
 
Loading...
Searching...
No Matches
testUnfold7c.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_unfold
3/// \notebook
4/// Test program for the classes TUnfoldDensity and TUnfoldBinning.
5///
6/// A toy test of the TUnfold package
7///
8///
9/// This example is documented in conference proceedings:
10///
11/// arXiv:1611.01927
12/// 12th Conference on Quark Confinement and the Hadron Spectrum (Confinement XII)
13///
14/// This is an example of unfolding a one-dimensional distribution. It compares
15/// various unfolding methods:
16///
17/// matrix inversion, template fit, L-curve scan,
18/// iterative unfolding, etc
19///
20/// Further details can be found in talk by S.Schmitt at:
21///
22/// XII Quark Confinement and the Hadron Spectrum
23/// 29.8. - 3.9.2016 Thessaloniki, Greece
24/// statictics session (+proceedings)
25///
26/// The example comprises several macros
27/// - testUnfold7a.C create root files with TTree objects for
28/// signal, background and data
29/// - write files testUnfold7_signal.root
30/// testUnfold7_background.root
31/// testUnfold7_data.root
32///
33/// - testUnfold7b.C loop over trees and fill histograms based on the
34/// TUnfoldBinning objects
35/// - read testUnfold7binning.xml
36/// testUnfold7_signal.root
37/// testUnfold7_background.root
38/// testUnfold7_data.root
39///
40/// - write testUnfold7_histograms.root
41///
42/// - testUnfold7c.C run the unfolding
43/// - read testUnfold7_histograms.root
44/// - write testUnfold7_result.root
45/// - write many histograms, to compare various unfolding methods
46///
47/// \macro_output
48/// \macro_image
49/// \macro_code
50///
51/// **Version 17.6, in parallel to changes in TUnfold**
52///
53/// This file is part of TUnfold.
54///
55/// TUnfold is free software: you can redistribute it and/or modify
56/// it under the terms of the GNU General Public License as published by
57/// the Free Software Foundation, either version 3 of the License, or
58/// (at your option) any later version.
59///
60/// TUnfold is distributed in the hope that it will be useful,
61/// but WITHOUT ANY WARRANTY; without even the implied warranty of
62/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
63/// GNU General Public License for more details.
64///
65/// You should have received a copy of the GNU General Public License
66/// along with TUnfold. If not, see <http://www.gnu.org/licenses/>.
67///
68/// \author Stefan Schmitt DESY, 14.10.2008
69
70#include <iostream>
71#include <cmath>
72#include <map>
73#include <TMath.h>
74#include <TCanvas.h>
75#include <TStyle.h>
76#include <TGraph.h>
77#include <TGraphErrors.h>
78#include <TFile.h>
79#include <TROOT.h>
80#include <TText.h>
81#include <TLine.h>
82#include <TLegend.h>
83#include <TH1.h>
84#include <TF1.h>
85#include <TFitter.h>
86#include <TMatrixD.h>
87#include <TMatrixDSym.h>
88#include <TVectorD.h>
89#include <TMatrixDSymEigen.h>
90#include <TFitResult.h>
91#include <TRandom3.h>
92#include "TUnfoldDensity.h"
93
94using std::vector, std::pair, std::cout;
95
96// #define PRINT_MATRIX_L
97
98#define TEST_INPUT_COVARIANCE
99
100void CreateHistogramCopies(TH1 *h[3],TUnfoldBinning const *binning);
101void CreateHistogramCopies(TH2 *h[3],TUnfoldBinning const *binningX);
102
103TH2 *AddOverflowXY(TH2 *h,double widthX,double widthY);
104TH1 *AddOverflowX(TH1 *h,double width);
105
106void DrawOverflowX(TH1 *h,double posy);
107void DrawOverflowY(TH1 *h,double posx);
108
109
110double const kLegendFontSize=0.05;
111int kNbinC=0;
112
113void DrawPadProbability(TH2 *h);
114void DrawPadEfficiency(TH1 *h);
115void DrawPadReco(TH1 *histMcRec,TH1 *histMcbgrRec,TH1 *histDataRec,
116 TH1 *histDataUnfold,TH2 *histProbability,TH2 *histRhoij);
117void DrawPadTruth(TH1 *histMcsigGen,TH1 *histDataGen,TH1 *histDataUnfold,
118 char const *text=nullptr,double tau=0.0,vector<double> const *r=nullptr,
119 TF1 *f=nullptr);
120void DrawPadCorrelations(TH2 *h,
121 vector<pair<TF1*,vector<double> > > const *table);
122
123TFitResultPtr DoFit(TH1 *h,TH2 *rho,TH1 *truth,char const *text,
124 vector<pair<TF1*,vector<double> > > &table,int niter=0);
125
126void GetNiterGraphs(int iFirst,int iLast,vector<pair<TF1*,
127 vector<double> > > const &table,int color,
128 TGraph *graph[4],int style);
129void GetNiterHist(int ifit,vector<pair<TF1*,vector<double> > > const &table,
130 TH1 *hist[4],int color,int style,int fillStyle);
131
132#ifdef WITH_IDS
133void IDSfirst(TVectorD *data, TVectorD *dataErr, TMatrixD *A_, Double_t lambdaL_, TVectorD* &unfres1IDS_,TVectorD *&soustr);
134
135void IDSiterate(TVectorD *data, TVectorD *dataErr, TMatrixD *A_,TMatrixD *Am_,
136 Double_t lambdaU_, Double_t lambdaM_, Double_t lambdaS_,
137 TVectorD* &unfres2IDS_ ,TVectorD *&soustr);
138#endif
139
140TRandom3 *g_rnd;
141
142void testUnfold7c()
143{
144 gErrorIgnoreLevel=kInfo;
145 // switch on histogram errors
146 TH1::SetDefaultSumw2();
147
148 gStyle->SetOptStat(0);
149
150 g_rnd=new TRandom3(4711);
151
152 //==============================================
153 // step 1 : open output file
154 TFile *outputFile=new TFile("testUnfold7_results.root","recreate");
155
156 //==============================================
157 // step 2 : read binning schemes and input histograms
158 TFile *inputFile=new TFile("testUnfold7_histograms.root");
159
160 outputFile->cd();
161
162 TUnfoldBinning *fineBinning,*coarseBinning;
163
164 inputFile->GetObject("fine",fineBinning);
165 inputFile->GetObject("coarse",coarseBinning);
166
167 if((!fineBinning)||(!coarseBinning)) {
168 cout<<"problem to read binning schemes\n";
169 }
170
171 // save binning schemes to output file
172 fineBinning->Write();
173 coarseBinning->Write();
174
175 // read histograms
176#define READ(TYPE,binning,name) \
177 TYPE *name[3]; inputFile->GetObject(#name,name[0]); \
178 name[0]->Write(); \
179 if(!name[0]) cout<<"Error reading " #name "\n"; \
180 CreateHistogramCopies(name,binning);
181
182 outputFile->cd();
183
184 READ(TH1,fineBinning,histDataRecF);
185 READ(TH1,coarseBinning,histDataRecC);
186 READ(TH1,fineBinning,histDataBgrF);
187 READ(TH1,coarseBinning,histDataBgrC);
188 READ(TH1,coarseBinning,histDataGen);
189
190 READ(TH2,fineBinning,histMcsigGenRecF);
191 READ(TH2,coarseBinning,histMcsigGenRecC);
192 READ(TH1,fineBinning,histMcsigRecF);
193 READ(TH1,coarseBinning,histMcsigRecC);
194 READ(TH1,coarseBinning,histMcsigGen);
195
196 READ(TH1,fineBinning,histMcbgrRecF);
197 READ(TH1,coarseBinning,histMcbgrRecC);
198
199 TH1 *histOutputCtau0[3];
200 TH2 *histRhoCtau0;
201 TH1 *histOutputCLCurve[3];
202 //TH2 *histRhoCLCurve;
203 TH2 *histProbC;
204 double tauMin=1.e-4;
205 double tauMax=1.e-1;
206 double fBgr=1.0; // 0.2/0.25;
207 double biasScale=1.0;
208 TUnfold::ERegMode mode=TUnfold::kRegModeSize; //Derivative;
209
210 //double tauC;
211 {
212 TUnfoldDensity *tunfoldC=
213 new TUnfoldDensity(histMcsigGenRecC[0],
214 TUnfold::kHistMapOutputHoriz,
215 mode,
216 TUnfold::kEConstraintNone,//Area,
217 TUnfoldDensity::kDensityModeNone,
218 coarseBinning,
219 coarseBinning);
220 tunfoldC->SetInput(histDataRecC[0],biasScale);
221 tunfoldC->SubtractBackground(histMcbgrRecC[0],"BGR",fBgr,0.0);
222 tunfoldC->DoUnfold(0.);
223 histOutputCtau0[0]=tunfoldC->GetOutput("histOutputCtau0");
224 histRhoCtau0=tunfoldC->GetRhoIJtotal("histRhoCtau0");
225 CreateHistogramCopies(histOutputCtau0,coarseBinning);
226 tunfoldC->ScanLcurve(50,tauMin,tauMax,nullptr);
227 /* tauC= */tunfoldC->GetTau();
228 //tunfoldC->ScanTau(50,1.E-7,1.E-1,0,TUnfoldDensity::kEScanTauRhoAvg);
229 histOutputCLCurve[0]=tunfoldC->GetOutput("histOutputCLCurve");
230 /* histRhoCLCurve= */tunfoldC->GetRhoIJtotal("histRhoCLCurve");
231 CreateHistogramCopies(histOutputCLCurve,coarseBinning);
232 histProbC=tunfoldC->GetProbabilityMatrix("histProbC",";P_T(gen);P_T(rec)");
233 }
234 TH1 *histOutputFtau0[3];
235 TH2 *histRhoFtau0;
236 TH1 *histOutputFLCurve[3];
237 //TH2 *histRhoFLCurve;
238 TH2 *histProbF;
239 TGraph *lCurve;
240 TSpline *logTauX,*logTauY;
241 tauMin=3.E-4;
242 tauMax=3.E-2;
243 //double tauF;
244 {
245 TUnfoldDensity *tunfoldF=
246 new TUnfoldDensity(histMcsigGenRecF[0],
247 TUnfold::kHistMapOutputHoriz,
248 mode,
249 TUnfold::kEConstraintNone,//Area,
250 TUnfoldDensity::kDensityModeNone,
251 coarseBinning,
252 fineBinning);
253 tunfoldF->SetInput(histDataRecF[0],biasScale);
254 tunfoldF->SubtractBackground(histMcbgrRecF[0],"BGR",fBgr,0.0);
255 tunfoldF->DoUnfold(0.);
256 histOutputFtau0[0]=tunfoldF->GetOutput("histOutputFtau0");
257 histRhoFtau0=tunfoldF->GetRhoIJtotal("histRhoFtau0");
258 CreateHistogramCopies(histOutputFtau0,coarseBinning);
259 tunfoldF->ScanLcurve(50,tauMin,tauMax,nullptr);
260 //tunfoldF->DoUnfold(tauC);
261 /* tauF= */tunfoldF->GetTau();
262 //tunfoldF->ScanTau(50,1.E-7,1.E-1,0,TUnfoldDensity::kEScanTauRhoAvg);
263 histOutputFLCurve[0]=tunfoldF->GetOutput("histOutputFLCurve");
264 /* histRhoFLCurve= */tunfoldF->GetRhoIJtotal("histRhoFLCurve");
265 CreateHistogramCopies(histOutputFLCurve,coarseBinning);
266 histProbF=tunfoldF->GetProbabilityMatrix("histProbF",";P_T(gen);P_T(rec)");
267 }
268 TH1 *histOutputFAtau0[3];
269 TH2 *histRhoFAtau0;
270 TH1 *histOutputFALCurve[3];
271 TH2 *histRhoFALCurve;
272 TH1 *histOutputFArho[3];
273 TH2 *histRhoFArho;
274 TSpline *rhoScan=nullptr;
275 TSpline *logTauCurvature=nullptr;
276
277 double tauFA,tauFArho;
278 {
279 TUnfoldDensity *tunfoldFA=
280 new TUnfoldDensity(histMcsigGenRecF[0],
281 TUnfold::kHistMapOutputHoriz,
282 mode,
283 TUnfold::kEConstraintArea,
284 TUnfoldDensity::kDensityModeNone,
285 coarseBinning,
286 fineBinning);
287 tunfoldFA->SetInput(histDataRecF[0],biasScale);
288 tunfoldFA->SubtractBackground(histMcbgrRecF[0],"BGR",fBgr,0.0);
289 tunfoldFA->DoUnfold(0.);
290 histOutputFAtau0[0]=tunfoldFA->GetOutput("histOutputFAtau0");
291 histRhoFAtau0=tunfoldFA->GetRhoIJtotal("histRhoFAtau0");
292 CreateHistogramCopies(histOutputFAtau0,coarseBinning);
293 tunfoldFA->ScanTau(50,tauMin,tauMax,&rhoScan,TUnfoldDensity::kEScanTauRhoAvg);
294 tauFArho=tunfoldFA->GetTau();
295 histOutputFArho[0]=tunfoldFA->GetOutput("histOutputFArho");
296 histRhoFArho=tunfoldFA->GetRhoIJtotal("histRhoFArho");
297 CreateHistogramCopies(histOutputFArho,coarseBinning);
298
299 tunfoldFA->ScanLcurve(50,tauMin,tauMax,&lCurve,&logTauX,&logTauY,&logTauCurvature);
300 tauFA=tunfoldFA->GetTau();
301 histOutputFALCurve[0]=tunfoldFA->GetOutput("histOutputFALCurve");
302 histRhoFALCurve=tunfoldFA->GetRhoIJtotal("histRhoFALCurve");
303 CreateHistogramCopies(histOutputFALCurve,coarseBinning);
304 }
305 lCurve->Write();
306 logTauX->Write();
307 logTauY->Write();
308
309
310 double widthC=coarseBinning->GetBinSize(histProbC->GetNbinsY()+1);
311 double widthF=fineBinning->GetBinSize(histProbF->GetNbinsY()+1);
312
313 TH2 *histProbCO=AddOverflowXY(histProbC,widthC,widthC);
314 TH2 *histProbFO=AddOverflowXY(histProbF,widthC,widthF);
315
316 // efficiency
317 TH1 *histEfficiencyC=histProbCO->ProjectionX("histEfficiencyC");
318 kNbinC=histProbCO->GetNbinsX();
319
320 // reconstructed quantities with overflow (coarse binning)
321 // MC: add signal and bgr
322 TH1 *histMcsigRecCO=AddOverflowX(histMcsigRecC[2],widthC);
323 TH1 *histMcbgrRecCO=AddOverflowX(histMcbgrRecC[2],widthC);
324 histMcbgrRecCO->Scale(fBgr);
325 TH1 *histMcRecCO=(TH1 *)histMcsigRecCO->Clone("histMcRecC0");
326 histMcRecCO->Add(histMcsigRecCO,histMcbgrRecCO);
327 TH1 *histDataRecCO=AddOverflowX(histDataRecC[2],widthC);
328 //TH1 *histDataRecCNO=AddOverflowX(histDataRecC[1],widthC);
329
330 TH1 *histMcsigRecFO=AddOverflowX(histMcsigRecF[2],widthF);
331 TH1 *histMcbgrRecFO=AddOverflowX(histMcbgrRecF[2],widthF);
332 histMcbgrRecFO->Scale(fBgr);
333 TH1 *histMcRecFO=(TH1 *)histMcsigRecFO->Clone("histMcRecF0");
334 histMcRecFO->Add(histMcsigRecFO,histMcbgrRecFO);
335 TH1 *histDataRecFO=AddOverflowX(histDataRecF[2],widthF);
336
337 // truth level with overflow
338 TH1 *histMcsigGenO=AddOverflowX(histMcsigGen[2],widthC);
339 TH1 *histDataGenO=AddOverflowX(histDataGen[2],widthC);
340
341 // unfolding result with overflow
342 TH1 *histOutputCtau0O=AddOverflowX(histOutputCtau0[2],widthC);
343 TH2 *histRhoCtau0O=AddOverflowXY(histRhoCtau0,widthC,widthC);
344 //TH1 *histOutputCLCurveO=AddOverflowX(histOutputCLCurve[2],widthC);
345 //TH2 *histRhoCLCurveO=AddOverflowXY(histRhoCLCurve,widthC,widthC);
346 TH1 *histOutputFtau0O=AddOverflowX(histOutputFtau0[2],widthC);
347 TH2 *histRhoFtau0O=AddOverflowXY(histRhoFtau0,widthC,widthC);
348 TH1 *histOutputFAtau0O=AddOverflowX(histOutputFAtau0[2],widthC);
349 TH2 *histRhoFAtau0O=AddOverflowXY(histRhoFAtau0,widthC,widthC);
350 //TH1 *histOutputFLCurveO=AddOverflowX(histOutputFLCurve[2],widthC);
351 //TH2 *histRhoFLCurveO=AddOverflowXY(histRhoFLCurve,widthC,widthC);
352 TH1 *histOutputFALCurveO=AddOverflowX(histOutputFALCurve[2],widthC);
353 TH2 *histRhoFALCurveO=AddOverflowXY(histRhoFALCurve,widthC,widthC);
354 TH1 *histOutputFArhoO=AddOverflowX(histOutputFArho[2],widthC);
355 TH2 *histRhoFArhoO=AddOverflowXY(histRhoFArho,widthC,widthC);
356
357 // bin-by-bin
358 TH2 *histRhoBBBO=(TH2 *)histRhoCtau0O->Clone("histRhoBBBO");
359 for(int i=1;i<=histRhoBBBO->GetNbinsX();i++) {
360 for(int j=1;j<=histRhoBBBO->GetNbinsX();j++) {
361 histRhoBBBO->SetBinContent(i,j,(i==j)?1.:0.);
362 }
363 }
364 TH1 *histDataBgrsub=(TH1 *)histDataRecCO->Clone("histDataBgrsub");
365 histDataBgrsub->Add(histMcbgrRecCO,-fBgr);
366 TH1 *histOutputBBBO=(TH1 *)histDataBgrsub->Clone("histOutputBBBO");
367 histOutputBBBO->Divide(histMcsigRecCO);
368 histOutputBBBO->Multiply(histMcsigGenO);
369
370 // iterative
371 int niter=1000;
372 cout<<"maximum number of iterations: "<<niter<<"\n";
373
374 vector <TH1 *>histOutputAgo,histOutputAgorep;
375 vector <TH2 *>histRhoAgo,histRhoAgorep;
376 vector<int> nIter;
377 histOutputAgo.push_back((TH1*)histMcsigGenO->Clone("histOutputAgo-1"));
378 histOutputAgorep.push_back((TH1*)histMcsigGenO->Clone("histOutputAgorep-1"));
379 histRhoAgo.push_back((TH2*)histRhoBBBO->Clone("histRhoAgo-1"));
380 histRhoAgorep.push_back((TH2*)histRhoBBBO->Clone("histRhoAgorep-1"));
381 nIter.push_back(-1);
382
383 int nx=histProbCO->GetNbinsX();
384 int ny=histProbCO->GetNbinsY();
385 TMatrixD covAgo(nx+ny,nx+ny);
386 TMatrixD A(ny,nx);
387 TMatrixD AToverEps(nx,ny);
388 for(int i=0;i<nx;i++) {
389 double epsilonI=0.;
390 for(int j=0;j<ny;j++) {
391 epsilonI+= histProbCO->GetBinContent(i+1,j+1);
392 }
393 for(int j=0;j<ny;j++) {
394 double aji=histProbCO->GetBinContent(i+1,j+1);
395 A(j,i)=aji;
396 AToverEps(i,j)=aji/epsilonI;
397 }
398 }
399 for(int i=0;i<nx;i++) {
400 covAgo(i,i)=TMath::Power
401 (histOutputAgo[0]->GetBinError(i+1)
402 *histOutputAgo[0]->GetXaxis()->GetBinWidth(i+1),2.);
403 }
404 for(int i=0;i<ny;i++) {
405 covAgo(i+nx,i+nx)=TMath::Power
406 (histDataRecCO->GetBinError(i+1)
407 *histDataRecCO->GetXaxis()->GetBinWidth(i+1),2.);
408 }
409#define NREPLICA 300
410 vector<TVectorD *> y(NREPLICA);
411 vector<TVectorD *> yMb(NREPLICA);
412 vector<TVectorD *> yErr(NREPLICA);
413 vector<TVectorD *> x(NREPLICA);
414 TVectorD b(nx);
415 for(int nr=0;nr<NREPLICA;nr++) {
416 x[nr]=new TVectorD(nx);
417 y[nr]=new TVectorD(ny);
418 yMb[nr]=new TVectorD(ny);
419 yErr[nr]=new TVectorD(ny);
420 }
421 for(int i=0;i<nx;i++) {
422 (*x[0])(i)=histOutputAgo[0]->GetBinContent(i+1)
423 *histOutputAgo[0]->GetXaxis()->GetBinWidth(i+1);
424 for(int nr=1;nr<NREPLICA;nr++) {
425 (*x[nr])(i)=(*x[0])(i);
426 }
427 }
428 for(int i=0;i<ny;i++) {
429 (*y[0])(i)=histDataRecCO->GetBinContent(i+1)
430 *histDataRecCO->GetXaxis()->GetBinWidth(i+1);
431 for(int nr=1;nr<NREPLICA;nr++) {
432 (*y[nr])(i)=g_rnd->Poisson((*y[0])(i));
433 (*yErr[nr])(i)=TMath::Sqrt((*y[nr])(i));
434 }
435 b(i)=histMcbgrRecCO->GetBinContent(i+1)*
436 histMcbgrRecCO->GetXaxis()->GetBinWidth(i+1);
437 for(int nr=0;nr<NREPLICA;nr++) {
438 (*yMb[nr])(i)=(*y[nr])(i)-b(i);
439 }
440 }
441 for(int iter=0;iter<=niter;iter++) {
442 if(!(iter %100)) cout<<iter<<"\n";
443 for(int nr=0;nr<NREPLICA;nr++) {
444 TVectorD yrec=A*(*x[nr])+b;
445 TVectorD yOverYrec(ny);
446 for(int j=0;j<ny;j++) {
447 yOverYrec(j)=(*y[nr])(j)/yrec(j);
448 }
449 TVectorD f=AToverEps * yOverYrec;
450 TVectorD xx(nx);
451 for(int i=0;i<nx;i++) {
452 xx(i) = (*x[nr])(i) * f(i);
453 }
454 if(nr==0) {
455 TMatrixD xdf_dr=AToverEps;
456 for(int i=0;i<nx;i++) {
457 for(int j=0;j<ny;j++) {
458 xdf_dr(i,j) *= (*x[nr])(i);
459 }
460 }
461 TMatrixD dr_dxdy(ny,nx+ny);
462 for(int j=0;j<ny;j++) {
463 dr_dxdy(j,nx+j)=1.0/yrec(j);
464 for(int i=0;i<nx;i++) {
465 dr_dxdy(j,i)= -yOverYrec(j)/yrec(j)*A(j,i);
466 }
467 }
468 TMatrixD dxy_dxy(nx+ny,nx+ny);
469 dxy_dxy.SetSub(0,0,xdf_dr*dr_dxdy);
470 for(int i=0;i<nx;i++) {
471 dxy_dxy(i,i) +=f(i);
472 }
473 for(int i=0;i<ny;i++) {
474 dxy_dxy(nx+i,nx+i) +=1.0;
475 }
476 TMatrixD VDT(covAgo,TMatrixD::kMultTranspose,dxy_dxy);
477 covAgo= dxy_dxy*VDT;
478 }
479 (*x[nr])=xx;
480 }
481 if((iter<=25)||
482 ((iter<=100)&&(iter %5==0))||
483 ((iter<=1000)&&(iter %50==0))||
484 (iter %1000==0)) {
485 nIter.push_back(iter);
486 TH1 * h=(TH1*)histOutputAgo[0]->Clone
487 (TString::Format("histOutputAgo%d",iter));
488 histOutputAgo.push_back(h);
489 for(int i=0;i<nx;i++) {
490 double bw=h->GetXaxis()->GetBinWidth(i+1);
491 h->SetBinContent(i+1,(*x[0])(i)/bw);
492 h->SetBinError(i+1,TMath::Sqrt(covAgo(i,i))/bw);
493 }
494 TH2 *h2=(TH2*)histRhoAgo[0]->Clone
495 (TString::Format("histRhoAgo%d",iter));
496 histRhoAgo.push_back(h2);
497 for(int i=0;i<nx;i++) {
498 for(int j=0;j<nx;j++) {
499 double rho= covAgo(i,j)/TMath::Sqrt(covAgo(i,i)*covAgo(j,j));
500 if((i!=j)&&(TMath::Abs(rho)>=1.0)) {
501 cout<<"bad error matrix: iter="<<iter<<"\n";
502 exit(0);
503 }
504 h2->SetBinContent(i+1,j+1,rho);
505 }
506 }
507 // error and correlations from replica analysis
508 h=(TH1*)histOutputAgo[0]->Clone
509 (TString::Format("histOutputAgorep%d",iter));
510 h2=(TH2*)histRhoAgo[0]->Clone
511 (TString::Format("histRhoAgorep%d",iter));
512 histOutputAgorep.push_back(h);
513 histRhoAgorep.push_back(h2);
514
515 TVectorD mean(nx);
516 double w=1./(NREPLICA-1.);
517 for(int nr=1;nr<NREPLICA;nr++) {
518 mean += w* *x[nr];
519 }
520 TMatrixD covAgorep(nx,nx);
521 for(int nr=1;nr<NREPLICA;nr++) {
522 //TMatrixD dx= (*x)-mean;
523 TMatrixD dx(nx,1);
524 for(int i=0;i<nx;i++) {
525 dx(i,0)= (*x[nr])(i)-(*x[0])(i);
526 }
527 covAgorep += w*TMatrixD(dx,TMatrixD::kMultTranspose,dx);
528 }
529
530 for(int i=0;i<nx;i++) {
531 double bw=h->GetXaxis()->GetBinWidth(i+1);
532 h->SetBinContent(i+1,(*x[0])(i)/bw);
533 h->SetBinError(i+1,TMath::Sqrt(covAgorep(i,i))/bw);
534 // cout<<i<<" "<<(*x[0])(i)/bw<<" +/-"<<TMath::Sqrt(covAgorep(i,i))/bw<<" "<<TMath::Sqrt(covAgo(i,i))/bw<<"\n";
535 }
536 for(int i=0;i<nx;i++) {
537 for(int j=0;j<nx;j++) {
538 double rho= covAgorep(i,j)/
539 TMath::Sqrt(covAgorep(i,i)*covAgorep(j,j));
540 if((i!=j)&&(TMath::Abs(rho)>=1.0)) {
541 cout<<"bad error matrix: iter="<<iter<<"\n";
542 exit(0);
543 }
544 h2->SetBinContent(i+1,j+1,rho);
545 }
546 }
547 }
548 }
549
550#ifdef WITH_IDS
551 // IDS Malaescu
552 int niterIDS=100;
553 vector<TVectorD*> unfresIDS(NREPLICA),soustr(NREPLICA);
554 cout<<"IDS number of iterations: "<<niterIDS<<"\n";
555 TMatrixD *Am_IDS[NREPLICA];
556 TMatrixD A_IDS(ny,nx);
557 for(int nr=0;nr<NREPLICA;nr++) {
558 Am_IDS[nr]=new TMatrixD(ny,nx);
559 }
560 for(int iy=0;iy<ny;iy++) {
561 for(int ix=0;ix<nx;ix++) {
562 A_IDS(iy,ix)=histMcsigGenRecC[0]->GetBinContent(ix+1,iy+1);
563 }
564 }
565 double lambdaL=0.;
566 Double_t lambdaUmin = 1.0000002;
567 Double_t lambdaMmin = 0.0000001;
568 Double_t lambdaS = 0.000001;
569 double lambdaU=lambdaUmin;
570 double lambdaM=lambdaMmin;
571 vector<TH1 *> histOutputIDS;
572 vector<TH2 *> histRhoIDS;
573 histOutputIDS.push_back((TH1*)histOutputAgo[0]->Clone("histOutputIDS-1"));
574 histRhoIDS.push_back((TH2*)histRhoAgo[0]->Clone("histRhoIDS-1"));
575 histOutputIDS.push_back((TH1*)histOutputAgo[0]->Clone("histOutputIDS0"));
576 histRhoIDS.push_back((TH2*)histRhoAgo[0]->Clone("histRhoIDS0"));
577 for(int iter=1;iter<=niterIDS;iter++) {
578 if(!(iter %10)) cout<<iter<<"\n";
579
580 for(int nr=0;nr<NREPLICA;nr++) {
581 if(iter==1) {
582 IDSfirst(yMb[nr],yErr[nr],&A_IDS,lambdaL,unfresIDS[nr],soustr[nr]);
583 } else {
584 IDSiterate(yMb[nr],yErr[nr],&A_IDS,Am_IDS[nr],
585 lambdaU,lambdaM,lambdaS,
586 unfresIDS[nr],soustr[nr]);
587 }
588 }
589 unsigned ix;
590 for(ix=0;ix<nIter.size();ix++) {
591 if(nIter[ix]==iter) break;
592 }
593 if(ix<nIter.size()) {
594 TH1 * h=(TH1*)histOutputIDS[0]->Clone
595 (TString::Format("histOutputIDS%d",iter));
596 TH2 *h2=(TH2*)histRhoIDS[0]->Clone
597 (TString::Format("histRhoIDS%d",iter));
598 histOutputIDS.push_back(h);
599 histRhoIDS.push_back(h2);
600 TVectorD mean(nx);
601 double w=1./(NREPLICA-1.);
602 for(int nr=1;nr<NREPLICA;nr++) {
603 mean += w* (*unfresIDS[nr]);
604 }
605 TMatrixD covIDSrep(nx,nx);
606 for(int nr=1;nr<NREPLICA;nr++) {
607 //TMatrixD dx= (*x)-mean;
608 TMatrixD dx(nx,1);
609 for(int i=0;i<nx;i++) {
610 dx(i,0)= (*unfresIDS[nr])(i)-(*unfresIDS[0])(i);
611 }
612 covIDSrep += w*TMatrixD(dx,TMatrixD::kMultTranspose,dx);
613 }
614 for(int i=0;i<nx;i++) {
615 double bw=h->GetXaxis()->GetBinWidth(i+1);
616 h->SetBinContent(i+1,(*unfresIDS[0])(i)/bw/
617 histEfficiencyC->GetBinContent(i+1));
618 h->SetBinError(i+1,TMath::Sqrt(covIDSrep(i,i))/bw/
619 histEfficiencyC->GetBinContent(i+1));
620 // cout<<i<<" "<<(*x[0])(i)/bw<<" +/-"<<TMath::Sqrt(covAgorep(i,i))/bw<<" "<<TMath::Sqrt(covAgo(i,i))/bw<<"\n";
621 }
622 for(int i=0;i<nx;i++) {
623 for(int j=0;j<nx;j++) {
624 double rho= covIDSrep(i,j)/
625 TMath::Sqrt(covIDSrep(i,i)*covIDSrep(j,j));
626 if((i!=j)&&(TMath::Abs(rho)>=1.0)) {
627 cout<<"bad error matrix: iter="<<iter<<"\n";
628 exit(0);
629 }
630 h2->SetBinContent(i+1,j+1,rho);
631 }
632 }
633 }
634 }
635#endif
636
637 //double NEdSmc=histDataBgrsub->GetSumOfWeights();
638
639 vector<pair<TF1 *,vector<double> > > table;
640
641 TCanvas *c1=new TCanvas("c1","",600,600);
642 TCanvas *c2sq=new TCanvas("c2sq","",600,600);
643 c2sq->Divide(1,2);
644 TCanvas *c2w=new TCanvas("c2w","",600,300);
645 c2w->Divide(2,1);
646 TCanvas *c4=new TCanvas("c4","",600,600);
647 c4->Divide(2,2);
648 //TCanvas *c3n=new TCanvas("c3n","",600,600);
649 TPad *subn[3];
650 //gROOT->SetStyle("xTimes2");
651 subn[0]= new TPad("subn0","",0.,0.5,1.,1.);
652 //gROOT->SetStyle("square");
653 subn[1]= new TPad("subn1","",0.,0.,0.5,0.5);
654 subn[2]= new TPad("subn2","",0.5,0.0,1.,0.5);
655 for(int i=0;i<3;i++) {
656 subn[i]->SetFillStyle(0);
657 subn[i]->Draw();
658 }
659 TCanvas *c3c=new TCanvas("c3c","",600,600);
660 TPad *subc[3];
661 //gROOT->SetStyle("xTimes2");
662 subc[0]= new TPad("sub0","",0.,0.5,1.,1.);
663 //gROOT->SetStyle("squareCOLZ");
664 subc[1]= new TPad("sub1","",0.,0.,0.5,0.5);
665 //gROOT->SetStyle("square");
666 subc[2]= new TPad("sub2","",0.5,0.0,1.,0.5);
667 for(int i=0;i<3;i++) {
668 subc[i]->SetFillStyle(0);
669 subc[i]->Draw();
670 }
671
672 //=========================== example ==================================
673
674 c2w->cd(1);
675 DrawPadTruth(histMcsigGenO,histDataGenO,nullptr);
676 c2w->cd(2);
677 DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,nullptr,nullptr,nullptr);
678 c2w->SaveAs("exampleTR.eps");
679
680 //=========================== example ==================================
681
682 c2w->cd(1);
683 DrawPadProbability(histProbCO);
684 c2w->cd(2);
685 DrawPadEfficiency(histEfficiencyC);
686 c2w->SaveAs("exampleAE.eps");
687
688 int iFitInversion=table.size();
689 DoFit(histOutputCtau0O,histRhoCtau0O,histDataGenO,"inversion",table);
690 //=========================== inversion ==================================
691
692 subc[0]->cd();
693 DrawPadTruth(histMcsigGenO,histDataGenO,histOutputCtau0O,"inversion",0.,
694 &table[table.size()-1].second);
695 subc[1]->cd();
696 DrawPadCorrelations(histRhoCtau0O,&table);
697 subc[2]->cd();
698 DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,
699 histOutputCtau0O,histProbCO,histRhoCtau0O);
700 c3c->SaveAs("inversion.eps");
701
702
703 DoFit(histOutputFtau0O,histRhoFtau0O,histDataGenO,"template",table);
704 //=========================== template ==================================
705
706 subc[0]->cd();
707 DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFtau0O,"fit",0.,
708 &table[table.size()-1].second);
709 subc[1]->cd();
710 DrawPadCorrelations(histRhoFtau0O,&table);
711 subc[2]->cd();
712 DrawPadReco(histMcRecFO,histMcbgrRecFO,histDataRecFO,
713 histOutputFtau0O,histProbFO,histRhoFtau0O);
714 c3c->SaveAs("template.eps");
715
716 DoFit(histOutputFAtau0O,histRhoFAtau0O,histDataGenO,"template+area",table);
717 //=========================== template+area ==================================
718
719 subc[0]->cd();
720 DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFAtau0O,"fit",0.,
721 &table[table.size()-1].second);
722 subc[1]->cd();
723 DrawPadCorrelations(histRhoFAtau0O,&table);
724 subc[2]->cd();
725 DrawPadReco(histMcRecFO,histMcbgrRecFO,histDataRecFO,
726 histOutputFAtau0O,histProbFO,histRhoFAtau0O);
727 c3c->SaveAs("templateA.eps");
728
729 int iFitFALCurve=table.size();
730 DoFit(histOutputFALCurveO,histRhoFALCurveO,histDataGenO,"Tikhonov+area",table);
731 //=========================== template+area+tikhonov =====================
732 subc[0]->cd();
733 DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFALCurveO,"Tikhonov",tauFA,
734 &table[table.size()-1].second);
735 subc[1]->cd();
736 DrawPadCorrelations(histRhoFALCurveO,&table);
737 subc[2]->cd();
738 DrawPadReco(histMcRecFO,histMcbgrRecFO,histDataRecFO,
739 histOutputFALCurveO,histProbFO,histRhoFALCurveO);
740 c3c->SaveAs("lcurveFA.eps");
741
742 int iFitFArho=table.size();
743 DoFit(histOutputFArhoO,histRhoFArhoO,histDataGenO,"min(rhomax)",table);
744 //=========================== template+area+tikhonov =====================
745 subc[0]->cd();
746 DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFArhoO,"Tikhonov",tauFArho,
747 &table[table.size()-1].second);
748 subc[1]->cd();
749 DrawPadCorrelations(histRhoFArho,&table);
750 subc[2]->cd();
751 DrawPadReco(histMcRecFO,histMcbgrRecFO,histDataRecFO,
752 histOutputFArhoO,histProbFO,histRhoFArhoO);
753 c3c->SaveAs("rhoscanFA.eps");
754
755 int iFitBinByBin=table.size();
756 DoFit(histOutputBBBO,histRhoBBBO,histDataGenO,"bin-by-bin",table);
757 //=========================== bin-by-bin =================================
758 //c->cd(1);
759 //DrawPadProbability(histProbCO);
760 //c->cd(2);
761 //DrawPadCorrelations(histRhoBBBO,&table);
762 c2sq->cd(1);
763 DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,
764 histOutputBBBO,histProbCO,histRhoBBBO);
765 c2sq->cd(2);
766 DrawPadTruth(histMcsigGenO,histDataGenO,histOutputBBBO,"bin-by-bin",0.,
767 &table[table.size()-1].second);
768 c2sq->SaveAs("binbybin.eps");
769
770
771 //=========================== iterative ===================================
772 int iAgoFirstFit=table.size();
773 for(size_t i=1;i<histRhoAgorep.size();i++) {
774 int n=nIter[i];
775 bool isFitted=false;
776 DoFit(histOutputAgorep[i],histRhoAgorep[i],histDataGenO,
777 TString::Format("iterative, N=%d",n),table,n);
778 isFitted=true;
779 subc[0]->cd();
780 DrawPadTruth(histMcsigGenO,histDataGenO,histOutputAgorep[i],
781 TString::Format("iterative N=%d",nIter[i]),0.,
782 isFitted ? &table[table.size()-1].second : nullptr);
783 subc[1]->cd();
784 DrawPadCorrelations(histRhoAgorep[i],&table);
785 subc[2]->cd();
786 DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,
787 histOutputAgorep[i],histProbCO,histRhoAgorep[i]);
788 c3c->SaveAs(TString::Format("iterative%d.eps",nIter[i]));
789 }
790 int iAgoLastFit=table.size();
791
792#ifdef WITH_IDS
793 int iIDSFirstFit=table.size();
794
795 //=========================== IDS ===================================
796
797 for(size_t i=2;i<histRhoIDS.size();i++) {
798 int n=nIter[i];
799 bool isFitted=false;
800 DoFit(histOutputIDS[i],histRhoIDS[i],histDataGenO,
801 TString::Format("IDS, N=%d",n),table,n);
802 isFitted=true;
803 subc[0]->cd();
804 DrawPadTruth(histMcsigGenO,histDataGenO,histOutputIDS[i],
805 TString::Format("IDS N=%d",nIter[i]),0.,
806 isFitted ? &table[table.size()-1].second : 0);
807 subc[1]->cd();
808 DrawPadCorrelations(histRhoIDS[i],&table);
809 subc[2]->cd();
810 DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,
811 histOutputIDS[i],histProbCO,histRhoIDS[i]);
812 c3c->SaveAs(TString::Format("ids%d.eps",nIter[i]));
813 }
814 int iIDSLastFit=table.size();
815#endif
816
817 int nfit=table.size();
818 TH1D *fitChindf=new TH1D("fitChindf",";algorithm;#chi^{2}/NDF",nfit,0,nfit);
819 TH1D *fitNorm=new TH1D("fitNorm",";algorithm;Landau amplitude [1/GeV]",nfit,0,nfit);
820 TH1D *fitMu=new TH1D("fitMu",";algorithm;Landau #mu [GeV]",nfit,0,nfit);
821 TH1D *fitSigma=new TH1D("fitSigma",";algorithm;Landau #sigma [GeV]",nfit,0,nfit);
822 for(int fit=0;fit<nfit;fit++) {
823 TF1 *f=table[fit].first;
824 vector<double> const &r=table[fit].second;
825 fitChindf->GetXaxis()->SetBinLabel(fit+1,f->GetName());
826 fitNorm->GetXaxis()->SetBinLabel(fit+1,f->GetName());
827 fitMu->GetXaxis()->SetBinLabel(fit+1,f->GetName());
828 fitSigma->GetXaxis()->SetBinLabel(fit+1,f->GetName());
829 double chi2=r[0];
830 double ndf=r[1];
831 fitChindf->SetBinContent(fit+1,chi2/ndf);
832 fitChindf->SetBinError(fit+1,TMath::Sqrt(2./ndf));
833 fitNorm->SetBinContent(fit+1,f->GetParameter(0));
834 fitNorm->SetBinError(fit+1,f->GetParError(0));
835 fitMu->SetBinContent(fit+1,f->GetParameter(1));
836 fitMu->SetBinError(fit+1,f->GetParError(1));
837 fitSigma->SetBinContent(fit+1,f->GetParameter(2));
838 fitSigma->SetBinError(fit+1,f->GetParError(2));
839 cout<<"\""<<f->GetName()<<"\","<<r[2]/r[3]<<","<<r[3]
840 <<","<<TMath::Prob(r[2],r[3])
841 <<","<<chi2/ndf
842 <<","<<ndf
843 <<","<<TMath::Prob(r[0],r[1])
844 <<","<<r[5];
845 for(int i=1;i<3;i++) {
846 cout<<","<<f->GetParameter(i)<<",\"\302\261\","<<f->GetParError(i);
847 }
848 cout<<"\n";
849 }
850
851 //=========================== L-curve ==========================
852 c4->cd(1);
853 lCurve->SetTitle("L curve;log_{10} L_{x};log_{10} L_{y}");
854 lCurve->SetLineColor(kRed);
855 lCurve->Draw("AL");
856 c4->cd(2);
857 gPad->Clear();
858 c4->cd(3);
859 logTauX->SetTitle(";log_{10} #tau;log_{10} L_{x}");
860 logTauX->SetLineColor(kBlue);
861 logTauX->Draw();
862 c4->cd(4);
863 logTauY->SetTitle(";log_{10} #tau;log_{10} L_{y}");
864 logTauY->SetLineColor(kBlue);
865 logTauY->Draw();
866 c4->SaveAs("lcurveL.eps");
867
868 //========================= rho and L-curve scan ===============
869 c4->cd(1);
870 logTauCurvature->SetTitle(";log_{10}(#tau);L curve curvature");
871 logTauCurvature->SetLineColor(kRed);
872 logTauCurvature->Draw();
873 c4->cd(2);
874 rhoScan->SetTitle(";log_{10}(#tau);average(#rho_{i})");
875 rhoScan->SetLineColor(kRed);
876 rhoScan->Draw();
877 c4->cd(3);
878 DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFALCurveO,"Tikhonov",tauFA,
879 &table[iFitFALCurve].second);
880 c4->cd(4);
881 DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFArhoO,"Tikhonov",tauFArho,
882 &table[iFitFArho].second);
883
884 c4->SaveAs("scanTau.eps");
885
886
887 TGraph *graphNiterAgoBay[4];
888 GetNiterGraphs(iAgoFirstFit,iAgoFirstFit+1,table,kRed-2,graphNiterAgoBay,20);
889 TGraph *graphNiterAgo[4];
890 GetNiterGraphs(iAgoFirstFit,iAgoLastFit,table,kRed,graphNiterAgo,24);
891#ifdef WITH_IDS
892 TGraph *graphNiterIDS[4];
893 GetNiterGraphs(iIDSFirstFit,iIDSLastFit,table,kMagenta,graphNiterIDS,21);
894#endif
895 TH1 *histNiterInversion[4];
896 GetNiterHist(iFitInversion,table,histNiterInversion,kCyan,1,1001);
897 TH1 *histNiterFALCurve[4];
898 GetNiterHist(iFitFALCurve,table,histNiterFALCurve,kBlue,2,3353);
899 TH1 *histNiterFArho[4];
900 GetNiterHist(iFitFArho,table,histNiterFArho,kAzure-4,3,3353);
901 TH1 *histNiterBinByBin[4];
902 GetNiterHist(iFitBinByBin,table,histNiterBinByBin,kOrange+1,3,3335);
903
904 histNiterInversion[0]->GetYaxis()->SetRangeUser(0.3,500.);
905 histNiterInversion[1]->GetYaxis()->SetRangeUser(-0.1,0.9);
906 histNiterInversion[2]->GetYaxis()->SetRangeUser(5.6,6.3);
907 histNiterInversion[3]->GetYaxis()->SetRangeUser(1.6,2.4);
908
909 TLine *line=nullptr;
910 c1->cd();
911 for(int i=0;i<2;i++) {
912 gPad->Clear();
913 gPad->SetLogx();
914 gPad->SetLogy((i<1));
915 if(! histNiterInversion[i]) continue;
916 histNiterInversion[i]->Draw("][");
917 histNiterFALCurve[i]->Draw("SAME ][");
918 histNiterFArho[i]->Draw("SAME ][");
919 histNiterBinByBin[i]->Draw("SAME ][");
920 graphNiterAgo[i]->Draw("LP");
921 graphNiterAgoBay[i]->SetMarkerStyle(20);
922 graphNiterAgoBay[i]->Draw("P");
923#ifdef WITH_IDS
924 graphNiterIDS[i]->Draw("LP");
925#endif
926 TLegend *legend;
927 if(i==1) {
928 legend=new TLegend(0.48,0.28,0.87,0.63);
929 } else {
930 legend=new TLegend(0.45,0.5,0.88,0.88);
931 }
932 legend->SetBorderSize(0);
933 legend->SetFillStyle(1001);
934 legend->SetFillColor(kWhite);
935 legend->SetTextSize(kLegendFontSize*0.7);
936 legend->AddEntry( histNiterInversion[0],"inversion","l");
937 legend->AddEntry( histNiterFALCurve[0],"Tikhonov L-curve","l");
938 legend->AddEntry( histNiterFArho[0],"Tikhonov global cor.","l");
939 legend->AddEntry( histNiterBinByBin[0],"bin-by-bin","l");
940 legend->AddEntry( graphNiterAgoBay[0],"\"Bayesian\"","p");
941 legend->AddEntry( graphNiterAgo[0],"iterative","p");
942#ifdef WITH_IDS
943 legend->AddEntry( graphNiterIDS[0],"IDS","p");
944#endif
945 legend->Draw();
946
947 c1->SaveAs(TString::Format("niter%d.eps",i));
948 }
949
950 //c4->cd(1);
951 //DrawPadCorrelations(histRhoFALCurveO,&table);
952 c2sq->cd(1);
953 DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFALCurveO,"Tikhonov",tauFA,
954 &table[iFitFALCurve].second,table[iFitFALCurve].first);
955
956 c2sq->cd(2);
957 gPad->SetLogx();
958 gPad->SetLogy(false);
959 histNiterInversion[3]->DrawClone("E2");
960 histNiterInversion[3]->SetFillStyle(0);
961 histNiterInversion[3]->Draw("SAME HIST ][");
962 histNiterFALCurve[3]->DrawClone("SAME E2");
963 histNiterFALCurve[3]->SetFillStyle(0);
964 histNiterFALCurve[3]->Draw("SAME HIST ][");
965 histNiterFArho[3]->DrawClone("SAME E2");
966 histNiterFArho[3]->SetFillStyle(0);
967 histNiterFArho[3]->Draw("SAME HIST ][");
968 histNiterBinByBin[3]->DrawClone("SAME E2");
969 histNiterBinByBin[3]->SetFillStyle(0);
970 histNiterBinByBin[3]->Draw("SAME HIST ][");
971 double yTrue=1.8;
972 line=new TLine(histNiterInversion[3]->GetXaxis()->GetXmin(),
973 yTrue,
974 histNiterInversion[3]->GetXaxis()->GetXmax(),
975 yTrue);
976 line->SetLineWidth(3);
977 line->Draw();
978
979 graphNiterAgo[3]->Draw("LP");
980 graphNiterAgoBay[3]->SetMarkerStyle(20);
981 graphNiterAgoBay[3]->Draw("P");
982#ifdef WITH_IDS
983 graphNiterIDS[3]->Draw("LP");
984#endif
985
986 TLegend *legend;
987 legend=new TLegend(0.55,0.53,0.95,0.97);
988 legend->SetBorderSize(0);
989 legend->SetFillStyle(1001);
990 legend->SetFillColor(kWhite);
991 legend->SetTextSize(kLegendFontSize);
992 legend->AddEntry( line,"truth","l");
993 legend->AddEntry( histNiterInversion[3],"inversion","l");
994 legend->AddEntry( histNiterFALCurve[3],"Tikhonov L-curve","l");
995 legend->AddEntry( histNiterFArho[3],"Tikhonov global cor.","l");
996 legend->AddEntry( histNiterBinByBin[3],"bin-by-bin","l");
997 legend->AddEntry( graphNiterAgoBay[3],"\"Bayesian\"","p");
998 legend->AddEntry( graphNiterAgo[3],"iterative","p");
999#ifdef WITH_IDS
1000 legend->AddEntry( graphNiterIDS[3],"IDS","p");
1001#endif
1002 legend->Draw();
1003
1004 c2sq->SaveAs("fitSigma.eps");
1005
1006 outputFile->Write();
1007
1008 delete outputFile;
1009
1010}
1011
1012void GetNiterGraphs(int iFirst,int iLast,vector<pair<TF1*,
1013 vector<double> > > const &table,int color,
1014 TGraph *graph[4],int style) {
1015 TVectorD niter(iLast-iFirst);
1016 TVectorD eniter(iLast-iFirst);
1017 TVectorD chi2(iLast-iFirst);
1018 TVectorD gcor(iLast-iFirst);
1019 TVectorD mean(iLast-iFirst);
1020 TVectorD emean(iLast-iFirst);
1021 TVectorD sigma(iLast-iFirst);
1022 TVectorD esigma(iLast-iFirst);
1023 for(int ifit=iFirst;ifit<iLast;ifit++) {
1024 vector<double> const &r=table[ifit].second;
1025 niter(ifit-iFirst)=r[4];
1026 chi2(ifit-iFirst)=r[2]/r[3];
1027 gcor(ifit-iFirst)=r[5];
1028 TF1 const *f=table[ifit].first;
1029 mean(ifit-iFirst)=f->GetParameter(1);
1030 emean(ifit-iFirst)=f->GetParError(1);
1031 sigma(ifit-iFirst)=f->GetParameter(2);
1032 esigma(ifit-iFirst)=f->GetParError(2);
1033 }
1034 graph[0]=new TGraph(niter,chi2);
1035 graph[1]=new TGraph(niter,gcor);
1036 graph[2]=new TGraphErrors(niter,mean,eniter,emean);
1037 graph[3]=new TGraphErrors(niter,sigma,eniter,esigma);
1038 for(int g=0;g<4;g++) {
1039 if(graph[g]) {
1040 graph[g]->SetLineColor(color);
1041 graph[g]->SetMarkerColor(color);
1042 graph[g]->SetMarkerStyle(style);
1043 }
1044 }
1045}
1046
1047void GetNiterHist(int ifit,vector<pair<TF1*,vector<double> > > const &table,
1048 TH1 *hist[4],int color,int style,int fillStyle) {
1049 vector<double> const &r=table[ifit].second;
1050 TF1 const *f=table[ifit].first;
1051 hist[0]=new TH1D(table[ifit].first->GetName()+TString("_chi2"),
1052 ";iteration;unfold-truth #chi^{2}/N_{D.F.}",1,0.2,1500.);
1053 hist[0]->SetBinContent(1,r[2]/r[3]);
1054
1055 hist[1]=new TH1D(table[ifit].first->GetName()+TString("_gcor"),
1056 ";iteration;avg(#rho_{i})",1,0.2,1500.);
1057 hist[1]->SetBinContent(1,r[5]);
1058 hist[2]=new TH1D(table[ifit].first->GetName()+TString("_mu"),
1059 ";iteration;parameter #mu",1,0.2,1500.);
1060 hist[2]->SetBinContent(1,f->GetParameter(1));
1061 hist[2]->SetBinError(1,f->GetParError(1));
1062 hist[3]=new TH1D(table[ifit].first->GetName()+TString("_sigma"),
1063 ";iteration;parameter #sigma",1,0.2,1500.);
1064 hist[3]->SetBinContent(1,f->GetParameter(2));
1065 hist[3]->SetBinError(1,f->GetParError(2));
1066 for(int h=0;h<4;h++) {
1067 if(hist[h]) {
1068 hist[h]->SetLineColor(color);
1069 hist[h]->SetLineStyle(style);
1070 if( hist[h]->GetBinError(1)>0.0) {
1071 hist[h]->SetFillColor(color-10);
1072 hist[h]->SetFillStyle(fillStyle);
1073 }
1074 hist[h]->SetMarkerStyle(0);
1075 }
1076 }
1077}
1078
1079void CreateHistogramCopies(TH1 *h[3],TUnfoldBinning const *binning) {
1080 TString baseName(h[0]->GetName());
1081 Int_t *binMap;
1082 h[1]=binning->CreateHistogram(baseName+"_axis",kTRUE,&binMap);
1083 h[2]=(TH1 *)h[1]->Clone(baseName+"_binw");
1084 Int_t nMax=binning->GetEndBin()+1;
1085 for(Int_t iSrc=0;iSrc<nMax;iSrc++) {
1086 Int_t iDest=binMap[iSrc];
1087 double c=h[0]->GetBinContent(iSrc)+h[1]->GetBinContent(iDest);
1088 double e=TMath::Hypot(h[0]->GetBinError(iSrc),h[1]->GetBinError(iDest));
1089 h[1]->SetBinContent(iDest,c);
1090 h[1]->SetBinError(iDest,e);
1091 h[2]->SetBinContent(iDest,c);
1092 h[2]->SetBinError(iDest,e);
1093 }
1094 for(int iDest=0;iDest<=h[2]->GetNbinsX()+1;iDest++) {
1095 double c=h[2]->GetBinContent(iDest);
1096 double e=h[2]->GetBinError(iDest);
1097 double bw=binning->GetBinSize(iDest);
1098 /* if(bw!=h[2]->GetBinWidth(iDest)) {
1099 cout<<"bin "<<iDest<<" width="<<bw<<" "<<h[2]->GetBinWidth(iDest)
1100 <<"\n";
1101 } */
1102 if(bw>0.0) {
1103 h[2]->SetBinContent(iDest,c/bw);
1104 h[2]->SetBinError(iDest,e/bw);
1105 } else {
1106 }
1107 }
1108}
1109
1110void CreateHistogramCopies(TH2 *h[3],TUnfoldBinning const *binningX) {
1111 h[1]=nullptr;
1112 h[2]=nullptr;
1113}
1114
1115TH2 *AddOverflowXY(TH2 *h,double widthX,double widthY) {
1116 // add overflow bin to X-axis
1117 int nx=h->GetNbinsX();
1118 int ny=h->GetNbinsY();
1119 double *xBins=new double[nx+2];
1120 double *yBins=new double[ny+2];
1121 for(int i=1;i<=nx;i++) {
1122 xBins[i-1]=h->GetXaxis()->GetBinLowEdge(i);
1123 }
1124 xBins[nx]=h->GetXaxis()->GetBinUpEdge(nx);
1125 xBins[nx+1]=xBins[nx]+widthX;
1126 for(int i=1;i<=ny;i++) {
1127 yBins[i-1]=h->GetYaxis()->GetBinLowEdge(i);
1128 }
1129 yBins[ny]=h->GetYaxis()->GetBinUpEdge(ny);
1130 yBins[ny+1]=yBins[ny]+widthY;
1131 TString name(h->GetName());
1132 name+="U";
1133 TH2 *r=new TH2D(name,h->GetTitle(),nx+1,xBins,ny+1,yBins);
1134 for(int ix=0;ix<=nx+1;ix++) {
1135 for(int iy=0;iy<=ny+1;iy++) {
1136 r->SetBinContent(ix,iy,h->GetBinContent(ix,iy));
1137 r->SetBinError(ix,iy,h->GetBinError(ix,iy));
1138 }
1139 }
1140 delete [] yBins;
1141 delete [] xBins;
1142 return r;
1143}
1144
1145TH1 *AddOverflowX(TH1 *h,double widthX) {
1146 // add overflow bin to X-axis
1147 int nx=h->GetNbinsX();
1148 double *xBins=new double[nx+2];
1149 for(int i=1;i<=nx;i++) {
1150 xBins[i-1]=h->GetXaxis()->GetBinLowEdge(i);
1151 }
1152 xBins[nx]=h->GetXaxis()->GetBinUpEdge(nx);
1153 xBins[nx+1]=xBins[nx]+widthX;
1154 TString name(h->GetName());
1155 name+="U";
1156 TH1 *r=new TH1D(name,h->GetTitle(),nx+1,xBins);
1157 for(int ix=0;ix<=nx+1;ix++) {
1158 r->SetBinContent(ix,h->GetBinContent(ix));
1159 r->SetBinError(ix,h->GetBinError(ix));
1160 }
1161 delete [] xBins;
1162 return r;
1163}
1164
1165void DrawOverflowX(TH1 *h,double posy) {
1166 double x1=h->GetXaxis()->GetBinLowEdge(h->GetNbinsX());
1167 double x2=h->GetXaxis()->GetBinUpEdge(h->GetNbinsX());
1168 double y0=h->GetYaxis()->GetBinLowEdge(1);
1169 double y2=h->GetYaxis()->GetBinUpEdge(h->GetNbinsY());;
1170 if(h->GetDimension()==1) {
1171 y0=h->GetMinimum();
1172 y2=h->GetMaximum();
1173 }
1174 double w1=-0.3;
1175 TText *textX=new TText((1.+w1)*x2-w1*x1,(1.-posy)*y0+posy*y2,"Overflow bin");
1176 textX->SetNDC(kFALSE);
1177 textX->SetTextSize(0.05);
1178 textX->SetTextAngle(90.);
1179 textX->Draw();
1180 TLine *lineX=new TLine(x1,y0,x1,y2);
1181 lineX->Draw();
1182}
1183
1184void DrawOverflowY(TH1 *h,double posx) {
1185 double x0=h->GetXaxis()->GetBinLowEdge(1);
1186 double x2=h->GetXaxis()->GetBinUpEdge(h->GetNbinsX());
1187 double y1=h->GetYaxis()->GetBinLowEdge(h->GetNbinsY());;
1188 double y2=h->GetYaxis()->GetBinUpEdge(h->GetNbinsY());;
1189 double w1=-0.3;
1190 TText *textY=new TText((1.-posx)*x0+posx*x2,(1.+w1)*y1-w1*y2,"Overflow bin");
1191 textY->SetNDC(kFALSE);
1192 textY->SetTextSize(0.05);
1193 textY->Draw();
1194 TLine *lineY=new TLine(x0,y1,x2,y1);
1195 lineY->Draw();
1196}
1197
1198void DrawPadProbability(TH2 *h) {
1199 h->Draw("COLZ");
1200 h->SetTitle("migration probabilities;P_{T}(gen) [GeV];P_{T}(rec) [GeV]");
1201 DrawOverflowX(h,0.05);
1202 DrawOverflowY(h,0.35);
1203}
1204
1205void DrawPadEfficiency(TH1 *h) {
1206 h->SetTitle("efficiency;P_{T}(gen) [GeV];#epsilon");
1207 h->SetLineColor(kBlue);
1208 h->SetMinimum(0.75);
1209 h->SetMaximum(1.0);
1210 h->Draw();
1211 DrawOverflowX(h,0.05);
1212 TLegend *legEfficiency=new TLegend(0.3,0.58,0.6,0.75);
1213 legEfficiency->SetBorderSize(0);
1214 legEfficiency->SetFillStyle(0);
1215 legEfficiency->SetTextSize(kLegendFontSize);
1216 legEfficiency->AddEntry(h,"reconstruction","l");
1217 legEfficiency->AddEntry((TObject*)nullptr," efficiency","");
1218 legEfficiency->Draw();
1219}
1220
1221void DrawPadReco(TH1 *histMcRec,TH1 *histMcbgrRec,TH1 *histDataRec,
1222 TH1 *histDataUnfold,TH2 *histProbability,TH2 *histRhoij) {
1223 //gPad->SetLogy(kTRUE);
1224 double amax=0.0;
1225 for(int i=1;i<=histMcRec->GetNbinsX();i++) {
1226 amax=TMath::Max(amax,histMcRec->GetBinContent(i)
1227 +5.0*histMcRec->GetBinError(i));
1228 amax=TMath::Max(amax,histDataRec->GetBinContent(i)
1229 +2.0*histDataRec->GetBinError(i));
1230 }
1231 histMcRec->SetTitle("Reconstructed;P_{T}(rec);Nevent / GeV");
1232 histMcRec->Draw("HIST");
1233 histMcRec->SetLineColor(kBlue);
1234 histMcRec->SetMinimum(1.0);
1235 histMcRec->SetMaximum(amax);
1236 //histMcbgrRec->SetFillMode(1);
1237 histMcbgrRec->SetLineColor(kBlue-6);
1238 histMcbgrRec->SetFillColor(kBlue-10);
1239 histMcbgrRec->Draw("SAME HIST");
1240
1241 TH1 * histFoldBack=nullptr;
1242 if(histDataUnfold && histProbability && histRhoij) {
1243 histFoldBack=(TH1 *)
1244 histMcRec->Clone(histDataUnfold->GetName()+TString("_folded"));
1245 int nrec=histFoldBack->GetNbinsX();
1246 if((nrec==histProbability->GetNbinsY())&&
1247 (nrec==histMcbgrRec->GetNbinsX())&&
1248 (nrec==histDataRec->GetNbinsX())
1249 ) {
1250 for(int ix=1;ix<=nrec;ix++) {
1251 double sum=0.0;
1252 double sume2=0.0;
1253 for(int iy=0;iy<=histProbability->GetNbinsX()+1;iy++) {
1254 sum += histDataUnfold->GetBinContent(iy)*
1255 histDataUnfold->GetBinWidth(iy)*
1256 histProbability->GetBinContent(iy,ix);
1257 for(int iy2=0;iy2<=histProbability->GetNbinsX()+1;iy2++) {
1258 sume2 += histDataUnfold->GetBinError(iy)*
1259 histDataUnfold->GetBinWidth(iy)*
1260 histProbability->GetBinContent(iy,ix)*
1261 histDataUnfold->GetBinError(iy2)*
1262 histDataUnfold->GetBinWidth(iy2)*
1263 histProbability->GetBinContent(iy2,ix)*
1264 histRhoij->GetBinContent(iy,iy2);
1265 }
1266 }
1267 sum /= histFoldBack->GetBinWidth(ix);
1268 sum += histMcbgrRec->GetBinContent(ix);
1269 histFoldBack->SetBinContent(ix,sum);
1270 histFoldBack->SetBinError(ix,TMath::Sqrt(sume2)
1271 /histFoldBack->GetBinWidth(ix));
1272 }
1273 } else {
1274 cout<<"can not fold back: "<<nrec
1275 <<" "<<histProbability->GetNbinsY()
1276 <<" "<<histMcbgrRec->GetNbinsX()
1277 <<" "<<histDataRec->GetNbinsX()
1278 <<"\n";
1279 exit(0);
1280 }
1281
1282 histFoldBack->SetLineColor(kBlack);
1283 histFoldBack->SetMarkerStyle(0);
1284 histFoldBack->Draw("SAME HIST");
1285 }
1286
1287 histDataRec->SetLineColor(kRed);
1288 histDataRec->SetMarkerColor(kRed);
1289 histDataRec->Draw("SAME");
1290 DrawOverflowX(histMcRec,0.5);
1291
1292 TLegend *legRec=new TLegend(0.4,0.5,0.68,0.85);
1293 legRec->SetBorderSize(0);
1294 legRec->SetFillStyle(0);
1295 legRec->SetTextSize(kLegendFontSize);
1296 legRec->AddEntry(histMcRec,"MC total","l");
1297 legRec->AddEntry(histMcbgrRec,"background","f");
1298 if(histFoldBack) {
1299 int ndf=-kNbinC;
1300 double sumD=0.,sumF=0.,chi2=0.;
1301 for(int i=1;i<=histDataRec->GetNbinsX();i++) {
1302 //cout<<histDataRec->GetBinContent(i)<<" "<<histFoldBack->GetBinContent(i)<<" "<<" w="<<histFoldBack->GetBinWidth(i)<<"\n";
1303 sumD+=histDataRec->GetBinContent(i)*histDataRec->GetBinWidth(i);
1304 sumF+=histFoldBack->GetBinContent(i)*histFoldBack->GetBinWidth(i);
1305 double pull=(histFoldBack->GetBinContent(i)-histDataRec->GetBinContent(i))/histDataRec->GetBinError(i);
1306 chi2+= pull*pull;
1307 ndf+=1;
1308 }
1309 legRec->AddEntry(histDataRec,TString::Format("data N_{evt}=%.0f",sumD),"lp");
1310 legRec->AddEntry(histFoldBack,TString::Format("folded N_{evt}=%.0f",sumF),"l");
1311 legRec->AddEntry((TObject*)nullptr,TString::Format("#chi^{2}=%.1f ndf=%d",chi2,ndf),"");
1312 //exit(0);
1313 } else {
1314 legRec->AddEntry(histDataRec,"data","lp");
1315 }
1316 legRec->Draw();
1317}
1318
1319void DrawPadTruth(TH1 *histMcsigGen,TH1 *histDataGen,TH1 *histDataUnfold,
1320 char const *text,double tau,vector<double> const *r,
1321 TF1 *f) {
1322 //gPad->SetLogy(kTRUE);
1323 double amin=0.;
1324 double amax=0.;
1325 for(int i=1;i<=histMcsigGen->GetNbinsX();i++) {
1326 if(histDataUnfold) {
1327 amin=TMath::Min(amin,histDataUnfold->GetBinContent(i)
1328 -1.1*histDataUnfold->GetBinError(i));
1329 amax=TMath::Max(amax,histDataUnfold->GetBinContent(i)
1330 +1.1*histDataUnfold->GetBinError(i));
1331 }
1332 amin=TMath::Min(amin,histMcsigGen->GetBinContent(i)
1333 -histMcsigGen->GetBinError(i));
1334 amin=TMath::Min(amin,histDataGen->GetBinContent(i)
1335 -histDataGen->GetBinError(i));
1336 amax=TMath::Max(amax,histMcsigGen->GetBinContent(i)
1337 +10.*histMcsigGen->GetBinError(i));
1338 amax=TMath::Max(amax,histDataGen->GetBinContent(i)
1339 +2.*histDataGen->GetBinError(i));
1340 }
1341 histMcsigGen->SetMinimum(amin);
1342 histMcsigGen->SetMaximum(amax);
1343
1344 histMcsigGen->SetTitle("Truth;P_{T};Nevent / GeV");
1345 histMcsigGen->SetLineColor(kBlue);
1346 histMcsigGen->Draw("HIST");
1347 histDataGen->SetLineColor(kRed);
1348 histDataGen->SetMarkerColor(kRed);
1349 histDataGen->SetMarkerSize(1.0);
1350 histDataGen->Draw("SAME HIST");
1351 if(histDataUnfold) {
1352 histDataUnfold->SetMarkerStyle(21);
1353 histDataUnfold->SetMarkerSize(0.7);
1354 histDataUnfold->Draw("SAME");
1355 }
1356 DrawOverflowX(histMcsigGen,0.5);
1357
1358 if(f) {
1359 f->SetLineStyle(1);
1360 f->Draw("SAME");
1361 }
1362
1363 TLegend *legTruth=new TLegend(0.32,0.65,0.6,0.9);
1364 legTruth->SetBorderSize(0);
1365 legTruth->SetFillStyle(0);
1366 legTruth->SetTextSize(kLegendFontSize);
1367 legTruth->AddEntry(histMcsigGen,"MC","l");
1368 if(!histDataUnfold) legTruth->AddEntry((TObject *)nullptr," Landau(5,2)","");
1369 legTruth->AddEntry(histDataGen,"data","l");
1370 if(!histDataUnfold) legTruth->AddEntry((TObject *)nullptr," Landau(6,1.8)","");
1371 if(histDataUnfold) {
1372 TString t;
1373 if(text) t=text;
1374 else t=histDataUnfold->GetName();
1375 if(tau>0) {
1376 t+=TString::Format(" #tau=%.2g",tau);
1377 }
1378 legTruth->AddEntry(histDataUnfold,t,"lp");
1379 if(r) {
1380 legTruth->AddEntry((TObject *)nullptr,"test wrt data:","");
1381 legTruth->AddEntry((TObject *)nullptr,TString::Format
1382 ("#chi^{2}/%d=%.1f prob=%.3f",
1383 (int)(*r)[3],(*r)[2]/(*r)[3],
1384 TMath::Prob((*r)[2],(*r)[3])),"");
1385 }
1386 }
1387 if(f) {
1388 legTruth->AddEntry(f,"fit","l");
1389 }
1390 legTruth->Draw();
1391 if(histDataUnfold ) {
1392 TPad *subpad = new TPad("subpad","",0.35,0.29,0.88,0.68);
1393 subpad->SetFillStyle(0);
1394 subpad->Draw();
1395 subpad->cd();
1396 amin=0.;
1397 amax=0.;
1398 int istart=11;
1399 for(int i=istart;i<=histMcsigGen->GetNbinsX();i++) {
1400 amin=TMath::Min(amin,histMcsigGen->GetBinContent(i)
1401 -histMcsigGen->GetBinError(i));
1402 amin=TMath::Min(amin,histDataGen->GetBinContent(i)
1403 -histDataGen->GetBinError(i));
1404 amin=TMath::Min(amin,histDataUnfold->GetBinContent(i)
1405 -histDataUnfold->GetBinError(i));
1406 amax=TMath::Max(amax,histMcsigGen->GetBinContent(i)
1407 +histMcsigGen->GetBinError(i));
1408 amax=TMath::Max(amax,histDataGen->GetBinContent(i)
1409 +histDataGen->GetBinError(i));
1410 amax=TMath::Max(amax,histDataUnfold->GetBinContent(i)
1411 +histDataUnfold->GetBinError(i));
1412 }
1413 TH1 *copyMcsigGen=(TH1*)histMcsigGen->Clone();
1414 TH1 *copyDataGen=(TH1*)histDataGen->Clone();
1415 TH1 *copyDataUnfold=(TH1*)histDataUnfold->Clone();
1416 copyMcsigGen->GetXaxis()->SetRangeUser
1417 (copyMcsigGen->GetXaxis()->GetBinLowEdge(istart),
1418 copyMcsigGen->GetXaxis()->GetBinUpEdge(copyMcsigGen->GetNbinsX()-1));
1419 copyMcsigGen->SetTitle(";;");
1420 copyMcsigGen->GetYaxis()->SetRangeUser(amin,amax);
1421 copyMcsigGen->Draw("HIST");
1422 copyDataGen->Draw("SAME HIST");
1423 copyDataUnfold->Draw("SAME");
1424 if(f) {
1425 ((TF1 *)f->Clone())->Draw("SAME");
1426 }
1427 }
1428}
1429
1430void DrawPadCorrelations(TH2 *h,
1431 vector<pair<TF1*,vector<double> > > const *table) {
1432 h->SetMinimum(-1.);
1433 h->SetMaximum(1.);
1434 h->SetTitle("correlation coefficients;P_{T}(gen) [GeV];P_{T}(gen) [GeV]");
1435 h->Draw("COLZ");
1436 DrawOverflowX(h,0.05);
1437 DrawOverflowY(h,0.05);
1438 if(table) {
1439 TLegend *legGCor=new TLegend(0.13,0.6,0.5,0.8);
1440 legGCor->SetBorderSize(0);
1441 legGCor->SetFillStyle(0);
1442 legGCor->SetTextSize(kLegendFontSize);
1443 vector<double> const &r=(*table)[table->size()-1].second;
1444 legGCor->AddEntry((TObject *)nullptr,TString::Format("min(#rho_{ij})=%5.2f",r[6]),"");
1445 legGCor->AddEntry((TObject *)nullptr,TString::Format("max(#rho_{ij})=%5.2f",r[7]),"");
1446 legGCor->AddEntry((TObject *)nullptr,TString::Format("avg(#rho_i)=%5.2f",r[5]),"");
1447 legGCor->Draw();
1448 }
1449}
1450
1451TH1 *g_fcnHist=nullptr;
1452TMatrixD *g_fcnMatrix=nullptr;
1453
1454void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag) {
1455 if(flag==0) {
1456 cout<<"fcn flag=0: npar="<<npar<<" gin="<<gin<<" par=[";
1457 for(int i=0;i<npar;i++) {
1458 cout<<" "<<u[i];
1459 }
1460 cout<<"]\n";
1461 }
1462 int n=g_fcnMatrix->GetNrows();
1463 TVectorD dy(n);
1464 double x0=0,y0=0.;
1465 for(int i=0;i<=n;i++) {
1466 double x1;
1467 if(i<1) x1=g_fcnHist->GetXaxis()->GetBinLowEdge(i+1);
1468 else x1=g_fcnHist->GetXaxis()->GetBinUpEdge(i);
1469 double y1=TMath::LandauI((x1-u[1])/u[2]);
1470 if(i>0) {
1471 double iy=u[0]*u[2]*(y1-y0)/(x1-x0);
1472 dy(i-1)=iy-g_fcnHist->GetBinContent(i);
1473 //cout<<"i="<<i<<" iy="<<iy<<" delta="<< dy(i-1)<<"\n";
1474 }
1475 x0=x1;
1476 y0=y1;
1477 //cout<<"i="<<i<<" y1="<<y1<<" x1="<<x1<<"\n";
1478 }
1479 TVectorD Hdy=(*g_fcnMatrix) * dy;
1480 //Hdy.Print();
1481 f=Hdy*dy;
1482 //exit(0);
1483}
1484
1485
1486TFitResultPtr DoFit(TH1 *h,TH2 *rho,TH1 *truth,const char *text,
1487 vector<pair<TF1 *,vector<double> > > &table,int niter) {
1488 TString option="IESN";
1489 cout<<h->GetName()<<"\n";
1490 double gcorAvg=0.;
1491 double rhoMin=0.;
1492 double rhoMax=0.;
1493 if(rho) {
1494 g_fcnHist=h;
1495 int n=h->GetNbinsX()-1; // overflow is included as extra bin, exclude in fit
1496 TMatrixDSym v(n);
1497 //g_fcnMatrix=new TMatrixD(n,n);
1498 for(int i=0;i<n;i++) {
1499 for(int j=0;j<n;j++) {
1500 v(i,j)=rho->GetBinContent(i+1,j+1)*
1501 (h->GetBinError(i+1)*h->GetBinError(j+1));
1502 }
1503 }
1504 TMatrixDSymEigen ev(v);
1505 TMatrixD d(n,n);
1506 TVectorD di(ev.GetEigenValues());
1507 for(int i=0;i<n;i++) {
1508 if(di(i)>0.0) {
1509 d(i,i)=1./di(i);
1510 } else {
1511 cout<<"bad eigenvalue i="<<i<<" di="<<di(i)<<"\n";
1512 exit(0);
1513 }
1514 }
1515 TMatrixD O(ev.GetEigenVectors());
1516 TMatrixD DOT(d,TMatrixD::kMultTranspose,O);
1517 g_fcnMatrix=new TMatrixD(O,TMatrixD::kMult,DOT);
1518 TMatrixD test(*g_fcnMatrix,TMatrixD::kMult,v);
1519 int error=0;
1520 for(int i=0;i<n;i++) {
1521 if(TMath::Abs(test(i,i)-1.0)>1.E-7) {
1522 error++;
1523 }
1524 for(int j=0;j<n;j++) {
1525 if(i==j) continue;
1526 if(TMath::Abs(test(i,j)>1.E-7)) error++;
1527 }
1528 }
1529 // calculate global correlation coefficient (all bins)
1530 TMatrixDSym v1(n+1);
1531 rhoMin=1.;
1532 rhoMax=-1.;
1533 for(int i=0;i<=n;i++) {
1534 for(int j=0;j<=n;j++) {
1535 double rho_ij=rho->GetBinContent(i+1,j+1);
1536 v1(i,j)=rho_ij*
1537 (h->GetBinError(i+1)*h->GetBinError(j+1));
1538 if(i!=j) {
1539 if(rho_ij<rhoMin) rhoMin=rho_ij;
1540 if(rho_ij>rhoMax) rhoMax=rho_ij;
1541 }
1542 }
1543 }
1544 TMatrixDSymEigen ev1(v1);
1545 TMatrixD d1(n+1,n+1);
1546 TVectorD di1(ev1.GetEigenValues());
1547 for(int i=0;i<=n;i++) {
1548 if(di1(i)>0.0) {
1549 d1(i,i)=1./di1(i);
1550 } else {
1551 cout<<"bad eigenvalue i="<<i<<" di1="<<di1(i)<<"\n";
1552 exit(0);
1553 }
1554 }
1555 TMatrixD O1(ev1.GetEigenVectors());
1556 TMatrixD DOT1(d1,TMatrixD::kMultTranspose,O1);
1557 TMatrixD vinv1(O1,TMatrixD::kMult,DOT1);
1558 for(int i=0;i<=n;i++) {
1559 double gcor2=1.-1./(vinv1(i,i)*v1(i,i));
1560 if(gcor2>=0.0) {
1561 double gcor=TMath::Sqrt(gcor2);
1562 gcorAvg += gcor;
1563 } else {
1564 cout<<"bad global correlation "<<i<<" "<<gcor2<<"\n";
1565 }
1566 }
1567 gcorAvg /=(n+1);
1568 /* if(error) {
1569 v.Print();
1570 g_fcnMatrix->Print();
1571 exit(0);
1572 } */
1573 //g_fcnMatrix->Invert();
1574 //from: HFitImpl.cxx
1575 // TVirtualFitter::FCNFunc_t userFcn = 0;
1576 // typedef void (* FCNFunc_t )(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
1577 // userFcn = (TVirtualFitter::GetFitter())->GetFCN();
1578 // (TVirtualFitter::GetFitter())->SetUserFunc(f1);
1579 //...
1580 //fitok = fitter->FitFCN( userFcn );
1581
1582 TVirtualFitter::Fitter(h)->SetFCN(fcn);
1583 option += "U";
1584
1585 }
1586 double xmax=h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()-1);
1587 TF1 *landau=new TF1(text,"[0]*TMath::Landau(x,[1],[2],0)",
1588 0.,xmax);
1589 landau->SetParameter(0,6000.);
1590 landau->SetParameter(1,5.);
1591 landau->SetParameter(2,2.);
1592 landau->SetParError(0,10.);
1593 landau->SetParError(1,0.5);
1594 landau->SetParError(2,0.1);
1595 TFitResultPtr s=h->Fit(landau,option,nullptr,0.,xmax);
1596 vector<double> r(8);
1597 int np=landau->GetNpar();
1598 fcn(np,nullptr,r[0],landau->GetParameters(),0);
1599 r[1]=h->GetNbinsX()-1-landau->GetNpar();
1600 for(int i=0;i<h->GetNbinsX()-1;i++) {
1601 double di=h->GetBinContent(i+1)-truth->GetBinContent(i+1);
1602 if(g_fcnMatrix) {
1603 for(int j=0;j<h->GetNbinsX()-1;j++) {
1604 double dj=h->GetBinContent(j+1)-truth->GetBinContent(j+1);
1605 r[2]+=di*dj*(*g_fcnMatrix)(i,j);
1606 }
1607 } else {
1608 double pull=di/h->GetBinError(i+1);
1609 r[2]+=pull*pull;
1610 }
1611 r[3]+=1.0;
1612 }
1613 r[4]=niter;
1614 if(!niter) r[4]=0.25;
1615 r[5]=gcorAvg;
1616 r[6]=rhoMin;
1617 r[7]=rhoMax;
1618 if(rho) {
1619 g_fcnHist=nullptr;
1620 delete g_fcnMatrix;
1621 g_fcnMatrix=nullptr;
1622 }
1623 table.push_back(make_pair(landau,r));
1624 return s;
1625}
1626
1627#ifdef WITH_IDS
1628
1629//===================== interface to IDS unfolding code follows here
1630// contact Bogdan Malescu to find it
1631
1632#include "ids_code.cc"
1633
1634void IDSfirst(TVectorD *data, TVectorD *dataErr, TMatrixD *A_, Double_t lambdaL_, TVectorD* &unfres1IDS_,TVectorD *&soustr){
1635
1636 int N_=data->GetNrows();
1637 soustr = new TVectorD(N_);
1638 for( Int_t i=0; i<N_; i++ ){ (*soustr)[i] = 0.; }
1639 unfres1IDS_ = Unfold( data, dataErr, A_, N_, lambdaL_, soustr );
1640}
1641
1642void IDSiterate(TVectorD *data, TVectorD *dataErr, TMatrixD *A_, TMatrixD *Am_, Double_t lambdaU_, Double_t lambdaM_, Double_t lambdaS_,TVectorD* &unfres2IDS_ ,TVectorD *&soustr) {
1643
1644 int N_=data->GetNrows();
1645 ModifyMatrix( Am_, A_, unfres2IDS_, dataErr, N_, lambdaM_, soustr, lambdaS_ );
1646 delete unfres2IDS_;
1647 unfres2IDS_ = Unfold( data, dataErr, Am_, N_, lambdaU_, soustr );
1648}
1649
1650#endif
d
#define d(i)
Definition RSha256.hxx:102
b
#define b(i)
Definition RSha256.hxx:100
f
#define f(i)
Definition RSha256.hxx:104
c
#define c(i)
Definition RSha256.hxx:101
g
#define g(i)
Definition RSha256.hxx:105
h
#define h(i)
Definition RSha256.hxx:106
e
#define e(i)
Definition RSha256.hxx:103
Int_t
int Int_t
Definition RtypesCore.h:45
kFALSE
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
Double_t
double Double_t
Definition RtypesCore.h:59
kTRUE
constexpr Bool_t kTRUE
Definition RtypesCore.h:100
kRed
@ kRed
Definition Rtypes.h:66
kOrange
@ kOrange
Definition Rtypes.h:67
kBlack
@ kBlack
Definition Rtypes.h:65
kMagenta
@ kMagenta
Definition Rtypes.h:66
kWhite
@ kWhite
Definition Rtypes.h:65
kCyan
@ kCyan
Definition Rtypes.h:66
kBlue
@ kBlue
Definition Rtypes.h:66
kAzure
@ kAzure
Definition Rtypes.h:67
kInfo
constexpr Int_t kInfo
Definition TError.h:44
gErrorIgnoreLevel
Int_t gErrorIgnoreLevel
Error handling routines.
Definition TError.cxx:31
option
Option_t Option_t option
Definition TGWin32VirtualXProxy.cxx:44
np
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t np
Definition TGWin32VirtualXProxy.cxx:222
data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Definition TGWin32VirtualXProxy.cxx:104
r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Definition TGWin32VirtualXProxy.cxx:168
x2
Option_t Option_t TPoint TPoint const char x2
Definition TGWin32VirtualXProxy.cxx:70
x1
Option_t Option_t TPoint TPoint const char x1
Definition TGWin32VirtualXProxy.cxx:70
mode
Option_t Option_t TPoint TPoint const char mode
Definition TGWin32VirtualXProxy.cxx:68
y2
Option_t Option_t TPoint TPoint const char y2
Definition TGWin32VirtualXProxy.cxx:70
width
Option_t Option_t width
Definition TGWin32VirtualXProxy.cxx:56
style
Option_t Option_t style
Definition TGWin32VirtualXProxy.cxx:46
text
Option_t Option_t TPoint TPoint const char text
Definition TGWin32VirtualXProxy.cxx:68
y1
Option_t Option_t TPoint TPoint const char y1
Definition TGWin32VirtualXProxy.cxx:70
name
char name[80]
Definition TGX11.cxx:110
xmax
float xmax
Definition THbookFile.cxx:95
TMatrixD
TMatrixT< Double_t > TMatrixD
Definition TMatrixDfwd.h:23
gStyle
R__EXTERN TStyle * gStyle
Definition TStyle.h:433
TVectorD
TVectorT< Double_t > TVectorD
Definition TVectorDfwd.h:23
gPad
#define gPad
Definition TVirtualPad.h:305
TAttFill::SetFillColor
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition TAttFill.h:37
TAttFill::SetFillStyle
virtual void SetFillStyle(Style_t fstyle)
Set the fill area style.
Definition TAttFill.h:39
TAttLine::SetLineStyle
virtual void SetLineStyle(Style_t lstyle)
Set the line style.
Definition TAttLine.h:42
TAttLine::SetLineWidth
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
Definition TAttLine.h:43
TAttLine::SetLineColor
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition TAttLine.h:40
TAttMarker::SetMarkerColor
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition TAttMarker.h:38
TAttMarker::SetMarkerStyle
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition TAttMarker.h:40
TAttMarker::SetMarkerSize
virtual void SetMarkerSize(Size_t msize=1)
Set the marker size.
Definition TAttMarker.h:45
TAttText::SetTextAngle
virtual void SetTextAngle(Float_t tangle=0)
Set the text angle.
Definition TAttText.h:43
TAttText::SetTextSize
virtual void SetTextSize(Float_t tsize=1)
Set the text size.
Definition TAttText.h:47
TAxis::SetBinLabel
virtual void SetBinLabel(Int_t bin, const char *label)
Set label for bin.
Definition TAxis.cxx:886
TAxis::GetBinLowEdge
virtual Double_t GetBinLowEdge(Int_t bin) const
Return low edge of bin.
Definition TAxis.cxx:518
TAxis::SetRangeUser
virtual void SetRangeUser(Double_t ufirst, Double_t ulast)
Set the viewing range for the axis from ufirst to ulast (in user coordinates, that is,...
Definition TAxis.cxx:1080
TAxis::GetBinWidth
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width.
Definition TAxis.cxx:540
TAxis::GetBinUpEdge
virtual Double_t GetBinUpEdge(Int_t bin) const
Return up edge of bin.
Definition TAxis.cxx:528
TCanvas
The Canvas class.
Definition TCanvas.h:23
TCanvas::cd
TVirtualPad * cd(Int_t subpadnumber=0) override
Set current canvas & pad.
Definition TCanvas.cxx:716
TDirectoryFile::cd
Bool_t cd() override
Change current directory to "this" directory.
Definition TDirectoryFile.cxx:334
TDirectory::GetObject
void GetObject(const char *namecycle, T *&ptr)
Get an object with proper type checking.
Definition TDirectory.h:212
TF1
1-Dim function class
Definition TF1.h:233
TF1::SetParError
virtual void SetParError(Int_t ipar, Double_t error)
Set error for parameter number ipar.
Definition TF1.cxx:3479
TF1::GetNpar
virtual Int_t GetNpar() const
Definition TF1.h:507
TF1::GetParameters
virtual Double_t * GetParameters() const
Definition TF1.h:546
TF1::SetParameter
virtual void SetParameter(Int_t param, Double_t value)
Definition TF1.h:660
TF1::GetParameter
virtual Double_t GetParameter(Int_t ipar) const
Definition TF1.h:538
TFile
A ROOT file is an on-disk file, usually with extension .root, that stores objects in a file-system-li...
Definition TFile.h:53
TFile::Write
Int_t Write(const char *name=nullptr, Int_t opt=0, Int_t bufsiz=0) override
Write memory objects to this file.
Definition TFile.cxx:2429
TFitResultPtr
Provides an indirection to the TFitResult class and with a semantics identical to a TFitResult pointe...
Definition TFitResultPtr.h:32
TGraphErrors
A TGraphErrors is a TGraph with error bars.
Definition TGraphErrors.h:26
TGraph
A TGraph is an object made of two arrays X and Y with npoints each.
Definition TGraph.h:41
TGraph::Draw
void Draw(Option_t *chopt="") override
Draw this graph with its current attributes.
Definition TGraph.cxx:809
TGraph::SetTitle
void SetTitle(const char *title="") override
Change (i.e.
Definition TGraph.cxx:2374
TH1D
1-D histogram with a double per channel (see TH1 documentation)
Definition TH1.h:669
TH1
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:59
TH1::Multiply
virtual Bool_t Multiply(TF1 *f1, Double_t c1=1)
Performs the operation:
Definition TH1.cxx:6017
TH1::SetTitle
void SetTitle(const char *title) override
Change/set the title.
Definition TH1.cxx:6686
TH1::GetNbinsY
virtual Int_t GetNbinsY() const
Definition TH1.h:298
TH1::GetBinError
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition TH1.cxx:9031
TH1::GetXaxis
TAxis * GetXaxis()
Definition TH1.h:324
TH1::GetNbinsX
virtual Int_t GetNbinsX() const
Definition TH1.h:297
TH1::SetMaximum
virtual void SetMaximum(Double_t maximum=-1111)
Definition TH1.h:403
TH1::Add
virtual Bool_t Add(TF1 *h1, Double_t c1=1, Option_t *option="")
Performs the operation: this = this + c1*f1 if errors are defined (see TH1::Sumw2),...
Definition TH1.cxx:826
TH1::SetBinError
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:9174
TH1::GetYaxis
TAxis * GetYaxis()
Definition TH1.h:325
TH1::Draw
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition TH1.cxx:3066
TH1::SetMinimum
virtual void SetMinimum(Double_t minimum=-1111)
Definition TH1.h:404
TH1::SetDefaultSumw2
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:6671
TH1::SetBinContent
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:9190
TH1::GetBinContent
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:5029
TH1::GetBinWidth
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width for 1D histogram.
Definition TH1.cxx:9131
TH1::Scale
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition TH1.cxx:6572
TH1::Clone
TObject * Clone(const char *newname="") const override
Make a complete copy of the underlying object.
Definition TH1.cxx:2752
TH1::Divide
virtual Bool_t Divide(TF1 *f1, Double_t c1=1)
Performs the operation: this = this/(c1*f1) if errors are defined (see TH1::Sumw2),...
Definition TH1.cxx:2840
TH2D
2-D histogram with a double per channel (see TH1 documentation)
Definition TH2.h:357
TH2
Service class for 2-D histogram classes.
Definition TH2.h:30
TH2::SetBinContent
void SetBinContent(Int_t bin, Double_t content) override
Set bin content.
Definition TH2.cxx:2616
TH2::ProjectionX
TH1D * ProjectionX(const char *name="_px", Int_t firstybin=0, Int_t lastybin=-1, Option_t *option="") const
Project a 2-D histogram into a 1-D histogram along X.
Definition TH2.cxx:2436
TH2::GetBinContent
Double_t GetBinContent(Int_t binx, Int_t biny) const override
Definition TH2.h:93
TLegend
This class displays a legend box (TPaveText) containing several legend entries.
Definition TLegend.h:23
TLegend::AddEntry
TLegendEntry * AddEntry(const TObject *obj, const char *label="", Option_t *option="lpf")
Add a new entry to this legend.
Definition TLegend.cxx:317
TLegend::Draw
void Draw(Option_t *option="") override
Draw this legend with its current attributes.
Definition TLegend.cxx:422
TLine
Use the TLine constructor to create a simple line.
Definition TLine.h:22
TMatrixDSymEigen
TMatrixDSymEigen.
Definition TMatrixDSymEigen.h:28
TMatrixTBase::GetNrows
Int_t GetNrows() const
Definition TMatrixTBase.h:123
TNamed::SetTitle
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition TNamed.cxx:164
TNamed::GetName
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
TObject
Mother of all ROOT objects.
Definition TObject.h:41
TObject::Clear
virtual void Clear(Option_t *="")
Definition TObject.h:119
TObject::DrawClone
virtual TObject * DrawClone(Option_t *option="") const
Draw a clone of this object in the current selected pad with: gROOT->SetSelectedPad(c1).
Definition TObject.cxx:299
TObject::Write
virtual Int_t Write(const char *name=nullptr, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
Definition TObject.cxx:880
TObject::Draw
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition TObject.cxx:274
TPad
The most important graphics class in the ROOT system.
Definition TPad.h:28
TPad::Divide
void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0) override
Automatic pad generation by division.
Definition TPad.cxx:1196
TPad::SetFillStyle
void SetFillStyle(Style_t fstyle) override
Override TAttFill::FillStyle for TPad because we want to handle style=0 as style 4000.
Definition TPad.cxx:6004
TPad::SaveAs
void SaveAs(const char *filename="", Option_t *option="") const override
Save the pad content in a file.
Definition TPad.cxx:5714
TPad::cd
TVirtualPad * cd(Int_t subpadnumber=0) override
Set Current pad.
Definition TPad.cxx:640
TPad::Draw
void Draw(Option_t *option="") override
Draw Pad in Current pad (re-parent pad if necessary).
Definition TPad.cxx:1311
TPave::SetBorderSize
virtual void SetBorderSize(Int_t bordersize=4)
Sets the border size of the TPave box and shadow.
Definition TPave.h:77
TRandom3
Random number generator class based on M.
Definition TRandom3.h:27
TRandom::Poisson
virtual ULong64_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Definition TRandom.cxx:404
TSpline
Base class for spline implementation containing the Draw/Paint methods.
Definition TSpline.h:31
TSpline::Draw
void Draw(Option_t *option="") override
Draw this function with its current attributes.
Definition TSpline.cxx:101
TString
Basic string class.
Definition TString.h:139
TString::Format
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:2378
TStyle::SetOptStat
void SetOptStat(Int_t stat=1)
The type of information printed in the histogram statistics box can be selected via the parameter mod...
Definition TStyle.cxx:1636
TText
Base class for several text objects.
Definition TText.h:22
TText::SetNDC
virtual void SetNDC(Bool_t isNDC=kTRUE)
Set NDC mode on if isNDC = kTRUE, off otherwise.
Definition TText.cxx:823
TUnfoldBinning
Binning schemes for use with the unfolding algorithm TUnfoldDensity.
Definition TUnfoldBinning.h:53
TUnfoldBinning::CreateHistogram
TH1 * CreateHistogram(const char *histogramName, Bool_t originalAxisBinning=kFALSE, Int_t **binMap=nullptr, const char *histogramTitle=nullptr, const char *axisSteering=nullptr) const
Create a THxx histogram capable to hold the bins of this binning node and its children.
Definition TUnfoldBinning.cxx:742
TUnfoldBinning::GetBinSize
Double_t GetBinSize(Int_t iBin) const
Get N-dimensional bin size.
Definition TUnfoldBinning.cxx:1723
TUnfoldBinning::GetEndBin
Int_t GetEndBin(void) const
last+1 bin of this node (includes children)
Definition TUnfoldBinning.h:141
TUnfoldDensity
An algorithm to unfold distributions from detector to truth level.
Definition TUnfoldDensity.h:48
TUnfoldDensity::kEScanTauRhoAvg
@ kEScanTauRhoAvg
average global correlation coefficient (from TUnfold::GetRhoI())
Definition TUnfoldDensity.h:107
TUnfoldDensity::GetRhoIJtotal
TH2 * GetRhoIJtotal(const char *histogramName, const char *histogramTitle=nullptr, const char *distributionName=nullptr, const char *projectionMode=nullptr, Bool_t useAxisBinning=kTRUE)
Retrieve correlation coefficients, including all uncertainties.
Definition TUnfoldDensity.cxx:1020
TUnfoldDensity::GetProbabilityMatrix
TH2 * GetProbabilityMatrix(const char *histogramName, const char *histogramTitle=nullptr, Bool_t useAxisBinning=kTRUE) const
Get matrix of probabilities in a new histogram.
Definition TUnfoldDensity.cxx:1161
TUnfoldDensity::GetOutput
TH1 * GetOutput(const char *histogramName, const char *histogramTitle=nullptr, const char *distributionName=nullptr, const char *projectionMode=nullptr, Bool_t useAxisBinning=kTRUE) const
retrieve unfolding result as a new histogram
Definition TUnfoldDensity.cxx:659
TUnfoldDensity::ScanTau
virtual Int_t ScanTau(Int_t nPoint, Double_t tauMin, Double_t tauMax, TSpline **scanResult, Int_t mode=kEScanTauRhoAvg, const char *distribution=nullptr, const char *projectionMode=nullptr, TGraph **lCurvePlot=nullptr, TSpline **logTauXPlot=nullptr, TSpline **logTauYPlot=nullptr)
Scan a function wrt tau and determine the minimum.
Definition TUnfoldDensity.cxx:1352
TUnfoldDensity::kDensityModeNone
@ kDensityModeNone
no scale factors, matrix L is similar to unity matrix
Definition TUnfoldDensity.h:65
TUnfoldSys::SubtractBackground
void SubtractBackground(const TH1 *hist_bgr, const char *name, Double_t scale=1.0, Double_t scale_error=0.0)
Specify a source of background.
Definition TUnfoldSys.cxx:483
TUnfoldSys::SetInput
Int_t SetInput(const TH1 *hist_y, Double_t scaleBias=0.0, Double_t oneOverZeroError=0.0, const TH2 *hist_vyy=nullptr, const TH2 *hist_vyy_inv=nullptr) override
Define the input data for subsequent calls to DoUnfold(Double_t).
Definition TUnfoldSys.cxx:444
TUnfold::DoUnfold
virtual Double_t DoUnfold(void)
Core unfolding algorithm.
Definition TUnfold.cxx:291
TUnfold::ScanLcurve
virtual Int_t ScanLcurve(Int_t nPoint, Double_t tauMin, Double_t tauMax, TGraph **lCurve, TSpline **logTauX=nullptr, TSpline **logTauY=nullptr, TSpline **logTauCurvature=nullptr)
Scan the L curve, determine tau and unfold at the final value of tau.
Definition TUnfold.cxx:2549
TUnfold::kEConstraintArea
@ kEConstraintArea
enforce preservation of the area
Definition TUnfold.h:115
TUnfold::kEConstraintNone
@ kEConstraintNone
use no extra constraint
Definition TUnfold.h:112
TUnfold::ERegMode
ERegMode
choice of regularisation scheme
Definition TUnfold.h:119
TUnfold::kRegModeSize
@ kRegModeSize
regularise the amplitude of the output distribution
Definition TUnfold.h:125
TUnfold::kHistMapOutputHoriz
@ kHistMapOutputHoriz
truth level on x-axis of the response matrix
Definition TUnfold.h:142
TUnfold::GetTau
Double_t GetTau(void) const
Return regularisation parameter.
Definition TUnfold.cxx:3186
TVirtualFitter::SetFCN
virtual void SetFCN(void(*fcn)(Int_t &, Double_t *, Double_t &f, Double_t *, Int_t))
To set the address of the minimization objective function called by the native compiler (see function...
Definition TVirtualFitter.cxx:267
TVirtualFitter::Fitter
static TVirtualFitter * Fitter(TObject *obj, Int_t maxpar=25)
Static function returning a pointer to the current fitter.
Definition TVirtualFitter.cxx:159
TVirtualPad::SetLogx
virtual void SetLogx(Int_t value=1)=0
TVirtualPad::Clear
void Clear(Option_t *option="") override=0
line
TLine * line
Definition entrylistblock_figure1.C:235
sigma
const Double_t sigma
Definition h1analysisProxy.h:11
y
Double_t y[n]
Definition legend1.C:17
c1
return c1
Definition legend1.C:41
x
Double_t x[n]
Definition legend1.C:17
n
const Int_t n
Definition legend1.C:16
PyTorch_Generate_CNN_Model.fit
fit(model, train_loader, val_loader, num_epochs, batch_size, optimizer, criterion, save_best, scheduler)
Definition PyTorch_Generate_CNN_Model.py:34
TMVA_SOFIE_GNN_Parser.h2
h2
Definition TMVA_SOFIE_GNN_Parser.py:188
TMath::Max
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:250
TMath::Prob
Double_t Prob(Double_t chi2, Int_t ndf)
Computation of the probability for a certain Chi-squared (chi2) and number of degrees of freedom (ndf...
Definition TMath.cxx:637
TMath::Sqrt
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:662
TMath::Power
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition TMath.h:721
TMath::Min
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition TMathBase.h:198
TMath::Hypot
Double_t Hypot(Double_t x, Double_t y)
Returns sqrt(x*x + y*y)
Definition TMath.cxx:59
TMath::LandauI
Double_t LandauI(Double_t x)
Returns the cumulative (lower tail integral) of the Landau distribution function at point x.
Definition TMath.cxx:2845
TMath::Abs
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123
graph
Definition graph.py:1
v
@ v
Definition rootcling_impl.cxx:3687
v1
@ v1
Definition rootcling_impl.cxx:3689
Draw
th1 Draw()
sum
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345