%%cpp -d
#include <iostream>
#include <cmath>
#include <TMath.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TGraph.h>
#include <TFile.h>
#include <TH1.h>
#include "TUnfoldDensity.h"

using std::cout;

 TH1::SetDefaultSumw2();

 TFile *outputFile=new TFile("testUnfold5_results.root","recreate");

 TFile *inputFile=new TFile("testUnfold5_histograms.root");

 outputFile->cd();

 TUnfoldBinning *detectorBinning,*generatorBinning;

 inputFile->GetObject("detector",detectorBinning);
 inputFile->GetObject("generator",generatorBinning);

 if((!detectorBinning)||(!generatorBinning)) {
    cout<<"problem to read binning schemes\n";
 }

 detectorBinning->Write();
 generatorBinning->Write();

 TH1 *histDataReco,*histDataTruth;
 TH2 *histMCGenRec;

 inputFile->GetObject("histDataReco",histDataReco);
 inputFile->GetObject("histDataTruth",histDataTruth);
 inputFile->GetObject("histMCGenRec",histMCGenRec);

#ifdef TEST_ZERO_UNCORR_ERROR

 for(int i=0;i<=histMCGenRec->GetNbinsX()+1;i++) {
    for(int j=0;j<=histMCGenRec->GetNbinsY()+1;j++) {
       histMCGenRec->SetBinError(i,j,0.0);
    }
 }
#endif

 histDataReco->Write();
 histDataTruth->Write();
 histMCGenRec->Write();

 if((!histDataReco)||(!histDataTruth)||(!histMCGenRec)) {
    cout<<"problem to read input histograms\n";
 }

 TUnfold::EConstraint constraintMode= TUnfold::kEConstraintArea;

 TUnfold::ERegMode regMode = TUnfold::kRegModeCurvature;

 TUnfoldDensity::EDensityMode densityFlags=
   TUnfoldDensity::kDensityModeBinWidth;

 const char *REGULARISATION_DISTRIBUTION=nullptr;
 const char *REGULARISATION_AXISSTEERING="*[B]";

 TUnfoldDensity unfold(histMCGenRec,TUnfold::kHistMapOutputHoriz,
                       regMode,constraintMode,densityFlags,
                       generatorBinning,detectorBinning,
           REGULARISATION_DISTRIBUTION,
           REGULARISATION_AXISSTEERING);

#ifdef TEST_INPUT_COVARIANCE

 TH2D *inputEmatrix=
    detectorBinning->CreateErrorMatrixHistogram("input_covar",true);
 for(int i=1;i<=inputEmatrix->GetNbinsX();i++) {
    Double_t e=histDataReco->GetBinError(i);
    inputEmatrix->SetBinContent(i,i,e*e);
    // test: non-zero covariance where variance is zero
    //   if(e<=0.) inputEmatrix->SetBinContent(i,i+1,1.0);
 }
 unfold.SetInput(histDataReco,0.0,0.0,inputEmatrix);
#else
 unfold.SetInput(histDataReco /* ,0.0,1.0 */);
#endif

#ifdef PRINT_MATRIX_L
 TH2 *histL= unfold.GetL("L");
 for(Int_t j=1;j<=histL->GetNbinsY();j++) {
    cout<<"L["<<unfold.GetLBinning()->GetBinName(j)<<"]";
    for(Int_t i=1;i<=histL->GetNbinsX();i++) {
       Double_t c=histL->GetBinContent(i,j);
       if(c!=0.0) cout<<" ["<<i<<"]="<<c;
    }
    cout<<"\n";
 }
#endif

 Int_t nScan=30;
 TSpline *rhoLogTau=nullptr;
 TGraph *lCurve=nullptr;

 const char *SCAN_DISTRIBUTION="signal";
 const char *SCAN_AXISSTEERING=nullptr;

 Int_t iBest=unfold.ScanTau(nScan,0.,0.,&rhoLogTau,
                            TUnfoldDensity::kEScanTauRhoMax,
                            SCAN_DISTRIBUTION,SCAN_AXISSTEERING,
                            &lCurve);

 Double_t t[1],rho[1],x[1],y[1];
 rhoLogTau->GetKnot(iBest,t[0],rho[0]);
 lCurve->GetPoint(iBest,x[0],y[0]);
 TGraph *bestRhoLogTau=new TGraph(1,t,rho);
 TGraph *bestLCurve=new TGraph(1,x,y);
 Double_t *tAll=new Double_t[nScan],*rhoAll=new Double_t[nScan];
 for(Int_t i=0;i<nScan;i++) {
    rhoLogTau->GetKnot(i,tAll[i],rhoAll[i]);
 }
 TGraph *knots=new TGraph(nScan,tAll,rhoAll);

 cout<<"chi**2="<<unfold.GetChi2A()<<"+"<<unfold.GetChi2L()
     <<" / "<<unfold.GetNdf()<<"\n";

 TH1 *histDataUnfold=unfold.GetOutput("unfolded signal",nullptr,nullptr,nullptr,kFALSE);

 TH1 *histMCReco=histMCGenRec->ProjectionY("histMCReco",0,-1,"e");
 TH1 *histMCTruth=histMCGenRec->ProjectionX("histMCTruth",0,-1,"e");
 Double_t scaleFactor=histDataTruth->GetSumOfWeights()/
    histMCTruth->GetSumOfWeights();
 histMCReco->Scale(scaleFactor);
 histMCTruth->Scale(scaleFactor);

 TH2 *histProbability=unfold.GetProbabilityMatrix("histProbability");

 /* TH1 *histGlobalCorr=*/ unfold.GetRhoItotal("histGlobalCorr",nullptr,nullptr,nullptr,kFALSE);
 TH1 *histGlobalCorrScan=unfold.GetRhoItotal
    ("histGlobalCorrScan",nullptr,SCAN_DISTRIBUTION,SCAN_AXISSTEERING,kFALSE);
 /* TH2 *histCorrCoeff=*/ unfold.GetRhoIJtotal("histCorrCoeff",nullptr,nullptr,nullptr,kFALSE);

 TCanvas canvas;
 canvas.Print("testUnfold5.ps[");

 canvas.Clear();
 canvas.Divide(3,2);

 canvas.cd(1);
 histDataReco->SetMinimum(0.0);
 histDataReco->Draw("E");
 histMCReco->SetLineColor(kBlue);
 histMCReco->Draw("SAME HIST");

 canvas.cd(2);
 histProbability->Draw("BOX");

 canvas.cd(3);
 gPad->SetLogy();
 histDataUnfold->Draw("E");
 histDataTruth->SetLineColor(kBlue);
 histDataTruth->Draw("SAME HIST");
 histMCTruth->SetLineColor(kRed);
 histMCTruth->Draw("SAME HIST");

 canvas.cd(4);
 rhoLogTau->Draw();
 knots->Draw("*");
 bestRhoLogTau->SetMarkerColor(kRed);
 bestRhoLogTau->Draw("*");

 canvas.cd(5);

 histGlobalCorrScan->Draw("HIST");

 canvas.cd(6);
 lCurve->Draw("AL");
 bestLCurve->SetMarkerColor(kRed);
 bestLCurve->Draw("*");


 canvas.Print("testUnfold5.ps");

 canvas.Print("testUnfold5.ps]");

gROOT->GetListOfCanvases()->Draw()