TUnfold: testUnfold7c.C Source File

00001 // Author: Stefan Schmitt
00002 // DESY, July 2016
00003 
00004 //  Version 17.6, example study of unfolding methods
00005 
00006 #include <iostream>
00007 #include <cmath>
00008 #include <map>
00009 #include <TMath.h>
00010 #include <TCanvas.h>
00011 #include <TStyle.h>
00012 #include <TGraph.h>
00013 #include <TGraphErrors.h>
00014 #include <TFile.h>
00015 #include <TROOT.h>
00016 #include <TText.h>
00017 #include <TLine.h>
00018 #include <TLegend.h>
00019 #include <TH1.h>
00020 #include <TF1.h>
00021 #include <TFitter.h>
00022 #include <TMatrixD.h>
00023 #include <TMatrixDSym.h>
00024 #include <TVectorD.h>
00025 #include <TMatrixDSymEigen.h>
00026 #include <TFitResult.h>
00027 #include <TRandom3.h>
00028 #include "TUnfoldDensity.h"
00029 
00030 using namespace std;
00031 
00032 /*
00033   This file is part of TUnfold.
00034 
00035   TUnfold is free software: you can redistribute it and/or modify
00036   it under the terms of the GNU General Public License as published by
00037   the Free Software Foundation, either version 3 of the License, or
00038   (at your option) any later version.
00039 
00040   TUnfold is distributed in the hope that it will be useful,
00041   but WITHOUT ANY WARRANTY; without even the implied warranty of
00042   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00043   GNU General Public License for more details.
00044 
00045   You should have received a copy of the GNU General Public License
00046   along with TUnfold.  If not, see <http://www.gnu.org/licenses/>.
00047 */
00048 
00049 ///////////////////////////////////////////////////////////////////////
00050 // 
00051 // Test program for the classes TUnfoldDensity and TUnfoldBinning
00052 //
00053 // A toy test of the TUnfold package
00054 //
00055 //
00056 // This example is documented in conference proceedings
00057 //   arXiv:1611.01927
00058 // 12th Conference on Quark Confinement and the Hadron Spectrum (Confinement XII)
00059 //
00060 // This is an example of unfolding a one-dimensional distribution
00061 //        this example compares various unfolding methods:
00062 //           matrix iunversion, template fit, L-curve scan,
00063 //           iterative unfolding, etc
00064 //
00065 //        Further details can  be found in:
00066 //           Talk by S.Schmitt at:
00067 //             XII Quark Confinement and the Hadron Spectrum
00068 //             29.8. - 3.9.2016  Thessaloniki, Greece
00069 //             statictics session (+proceedings)
00070 //
00071 // The example comprizes several macros
00072 //   testUnfold7a.C   create root files with TTree objects for
00073 //                      signal, background and data
00074 //            -> write files  testUnfold7_signal.root
00075 //                            testUnfold7_background.root
00076 //                            testUnfold7_data.root
00077 //
00078 //   testUnfold7b.C   loop over trees and fill histograms based on the
00079 //                      TUnfoldBinning objects
00080 //            -> read  testUnfold7binning.xml
00081 //                     testUnfold7_signal.root
00082 //                     testUnfold7_background.root
00083 //                     testUnfold7_data.root
00084 //
00085 //            -> write testUnfold7_histograms.root
00086 //
00087 //   testUnfold7c.C   run the unfolding
00088 //            -> read  testUnfold7_histograms.root
00089 //            -> write testUnfold7_result.root
00090 //            -> write many histograms, to compare various unfolding methods
00091 // 
00092 ///////////////////////////////////////////////////////////////////////
00093 
00094 // #define PRINT_MATRIX_L
00095 
00096 #define TEST_INPUT_COVARIANCE
00097 
00098 void CreateHistogramCopies(TH1 *h[3],TUnfoldBinning const *binning);
00099 void CreateHistogramCopies(TH2 *h[3],TUnfoldBinning const *binningX);
00100 
00101 TH2 *AddOverflowXY(TH2 *h,double widthX,double widthY);
00102 TH1 *AddOverflowX(TH1 *h,double width);
00103 
00104 void DrawOverflowX(TH1 *h,double posy);
00105 void DrawOverflowY(TH1 *h,double posx);
00106 
00107 
00108 double const kLegendFontSize=0.05;
00109 int kNbinC=0;
00110 
00111 void DrawPadProbability(TH2 *h);
00112 void DrawPadEfficiency(TH1 *h);
00113 void DrawPadReco(TH1 *histMcRec,TH1 *histMcbgrRec,TH1 *histDataRec,
00114                  TH1 *histDataUnfold,TH2 *histProbability,TH2 *histRhoij);
00115 void DrawPadTruth(TH1 *histMcsigGen,TH1 *histDataGen,TH1 *histDataUnfold,
00116                   char const *text=0,double tau=0.0,vector<double> const *r=0,
00117                   TF1 *f=0);
00118 void DrawPadCorrelations(TH2 *h,
00119                          vector<pair<TF1*,vector<double> > > const *table);
00120 
00121 TFitResultPtr DoFit(TH1 *h,TH2 *rho,TH1 *truth,char const *text,
00122                     vector<pair<TF1*,vector<double> > > &table,int niter=0);
00123 
00124 void GetNiterGraphs(int iFirst,int iLast,vector<pair<TF1*,
00125                     vector<double> > > const &table,int color,
00126                     TGraph *graph[4],int style);
00127 void GetNiterHist(int ifit,vector<pair<TF1*,vector<double> > > const &table,
00128                   TH1 *hist[4],int color,int style,int fillStyle);
00129 
00130 #ifdef WITH_IDS
00131 void IDSfirst(TVectorD *data, TVectorD *dataErr, TMatrixD *A_, Double_t lambdaL_, TVectorD* &unfres1IDS_,TVectorD *&soustr);
00132 
00133 void IDSiterate(TVectorD *data, TVectorD *dataErr, TMatrixD *A_,TMatrixD *Am_,
00134                 Double_t lambdaU_, Double_t lambdaM_, Double_t lambdaS_,
00135                 TVectorD* &unfres2IDS_ ,TVectorD *&soustr);
00136 #endif
00137 
00138 TRandom3 *g_rnd;
00139 
00140 void testUnfold7c()
00141 {
00142    gErrorIgnoreLevel=kInfo;
00143   // switch on histogram errors
00144   TH1::SetDefaultSumw2();
00145 
00146   gStyle->SetOptStat(0);
00147 
00148   g_rnd=new TRandom3(4711);
00149 
00150   //==============================================
00151   // step 1 : open output file
00152   TFile *outputFile=new TFile("testUnfold7_results.root","recreate");
00153 
00154   //==============================================
00155   // step 2 : read binning schemes and input histograms
00156   TFile *inputFile=new TFile("testUnfold7_histograms.root");
00157 
00158   outputFile->cd();
00159 
00160   TUnfoldBinning *fineBinning,*coarseBinning;
00161 
00162   inputFile->GetObject("fine",fineBinning);
00163   inputFile->GetObject("coarse",coarseBinning);
00164 
00165   if((!fineBinning)||(!coarseBinning)) {
00166      cout<<"problem to read binning schemes\n";
00167   }
00168 
00169   // save binning schemes to output file
00170   fineBinning->Write();
00171   coarseBinning->Write();
00172 
00173   // read histograms
00174 #define READ(TYPE,binning,name)                       \
00175   TYPE *name[3]; inputFile->GetObject(#name,name[0]); \
00176   name[0]->Write();                                   \
00177   if(!name[0]) cout<<"Error reading " #name "\n";     \
00178   CreateHistogramCopies(name,binning);
00179 
00180   outputFile->cd();
00181 
00182   READ(TH1,fineBinning,histDataRecF);
00183   READ(TH1,coarseBinning,histDataRecC);
00184   READ(TH1,fineBinning,histDataBgrF);
00185   READ(TH1,coarseBinning,histDataBgrC);
00186   READ(TH1,coarseBinning,histDataGen);
00187 
00188   READ(TH2,fineBinning,histMcsigGenRecF);
00189   READ(TH2,coarseBinning,histMcsigGenRecC);
00190   READ(TH1,fineBinning,histMcsigRecF);
00191   READ(TH1,coarseBinning,histMcsigRecC);
00192   READ(TH1,coarseBinning,histMcsigGen);
00193 
00194   READ(TH1,fineBinning,histMcbgrRecF);
00195   READ(TH1,coarseBinning,histMcbgrRecC);
00196 
00197   TH1 *histOutputCtau0[3];
00198   TH2 *histRhoCtau0;
00199   TH1 *histOutputCLCurve[3];
00200   TH2 *histRhoCLCurve;
00201   TH2 *histProbC;
00202   double tauMin=1.e-4;
00203   double tauMax=1.e-1;
00204   double fBgr=1.0; // 0.2/0.25;
00205   double biasScale=1.0;
00206   TUnfold::ERegMode mode=TUnfold::kRegModeSize; //Derivative;
00207 
00208   double tauC;
00209   {
00210   TUnfoldDensity *tunfoldC=
00211      new TUnfoldDensity(histMcsigGenRecC[0],
00212                         TUnfold::kHistMapOutputHoriz,
00213                         mode,
00214                         TUnfold::kEConstraintNone,//Area,
00215                         TUnfoldDensity::kDensityModeNone,
00216                         coarseBinning,
00217                         coarseBinning);
00218   tunfoldC->SetInput(histDataRecC[0],biasScale);
00219   tunfoldC->SubtractBackground(histMcbgrRecC[0],"BGR",fBgr,0.0);
00220   tunfoldC->DoUnfold(0.);
00221   histOutputCtau0[0]=tunfoldC->GetOutput("histOutputCtau0");
00222   histRhoCtau0=tunfoldC->GetRhoIJtotal("histRhoCtau0");
00223   CreateHistogramCopies(histOutputCtau0,coarseBinning);
00224   tunfoldC->ScanLcurve(50,tauMin,tauMax,0);
00225   tauC=tunfoldC->GetTau();
00226   //tunfoldC->ScanTau(50,1.E-7,1.E-1,0,TUnfoldDensity::kEScanTauRhoAvg);
00227   histOutputCLCurve[0]=tunfoldC->GetOutput("histOutputCLCurve");
00228   histRhoCLCurve=tunfoldC->GetRhoIJtotal("histRhoCLCurve");
00229   CreateHistogramCopies(histOutputCLCurve,coarseBinning);
00230   histProbC=tunfoldC->GetProbabilityMatrix("histProbC",";P_T(gen);P_T(rec)");
00231   }
00232   TH1 *histOutputFtau0[3];
00233   TH2 *histRhoFtau0;
00234   TH1 *histOutputFLCurve[3];
00235   TH2 *histRhoFLCurve;
00236   TH2 *histProbF;
00237   TGraph *lCurve;
00238   TSpline *logTauX,*logTauY;
00239   tauMin=3.E-4;
00240   tauMax=3.E-2;
00241   double tauF;
00242   {
00243   TUnfoldDensity *tunfoldF=
00244      new TUnfoldDensity(histMcsigGenRecF[0],
00245                         TUnfold::kHistMapOutputHoriz,
00246                         mode,
00247                         TUnfold::kEConstraintNone,//Area,
00248                         TUnfoldDensity::kDensityModeNone,
00249                         coarseBinning,
00250                         fineBinning);
00251   tunfoldF->SetInput(histDataRecF[0],biasScale);
00252   tunfoldF->SubtractBackground(histMcbgrRecF[0],"BGR",fBgr,0.0);
00253   tunfoldF->DoUnfold(0.);
00254   histOutputFtau0[0]=tunfoldF->GetOutput("histOutputFtau0");
00255   histRhoFtau0=tunfoldF->GetRhoIJtotal("histRhoFtau0");
00256   CreateHistogramCopies(histOutputFtau0,coarseBinning);
00257   tunfoldF->ScanLcurve(50,tauMin,tauMax,0);
00258   //tunfoldF->DoUnfold(tauC);
00259   tauF=tunfoldF->GetTau();
00260   //tunfoldF->ScanTau(50,1.E-7,1.E-1,0,TUnfoldDensity::kEScanTauRhoAvg);
00261   histOutputFLCurve[0]=tunfoldF->GetOutput("histOutputFLCurve");
00262   histRhoFLCurve=tunfoldF->GetRhoIJtotal("histRhoFLCurve");
00263   CreateHistogramCopies(histOutputFLCurve,coarseBinning);
00264   histProbF=tunfoldF->GetProbabilityMatrix("histProbF",";P_T(gen);P_T(rec)");
00265   }
00266   TH1 *histOutputFAtau0[3];
00267   TH2 *histRhoFAtau0;
00268   TH1 *histOutputFALCurve[3];
00269   TH2 *histRhoFALCurve;
00270   TH1 *histOutputFArho[3];
00271   TH2 *histRhoFArho;
00272   TSpline *rhoScan=0;
00273   TSpline *logTauCurvature=0;
00274   
00275   double tauFA,tauFArho;
00276   {
00277   TUnfoldDensity *tunfoldFA=
00278      new TUnfoldDensity(histMcsigGenRecF[0],
00279                         TUnfold::kHistMapOutputHoriz,
00280                         mode,
00281                         TUnfold::kEConstraintArea,
00282                         TUnfoldDensity::kDensityModeNone,
00283                         coarseBinning,
00284                         fineBinning);
00285   tunfoldFA->SetInput(histDataRecF[0],biasScale);
00286   tunfoldFA->SubtractBackground(histMcbgrRecF[0],"BGR",fBgr,0.0);
00287   tunfoldFA->DoUnfold(0.);
00288   histOutputFAtau0[0]=tunfoldFA->GetOutput("histOutputFAtau0");
00289   histRhoFAtau0=tunfoldFA->GetRhoIJtotal("histRhoFAtau0");
00290   CreateHistogramCopies(histOutputFAtau0,coarseBinning);
00291   tunfoldFA->ScanTau(50,tauMin,tauMax,&rhoScan,TUnfoldDensity::kEScanTauRhoAvg);
00292   tauFArho=tunfoldFA->GetTau();
00293   histOutputFArho[0]=tunfoldFA->GetOutput("histOutputFArho");
00294   histRhoFArho=tunfoldFA->GetRhoIJtotal("histRhoFArho");
00295   CreateHistogramCopies(histOutputFArho,coarseBinning);
00296 
00297   tunfoldFA->ScanLcurve(50,tauMin,tauMax,&lCurve,&logTauX,&logTauY,&logTauCurvature);
00298   tauFA=tunfoldFA->GetTau();
00299   histOutputFALCurve[0]=tunfoldFA->GetOutput("histOutputFALCurve");
00300   histRhoFALCurve=tunfoldFA->GetRhoIJtotal("histRhoFALCurve");
00301   CreateHistogramCopies(histOutputFALCurve,coarseBinning);
00302   }
00303   lCurve->Write();
00304   logTauX->Write();
00305   logTauY->Write();
00306 
00307 
00308   double widthC=coarseBinning->GetBinSize(histProbC->GetNbinsY()+1);
00309   double widthF=fineBinning->GetBinSize(histProbF->GetNbinsY()+1);
00310 
00311   TH2 *histProbCO=AddOverflowXY(histProbC,widthC,widthC);
00312   TH2 *histProbFO=AddOverflowXY(histProbF,widthC,widthF);
00313 
00314   // efficiency
00315   TH1 *histEfficiencyC=histProbCO->ProjectionX("histEfficiencyC");
00316   kNbinC=histProbCO->GetNbinsX();
00317 
00318   // reconstructed quantities with overflow (coarse binning)
00319   // MC: add signal and bgr
00320   TH1 *histMcsigRecCO=AddOverflowX(histMcsigRecC[2],widthC);
00321   TH1 *histMcbgrRecCO=AddOverflowX(histMcbgrRecC[2],widthC);
00322   histMcbgrRecCO->Scale(fBgr);
00323   TH1 *histMcRecCO=(TH1 *)histMcsigRecCO->Clone("histMcRecC0");
00324   histMcRecCO->Add(histMcsigRecCO,histMcbgrRecCO);
00325   TH1 *histDataRecCO=AddOverflowX(histDataRecC[2],widthC);
00326   //TH1 *histDataRecCNO=AddOverflowX(histDataRecC[1],widthC);
00327 
00328   TH1 *histMcsigRecFO=AddOverflowX(histMcsigRecF[2],widthF);
00329   TH1 *histMcbgrRecFO=AddOverflowX(histMcbgrRecF[2],widthF);
00330   histMcbgrRecFO->Scale(fBgr);
00331   TH1 *histMcRecFO=(TH1 *)histMcsigRecFO->Clone("histMcRecF0");
00332   histMcRecFO->Add(histMcsigRecFO,histMcbgrRecFO);
00333   TH1 *histDataRecFO=AddOverflowX(histDataRecF[2],widthF);
00334 
00335   // truth level with overflow
00336   TH1 *histMcsigGenO=AddOverflowX(histMcsigGen[2],widthC);
00337   TH1 *histDataGenO=AddOverflowX(histDataGen[2],widthC);
00338 
00339   // unfolding result with overflow
00340   TH1 *histOutputCtau0O=AddOverflowX(histOutputCtau0[2],widthC);
00341   TH2 *histRhoCtau0O=AddOverflowXY(histRhoCtau0,widthC,widthC);
00342   //TH1 *histOutputCLCurveO=AddOverflowX(histOutputCLCurve[2],widthC);
00343   //TH2 *histRhoCLCurveO=AddOverflowXY(histRhoCLCurve,widthC,widthC);
00344   TH1 *histOutputFtau0O=AddOverflowX(histOutputFtau0[2],widthC);
00345   TH2 *histRhoFtau0O=AddOverflowXY(histRhoFtau0,widthC,widthC);
00346   TH1 *histOutputFAtau0O=AddOverflowX(histOutputFAtau0[2],widthC);
00347   TH2 *histRhoFAtau0O=AddOverflowXY(histRhoFAtau0,widthC,widthC);
00348   //TH1 *histOutputFLCurveO=AddOverflowX(histOutputFLCurve[2],widthC);
00349   //TH2 *histRhoFLCurveO=AddOverflowXY(histRhoFLCurve,widthC,widthC);
00350   TH1 *histOutputFALCurveO=AddOverflowX(histOutputFALCurve[2],widthC);
00351   TH2 *histRhoFALCurveO=AddOverflowXY(histRhoFALCurve,widthC,widthC);
00352   TH1 *histOutputFArhoO=AddOverflowX(histOutputFArho[2],widthC);
00353   TH2 *histRhoFArhoO=AddOverflowXY(histRhoFArho,widthC,widthC);
00354 
00355   // bin-by-bin
00356   TH2 *histRhoBBBO=(TH2 *)histRhoCtau0O->Clone("histRhoBBBO");
00357   for(int i=1;i<=histRhoBBBO->GetNbinsX();i++) {
00358      for(int j=1;j<=histRhoBBBO->GetNbinsX();j++) {
00359         histRhoBBBO->SetBinContent(i,j,(i==j)?1.:0.);
00360      }
00361   }
00362   TH1 *histDataBgrsub=(TH1 *)histDataRecCO->Clone("histDataBgrsub");
00363   histDataBgrsub->Add(histMcbgrRecCO,-fBgr);
00364   TH1 *histOutputBBBO=(TH1 *)histDataBgrsub->Clone("histOutputBBBO");
00365   histOutputBBBO->Divide(histMcsigRecCO);
00366   histOutputBBBO->Multiply(histMcsigGenO);
00367 
00368   // iterative
00369   int niter=1000;
00370   cout<<"maximum number of iterations: "<<niter<<"\n";
00371 
00372   vector <TH1 *>histOutputAgo,histOutputAgorep;
00373   vector <TH2 *>histRhoAgo,histRhoAgorep;
00374   vector<int> nIter;
00375   histOutputAgo.push_back((TH1*)histMcsigGenO->Clone("histOutputAgo-1"));
00376   histOutputAgorep.push_back((TH1*)histMcsigGenO->Clone("histOutputAgorep-1"));
00377   histRhoAgo.push_back((TH2*)histRhoBBBO->Clone("histRhoAgo-1"));
00378   histRhoAgorep.push_back((TH2*)histRhoBBBO->Clone("histRhoAgorep-1"));
00379   nIter.push_back(-1);
00380 
00381   int nx=histProbCO->GetNbinsX();
00382   int ny=histProbCO->GetNbinsY();
00383   TMatrixD covAgo(nx+ny,nx+ny);
00384   TMatrixD A(ny,nx);
00385   TMatrixD AToverEps(nx,ny);
00386   for(int i=0;i<nx;i++) {
00387      double epsilonI=0.;
00388      for(int j=0;j<ny;j++) {
00389         epsilonI+= histProbCO->GetBinContent(i+1,j+1);
00390      }
00391      for(int j=0;j<ny;j++) {
00392         double aji=histProbCO->GetBinContent(i+1,j+1);
00393         A(j,i)=aji;
00394         AToverEps(i,j)=aji/epsilonI;
00395      }
00396   }     
00397   for(int i=0;i<nx;i++) {
00398      covAgo(i,i)=TMath::Power
00399         (histOutputAgo[0]->GetBinError(i+1)
00400          *histOutputAgo[0]->GetXaxis()->GetBinWidth(i+1),2.);
00401   }
00402   for(int i=0;i<ny;i++) {
00403      covAgo(i+nx,i+nx)=TMath::Power
00404         (histDataRecCO->GetBinError(i+1)
00405          *histDataRecCO->GetXaxis()->GetBinWidth(i+1),2.);
00406   }
00407 #define NREPLICA 300
00408   vector<TVectorD *> y(NREPLICA);
00409   vector<TVectorD *> yMb(NREPLICA);
00410   vector<TVectorD *> yErr(NREPLICA);
00411   vector<TVectorD *> x(NREPLICA);
00412   TVectorD b(nx);
00413   for(int nr=0;nr<NREPLICA;nr++) {
00414      x[nr]=new TVectorD(nx);
00415      y[nr]=new TVectorD(ny);
00416      yMb[nr]=new TVectorD(ny);
00417      yErr[nr]=new TVectorD(ny);
00418   }
00419   for(int i=0;i<nx;i++) {
00420      (*x[0])(i)=histOutputAgo[0]->GetBinContent(i+1)
00421         *histOutputAgo[0]->GetXaxis()->GetBinWidth(i+1);
00422      for(int nr=1;nr<NREPLICA;nr++) {
00423         (*x[nr])(i)=(*x[0])(i);
00424      }
00425   }
00426   for(int i=0;i<ny;i++) {
00427      (*y[0])(i)=histDataRecCO->GetBinContent(i+1)
00428         *histDataRecCO->GetXaxis()->GetBinWidth(i+1);
00429      for(int nr=1;nr<NREPLICA;nr++) {
00430         (*y[nr])(i)=g_rnd->Poisson((*y[0])(i));
00431         (*yErr[nr])(i)=TMath::Sqrt((*y[nr])(i));
00432      }
00433      b(i)=histMcbgrRecCO->GetBinContent(i+1)*
00434         histMcbgrRecCO->GetXaxis()->GetBinWidth(i+1);
00435      for(int nr=0;nr<NREPLICA;nr++) {
00436         (*yMb[nr])(i)=(*y[nr])(i)-b(i);
00437      }
00438   }
00439   for(int iter=0;iter<=niter;iter++) {
00440      if(!(iter %100)) cout<<iter<<"\n";
00441      for(int nr=0;nr<NREPLICA;nr++) {
00442         TVectorD yrec=A*(*x[nr])+b;
00443         TVectorD yOverYrec(ny);
00444         for(int j=0;j<ny;j++) {
00445            yOverYrec(j)=(*y[nr])(j)/yrec(j);
00446         }
00447         TVectorD f=AToverEps * yOverYrec;
00448         TVectorD xx(nx);
00449         for(int i=0;i<nx;i++) {
00450            xx(i) = (*x[nr])(i) * f(i);
00451         }
00452         if(nr==0) {
00453            TMatrixD xdf_dr=AToverEps;
00454            for(int i=0;i<nx;i++) {
00455               for(int j=0;j<ny;j++) {
00456                  xdf_dr(i,j) *= (*x[nr])(i);
00457               }
00458            }
00459            TMatrixD dr_dxdy(ny,nx+ny);
00460            for(int j=0;j<ny;j++) {
00461               dr_dxdy(j,nx+j)=1.0/yrec(j);
00462               for(int i=0;i<nx;i++) {
00463                  dr_dxdy(j,i)= -yOverYrec(j)/yrec(j)*A(j,i);
00464               }
00465            }
00466            TMatrixD dxy_dxy(nx+ny,nx+ny);
00467            dxy_dxy.SetSub(0,0,xdf_dr*dr_dxdy);
00468            for(int i=0;i<nx;i++) {
00469               dxy_dxy(i,i) +=f(i);
00470            }
00471            for(int i=0;i<ny;i++) {
00472               dxy_dxy(nx+i,nx+i) +=1.0;
00473            }
00474            TMatrixD VDT(covAgo,TMatrixD::kMultTranspose,dxy_dxy);
00475            covAgo= dxy_dxy*VDT;
00476         }
00477         (*x[nr])=xx;
00478      }
00479      if((iter<=25)||
00480         ((iter<=100)&&(iter %5==0))||
00481         ((iter<=1000)&&(iter %50==0))||
00482         (iter %1000==0)) {
00483         nIter.push_back(iter);
00484         TH1 * h=(TH1*)histOutputAgo[0]->Clone
00485            (TString::Format("histOutputAgo%d",iter));
00486         histOutputAgo.push_back(h);
00487         for(int i=0;i<nx;i++) {
00488            double bw=h->GetXaxis()->GetBinWidth(i+1);
00489            h->SetBinContent(i+1,(*x[0])(i)/bw);
00490            h->SetBinError(i+1,TMath::Sqrt(covAgo(i,i))/bw);
00491         }
00492         TH2 *h2=(TH2*)histRhoAgo[0]->Clone
00493            (TString::Format("histRhoAgo%d",iter));
00494         histRhoAgo.push_back(h2);
00495         for(int i=0;i<nx;i++) {
00496            for(int j=0;j<nx;j++) {
00497               double rho= covAgo(i,j)/TMath::Sqrt(covAgo(i,i)*covAgo(j,j));
00498               if((i!=j)&&(TMath::Abs(rho)>=1.0)) {
00499                  cout<<"bad error matrix: iter="<<iter<<"\n";
00500                  exit(0);
00501               }
00502               h2->SetBinContent(i+1,j+1,rho);
00503            }
00504         }
00505         // error and correlations from replica analysis
00506         h=(TH1*)histOutputAgo[0]->Clone
00507            (TString::Format("histOutputAgorep%d",iter));
00508         h2=(TH2*)histRhoAgo[0]->Clone
00509            (TString::Format("histRhoAgorep%d",iter));
00510         histOutputAgorep.push_back(h);
00511         histRhoAgorep.push_back(h2);
00512         
00513         TVectorD mean(nx);
00514         double w=1./(NREPLICA-1.);
00515         for(int nr=1;nr<NREPLICA;nr++) {
00516            mean += w* *x[nr];
00517         }
00518         TMatrixD covAgorep(nx,nx);
00519         for(int nr=1;nr<NREPLICA;nr++) {
00520            //TMatrixD dx= (*x)-mean;
00521            TMatrixD dx(nx,1);
00522            for(int i=0;i<nx;i++) {
00523               dx(i,0)= (*x[nr])(i)-(*x[0])(i);
00524            }
00525            covAgorep += w*TMatrixD(dx,TMatrixD::kMultTranspose,dx);
00526         }
00527         
00528         for(int i=0;i<nx;i++) {
00529            double bw=h->GetXaxis()->GetBinWidth(i+1);
00530            h->SetBinContent(i+1,(*x[0])(i)/bw);
00531            h->SetBinError(i+1,TMath::Sqrt(covAgorep(i,i))/bw);
00532            // cout<<i<<" "<<(*x[0])(i)/bw<<" +/-"<<TMath::Sqrt(covAgorep(i,i))/bw<<" "<<TMath::Sqrt(covAgo(i,i))/bw<<"\n";
00533         }
00534         for(int i=0;i<nx;i++) {
00535            for(int j=0;j<nx;j++) {
00536               double rho= covAgorep(i,j)/
00537                  TMath::Sqrt(covAgorep(i,i)*covAgorep(j,j));
00538               if((i!=j)&&(TMath::Abs(rho)>=1.0)) {
00539                  cout<<"bad error matrix: iter="<<iter<<"\n";
00540                  exit(0);
00541               }
00542               h2->SetBinContent(i+1,j+1,rho);
00543            }
00544         }
00545      }
00546   }
00547 
00548 #ifdef WITH_IDS
00549   // IDS Malaescu
00550   int niterIDS=100;
00551   vector<TVectorD*> unfresIDS(NREPLICA),soustr(NREPLICA);
00552   cout<<"IDS number of iterations: "<<niterIDS<<"\n";
00553   TMatrixD *Am_IDS[NREPLICA];
00554   TMatrixD A_IDS(ny,nx);
00555   for(int nr=0;nr<NREPLICA;nr++) {
00556      Am_IDS[nr]=new TMatrixD(ny,nx);
00557   }
00558   for(int iy=0;iy<ny;iy++) {
00559      for(int ix=0;ix<nx;ix++) {
00560         A_IDS(iy,ix)=histMcsigGenRecC[0]->GetBinContent(ix+1,iy+1);
00561      }
00562   }
00563   double lambdaL=0.;
00564   Double_t lambdaUmin = 1.0000002; 
00565   Double_t lambdaMmin = 0.0000001; 
00566   Double_t lambdaS = 0.000001; 
00567   double lambdaU=lambdaUmin;
00568   double lambdaM=lambdaMmin;
00569   vector<TH1 *> histOutputIDS;
00570   vector<TH2 *> histRhoIDS;
00571   histOutputIDS.push_back((TH1*)histOutputAgo[0]->Clone("histOutputIDS-1"));
00572   histRhoIDS.push_back((TH2*)histRhoAgo[0]->Clone("histRhoIDS-1"));
00573   histOutputIDS.push_back((TH1*)histOutputAgo[0]->Clone("histOutputIDS0"));
00574   histRhoIDS.push_back((TH2*)histRhoAgo[0]->Clone("histRhoIDS0"));
00575   for(int iter=1;iter<=niterIDS;iter++) {
00576      if(!(iter %10)) cout<<iter<<"\n";
00577      
00578      for(int nr=0;nr<NREPLICA;nr++) {
00579         if(iter==1) {
00580            IDSfirst(yMb[nr],yErr[nr],&A_IDS,lambdaL,unfresIDS[nr],soustr[nr]);
00581         } else {
00582            IDSiterate(yMb[nr],yErr[nr],&A_IDS,Am_IDS[nr],
00583                       lambdaU,lambdaM,lambdaS,
00584                       unfresIDS[nr],soustr[nr]);
00585         }
00586      }
00587      unsigned ix;
00588      for(ix=0;ix<nIter.size();ix++) {
00589         if(nIter[ix]==iter) break;
00590      }
00591      if(ix<nIter.size()) {
00592         TH1 * h=(TH1*)histOutputIDS[0]->Clone
00593            (TString::Format("histOutputIDS%d",iter));
00594         TH2 *h2=(TH2*)histRhoIDS[0]->Clone
00595            (TString::Format("histRhoIDS%d",iter));
00596         histOutputIDS.push_back(h);
00597         histRhoIDS.push_back(h2);
00598         TVectorD mean(nx);
00599         double w=1./(NREPLICA-1.);
00600         for(int nr=1;nr<NREPLICA;nr++) {
00601            mean += w* (*unfresIDS[nr]);
00602         }
00603         TMatrixD covIDSrep(nx,nx);
00604         for(int nr=1;nr<NREPLICA;nr++) {
00605            //TMatrixD dx= (*x)-mean;
00606            TMatrixD dx(nx,1);
00607            for(int i=0;i<nx;i++) {
00608               dx(i,0)= (*unfresIDS[nr])(i)-(*unfresIDS[0])(i);
00609            }
00610            covIDSrep += w*TMatrixD(dx,TMatrixD::kMultTranspose,dx);
00611         }
00612         for(int i=0;i<nx;i++) {
00613            double bw=h->GetXaxis()->GetBinWidth(i+1);
00614            h->SetBinContent(i+1,(*unfresIDS[0])(i)/bw/
00615                             histEfficiencyC->GetBinContent(i+1));
00616            h->SetBinError(i+1,TMath::Sqrt(covIDSrep(i,i))/bw/
00617                           histEfficiencyC->GetBinContent(i+1));
00618            // cout<<i<<" "<<(*x[0])(i)/bw<<" +/-"<<TMath::Sqrt(covAgorep(i,i))/bw<<" "<<TMath::Sqrt(covAgo(i,i))/bw<<"\n";
00619         }
00620         for(int i=0;i<nx;i++) {
00621            for(int j=0;j<nx;j++) {
00622               double rho= covIDSrep(i,j)/
00623                  TMath::Sqrt(covIDSrep(i,i)*covIDSrep(j,j));
00624               if((i!=j)&&(TMath::Abs(rho)>=1.0)) {
00625                  cout<<"bad error matrix: iter="<<iter<<"\n";
00626                  exit(0);
00627               }
00628               h2->SetBinContent(i+1,j+1,rho);
00629            }
00630         }
00631      }
00632   }
00633 #endif
00634 
00635   //double NEdSmc=histDataBgrsub->GetSumOfWeights();
00636 
00637   vector<pair<TF1 *,vector<double> > > table;
00638 
00639   TCanvas *c1=new TCanvas("c1","",600,600);
00640   TCanvas *c2sq=new TCanvas("c2sq","",600,600);
00641   c2sq->Divide(1,2);
00642   TCanvas *c2w=new TCanvas("c2w","",600,300);
00643   c2w->Divide(2,1);
00644   TCanvas *c4=new TCanvas("c4","",600,600);
00645   c4->Divide(2,2);
00646   //TCanvas *c3n=new TCanvas("c3n","",600,600);
00647   TPad *subn[3];
00648   //gROOT->SetStyle("xTimes2");
00649   subn[0]= new TPad("subn0","",0.,0.5,1.,1.);
00650   //gROOT->SetStyle("square");
00651   subn[1]= new TPad("subn1","",0.,0.,0.5,0.5);
00652   subn[2]= new TPad("subn2","",0.5,0.0,1.,0.5);
00653   for(int i=0;i<3;i++) {
00654      subn[i]->SetFillStyle(0);
00655      subn[i]->Draw();
00656   }
00657   TCanvas *c3c=new TCanvas("c3c","",600,600);
00658   TPad *subc[3];
00659   //gROOT->SetStyle("xTimes2");
00660   subc[0]= new TPad("sub0","",0.,0.5,1.,1.);
00661   //gROOT->SetStyle("squareCOLZ");
00662   subc[1]= new TPad("sub1","",0.,0.,0.5,0.5);
00663   //gROOT->SetStyle("square");
00664   subc[2]= new TPad("sub2","",0.5,0.0,1.,0.5);
00665   for(int i=0;i<3;i++) {
00666      subc[i]->SetFillStyle(0);
00667      subc[i]->Draw();
00668   }
00669 
00670   //=========================== example ==================================
00671 
00672   c2w->cd(1);
00673   DrawPadTruth(histMcsigGenO,histDataGenO,0);
00674   c2w->cd(2);
00675   DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,0,0,0);
00676   c2w->SaveAs("exampleTR.eps");
00677 
00678   //=========================== example ==================================
00679 
00680   c2w->cd(1);
00681   DrawPadProbability(histProbCO);
00682   c2w->cd(2);
00683   DrawPadEfficiency(histEfficiencyC);
00684   c2w->SaveAs("exampleAE.eps");
00685 
00686   int iFitInversion=table.size();
00687   DoFit(histOutputCtau0O,histRhoCtau0O,histDataGenO,"inversion",table);
00688   //=========================== inversion ==================================
00689 
00690   subc[0]->cd();
00691   DrawPadTruth(histMcsigGenO,histDataGenO,histOutputCtau0O,"inversion",0.,
00692                &table[table.size()-1].second);
00693   subc[1]->cd();
00694   DrawPadCorrelations(histRhoCtau0O,&table);
00695   subc[2]->cd();
00696   DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,
00697               histOutputCtau0O,histProbCO,histRhoCtau0O);
00698   c3c->SaveAs("inversion.eps");
00699 
00700 
00701   DoFit(histOutputFtau0O,histRhoFtau0O,histDataGenO,"template",table);
00702   //=========================== template ==================================
00703 
00704   subc[0]->cd();
00705   DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFtau0O,"fit",0.,
00706                &table[table.size()-1].second);
00707   subc[1]->cd();
00708   DrawPadCorrelations(histRhoFtau0O,&table);
00709   subc[2]->cd();
00710   DrawPadReco(histMcRecFO,histMcbgrRecFO,histDataRecFO,
00711               histOutputFtau0O,histProbFO,histRhoFtau0O);
00712   c3c->SaveAs("template.eps");
00713 
00714   DoFit(histOutputFAtau0O,histRhoFAtau0O,histDataGenO,"template+area",table);
00715   //=========================== template+area ==================================
00716 
00717   subc[0]->cd();
00718   DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFAtau0O,"fit",0.,
00719                &table[table.size()-1].second);
00720   subc[1]->cd();
00721   DrawPadCorrelations(histRhoFAtau0O,&table);
00722   subc[2]->cd();
00723   DrawPadReco(histMcRecFO,histMcbgrRecFO,histDataRecFO,
00724               histOutputFAtau0O,histProbFO,histRhoFAtau0O);
00725   c3c->SaveAs("templateA.eps");
00726 
00727   int iFitFALCurve=table.size();
00728   DoFit(histOutputFALCurveO,histRhoFALCurveO,histDataGenO,"Tikhonov+area",table);
00729   //=========================== template+area+tikhonov =====================
00730   subc[0]->cd();
00731   DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFALCurveO,"Tikhonov",tauFA,
00732                 &table[table.size()-1].second);
00733   subc[1]->cd();
00734   DrawPadCorrelations(histRhoFALCurveO,&table);
00735   subc[2]->cd();
00736   DrawPadReco(histMcRecFO,histMcbgrRecFO,histDataRecFO,
00737               histOutputFALCurveO,histProbFO,histRhoFALCurveO);
00738   c3c->SaveAs("lcurveFA.eps");
00739 
00740   int iFitFArho=table.size();
00741   DoFit(histOutputFArhoO,histRhoFArhoO,histDataGenO,"min(rhomax)",table);
00742   //=========================== template+area+tikhonov =====================
00743   subc[0]->cd();
00744   DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFArhoO,"Tikhonov",tauFArho,
00745                 &table[table.size()-1].second);
00746   subc[1]->cd();
00747   DrawPadCorrelations(histRhoFArho,&table);
00748   subc[2]->cd();
00749   DrawPadReco(histMcRecFO,histMcbgrRecFO,histDataRecFO,
00750               histOutputFArhoO,histProbFO,histRhoFArhoO);
00751   c3c->SaveAs("rhoscanFA.eps");
00752 
00753   int iFitBinByBin=table.size();
00754   DoFit(histOutputBBBO,histRhoBBBO,histDataGenO,"bin-by-bin",table);
00755   //=========================== bin-by-bin =================================
00756   //c->cd(1);
00757   //DrawPadProbability(histProbCO);
00758   //c->cd(2);
00759   //DrawPadCorrelations(histRhoBBBO,&table);
00760   c2sq->cd(1);
00761   DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,
00762               histOutputBBBO,histProbCO,histRhoBBBO);
00763   c2sq->cd(2);
00764   DrawPadTruth(histMcsigGenO,histDataGenO,histOutputBBBO,"bin-by-bin",0.,
00765                &table[table.size()-1].second);
00766   c2sq->SaveAs("binbybin.eps");
00767 
00768 
00769   //=========================== iterative ===================================
00770   int iAgoFirstFit=table.size();
00771   for(size_t i=1;i<histRhoAgorep.size();i++) {
00772      int n=nIter[i];
00773      bool isFitted=false;
00774      DoFit(histOutputAgorep[i],histRhoAgorep[i],histDataGenO,
00775            TString::Format("iterative, N=%d",n),table,n);
00776      isFitted=true;
00777      subc[0]->cd();
00778      DrawPadTruth(histMcsigGenO,histDataGenO,histOutputAgorep[i],
00779                   TString::Format("iterative N=%d",nIter[i]),0.,
00780                   isFitted ? &table[table.size()-1].second : 0);
00781      subc[1]->cd();
00782      DrawPadCorrelations(histRhoAgorep[i],&table);
00783      subc[2]->cd();
00784      DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,
00785                  histOutputAgorep[i],histProbCO,histRhoAgorep[i]);
00786      c3c->SaveAs(TString::Format("iterative%d.eps",nIter[i]));
00787   }
00788   int iAgoLastFit=table.size();
00789 
00790 #ifdef WITH_IDS
00791   int iIDSFirstFit=table.size();
00792 
00793   //=========================== IDS ===================================
00794 
00795   for(size_t i=2;i<histRhoIDS.size();i++) {
00796      int n=nIter[i];
00797      bool isFitted=false;
00798      DoFit(histOutputIDS[i],histRhoIDS[i],histDataGenO,
00799            TString::Format("IDS, N=%d",n),table,n);
00800      isFitted=true;
00801      subc[0]->cd();
00802      DrawPadTruth(histMcsigGenO,histDataGenO,histOutputIDS[i],
00803                   TString::Format("IDS N=%d",nIter[i]),0.,
00804                   isFitted ? &table[table.size()-1].second : 0);
00805      subc[1]->cd();
00806      DrawPadCorrelations(histRhoIDS[i],&table);
00807      subc[2]->cd();
00808      DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,
00809                  histOutputIDS[i],histProbCO,histRhoIDS[i]);
00810      c3c->SaveAs(TString::Format("ids%d.eps",nIter[i]));
00811   }
00812   int iIDSLastFit=table.size();
00813 #endif
00814 
00815   int nfit=table.size();
00816   TH1D *fitChindf=new TH1D("fitChindf",";algorithm;#chi^{2}/NDF",nfit,0,nfit);
00817   TH1D *fitNorm=new TH1D("fitNorm",";algorithm;Landau amplitude [1/GeV]",nfit,0,nfit);
00818   TH1D *fitMu=new TH1D("fitMu",";algorithm;Landau #mu [GeV]",nfit,0,nfit);
00819   TH1D *fitSigma=new TH1D("fitSigma",";algorithm;Landau #sigma [GeV]",nfit,0,nfit);
00820   for(int fit=0;fit<nfit;fit++) {
00821      TF1 *f=table[fit].first;
00822      vector<double> const &r=table[fit].second;
00823      fitChindf->GetXaxis()->SetBinLabel(fit+1,f->GetName());
00824      fitNorm->GetXaxis()->SetBinLabel(fit+1,f->GetName());
00825      fitMu->GetXaxis()->SetBinLabel(fit+1,f->GetName());
00826      fitSigma->GetXaxis()->SetBinLabel(fit+1,f->GetName());
00827      double chi2=r[0];
00828      double ndf=r[1];
00829      fitChindf->SetBinContent(fit+1,chi2/ndf);
00830      fitChindf->SetBinError(fit+1,TMath::Sqrt(2./ndf));
00831      fitNorm->SetBinContent(fit+1,f->GetParameter(0));
00832      fitNorm->SetBinError(fit+1,f->GetParError(0));
00833      fitMu->SetBinContent(fit+1,f->GetParameter(1));
00834      fitMu->SetBinError(fit+1,f->GetParError(1));
00835      fitSigma->SetBinContent(fit+1,f->GetParameter(2));
00836      fitSigma->SetBinError(fit+1,f->GetParError(2));
00837      cout<<"\""<<f->GetName()<<"\","<<r[2]/r[3]<<","<<r[3]
00838          <<","<<TMath::Prob(r[2],r[3])
00839          <<","<<chi2/ndf
00840          <<","<<ndf
00841          <<","<<TMath::Prob(r[0],r[1])
00842          <<","<<r[5];
00843      for(int i=1;i<3;i++) {
00844         cout<<","<<f->GetParameter(i)<<",\"\302\261\","<<f->GetParError(i);
00845      }
00846      cout<<"\n";
00847   }
00848 
00849   //=========================== L-curve ==========================
00850   c4->cd(1);
00851   lCurve->SetTitle("L curve;log_{10} L_{x};log_{10} L_{y}");
00852   lCurve->SetLineColor(kRed);
00853   lCurve->Draw("AL");
00854   c4->cd(2);
00855   gPad->Clear();
00856   c4->cd(3);
00857   logTauX->SetTitle(";log_{10} #tau;log_{10} L_{x}");
00858   logTauX->SetLineColor(kBlue);
00859   logTauX->Draw();
00860   c4->cd(4);
00861   logTauY->SetTitle(";log_{10} #tau;log_{10} L_{y}");
00862   logTauY->SetLineColor(kBlue);
00863   logTauY->Draw();
00864   c4->SaveAs("lcurveL.eps");
00865 
00866   //========================= rho and L-curve scan =============== 
00867   c4->cd(1);
00868   logTauCurvature->SetTitle(";log_{10}(#tau);L curve curvature");
00869   logTauCurvature->SetLineColor(kRed);
00870   logTauCurvature->Draw();
00871   c4->cd(2);
00872   rhoScan->SetTitle(";log_{10}(#tau);average(#rho_{i})");
00873   rhoScan->SetLineColor(kRed);
00874   rhoScan->Draw();
00875   c4->cd(3);
00876   DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFALCurveO,"Tikhonov",tauFA,
00877                 &table[iFitFALCurve].second);
00878   c4->cd(4);
00879   DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFArhoO,"Tikhonov",tauFArho,
00880                 &table[iFitFArho].second);
00881 
00882   c4->SaveAs("scanTau.eps");
00883 
00884 
00885   TGraph *graphNiterAgoBay[4];
00886   GetNiterGraphs(iAgoFirstFit,iAgoFirstFit+1,table,kRed-2,graphNiterAgoBay,20);
00887   TGraph *graphNiterAgo[4];
00888   GetNiterGraphs(iAgoFirstFit,iAgoLastFit,table,kRed,graphNiterAgo,24);
00889 #ifdef WITH_IDS
00890   TGraph *graphNiterIDS[4];
00891   GetNiterGraphs(iIDSFirstFit,iIDSLastFit,table,kMagenta,graphNiterIDS,21);
00892 #endif
00893   TH1 *histNiterInversion[4];
00894   GetNiterHist(iFitInversion,table,histNiterInversion,kCyan,1,1001);
00895   TH1 *histNiterFALCurve[4];
00896   GetNiterHist(iFitFALCurve,table,histNiterFALCurve,kBlue,2,3353);
00897   TH1 *histNiterFArho[4];
00898   GetNiterHist(iFitFArho,table,histNiterFArho,kAzure-4,3,3353);
00899   TH1 *histNiterBinByBin[4];
00900   GetNiterHist(iFitBinByBin,table,histNiterBinByBin,kOrange+1,3,3335);
00901 
00902   histNiterInversion[0]->GetYaxis()->SetRangeUser(0.3,500.);
00903   histNiterInversion[1]->GetYaxis()->SetRangeUser(-0.1,0.9);
00904   histNiterInversion[2]->GetYaxis()->SetRangeUser(5.6,6.3);
00905   histNiterInversion[3]->GetYaxis()->SetRangeUser(1.6,2.4);
00906 
00907   TLine *line=0;
00908   c1->cd();
00909   for(int i=0;i<2;i++) {
00910      gPad->Clear();
00911      gPad->SetLogx();
00912      gPad->SetLogy((i<1));
00913      if(! histNiterInversion[i]) continue;
00914      histNiterInversion[i]->Draw("][");
00915      histNiterFALCurve[i]->Draw("SAME ][");
00916      histNiterFArho[i]->Draw("SAME ][");
00917      histNiterBinByBin[i]->Draw("SAME ][");
00918      graphNiterAgo[i]->Draw("LP");
00919      graphNiterAgoBay[i]->SetMarkerStyle(20);
00920      graphNiterAgoBay[i]->Draw("P");
00921 #ifdef WITH_IDS
00922      graphNiterIDS[i]->Draw("LP");
00923 #endif
00924      TLegend *legend;
00925      if(i==1) {
00926         legend=new TLegend(0.48,0.28,0.87,0.63);
00927      } else {
00928         legend=new TLegend(0.45,0.5,0.88,0.88);
00929      }
00930      legend->SetBorderSize(0);
00931      legend->SetFillStyle(1001);
00932      legend->SetFillColor(kWhite);
00933      legend->SetTextSize(kLegendFontSize*0.7);
00934      legend->AddEntry( histNiterInversion[0],"inversion","l");
00935      legend->AddEntry( histNiterFALCurve[0],"Tikhonov L-curve","l");
00936      legend->AddEntry( histNiterFArho[0],"Tikhonov global cor.","l");
00937      legend->AddEntry( histNiterBinByBin[0],"bin-by-bin","l");
00938      legend->AddEntry( graphNiterAgoBay[0],"\"Bayesian\"","p");
00939      legend->AddEntry( graphNiterAgo[0],"iterative","p");
00940 #ifdef WITH_IDS
00941      legend->AddEntry( graphNiterIDS[0],"IDS","p");
00942 #endif
00943      legend->Draw();
00944 
00945      c1->SaveAs(TString::Format("niter%d.eps",i));
00946   }
00947 
00948   //c4->cd(1);
00949   //DrawPadCorrelations(histRhoFALCurveO,&table);
00950   c2sq->cd(1);
00951   DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFALCurveO,"Tikhonov",tauFA,
00952                &table[iFitFALCurve].second,table[iFitFALCurve].first);
00953 
00954   c2sq->cd(2);
00955   gPad->SetLogx();
00956   gPad->SetLogy(false);
00957   histNiterInversion[3]->DrawClone("E2");
00958   histNiterInversion[3]->SetFillStyle(0);
00959   histNiterInversion[3]->Draw("SAME HIST ][");
00960   histNiterFALCurve[3]->DrawClone("SAME E2");
00961   histNiterFALCurve[3]->SetFillStyle(0);
00962   histNiterFALCurve[3]->Draw("SAME HIST ][");
00963   histNiterFArho[3]->DrawClone("SAME E2");
00964   histNiterFArho[3]->SetFillStyle(0);
00965   histNiterFArho[3]->Draw("SAME HIST ][");
00966   histNiterBinByBin[3]->DrawClone("SAME E2");
00967   histNiterBinByBin[3]->SetFillStyle(0);
00968   histNiterBinByBin[3]->Draw("SAME HIST ][");
00969   double yTrue=1.8;
00970   line=new TLine(histNiterInversion[3]->GetXaxis()->GetXmin(),
00971                  yTrue,
00972                  histNiterInversion[3]->GetXaxis()->GetXmax(),
00973                  yTrue);
00974   line->SetLineWidth(3);
00975   line->Draw();
00976 
00977   graphNiterAgo[3]->Draw("LP");
00978   graphNiterAgoBay[3]->SetMarkerStyle(20);
00979   graphNiterAgoBay[3]->Draw("P");
00980 #ifdef WITH_IDS
00981   graphNiterIDS[3]->Draw("LP");
00982 #endif
00983 
00984   TLegend *legend;
00985   legend=new TLegend(0.55,0.53,0.95,0.97);
00986   legend->SetBorderSize(0);
00987   legend->SetFillStyle(1001);
00988   legend->SetFillColor(kWhite);
00989   legend->SetTextSize(kLegendFontSize);
00990   legend->AddEntry( line,"truth","l");
00991   legend->AddEntry( histNiterInversion[3],"inversion","l");
00992   legend->AddEntry( histNiterFALCurve[3],"Tikhonov L-curve","l");
00993   legend->AddEntry( histNiterFArho[3],"Tikhonov global cor.","l");
00994   legend->AddEntry( histNiterBinByBin[3],"bin-by-bin","l");
00995   legend->AddEntry( graphNiterAgoBay[3],"\"Bayesian\"","p");
00996   legend->AddEntry( graphNiterAgo[3],"iterative","p");
00997 #ifdef WITH_IDS
00998   legend->AddEntry( graphNiterIDS[3],"IDS","p");
00999 #endif
01000   legend->Draw();
01001 
01002   c2sq->SaveAs("fitSigma.eps");
01003 
01004   outputFile->Write();
01005 
01006   delete outputFile;
01007 
01008 }
01009 
01010 void GetNiterGraphs(int iFirst,int iLast,vector<pair<TF1*,
01011                     vector<double> > > const &table,int color,
01012                     TGraph *graph[4],int style) {
01013    TVectorD niter(iLast-iFirst);
01014    TVectorD eniter(iLast-iFirst);
01015    TVectorD chi2(iLast-iFirst);
01016    TVectorD gcor(iLast-iFirst);
01017    TVectorD mean(iLast-iFirst);
01018    TVectorD emean(iLast-iFirst);
01019    TVectorD sigma(iLast-iFirst);
01020    TVectorD esigma(iLast-iFirst);
01021    for(int ifit=iFirst;ifit<iLast;ifit++) {
01022       vector<double> const &r=table[ifit].second;
01023       niter(ifit-iFirst)=r[4];
01024       chi2(ifit-iFirst)=r[2]/r[3];
01025       gcor(ifit-iFirst)=r[5];
01026       TF1 const *f=table[ifit].first;
01027       mean(ifit-iFirst)=f->GetParameter(1);
01028       emean(ifit-iFirst)=f->GetParError(1);
01029       sigma(ifit-iFirst)=f->GetParameter(2);
01030       esigma(ifit-iFirst)=f->GetParError(2);
01031    }
01032    graph[0]=new TGraph(niter,chi2);
01033    graph[1]=new TGraph(niter,gcor);
01034    graph[2]=new TGraphErrors(niter,mean,eniter,emean);
01035    graph[3]=new TGraphErrors(niter,sigma,eniter,esigma);
01036    for(int g=0;g<4;g++) {
01037       if(graph[g]) {
01038          graph[g]->SetLineColor(color);
01039          graph[g]->SetMarkerColor(color);
01040          graph[g]->SetMarkerStyle(style);
01041       }
01042    }
01043 }
01044 
01045 void GetNiterHist(int ifit,vector<pair<TF1*,vector<double> > > const &table,
01046                   TH1 *hist[4],int color,int style,int fillStyle) {
01047    vector<double> const &r=table[ifit].second;
01048    TF1 const *f=table[ifit].first;
01049    hist[0]=new TH1D(table[ifit].first->GetName()+TString("_chi2"),
01050                     ";iteration;unfold-truth #chi^{2}/N_{D.F.}",1,0.2,1500.);
01051    hist[0]->SetBinContent(1,r[2]/r[3]);
01052 
01053    hist[1]=new TH1D(table[ifit].first->GetName()+TString("_gcor"),
01054                     ";iteration;avg(#rho_{i})",1,0.2,1500.);
01055    hist[1]->SetBinContent(1,r[5]);
01056    hist[2]=new TH1D(table[ifit].first->GetName()+TString("_mu"),
01057                      ";iteration;parameter #mu",1,0.2,1500.);
01058    hist[2]->SetBinContent(1,f->GetParameter(1));
01059    hist[2]->SetBinError(1,f->GetParError(1));
01060    hist[3]=new TH1D(table[ifit].first->GetName()+TString("_sigma"),
01061                     ";iteration;parameter #sigma",1,0.2,1500.);
01062    hist[3]->SetBinContent(1,f->GetParameter(2));
01063    hist[3]->SetBinError(1,f->GetParError(2));
01064    for(int h=0;h<4;h++) {
01065       if(hist[h]) {
01066          hist[h]->SetLineColor(color);
01067          hist[h]->SetLineStyle(style);
01068          if( hist[h]->GetBinError(1)>0.0) {
01069             hist[h]->SetFillColor(color-10);
01070             hist[h]->SetFillStyle(fillStyle);
01071          }
01072          hist[h]->SetMarkerStyle(0);
01073       }
01074    }
01075 }
01076 
01077 void CreateHistogramCopies(TH1 *h[3],TUnfoldBinning const *binning) {
01078    TString baseName(h[0]->GetName());
01079    Int_t *binMap;
01080    h[1]=binning->CreateHistogram(baseName+"_axis",kTRUE,&binMap);
01081    h[2]=(TH1 *)h[1]->Clone(baseName+"_binw");
01082    Int_t nMax=binning->GetEndBin()+1;
01083    for(Int_t iSrc=0;iSrc<nMax;iSrc++) {
01084       Int_t iDest=binMap[iSrc];
01085       double c=h[0]->GetBinContent(iSrc)+h[1]->GetBinContent(iDest);
01086       double e=TMath::Hypot(h[0]->GetBinError(iSrc),h[1]->GetBinError(iDest));
01087       h[1]->SetBinContent(iDest,c);
01088       h[1]->SetBinError(iDest,e);
01089       h[2]->SetBinContent(iDest,c);
01090       h[2]->SetBinError(iDest,e);
01091    }
01092    for(int iDest=0;iDest<=h[2]->GetNbinsX()+1;iDest++) {
01093       double c=h[2]->GetBinContent(iDest);
01094       double e=h[2]->GetBinError(iDest);
01095       double bw=binning->GetBinSize(iDest);
01096       /* if(bw!=h[2]->GetBinWidth(iDest)) {
01097          cout<<"bin "<<iDest<<" width="<<bw<<" "<<h[2]->GetBinWidth(iDest)
01098              <<"\n";
01099              } */
01100       if(bw>0.0) {
01101          h[2]->SetBinContent(iDest,c/bw);
01102          h[2]->SetBinError(iDest,e/bw);
01103       } else {
01104       }
01105    }
01106 }
01107 
01108 void CreateHistogramCopies(TH2 *h[3],TUnfoldBinning const *binningX) {
01109    h[2]=0;
01110    h[3]=0;
01111 }
01112 
01113 TH2 *AddOverflowXY(TH2 *h,double widthX,double widthY) {
01114    // add overflow bin to X-axis
01115    int nx=h->GetNbinsX();
01116    int ny=h->GetNbinsY();
01117    double *xBins=new double[nx+2];
01118    double *yBins=new double[ny+2];
01119    for(int i=1;i<=nx;i++) {
01120       xBins[i-1]=h->GetXaxis()->GetBinLowEdge(i);
01121    }
01122    xBins[nx]=h->GetXaxis()->GetBinUpEdge(nx);
01123    xBins[nx+1]=xBins[nx]+widthX;
01124    for(int i=1;i<=ny;i++) {
01125       yBins[i-1]=h->GetYaxis()->GetBinLowEdge(i);
01126    }
01127    yBins[ny]=h->GetYaxis()->GetBinUpEdge(ny);
01128    yBins[ny+1]=yBins[ny]+widthY;
01129    TString name(h->GetName());
01130    name+="U";
01131    TH2 *r=new TH2D(name,h->GetTitle(),nx+1,xBins,ny+1,yBins);
01132    for(int ix=0;ix<=nx+1;ix++) {
01133      for(int iy=0;iy<=ny+1;iy++) {
01134         r->SetBinContent(ix,iy,h->GetBinContent(ix,iy));
01135         r->SetBinError(ix,iy,h->GetBinError(ix,iy));
01136      }
01137    }
01138    delete [] yBins;
01139    delete [] xBins;
01140    return r;
01141 }
01142 
01143 TH1 *AddOverflowX(TH1 *h,double widthX) {
01144    // add overflow bin to X-axis
01145    int nx=h->GetNbinsX();
01146    double *xBins=new double[nx+2];
01147    for(int i=1;i<=nx;i++) {
01148       xBins[i-1]=h->GetXaxis()->GetBinLowEdge(i);
01149    }
01150    xBins[nx]=h->GetXaxis()->GetBinUpEdge(nx);
01151    xBins[nx+1]=xBins[nx]+widthX;
01152    TString name(h->GetName());
01153    name+="U";
01154    TH1 *r=new TH1D(name,h->GetTitle(),nx+1,xBins);
01155    for(int ix=0;ix<=nx+1;ix++) {
01156       r->SetBinContent(ix,h->GetBinContent(ix));
01157       r->SetBinError(ix,h->GetBinError(ix));
01158    }
01159    delete [] xBins;
01160    return r;
01161 }
01162 
01163 void DrawOverflowX(TH1 *h,double posy) {
01164   double x1=h->GetXaxis()->GetBinLowEdge(h->GetNbinsX());
01165   double x2=h->GetXaxis()->GetBinUpEdge(h->GetNbinsX());
01166   double y0=h->GetYaxis()->GetBinLowEdge(1);
01167   double y2=h->GetYaxis()->GetBinUpEdge(h->GetNbinsY());;
01168   if(h->GetDimension()==1) {
01169      y0=h->GetMinimum();
01170      y2=h->GetMaximum();
01171   }
01172   double w1=-0.3;
01173   TText *textX=new TText((1.+w1)*x2-w1*x1,(1.-posy)*y0+posy*y2,"Overflow bin");
01174   textX->SetNDC(kFALSE);
01175   textX->SetTextSize(0.05);
01176   textX->SetTextAngle(90.);
01177   textX->Draw();
01178   TLine *lineX=new TLine(x1,y0,x1,y2);
01179   lineX->Draw();
01180 }
01181 
01182 void DrawOverflowY(TH1 *h,double posx) {
01183    double x0=h->GetXaxis()->GetBinLowEdge(1);
01184   double x2=h->GetXaxis()->GetBinUpEdge(h->GetNbinsX());
01185   double y1=h->GetYaxis()->GetBinLowEdge(h->GetNbinsY());;
01186   double y2=h->GetYaxis()->GetBinUpEdge(h->GetNbinsY());;
01187   double w1=-0.3;
01188   TText *textY=new TText((1.-posx)*x0+posx*x2,(1.+w1)*y1-w1*y2,"Overflow bin");
01189   textY->SetNDC(kFALSE);
01190   textY->SetTextSize(0.05);
01191   textY->Draw();
01192   TLine *lineY=new TLine(x0,y1,x2,y1);
01193   lineY->Draw();
01194 }
01195 
01196 void DrawPadProbability(TH2 *h) {
01197   h->Draw("COLZ");
01198   h->SetTitle("migration probabilities;P_{T}(gen) [GeV];P_{T}(rec) [GeV]");
01199   DrawOverflowX(h,0.05);
01200   DrawOverflowY(h,0.35);
01201 }
01202 
01203 void DrawPadEfficiency(TH1 *h) {
01204   h->SetTitle("efficiency;P_{T}(gen) [GeV];#epsilon");
01205   h->SetLineColor(kBlue);
01206   h->SetMinimum(0.75);
01207   h->SetMaximum(1.0);
01208   h->Draw();
01209   DrawOverflowX(h,0.05);
01210   TLegend *legEfficiency=new TLegend(0.3,0.58,0.6,0.75);
01211   legEfficiency->SetBorderSize(0);
01212   legEfficiency->SetFillStyle(0);
01213   legEfficiency->SetTextSize(kLegendFontSize);
01214   legEfficiency->AddEntry(h,"reconstruction","l");
01215   legEfficiency->AddEntry((TObject*)0,"   efficiency","");
01216   legEfficiency->Draw();
01217 }
01218 
01219 void DrawPadReco(TH1 *histMcRec,TH1 *histMcbgrRec,TH1 *histDataRec,
01220                  TH1 *histDataUnfold,TH2 *histProbability,TH2 *histRhoij) {
01221    //gPad->SetLogy(kTRUE);
01222    double amax=0.0;
01223    for(int i=1;i<=histMcRec->GetNbinsX();i++) {
01224       amax=TMath::Max(amax,histMcRec->GetBinContent(i)
01225                       +5.0*histMcRec->GetBinError(i));
01226       amax=TMath::Max(amax,histDataRec->GetBinContent(i)
01227                       +2.0*histDataRec->GetBinError(i));
01228    }
01229    histMcRec->SetTitle("Reconstructed;P_{T}(rec);Nevent / GeV");
01230    histMcRec->Draw("HIST");
01231    histMcRec->SetLineColor(kBlue);
01232    histMcRec->SetMinimum(1.0);
01233    histMcRec->SetMaximum(amax);
01234    //histMcbgrRec->SetFillMode(1);
01235    histMcbgrRec->SetLineColor(kBlue-6);
01236    histMcbgrRec->SetFillColor(kBlue-10);
01237    histMcbgrRec->Draw("SAME HIST");
01238 
01239    TH1 * histFoldBack=0;
01240    if(histDataUnfold && histProbability && histRhoij) {
01241       histFoldBack=(TH1 *)
01242          histMcRec->Clone(histDataUnfold->GetName()+TString("_folded"));
01243       int nrec=histFoldBack->GetNbinsX();
01244       if((nrec==histProbability->GetNbinsY())&&
01245          (nrec==histMcbgrRec->GetNbinsX())&&
01246          (nrec==histDataRec->GetNbinsX())
01247          ) {
01248          for(int ix=1;ix<=nrec;ix++) {
01249             double sum=0.0;
01250             double sume2=0.0;
01251             for(int iy=0;iy<=histProbability->GetNbinsX()+1;iy++) {
01252                sum += histDataUnfold->GetBinContent(iy)*
01253                   histDataUnfold->GetBinWidth(iy)*
01254                   histProbability->GetBinContent(iy,ix);
01255                for(int iy2=0;iy2<=histProbability->GetNbinsX()+1;iy2++) {
01256                   sume2 += histDataUnfold->GetBinError(iy)*
01257                      histDataUnfold->GetBinWidth(iy)*
01258                      histProbability->GetBinContent(iy,ix)*
01259                      histDataUnfold->GetBinError(iy2)*
01260                      histDataUnfold->GetBinWidth(iy2)*
01261                      histProbability->GetBinContent(iy2,ix)*
01262                      histRhoij->GetBinContent(iy,iy2);
01263                }
01264             }
01265             sum /= histFoldBack->GetBinWidth(ix);
01266             sum += histMcbgrRec->GetBinContent(ix);
01267             histFoldBack->SetBinContent(ix,sum);
01268             histFoldBack->SetBinError(ix,TMath::Sqrt(sume2)
01269                                       /histFoldBack->GetBinWidth(ix));
01270          }
01271       } else {
01272          cout<<"can not fold back: "<<nrec
01273              <<" "<<histProbability->GetNbinsY()
01274              <<" "<<histMcbgrRec->GetNbinsX()
01275              <<" "<<histDataRec->GetNbinsX()
01276              <<"\n";
01277          exit(0);
01278       }
01279 
01280       histFoldBack->SetLineColor(kBlack);
01281       histFoldBack->SetMarkerStyle(0);
01282       histFoldBack->Draw("SAME HIST");
01283    }
01284 
01285    histDataRec->SetLineColor(kRed);
01286    histDataRec->SetMarkerColor(kRed);
01287    histDataRec->Draw("SAME");
01288    DrawOverflowX(histMcRec,0.5);
01289    
01290    TLegend *legRec=new TLegend(0.4,0.5,0.68,0.85);
01291    legRec->SetBorderSize(0);
01292    legRec->SetFillStyle(0);
01293    legRec->SetTextSize(kLegendFontSize);
01294    legRec->AddEntry(histMcRec,"MC total","l");
01295    legRec->AddEntry(histMcbgrRec,"background","f");
01296    if(histFoldBack) {
01297       int ndf=-kNbinC;
01298       double sumD=0.,sumF=0.,chi2=0.;
01299       for(int i=1;i<=histDataRec->GetNbinsX();i++) {
01300          //cout<<histDataRec->GetBinContent(i)<<" "<<histFoldBack->GetBinContent(i)<<" "<<" w="<<histFoldBack->GetBinWidth(i)<<"\n";
01301          sumD+=histDataRec->GetBinContent(i)*histDataRec->GetBinWidth(i);
01302          sumF+=histFoldBack->GetBinContent(i)*histFoldBack->GetBinWidth(i);
01303          double pull=(histFoldBack->GetBinContent(i)-histDataRec->GetBinContent(i))/histDataRec->GetBinError(i);
01304          chi2+= pull*pull;
01305          ndf+=1;
01306       }
01307       legRec->AddEntry(histDataRec,TString::Format("data N_{evt}=%.0f",sumD),"lp");
01308       legRec->AddEntry(histFoldBack,TString::Format("folded N_{evt}=%.0f",sumF),"l");
01309       legRec->AddEntry((TObject*)0,TString::Format("#chi^{2}=%.1f ndf=%d",chi2,ndf),"");
01310       //exit(0);
01311    } else {
01312       legRec->AddEntry(histDataRec,"data","lp");
01313    }
01314    legRec->Draw();
01315 }
01316 
01317 void DrawPadTruth(TH1 *histMcsigGen,TH1 *histDataGen,TH1 *histDataUnfold,
01318                   char const *text,double tau,vector<double> const *r,
01319                   TF1 *f) {
01320   //gPad->SetLogy(kTRUE);
01321    double amin=0.;
01322    double amax=0.;
01323    for(int i=1;i<=histMcsigGen->GetNbinsX();i++) {
01324       if(histDataUnfold) {
01325          amin=TMath::Min(amin,histDataUnfold->GetBinContent(i)
01326                          -1.1*histDataUnfold->GetBinError(i));
01327          amax=TMath::Max(amax,histDataUnfold->GetBinContent(i)
01328                          +1.1*histDataUnfold->GetBinError(i));
01329       }
01330       amin=TMath::Min(amin,histMcsigGen->GetBinContent(i)
01331                       -histMcsigGen->GetBinError(i));
01332       amin=TMath::Min(amin,histDataGen->GetBinContent(i)
01333                       -histDataGen->GetBinError(i));
01334       amax=TMath::Max(amax,histMcsigGen->GetBinContent(i)
01335                       +10.*histMcsigGen->GetBinError(i));
01336       amax=TMath::Max(amax,histDataGen->GetBinContent(i)
01337                       +2.*histDataGen->GetBinError(i));
01338    }
01339    histMcsigGen->SetMinimum(amin);
01340    histMcsigGen->SetMaximum(amax);
01341 
01342    histMcsigGen->SetTitle("Truth;P_{T};Nevent / GeV");
01343    histMcsigGen->SetLineColor(kBlue);
01344    histMcsigGen->Draw("HIST");
01345    histDataGen->SetLineColor(kRed);
01346    histDataGen->SetMarkerColor(kRed);
01347    histDataGen->SetMarkerSize(1.0);
01348    histDataGen->Draw("SAME HIST");
01349    if(histDataUnfold) {
01350       histDataUnfold->SetMarkerStyle(21);
01351       histDataUnfold->SetMarkerSize(0.7);
01352       histDataUnfold->Draw("SAME");
01353    }
01354    DrawOverflowX(histMcsigGen,0.5);
01355 
01356    if(f) {
01357       f->SetLineStyle(1);
01358       f->Draw("SAME");
01359    }
01360 
01361    TLegend *legTruth=new TLegend(0.32,0.65,0.6,0.9);
01362    legTruth->SetBorderSize(0);
01363    legTruth->SetFillStyle(0);
01364    legTruth->SetTextSize(kLegendFontSize);
01365    legTruth->AddEntry(histMcsigGen,"MC","l");
01366    if(!histDataUnfold) legTruth->AddEntry((TObject *)0,"  Landau(5,2)","");
01367    legTruth->AddEntry(histDataGen,"data","l");
01368    if(!histDataUnfold) legTruth->AddEntry((TObject *)0,"  Landau(6,1.8)","");
01369    if(histDataUnfold) {
01370       TString t;
01371       if(text) t=text;
01372       else t=histDataUnfold->GetName();
01373       if(tau>0) {
01374          t+=TString::Format(" #tau=%.2g",tau);
01375       }
01376       legTruth->AddEntry(histDataUnfold,t,"lp");
01377       if(r) {
01378          legTruth->AddEntry((TObject *)0,"test wrt data:","");
01379          legTruth->AddEntry((TObject *)0,TString::Format
01380                             ("#chi^{2}/%d=%.1f prob=%.3f",
01381                              (int)(*r)[3],(*r)[2]/(*r)[3],
01382                              TMath::Prob((*r)[2],(*r)[3])),"");
01383       }
01384    }
01385    if(f) {
01386       legTruth->AddEntry(f,"fit","l");
01387    }
01388    legTruth->Draw();
01389    if(histDataUnfold ) {
01390       TPad *subpad = new TPad("subpad","",0.35,0.29,0.88,0.68);
01391       subpad->SetFillStyle(0);
01392       subpad->Draw();
01393       subpad->cd();
01394       double amin=0.,amax=0.;
01395       int istart=11;
01396       for(int i=istart;i<=histMcsigGen->GetNbinsX();i++) {
01397          amin=TMath::Min(amin,histMcsigGen->GetBinContent(i)
01398                          -histMcsigGen->GetBinError(i));
01399          amin=TMath::Min(amin,histDataGen->GetBinContent(i)
01400                          -histDataGen->GetBinError(i));
01401          amin=TMath::Min(amin,histDataUnfold->GetBinContent(i)
01402                          -histDataUnfold->GetBinError(i));
01403          amax=TMath::Max(amax,histMcsigGen->GetBinContent(i)
01404                          +histMcsigGen->GetBinError(i));
01405          amax=TMath::Max(amax,histDataGen->GetBinContent(i)
01406                          +histDataGen->GetBinError(i));
01407          amax=TMath::Max(amax,histDataUnfold->GetBinContent(i)
01408                          +histDataUnfold->GetBinError(i));
01409       }
01410       TH1 *copyMcsigGen=(TH1*)histMcsigGen->Clone();
01411       TH1 *copyDataGen=(TH1*)histDataGen->Clone();
01412       TH1 *copyDataUnfold=(TH1*)histDataUnfold->Clone();
01413       copyMcsigGen->GetXaxis()->SetRangeUser
01414          (copyMcsigGen->GetXaxis()->GetBinLowEdge(istart),
01415           copyMcsigGen->GetXaxis()->GetBinUpEdge(copyMcsigGen->GetNbinsX()-1));
01416       copyMcsigGen->SetTitle(";;");
01417       copyMcsigGen->GetYaxis()->SetRangeUser(amin,amax);
01418       copyMcsigGen->Draw("HIST");
01419       copyDataGen->Draw("SAME HIST");
01420       copyDataUnfold->Draw("SAME");
01421       if(f) {
01422          ((TF1 *)f->Clone())->Draw("SAME");
01423       }
01424    }
01425 }
01426 
01427 void DrawPadCorrelations(TH2 *h,
01428                          vector<pair<TF1*,vector<double> > > const *table) {
01429    h->SetMinimum(-1.);
01430    h->SetMaximum(1.);
01431    h->SetTitle("correlation coefficients;P_{T}(gen) [GeV];P_{T}(gen) [GeV]");
01432    h->Draw("COLZ");
01433    DrawOverflowX(h,0.05);
01434    DrawOverflowY(h,0.05);
01435    if(table) {
01436       TLegend *legGCor=new TLegend(0.13,0.6,0.5,0.8);
01437       legGCor->SetBorderSize(0);
01438       legGCor->SetFillStyle(0);
01439       legGCor->SetTextSize(kLegendFontSize);
01440       vector<double> const &r=(*table)[table->size()-1].second;
01441       legGCor->AddEntry((TObject *)0,TString::Format("min(#rho_{ij})=%5.2f",r[6]),"");
01442       legGCor->AddEntry((TObject *)0,TString::Format("max(#rho_{ij})=%5.2f",r[7]),"");
01443       legGCor->AddEntry((TObject *)0,TString::Format("avg(#rho_i)=%5.2f",r[5]),"");
01444       legGCor->Draw();
01445    }
01446 }
01447 
01448 TH1 *g_fcnHist=0;
01449 TMatrixD *g_fcnMatrix=0;
01450 
01451 void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag) {
01452    int n=g_fcnMatrix->GetNrows();
01453    TVectorD dy(n);
01454    double x0=0,y0=0.;
01455    for(int i=0;i<=n;i++) {
01456       double x1;
01457       if(i<1) x1=g_fcnHist->GetXaxis()->GetBinLowEdge(i+1);
01458       else x1=g_fcnHist->GetXaxis()->GetBinUpEdge(i);
01459       double y1=TMath::LandauI((x1-u[1])/u[2]);
01460       if(i>0) {
01461          double iy=u[0]*u[2]*(y1-y0)/(x1-x0);
01462          dy(i-1)=iy-g_fcnHist->GetBinContent(i);
01463          //cout<<"i="<<i<<" iy="<<iy<<" delta="<< dy(i-1)<<"\n";
01464       }
01465       x0=x1;
01466       y0=y1;
01467       //cout<<"i="<<i<<" y1="<<y1<<" x1="<<x1<<"\n";
01468    }
01469    TVectorD Hdy=(*g_fcnMatrix) * dy;
01470    //Hdy.Print();
01471    f=Hdy*dy;
01472    //exit(0);
01473 }
01474 
01475 
01476 TFitResultPtr DoFit(TH1 *h,TH2 *rho,TH1 *truth,const char *text,
01477                     vector<pair<TF1 *,vector<double> > > &table,int niter) {
01478    TString option="IESN";
01479    cout<<h->GetName()<<"\n";
01480    double gcorAvg=0.;
01481    double rhoMin=0.;
01482    double rhoMax=0.;
01483    if(rho) {
01484       g_fcnHist=h;
01485       int n=h->GetNbinsX()-1; // overflow is included as extra bin, exclude in fit
01486       TMatrixDSym v(n);
01487       //g_fcnMatrix=new TMatrixD(n,n);
01488       for(int i=0;i<n;i++) {
01489          for(int j=0;j<n;j++) {
01490             v(i,j)=rho->GetBinContent(i+1,j+1)*
01491                (h->GetBinError(i+1)*h->GetBinError(j+1));
01492          }
01493       }
01494       TMatrixDSymEigen ev(v);
01495       TMatrixD d(n,n);
01496       TVectorD di(ev.GetEigenValues());
01497       for(int i=0;i<n;i++) {
01498          if(di(i)>0.0) {
01499             d(i,i)=1./di(i);
01500          } else {
01501             cout<<"bad eigenvalue i="<<i<<" di="<<di(i)<<"\n";
01502             exit(0);
01503          }
01504       }
01505       TMatrixD O(ev.GetEigenVectors());
01506       TMatrixD DOT(d,TMatrixD::kMultTranspose,O);
01507       g_fcnMatrix=new TMatrixD(O,TMatrixD::kMult,DOT);
01508       TMatrixD test(*g_fcnMatrix,TMatrixD::kMult,v);
01509       int error=0;
01510       for(int i=0;i<n;i++) {
01511          if(TMath::Abs(test(i,i)-1.0)>1.E-7) {
01512             error++;
01513          }
01514          for(int j=0;j<n;j++) {
01515             if(i==j) continue;
01516             if(TMath::Abs(test(i,j)>1.E-7)) error++;
01517          }
01518       }
01519       // calculate global correlation coefficient (all bins)
01520       TMatrixDSym v1(n+1);
01521       rhoMin=1.;
01522       rhoMax=-1.;
01523       for(int i=0;i<=n;i++) {
01524          for(int j=0;j<=n;j++) {
01525             double rho_ij=rho->GetBinContent(i+1,j+1);
01526             v1(i,j)=rho_ij*
01527                (h->GetBinError(i+1)*h->GetBinError(j+1));
01528             if(i!=j) {
01529                if(rho_ij<rhoMin) rhoMin=rho_ij;
01530                if(rho_ij>rhoMax) rhoMax=rho_ij;
01531             }
01532          }
01533       }
01534       TMatrixDSymEigen ev1(v1);
01535       TMatrixD d1(n+1,n+1);
01536       TVectorD di1(ev1.GetEigenValues());
01537       for(int i=0;i<=n;i++) {
01538          if(di1(i)>0.0) {
01539             d1(i,i)=1./di1(i);
01540          } else {
01541             cout<<"bad eigenvalue i="<<i<<" di1="<<di1(i)<<"\n";
01542             exit(0);
01543          }
01544       }
01545       TMatrixD O1(ev1.GetEigenVectors());
01546       TMatrixD DOT1(d1,TMatrixD::kMultTranspose,O1);
01547       TMatrixD vinv1(O1,TMatrixD::kMult,DOT1);
01548       for(int i=0;i<=n;i++) {
01549          double gcor2=1.-1./(vinv1(i,i)*v1(i,i));
01550          if(gcor2>=0.0) {
01551             double gcor=TMath::Sqrt(gcor2);
01552             gcorAvg += gcor;
01553          } else {
01554             cout<<"bad global correlation "<<i<<" "<<gcor2<<"\n";
01555          }
01556       }
01557       gcorAvg /=(n+1);
01558       /* if(error) {
01559          v.Print();
01560          g_fcnMatrix->Print();
01561          exit(0);
01562          } */
01563       //g_fcnMatrix->Invert();
01564       //from: HFitImpl.cxx
01565       // TVirtualFitter::FCNFunc_t  userFcn = 0;
01566       //    typedef void   (* FCNFunc_t )(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
01567       // userFcn = (TVirtualFitter::GetFitter())->GetFCN();
01568       //  (TVirtualFitter::GetFitter())->SetUserFunc(f1); 
01569       //...
01570       //fitok = fitter->FitFCN( userFcn );
01571 
01572       TVirtualFitter::Fitter(h)->SetFCN(fcn);
01573       option += "U";
01574 
01575    }
01576    double xmax=h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()-1);
01577    TF1 *landau=new TF1(text,"[0]*TMath::Landau(x,[1],[2],0)",
01578                        0.,xmax);
01579    landau->SetParameter(0,6000.);
01580    landau->SetParameter(1,5.);
01581    landau->SetParameter(2,2.);
01582    landau->SetParError(0,10.);
01583    landau->SetParError(1,0.5);
01584    landau->SetParError(2,0.1);
01585    TFitResultPtr s=h->Fit(landau,option,0,0.,xmax);
01586    vector<double> r(8);
01587    int np=landau->GetNpar();
01588    fcn(np,0,r[0],landau->GetParameters(),0);
01589    r[1]=h->GetNbinsX()-1-landau->GetNpar();
01590    for(int i=0;i<h->GetNbinsX()-1;i++) {
01591       double di=h->GetBinContent(i+1)-truth->GetBinContent(i+1);
01592       if(g_fcnMatrix) {
01593          for(int j=0;j<h->GetNbinsX()-1;j++) {
01594             double dj=h->GetBinContent(j+1)-truth->GetBinContent(j+1);
01595             r[2]+=di*dj*(*g_fcnMatrix)(i,j);
01596          }
01597       } else {
01598          double pull=di/h->GetBinError(i+1);
01599          r[2]+=pull*pull;
01600       }
01601       r[3]+=1.0;
01602    }
01603    r[4]=niter;
01604    if(!niter) r[4]=0.25;
01605    r[5]=gcorAvg;
01606    r[6]=rhoMin;
01607    r[7]=rhoMax;
01608    if(rho) {
01609       g_fcnHist=0;
01610       delete g_fcnMatrix;
01611       g_fcnMatrix=0;
01612    }
01613    table.push_back(make_pair(landau,r));
01614    return s;
01615 }
01616 
01617 #ifdef WITH_IDS
01618 
01619 //===================== interface to IDS unfolding code follows here 
01620 // contact Bogdan Malescu to find it
01621 
01622 #include "ids_code.cc"
01623 
01624 void IDSfirst(TVectorD *data, TVectorD *dataErr, TMatrixD *A_, Double_t lambdaL_, TVectorD* &unfres1IDS_,TVectorD *&soustr){
01625 
01626    int N_=data->GetNrows();
01627    soustr = new TVectorD(N_);
01628    for( Int_t i=0; i<N_; i++ ){ (*soustr)[i] = 0.; }
01629    unfres1IDS_ = Unfold( data, dataErr, A_, N_, lambdaL_, soustr );
01630 }
01631 
01632 void IDSiterate(TVectorD *data, TVectorD *dataErr, TMatrixD *A_, TMatrixD *Am_, Double_t lambdaU_, Double_t lambdaM_, Double_t lambdaS_,TVectorD* &unfres2IDS_ ,TVectorD *&soustr) {
01633 
01634    int N_=data->GetNrows();
01635    ModifyMatrix( Am_, A_, unfres2IDS_, dataErr, N_, lambdaM_, soustr, lambdaS_ );
01636    delete unfres2IDS_;
01637    unfres2IDS_ = Unfold( data, dataErr, Am_, N_, lambdaU_, soustr );
01638 }
01639 
01640 #endif