98#define TEST_INPUT_COVARIANCE
103TH2 *AddOverflowXY(
TH2 *
h,
double widthX,
double widthY);
106void DrawOverflowX(
TH1 *
h,
double posy);
107void DrawOverflowY(
TH1 *
h,
double posx);
110double const kLegendFontSize=0.05;
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=0,
double tau=0.0,vector<double>
const *
r=0,
120void DrawPadCorrelations(
TH2 *
h,
121 vector<pair<
TF1*,vector<double> > >
const *table);
124 vector<pair<
TF1*,vector<double> > > &table,
int niter=0);
126void GetNiterGraphs(
int iFirst,
int iLast,vector<pair<
TF1*,
127 vector<double> > >
const &table,
int color,
129void GetNiterHist(
int ifit,vector<pair<
TF1*,vector<double> > >
const &table,
130 TH1 *hist[4],
int color,
int style,
int fillStyle);
154 TFile *outputFile=
new TFile(
"testUnfold7_results.root",
"recreate");
158 TFile *inputFile=
new TFile(
"testUnfold7_histograms.root");
164 inputFile->
GetObject(
"fine",fineBinning);
165 inputFile->
GetObject(
"coarse",coarseBinning);
167 if((!fineBinning)||(!coarseBinning)) {
168 cout<<
"problem to read binning schemes\n";
172 fineBinning->
Write();
173 coarseBinning->
Write();
176#define READ(TYPE,binning,name) \
177 TYPE *name[3]; inputFile->GetObject(#name,name[0]); \
179 if(!name[0]) cout<<"Error reading " #name "\n"; \
180 CreateHistogramCopies(name,binning);
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);
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);
196 READ(
TH1,fineBinning,histMcbgrRecF);
197 READ(
TH1,coarseBinning,histMcbgrRecC);
199 TH1 *histOutputCtau0[3];
201 TH1 *histOutputCLCurve[3];
207 double biasScale=1.0;
220 tunfoldC->
SetInput(histDataRecC[0],biasScale);
223 histOutputCtau0[0]=tunfoldC->
GetOutput(
"histOutputCtau0");
225 CreateHistogramCopies(histOutputCtau0,coarseBinning);
229 histOutputCLCurve[0]=tunfoldC->
GetOutput(
"histOutputCLCurve");
231 CreateHistogramCopies(histOutputCLCurve,coarseBinning);
234 TH1 *histOutputFtau0[3];
236 TH1 *histOutputFLCurve[3];
253 tunfoldF->
SetInput(histDataRecF[0],biasScale);
256 histOutputFtau0[0]=tunfoldF->
GetOutput(
"histOutputFtau0");
258 CreateHistogramCopies(histOutputFtau0,coarseBinning);
263 histOutputFLCurve[0]=tunfoldF->
GetOutput(
"histOutputFLCurve");
265 CreateHistogramCopies(histOutputFLCurve,coarseBinning);
268 TH1 *histOutputFAtau0[3];
270 TH1 *histOutputFALCurve[3];
271 TH2 *histRhoFALCurve;
272 TH1 *histOutputFArho[3];
277 double tauFA,tauFArho;
287 tunfoldFA->
SetInput(histDataRecF[0],biasScale);
290 histOutputFAtau0[0]=tunfoldFA->
GetOutput(
"histOutputFAtau0");
292 CreateHistogramCopies(histOutputFAtau0,coarseBinning);
294 tauFArho=tunfoldFA->
GetTau();
295 histOutputFArho[0]=tunfoldFA->
GetOutput(
"histOutputFArho");
297 CreateHistogramCopies(histOutputFArho,coarseBinning);
299 tunfoldFA->
ScanLcurve(50,tauMin,tauMax,&lCurve,&logTauX,&logTauY,&logTauCurvature);
300 tauFA=tunfoldFA->
GetTau();
301 histOutputFALCurve[0]=tunfoldFA->
GetOutput(
"histOutputFALCurve");
303 CreateHistogramCopies(histOutputFALCurve,coarseBinning);
313 TH2 *histProbCO=AddOverflowXY(histProbC,widthC,widthC);
314 TH2 *histProbFO=AddOverflowXY(histProbF,widthC,widthF);
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);
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);
338 TH1 *histMcsigGenO=AddOverflowX(histMcsigGen[2],widthC);
339 TH1 *histDataGenO=AddOverflowX(histDataGen[2],widthC);
342 TH1 *histOutputCtau0O=AddOverflowX(histOutputCtau0[2],widthC);
343 TH2 *histRhoCtau0O=AddOverflowXY(histRhoCtau0,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);
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);
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++) {
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);
372 cout<<
"maximum number of iterations: "<<niter<<
"\n";
374 vector <TH1 *>histOutputAgo,histOutputAgorep;
375 vector <TH2 *>histRhoAgo,histRhoAgorep;
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"));
388 for(
int i=0;i<nx;i++) {
390 for(
int j=0;j<ny;j++) {
393 for(
int j=0;j<ny;j++) {
396 AToverEps(i,j)=aji/epsilonI;
399 for(
int i=0;i<nx;i++) {
401 (histOutputAgo[0]->GetBinError(i+1)
402 *histOutputAgo[0]->GetXaxis()->GetBinWidth(i+1),2.);
404 for(
int i=0;i<ny;i++) {
410 vector<TVectorD *>
y(NREPLICA);
411 vector<TVectorD *> yMb(NREPLICA);
412 vector<TVectorD *> yErr(NREPLICA);
413 vector<TVectorD *>
x(NREPLICA);
415 for(
int nr=0;nr<NREPLICA;nr++) {
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);
428 for(
int i=0;i<ny;i++) {
431 for(
int nr=1;nr<NREPLICA;nr++) {
437 for(
int nr=0;nr<NREPLICA;nr++) {
438 (*yMb[nr])(i)=(*
y[nr])(i)-
b(i);
441 for(
int iter=0;iter<=niter;iter++) {
442 if(!(iter %100)) cout<<iter<<
"\n";
443 for(
int nr=0;nr<NREPLICA;nr++) {
446 for(
int j=0;j<ny;j++) {
447 yOverYrec(j)=(*
y[nr])(j)/yrec(j);
451 for(
int i=0;i<nx;i++) {
452 xx(i) = (*
x[nr])(i) *
f(i);
456 for(
int i=0;i<nx;i++) {
457 for(
int j=0;j<ny;j++) {
458 xdf_dr(i,j) *= (*
x[nr])(i);
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);
469 dxy_dxy.SetSub(0,0,xdf_dr*dr_dxdy);
470 for(
int i=0;i<nx;i++) {
473 for(
int i=0;i<ny;i++) {
474 dxy_dxy(nx+i,nx+i) +=1.0;
482 ((iter<=100)&&(iter %5==0))||
483 ((iter<=1000)&&(iter %50==0))||
485 nIter.push_back(iter);
486 TH1 *
h=(
TH1*)histOutputAgo[0]->Clone
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);
494 TH2 *h2=(
TH2*)histRhoAgo[0]->Clone
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));
501 cout<<
"bad error matrix: iter="<<iter<<
"\n";
508 h=(
TH1*)histOutputAgo[0]->Clone
510 h2=(
TH2*)histRhoAgo[0]->Clone
512 histOutputAgorep.push_back(
h);
513 histRhoAgorep.push_back(h2);
516 double w=1./(NREPLICA-1.);
517 for(
int nr=1;nr<NREPLICA;nr++) {
521 for(
int nr=1;nr<NREPLICA;nr++) {
524 for(
int i=0;i<nx;i++) {
525 dx(i,0)= (*
x[nr])(i)-(*
x[0])(i);
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);
536 for(
int i=0;i<nx;i++) {
537 for(
int j=0;j<nx;j++) {
538 double rho= covAgorep(i,j)/
541 cout<<
"bad error matrix: iter="<<iter<<
"\n";
553 vector<TVectorD*> unfresIDS(NREPLICA),soustr(NREPLICA);
554 cout<<
"IDS number of iterations: "<<niterIDS<<
"\n";
557 for(
int nr=0;nr<NREPLICA;nr++) {
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);
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";
580 for(
int nr=0;nr<NREPLICA;nr++) {
582 IDSfirst(yMb[nr],yErr[nr],&A_IDS,lambdaL,unfresIDS[nr],soustr[nr]);
584 IDSiterate(yMb[nr],yErr[nr],&A_IDS,Am_IDS[nr],
585 lambdaU,lambdaM,lambdaS,
586 unfresIDS[nr],soustr[nr]);
590 for(ix=0;ix<nIter.size();ix++) {
591 if(nIter[ix]==iter)
break;
593 if(ix<nIter.size()) {
594 TH1 *
h=(
TH1*)histOutputIDS[0]->Clone
596 TH2 *h2=(
TH2*)histRhoIDS[0]->Clone
598 histOutputIDS.push_back(
h);
599 histRhoIDS.push_back(h2);
601 double w=1./(NREPLICA-1.);
602 for(
int nr=1;nr<NREPLICA;nr++) {
603 mean += w* (*unfresIDS[nr]);
606 for(
int nr=1;nr<NREPLICA;nr++) {
609 for(
int i=0;i<nx;i++) {
610 dx(i,0)= (*unfresIDS[nr])(i)-(*unfresIDS[0])(i);
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/
622 for(
int i=0;i<nx;i++) {
623 for(
int j=0;j<nx;j++) {
624 double rho= covIDSrep(i,j)/
627 cout<<
"bad error matrix: iter="<<iter<<
"\n";
639 vector<pair<TF1 *,vector<double> > > table;
651 subn[0]=
new TPad(
"subn0",
"",0.,0.5,1.,1.);
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++) {
662 subc[0]=
new TPad(
"sub0",
"",0.,0.5,1.,1.);
664 subc[1]=
new TPad(
"sub1",
"",0.,0.,0.5,0.5);
666 subc[2]=
new TPad(
"sub2",
"",0.5,0.0,1.,0.5);
667 for(
int i=0;i<3;i++) {
675 DrawPadTruth(histMcsigGenO,histDataGenO,0);
677 DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,0,0,0);
678 c2w->
SaveAs(
"exampleTR.eps");
683 DrawPadProbability(histProbCO);
685 DrawPadEfficiency(histEfficiencyC);
686 c2w->
SaveAs(
"exampleAE.eps");
688 int iFitInversion=table.size();
689 DoFit(histOutputCtau0O,histRhoCtau0O,histDataGenO,
"inversion",table);
693 DrawPadTruth(histMcsigGenO,histDataGenO,histOutputCtau0O,
"inversion",0.,
694 &table[table.size()-1].second);
696 DrawPadCorrelations(histRhoCtau0O,&table);
698 DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,
699 histOutputCtau0O,histProbCO,histRhoCtau0O);
700 c3c->
SaveAs(
"inversion.eps");
703 DoFit(histOutputFtau0O,histRhoFtau0O,histDataGenO,
"template",table);
707 DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFtau0O,
"fit",0.,
708 &table[table.size()-1].second);
710 DrawPadCorrelations(histRhoFtau0O,&table);
712 DrawPadReco(histMcRecFO,histMcbgrRecFO,histDataRecFO,
713 histOutputFtau0O,histProbFO,histRhoFtau0O);
714 c3c->
SaveAs(
"template.eps");
716 DoFit(histOutputFAtau0O,histRhoFAtau0O,histDataGenO,
"template+area",table);
720 DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFAtau0O,
"fit",0.,
721 &table[table.size()-1].second);
723 DrawPadCorrelations(histRhoFAtau0O,&table);
725 DrawPadReco(histMcRecFO,histMcbgrRecFO,histDataRecFO,
726 histOutputFAtau0O,histProbFO,histRhoFAtau0O);
727 c3c->
SaveAs(
"templateA.eps");
729 int iFitFALCurve=table.size();
730 DoFit(histOutputFALCurveO,histRhoFALCurveO,histDataGenO,
"Tikhonov+area",table);
733 DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFALCurveO,
"Tikhonov",tauFA,
734 &table[table.size()-1].second);
736 DrawPadCorrelations(histRhoFALCurveO,&table);
738 DrawPadReco(histMcRecFO,histMcbgrRecFO,histDataRecFO,
739 histOutputFALCurveO,histProbFO,histRhoFALCurveO);
740 c3c->
SaveAs(
"lcurveFA.eps");
742 int iFitFArho=table.size();
743 DoFit(histOutputFArhoO,histRhoFArhoO,histDataGenO,
"min(rhomax)",table);
746 DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFArhoO,
"Tikhonov",tauFArho,
747 &table[table.size()-1].second);
749 DrawPadCorrelations(histRhoFArho,&table);
751 DrawPadReco(histMcRecFO,histMcbgrRecFO,histDataRecFO,
752 histOutputFArhoO,histProbFO,histRhoFArhoO);
753 c3c->
SaveAs(
"rhoscanFA.eps");
755 int iFitBinByBin=table.size();
756 DoFit(histOutputBBBO,histRhoBBBO,histDataGenO,
"bin-by-bin",table);
763 DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,
764 histOutputBBBO,histProbCO,histRhoBBBO);
766 DrawPadTruth(histMcsigGenO,histDataGenO,histOutputBBBO,
"bin-by-bin",0.,
767 &table[table.size()-1].second);
768 c2sq->
SaveAs(
"binbybin.eps");
772 int iAgoFirstFit=table.size();
773 for(
size_t i=1;i<histRhoAgorep.size();i++) {
776 DoFit(histOutputAgorep[i],histRhoAgorep[i],histDataGenO,
780 DrawPadTruth(histMcsigGenO,histDataGenO,histOutputAgorep[i],
782 isFitted ? &table[table.size()-1].second : 0);
784 DrawPadCorrelations(histRhoAgorep[i],&table);
786 DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,
787 histOutputAgorep[i],histProbCO,histRhoAgorep[i]);
790 int iAgoLastFit=table.size();
793 int iIDSFirstFit=table.size();
797 for(
size_t i=2;i<histRhoIDS.size();i++) {
800 DoFit(histOutputIDS[i],histRhoIDS[i],histDataGenO,
804 DrawPadTruth(histMcsigGenO,histDataGenO,histOutputIDS[i],
806 isFitted ? &table[table.size()-1].second : 0);
808 DrawPadCorrelations(histRhoIDS[i],&table);
810 DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,
811 histOutputIDS[i],histProbCO,histRhoIDS[i]);
814 int iIDSLastFit=table.size();
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;
839 cout<<
"\""<<
f->GetName()<<
"\","<<
r[2]/
r[3]<<
","<<
r[3]
845 for(
int i=1;i<3;i++) {
846 cout<<
","<<
f->GetParameter(i)<<
",\"\302\261\","<<
f->GetParError(i);
853 lCurve->
SetTitle(
"L curve;log_{10} L_{x};log_{10} L_{y}");
859 logTauX->
SetTitle(
";log_{10} #tau;log_{10} L_{x}");
863 logTauY->
SetTitle(
";log_{10} #tau;log_{10} L_{y}");
866 c4->
SaveAs(
"lcurveL.eps");
870 logTauCurvature->
SetTitle(
";log_{10}(#tau);L curve curvature");
872 logTauCurvature->
Draw();
874 rhoScan->
SetTitle(
";log_{10}(#tau);average(#rho_{i})");
878 DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFALCurveO,
"Tikhonov",tauFA,
879 &table[iFitFALCurve].
second);
881 DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFArhoO,
"Tikhonov",tauFArho,
882 &table[iFitFArho].
second);
884 c4->
SaveAs(
"scanTau.eps");
887 TGraph *graphNiterAgoBay[4];
888 GetNiterGraphs(iAgoFirstFit,iAgoFirstFit+1,table,
kRed-2,graphNiterAgoBay,20);
890 GetNiterGraphs(iAgoFirstFit,iAgoLastFit,table,
kRed,graphNiterAgo,24);
893 GetNiterGraphs(iIDSFirstFit,iIDSLastFit,table,
kMagenta,graphNiterIDS,21);
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);
911 for(
int i=0;i<2;i++) {
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");
922 graphNiterAgoBay[i]->
Draw(
"P");
924 graphNiterIDS[i]->
Draw(
"LP");
928 legend=
new TLegend(0.48,0.28,0.87,0.63);
930 legend=
new TLegend(0.45,0.5,0.88,0.88);
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");
943 legend->
AddEntry( graphNiterIDS[0],
"IDS",
"p");
953 DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFALCurveO,
"Tikhonov",tauFA,
954 &table[iFitFALCurve].
second,table[iFitFALCurve].
first);
958 gPad->SetLogy(
false);
961 histNiterInversion[3]->
Draw(
"SAME HIST ][");
962 histNiterFALCurve[3]->
DrawClone(
"SAME E2");
964 histNiterFALCurve[3]->
Draw(
"SAME HIST ][");
967 histNiterFArho[3]->
Draw(
"SAME HIST ][");
968 histNiterBinByBin[3]->
DrawClone(
"SAME E2");
970 histNiterBinByBin[3]->
Draw(
"SAME HIST ][");
972 line=
new TLine(histNiterInversion[3]->GetXaxis()->GetXmin(),
974 histNiterInversion[3]->GetXaxis()->GetXmax(),
979 graphNiterAgo[3]->
Draw(
"LP");
981 graphNiterAgoBay[3]->
Draw(
"P");
983 graphNiterIDS[3]->
Draw(
"LP");
987 legend=
new TLegend(0.55,0.53,0.95,0.97);
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");
1000 legend->
AddEntry( graphNiterIDS[3],
"IDS",
"p");
1004 c2sq->
SaveAs(
"fitSigma.eps");
1006 outputFile->
Write();
1012void GetNiterGraphs(
int iFirst,
int iLast,vector<pair<
TF1*,
1013 vector<double> > >
const &table,
int color,
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);
1038 for(
int g=0;
g<4;
g++) {
1040 graph[
g]->SetLineColor(color);
1041 graph[
g]->SetMarkerColor(color);
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;
1052 ";iteration;unfold-truth #chi^{2}/N_{D.F.}",1,0.2,1500.);
1056 ";iteration;avg(#rho_{i})",1,0.2,1500.);
1059 ";iteration;parameter #mu",1,0.2,1500.);
1063 ";iteration;parameter #sigma",1,0.2,1500.);
1066 for(
int h=0;
h<4;
h++) {
1070 if( hist[
h]->GetBinError(1)>0.0) {
1083 h[2]=(
TH1 *)
h[1]->Clone(baseName+
"_binw");
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);
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);
1103 h[2]->SetBinContent(iDest,
c/bw);
1104 h[2]->SetBinError(iDest,
e/bw);
1115TH2 *AddOverflowXY(
TH2 *
h,
double widthX,
double widthY) {
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);
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);
1129 yBins[ny]=
h->GetYaxis()->GetBinUpEdge(ny);
1130 yBins[ny+1]=yBins[ny]+widthY;
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));
1145TH1 *AddOverflowX(
TH1 *
h,
double widthX) {
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);
1152 xBins[nx]=
h->GetXaxis()->GetBinUpEdge(nx);
1153 xBins[nx+1]=xBins[nx]+widthX;
1157 for(
int ix=0;ix<=nx+1;ix++) {
1158 r->SetBinContent(ix,
h->GetBinContent(ix));
1159 r->SetBinError(ix,
h->GetBinError(ix));
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) {
1175 TText *textX=
new TText((1.+w1)*
x2-w1*
x1,(1.-posy)*y0+posy*y2,
"Overflow bin");
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());;
1190 TText *textY=
new TText((1.-posx)*x0+posx*
x2,(1.+w1)*y1-w1*y2,
"Overflow bin");
1198void DrawPadProbability(
TH2 *
h) {
1200 h->SetTitle(
"migration probabilities;P_{T}(gen) [GeV];P_{T}(rec) [GeV]");
1201 DrawOverflowX(
h,0.05);
1202 DrawOverflowY(
h,0.35);
1205void DrawPadEfficiency(
TH1 *
h) {
1206 h->SetTitle(
"efficiency;P_{T}(gen) [GeV];#epsilon");
1208 h->SetMinimum(0.75);
1211 DrawOverflowX(
h,0.05);
1216 legEfficiency->
AddEntry(
h,
"reconstruction",
"l");
1218 legEfficiency->
Draw();
1221void DrawPadReco(
TH1 *histMcRec,
TH1 *histMcbgrRec,
TH1 *histDataRec,
1222 TH1 *histDataUnfold,
TH2 *histProbability,
TH2 *histRhoij) {
1225 for(
int i=1;i<=histMcRec->
GetNbinsX();i++) {
1231 histMcRec->
SetTitle(
"Reconstructed;P_{T}(rec);Nevent / GeV");
1232 histMcRec->
Draw(
"HIST");
1239 histMcbgrRec->
Draw(
"SAME HIST");
1241 TH1 * histFoldBack=0;
1242 if(histDataUnfold && histProbability && histRhoij) {
1243 histFoldBack=(
TH1 *)
1246 if((nrec==histProbability->
GetNbinsY())&&
1250 for(
int ix=1;ix<=nrec;ix++) {
1253 for(
int iy=0;iy<=histProbability->
GetNbinsX()+1;iy++) {
1257 for(
int iy2=0;iy2<=histProbability->
GetNbinsX()+1;iy2++) {
1274 cout<<
"can not fold back: "<<nrec
1284 histFoldBack->
Draw(
"SAME HIST");
1289 histDataRec->
Draw(
"SAME");
1290 DrawOverflowX(histMcRec,0.5);
1296 legRec->
AddEntry(histMcRec,
"MC total",
"l");
1297 legRec->
AddEntry(histMcbgrRec,
"background",
"f");
1300 double sumD=0.,sumF=0.,chi2=0.;
1301 for(
int i=1;i<=histDataRec->
GetNbinsX();i++) {
1314 legRec->
AddEntry(histDataRec,
"data",
"lp");
1319void DrawPadTruth(
TH1 *histMcsigGen,
TH1 *histDataGen,
TH1 *histDataUnfold,
1320 char const *
text,
double tau,vector<double>
const *
r,
1325 for(
int i=1;i<=histMcsigGen->
GetNbinsX();i++) {
1326 if(histDataUnfold) {
1344 histMcsigGen->
SetTitle(
"Truth;P_{T};Nevent / GeV");
1346 histMcsigGen->
Draw(
"HIST");
1350 histDataGen->
Draw(
"SAME HIST");
1351 if(histDataUnfold) {
1354 histDataUnfold->
Draw(
"SAME");
1356 DrawOverflowX(histMcsigGen,0.5);
1367 legTruth->
AddEntry(histMcsigGen,
"MC",
"l");
1368 if(!histDataUnfold) legTruth->
AddEntry((
TObject *)0,
" Landau(5,2)",
"");
1369 legTruth->
AddEntry(histDataGen,
"data",
"l");
1370 if(!histDataUnfold) legTruth->
AddEntry((
TObject *)0,
" Landau(6,1.8)",
"");
1371 if(histDataUnfold) {
1374 else t=histDataUnfold->
GetName();
1378 legTruth->
AddEntry(histDataUnfold,t,
"lp");
1382 (
"#chi^{2}/%d=%.1f prob=%.3f",
1383 (
int)(*
r)[3],(*
r)[2]/(*
r)[3],
1391 if(histDataUnfold ) {
1392 TPad *subpad =
new TPad(
"subpad",
"",0.35,0.29,0.88,0.68);
1399 for(
int i=istart;i<=histMcsigGen->
GetNbinsX();i++) {
1415 TH1 *copyDataUnfold=(
TH1*)histDataUnfold->
Clone();
1421 copyMcsigGen->
Draw(
"HIST");
1422 copyDataGen->
Draw(
"SAME HIST");
1423 copyDataUnfold->
Draw(
"SAME");
1430void DrawPadCorrelations(
TH2 *
h,
1431 vector<pair<
TF1*,vector<double> > >
const *table) {
1434 h->SetTitle(
"correlation coefficients;P_{T}(gen) [GeV];P_{T}(gen) [GeV]");
1436 DrawOverflowX(
h,0.05);
1437 DrawOverflowY(
h,0.05);
1443 vector<double>
const &
r=(*table)[table->size()-1].second;
1456 cout<<
"fcn flag=0: npar="<<npar<<
" gin="<<gin<<
" par=[";
1457 for(
int i=0;i<npar;i++) {
1465 for(
int i=0;i<=
n;i++) {
1471 double iy=u[0]*u[2]*(y1-y0)/(
x1-x0);
1487 vector<pair<
TF1 *,vector<double> > > &table,
int niter) {
1489 cout<<
h->GetName()<<
"\n";
1495 int n=
h->GetNbinsX()-1;
1498 for(
int i=0;i<
n;i++) {
1499 for(
int j=0;j<
n;j++) {
1501 (
h->GetBinError(i+1)*
h->GetBinError(j+1));
1507 for(
int i=0;i<
n;i++) {
1511 cout<<
"bad eigenvalue i="<<i<<
" di="<<di(i)<<
"\n";
1520 for(
int i=0;i<
n;i++) {
1524 for(
int j=0;j<
n;j++) {
1533 for(
int i=0;i<=
n;i++) {
1534 for(
int j=0;j<=
n;j++) {
1537 (
h->GetBinError(i+1)*
h->GetBinError(j+1));
1539 if(rho_ij<rhoMin) rhoMin=rho_ij;
1540 if(rho_ij>rhoMax) rhoMax=rho_ij;
1546 TVectorD di1(ev1.GetEigenValues());
1547 for(
int i=0;i<=
n;i++) {
1551 cout<<
"bad eigenvalue i="<<i<<
" di1="<<di1(i)<<
"\n";
1555 TMatrixD O1(ev1.GetEigenVectors());
1558 for(
int i=0;i<=
n;i++) {
1559 double gcor2=1.-1./(vinv1(i,i)*
v1(i,i));
1564 cout<<
"bad global correlation "<<i<<
" "<<gcor2<<
"\n";
1586 double xmax=
h->GetXaxis()->GetBinUpEdge(
h->GetNbinsX()-1);
1587 TF1 *landau=
new TF1(
text,
"[0]*TMath::Landau(x,[1],[2],0)",
1596 vector<double>
r(8);
1599 r[1]=
h->GetNbinsX()-1-landau->
GetNpar();
1600 for(
int i=0;i<
h->GetNbinsX()-1;i++) {
1603 for(
int j=0;j<
h->GetNbinsX()-1;j++) {
1605 r[2]+=di*dj*(*g_fcnMatrix)(i,j);
1608 double pull=di/
h->GetBinError(i+1);
1614 if(!niter)
r[4]=0.25;
1623 table.push_back(make_pair(landau,
r));
1632#include "ids_code.cc"
1638 for(
Int_t i=0; i<N_; i++ ){ (*soustr)[i] = 0.; }
1639 unfres1IDS_ = Unfold( data, dataErr, A_, N_, lambdaL_, soustr );
1645 ModifyMatrix( Am_, A_, unfres2IDS_, dataErr, N_, lambdaM_, soustr, lambdaS_ );
1647 unfres2IDS_ = Unfold( data, dataErr, Am_, N_, lambdaU_, soustr );
static const double x2[5]
static const double x1[5]
include TDocParser_001 C image html pict1_TDocParser_001 png width
R__EXTERN Int_t gErrorIgnoreLevel
TMatrixT< Double_t > TMatrixD
R__EXTERN TStyle * gStyle
TVectorT< Double_t > TVectorD
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
virtual void SetFillStyle(Style_t fstyle)
Set the fill area style.
virtual void SetLineStyle(Style_t lstyle)
Set the line style.
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
virtual void SetLineColor(Color_t lcolor)
Set the line color.
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
virtual void SetMarkerSize(Size_t msize=1)
Set the marker size.
virtual void SetTextAngle(Float_t tangle=0)
Set the text angle.
virtual void SetTextSize(Float_t tsize=1)
Set the text size.
virtual void SetBinLabel(Int_t bin, const char *label)
Set label for bin.
virtual Double_t GetBinLowEdge(Int_t bin) const
Return low edge of bin.
virtual void SetRangeUser(Double_t ufirst, Double_t ulast)
Set the viewing range for the axis from ufirst to ulast (in user coordinates).
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width.
virtual Double_t GetBinUpEdge(Int_t bin) const
Return up edge of bin.
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Bool_t cd(const char *path=nullptr) override
Change current directory to "this" directory.
void GetObject(const char *namecycle, T *&ptr)
virtual void SetParError(Int_t ipar, Double_t error)
Set error for parameter number ipar.
virtual Int_t GetNpar() const
virtual Double_t * GetParameters() const
virtual void SetParameter(Int_t param, Double_t value)
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Int_t Write(const char *name=nullptr, Int_t opt=0, Int_t bufsiz=0) override
Write memory objects to this file.
Provides an indirection to the TFitResult class and with a semantics identical to a TFitResult pointe...
A TGraphErrors is a TGraph with error bars.
A TGraph is an object made of two arrays X and Y with npoints each.
virtual void SetTitle(const char *title="")
Change (i.e.
virtual void Draw(Option_t *chopt="")
Draw this graph with its current attributes.
1-D histogram with a double per channel (see TH1 documentation)}
virtual void SetTitle(const char *title)
See GetStatOverflows for more information.
virtual Bool_t Multiply(TF1 *f1, Double_t c1=1)
Performs the operation:
virtual Int_t GetNbinsY() const
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
virtual Int_t GetNbinsX() const
virtual void SetMaximum(Double_t maximum=-1111)
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),...
virtual void SetBinError(Int_t bin, Double_t error)
Set the bin Error Note that this resets the bin eror option to be of Normal Type and for the non-empt...
virtual void SetMinimum(Double_t minimum=-1111)
static void SetDefaultSumw2(Bool_t sumw2=kTRUE)
When this static function is called with sumw2=kTRUE, all new histograms will automatically activate ...
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
virtual void Draw(Option_t *option="")
Draw this histogram with options.
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width for 1D histogram.
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
virtual Bool_t Divide(TF1 *f1, Double_t c1=1)
Performs the operation: this = this/(c1*f1) if errors are defined (see TH1::Sumw2),...
2-D histogram with a double per channel (see TH1 documentation)}
Service class for 2-Dim histogram classes.
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.
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content.
This class displays a legend box (TPaveText) containing several legend entries.
TLegendEntry * AddEntry(const TObject *obj, const char *label="", Option_t *option="lpf")
Add a new entry to this legend.
virtual void Draw(Option_t *option="")
Draw this legend with its current attributes.
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
virtual const char * GetName() const
Returns name of object.
Mother of all ROOT objects.
virtual Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
virtual TObject * DrawClone(Option_t *option="") const
Draw a clone of this object in the current selected pad for instance with: gROOT->SetSelectedPad(gPad...
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
The most important graphics class in the ROOT system.
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)
Automatic pad generation by division.
virtual void SaveAs(const char *filename="", Option_t *option="") const
Save Pad contents in a file in one of various formats.
virtual void Draw(Option_t *option="")
Draw Pad in Current pad (re-parent pad if necessary).
TVirtualPad * cd(Int_t subpadnumber=0)
Set Current pad.
virtual void SetFillStyle(Style_t fstyle)
Override TAttFill::FillStyle for TPad because we want to handle style=0 as style 4000.
virtual void SetBorderSize(Int_t bordersize=4)
Random number generator class based on M.
virtual Int_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Base class for spline implementation containing the Draw/Paint methods.
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString.
void SetOptStat(Int_t stat=1)
The type of information printed in the histogram statistics box can be selected via the parameter mod...
Base class for several text objects.
virtual void SetNDC(Bool_t isNDC=kTRUE)
Set NDC mode on if isNDC = kTRUE, off otherwise.
Binning schemes for use with the unfolding algorithm TUnfoldDensity.
TH1 * CreateHistogram(const char *histogramName, Bool_t originalAxisBinning=kFALSE, Int_t **binMap=0, const char *histogramTitle=0, const char *axisSteering=0) const
Create a THxx histogram capable to hold the bins of this binning node and its children.
Double_t GetBinSize(Int_t iBin) const
Get N-dimensional bin size.
Int_t GetEndBin(void) const
last+1 bin of this node (includes children)
An algorithm to unfold distributions from detector to truth level.
TH2 * GetRhoIJtotal(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE)
Retrieve correlation coefficients, including all uncertainties.
@ kEScanTauRhoAvg
average global correlation coefficient (from TUnfold::GetRhoI())
TH2 * GetProbabilityMatrix(const char *histogramName, const char *histogramTitle=0, Bool_t useAxisBinning=kTRUE) const
Get matrix of probabilities in a new histogram.
TH1 * GetOutput(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE) const
retrieve unfolding result as a new histogram
virtual Int_t ScanTau(Int_t nPoint, Double_t tauMin, Double_t tauMax, TSpline **scanResult, Int_t mode=kEScanTauRhoAvg, const char *distribution=0, const char *projectionMode=0, TGraph **lCurvePlot=0, TSpline **logTauXPlot=0, TSpline **logTauYPlot=0)
Scan a function wrt tau and determine the minimum.
@ kDensityModeNone
no scale factors, matrix L is similar to unity matrix
virtual Int_t SetInput(const TH1 *hist_y, Double_t scaleBias=0.0, Double_t oneOverZeroError=0.0, const TH2 *hist_vyy=0, const TH2 *hist_vyy_inv=0)
Define the input data for subsequent calls to DoUnfold(Double_t).
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.
virtual Double_t DoUnfold(void)
Core unfolding algorithm.
virtual Int_t ScanLcurve(Int_t nPoint, Double_t tauMin, Double_t tauMax, TGraph **lCurve, TSpline **logTauX=0, TSpline **logTauY=0, TSpline **logTauCurvature=0)
Scan the L curve, determine tau and unfold at the final value of tau.
@ kEConstraintArea
enforce preservation of the area
@ kEConstraintNone
use no extra constraint
ERegMode
choice of regularisation scheme
@ kRegModeSize
regularise the amplitude of the output distribution
@ kHistMapOutputHoriz
truth level on x-axis of the response matrix
Double_t GetTau(void) const
Return regularisation parameter.
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...
static TVirtualFitter * Fitter(TObject *obj, Int_t maxpar=25)
Static function returning a pointer to the current fitter.
static constexpr double s
static constexpr double second
Short_t Max(Short_t a, Short_t b)
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...
Double_t Sqrt(Double_t x)
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Short_t Min(Short_t a, Short_t b)
Double_t Hypot(Double_t x, Double_t y)
Double_t LandauI(Double_t x)
Returns the value of the Landau distribution function at point x.
static long int sum(long int i)