dagostini.C

// Purpose: d'Agostini ("Bayesian" or Richardson-Lucy) unfolding, including
//          response matrix generation from NLO theory and parameterized JER
// Author:  mikko.voutilainen@cern.ch
// Created: September 2, 2012
// Updated: June 5, 2015
#include "TFile.h"
#include "TDirectory.h"
#include "TList.h"
#include "TObject.h"
#include "TKey.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TGraphErrors.h"
#include "TF1.h"
#include "TMath.h"
#include "TMatrixD.h"
#include "TCanvas.h"
#include "TLine.h"
#include "TStyle.h"
#include "TLegend.h"
#include "TDecompSVD.h"

#include "RooUnfold/src/RooUnfold.h"
#include "RooUnfold/src/RooUnfoldBayes.h"
#include "RooUnfold/src/RooUnfoldBinByBin.h"
#include "RooUnfold/src/RooUnfoldSvd.h"
#include "RooUnfold/src/RooUnfoldResponse.h"
//#include "RooUnfold.h"

#include "TUnfold.h"
#include "TUnfoldDensity.h"

#include "tdrstyle_mod18.C"
#include "ptresolution.h"
#include "settings.h"
#include "tools.h"

#include <iostream>

using namespace std;

bool _jet = false;

// Resolution function
Double_t fPtRes(Double_t *x, Double_t *p) { return ptresolution(x[0], p[0]);}

// Ansatz Kernel 
int cnt_a = 0;
const int nk = 3; // number of kernel parameters (excluding pt, eta)

//quark and gluon response fit which is normilized to inc jets
TF1 *RNF_g=0;
TF1 *RNF_q=0;

Double_t smearedAnsatzKernel(Double_t *x, Double_t *p) {

  if (++cnt_a%1000000==0) cout << "+" << flush;

  const double pt = x[0]; // true pT
  const double ptmeas = p[0]; // measured pT
  const double eta = p[1]; // rapidity

  
  double mean=0;
  if(jp::isgluon && !jp::isquark){
      if (!RNF_g) RNF_g=new TF1("RNF_g","([0]+[1]*pow(x,[2]))/(0.997212+0.568271*pow(x,-1.16785))", 1,5000);
      RNF_g->SetParameters(0.984455,0.0953105,-0.493221);
      mean = RNF_g->Eval(pt)*pt;
  }else if(jp::isquark && !jp::isgluon){
      if (!RNF_q) RNF_q=new TF1("RNF_q","([0]+[1]*pow(x,[2]))/(0.997212+0.568271*pow(x,-1.16785))", 1,5000);
      RNF_q->SetParameters(0.989670,0.150814,-0.331038);
      mean = RNF_q->Eval(pt)*pt;
  }else{
      mean = pt;
}

  
  double res = ptresolution(pt, eta+1e-3) * pt;
  const double s = TMath::Gaus(ptmeas, mean, res, kTRUE);
  double scale = 1;
  //const double s = TMath::Gaus(ptmeas, pt*scale, res, kTRUE);
  const double f = p[2] * pow(pt, p[3]) * pow(1 - pt*cosh(eta) / jp::emax, p[4]);
  
  return (f * s);
}

// Smeared Ansatz
double _epsilon = 1e-12;
TF1 *_kernel = 0; // global variable, not pretty but works
Double_t smearedAnsatz(Double_t *x, Double_t *p) {

  const double pt = x[0];
  const double eta = p[0];

  if (!_kernel) _kernel = new TF1("_kernel", smearedAnsatzKernel, 1., jp::emax/cosh(eta), nk+2);

  double res = ptresolution(pt, eta+1e-3) * pt;
  const double sigma = max(0.10, min(res/pt, 0.30));
  double ptmin = pt / (1. + 4.*sigma); // xmin*(1+4*sigma)=x
  ptmin = max(1.,ptmin); // safety check
  double ptmax = pt / (1. - 3.*sigma); // xmax*(1-3*sigma)=x
  //  cout << Form("1pt %10.5f sigma %10.5f ptmin %10.5f ptmax %10.5f eta %10.5f",pt, sigma, ptmin, ptmax, eta) << endl << flush;
  ptmax = min(jp::emax/cosh(eta), ptmax); // safety check
  //  cout << Form("2pt %10.5f sigma %10.5f ptmin %10.5f ptmax %10.5f eta %10.5f",pt, sigma, ptmin, ptmax, eta) << endl << flush;

  const double par[nk+2] = {pt, eta, p[1], p[2], p[3]};
  _kernel->SetParameters(&par[0]);

  // Set pT bin limits needed in smearing matrix generation
  if (p[4]>0 && p[4]<jp::emax/cosh(eta)) ptmin = p[4];
  if (p[5]>0 && p[5]<jp::emax/cosh(eta)) ptmax = p[5];

  return ( _kernel->Integral(ptmin, ptmax, _epsilon) );
  //  return ( 1.0); // integral fails due to nan ptmin ptmax
}
void recurseFile(TDirectory *indir, TDirectory *indir2, TDirectory *outdir,
                 bool ismc);
void dagostiniUnfold_histo(TH1D *hpt, TH1D *hpt2, TDirectory *outdir,
			   bool ismc, bool kscale = false, string id = "");


void dagostiniUnfold(string type) {

  //TFile *fin = new TFile(Form("../ROOTFiles_of_Laura_for_test/output-%s-2b.root",type.c_str()),"READ");
  //TFile *fin = new TFile(Form("output-%s-1.root","DATA"),"READ");
  TFile *fin = new TFile(Form("output-%s-QGUL.root",type.c_str()),"READ");
  assert(fin && !fin->IsZombie());

  // TFile *fin2 = new TFile(Form("output-%s-2c.root",type.c_str()),"READ");
  // TFile *fin2 = new TFile(Form("output-%s-2b.root",type.c_str()),"READ");
  TFile *fin2 = new TFile(Form("output-%s-QGUL.root","MC"),"READ");
  //TFile *fin2 = new TFile("../ROOTFiles_of_Laura_for_test/output-MC-2b.root","READ");
  //TFile *fin2 = new TFile(jp::dagfile1 ? "output-MC-1.root" : "../ROOTFiles_of_Laura_for_test/output-MC-2b.root","READ");
  assert(fin2 && !fin2->IsZombie());
  
  TFile *fout;
  if(jp::isgluon && !jp::isquark){ fout = new TFile(Form("output-%s-3_gluon.root",type.c_str()),"RECREATE"); assert(fout && !fout->IsZombie());}
  else if(jp::isquark && !jp::isgluon){ fout = new TFile(Form("output-%s-3_quark.root",type.c_str()),"RECREATE"); assert(fout && !fout->IsZombie());}
  else{ fout = new TFile(Form("output-%s-3.root",type.c_str()),"RECREATE"); assert(fout && !fout->IsZombie());}

  bool ismc = jp::ismc;

  recurseFile(fin, fin2, fout, ismc);

  cout << "Output stored in " << fout->GetName() << endl;
  fout->Close();
  fout->Delete();

  fin->Close();
  fin->Delete();
}


void recurseFile(TDirectory *indir, TDirectory *indir2, TDirectory *outdir,
                 bool ismc) {

  TDirectory *curdir = gDirectory;

  // Automatically go through the list of keys (directories)
  TList *keys = indir->GetListOfKeys();
  TListIter itkey(keys);
  TObject *key, *obj;

  while ( (key = itkey.Next()) ) {

    obj = ((TKey*)key)->ReadObj(); assert(obj);

    // Found a subdirectory: copy it to output and go deeper
    if (obj->InheritsFrom("TDirectory")) {

      if (jp::debug) cout << key->GetName() << endl;

      assert(outdir->mkdir(obj->GetName()));
      outdir->mkdir(obj->GetName());
      assert(outdir->cd(obj->GetName()));
      TDirectory *outdir2 = outdir->GetDirectory(obj->GetName()); assert(outdir2);
      outdir2->cd();

      assert(indir->cd(obj->GetName()));
      TDirectory *indir2a = indir->GetDirectory(obj->GetName()); assert(indir2a);
      indir2a->cd();

      if (indir2->cd(obj->GetName())) {
        TDirectory *indir2b = indir2->GetDirectory(obj->GetName()); assert(indir2b);

        recurseFile(indir2a, indir2b, outdir2, ismc);
      }
    } // inherits from TDirectory

    // Found hpt plots: call unfolding routine
    if(jp::isgluon && !jp::isquark){
        cout << "gluon unfolding routine"<<endl;
        if (obj->InheritsFrom("TH1") && (string(obj->GetName())=="hgpt" || string(obj->GetName())=="hpt_jet" )) {
            cout << "+" << flush;
            
            _jet = TString(obj->GetName()).Contains("hpt_jet");
            TH1D *hpt = (TH1D*)obj;
            //TH1D *hpt2 = (TH1D*)indir2->Get("hnlo"); assert(hpt2);
            //    TH1D *hpt2 = (TH1D*)indir2->Get(jp::dagfile1 ? "mc/hpt_g" : "hgpt"); assert(hpt2);
            TH1D *hpt2 = (TH1D*)indir2->Get(jp::dagfile1 ? "mc/hgpt_g" : "hgpt"); assert(hpt2);
            if (hpt2)
                dagostiniUnfold_histo(hpt, hpt2, outdir, ismc);
        } //hgpt plots
    }else if(jp::isquark && !jp::isgluon){
        cout << "quark unfolding routine"<<endl;
        if (obj->InheritsFrom("TH1") && (string(obj->GetName())=="hqpt" || string(obj->GetName())=="hpt_jet" )) {
            cout << "+" << flush;
            
            _jet = TString(obj->GetName()).Contains("hpt_jet");
            TH1D *hpt = (TH1D*)obj;
            //TH1D *hpt2 = (TH1D*)indir2->Get("hnlo"); assert(hpt2);
            //    TH1D *hpt2 = (TH1D*)indir2->Get(jp::dagfile1 ? "mc/hpt_g" : "hgpt"); assert(hpt2);
            TH1D *hpt2 = (TH1D*)indir2->Get(jp::dagfile1 ? "mc/hqpt_g" : "hqpt"); assert(hpt2);
            if (hpt2)
                dagostiniUnfold_histo(hpt, hpt2, outdir, ismc);
        } //hqpt plots
    }else{
        cout << "all unfolding routine"<<endl;
        if (obj->InheritsFrom("TH1") && (string(obj->GetName())=="hpt" || string(obj->GetName())=="hpt_jet" )) {
            cout << "+" << flush;
            
            _jet = TString(obj->GetName()).Contains("hpt_jet");
            TH1D *hpt = (TH1D*)obj;
            //TH1D *hpt2 = (TH1D*)indir2->Get("hnlo"); assert(hpt2);
            //    TH1D *hpt2 = (TH1D*)indir2->Get(jp::dagfile1 ? "mc/hpt_g" : "hgpt"); assert(hpt2);
            TH1D *hpt2 = (TH1D*)indir2->Get(jp::dagfile1 ? "mc/hpt_g" : "hpt"); assert(hpt2);
            if (hpt2)
                dagostiniUnfold_histo(hpt, hpt2, outdir, ismc);
            } //hpt for inclusive for other cases 
        
    } // hpt plots

  } // while key

  curdir->cd();
} // recurseFile


void dagostiniUnfold_histo(TH1D *hpt, TH1D *hnlo, TDirectory *outdir,
			   bool ismc, bool kscale, string id) {

  float y1, y2;
  assert(sscanf(outdir->GetName(),"Eta_%f-%f",&y1,&y2)==2);
  sscanf(outdir->GetName(),"Eta_%f-%f",&y1,&y2);
  cout << outdir->GetName() << " y1:" << y1 << " y2: " << y2 << endl;
  const char *c = id.c_str();
  if (_jet) c = "_jet";

  _ismcjer = ismc;

  // initial fit of the NLO curve to a histogram
  TF1 *fnlo = new TF1(Form("fus%s",c),
		      //   "[0]*exp([1]/x)*pow(x,[2])"
		           "[0]*pow(x,[1])"
                      "*pow(1-x*cosh([3])/[4],[2])", //10., 1000.);
		      jp::unfptminnlo, min(jp::xmax, jp::emax/cosh(y1)));

  fnlo->SetParameters(5e10,-5.2,8.9,y1,jp::emax);
  //  fnlo->SetParameters(2e14*2e-10,-18,-5,10,y1,jp::emax);
  fnlo->FixParameter(3,y1);
  fnlo->FixParameter(4,jp::emax);

  //hnlo->Fit(fnlo,"QRN");
  //hnlo->Scale(2e-10); // TEMP PATCH
  fnlo->SetRange(max(60.,jp::unfptminnlo), min(jp::xmax, jp::emax/cosh(y1)));
  cout << "fit hnlo" << endl;
  hnlo->Fit(fnlo,"RN");  // There seems to be abnormal terminations
  fnlo->SetRange(jp::unfptminnlo, min(jp::xmax, jp::emax/cosh(y1)));

  // Graph of theory points with centered bins
  const double minerr = 0.02;
  TGraphErrors *gnlo = new TGraphErrors(0);
  TGraphErrors *gnlo2 = new TGraphErrors(0); // above + minerr
  gnlo->SetName("gnlo");
  gnlo2->SetName("gnlo2");
  for (int i = 1; i != hnlo->GetNbinsX()+1; ++i) {

    double y = hnlo->GetBinContent(i);
    double dy = hnlo->GetBinError(i);
    double ptmin = hnlo->GetBinLowEdge(i);
    double ptmax = hnlo->GetBinLowEdge(i+1);

    double y0 = fnlo->Integral(ptmin, ptmax) / (ptmax - ptmin);
    double x = fnlo->GetX(y0, ptmin, ptmax);

    int n = gnlo->GetN();
    tools::SetPoint(gnlo, n, x, y, 0, dy);
    tools::SetPoint(gnlo2, n, x, y, 0, tools::oplus(dy, minerr*y));
  }

  // Second fit to properly centered graph
  //gnlo2->Fit(fnlo,"QRN");
  fnlo->SetRange(max(60.,jp::unfptminnlo), min(jp::xmax, jp::emax/cosh(y1)));   // Fit ranges here...
  cout << "fit to gnlo2" << endl;
  gnlo2->Fit(fnlo,"RN");
  fnlo->SetRange(jp::unfptminnlo, min(jp::xmax, jp::emax/cosh(y1)));

  // Bin-centered data points
  TGraphErrors *gpt = new TGraphErrors(0);
  gpt->SetName(Form("gpt%s",c));
    for (int i = 1; i != hpt->GetNbinsX()+1; ++i) {
 
      double ptmin = hpt->GetBinLowEdge(i);
      double ptmax = hpt->GetBinLowEdge(i+1);
      double y = fnlo->Integral(ptmin, ptmax) / (ptmax - ptmin);
      double x = fnlo->GetX(y, ptmin, ptmax);
      double ym = hpt->GetBinContent(i);
      double ym_err = hpt->GetBinError(i);
      if (ym>0) {
	tools::SetPoint(gpt, gpt->GetN(), x, ym, 0., ym_err);
      }
  } // for i

  // Create smeared theory curve
  double maxpt = jp::emax/cosh(y1);
  cout << "y1 "<< y1 << " c "<< c <<endl<<flush;
  TF1 *fnlos = new TF1(Form("fs%s",c),smearedAnsatz,jp::unfptminnlo,maxpt,nk+3); 
  fnlos->SetParameters(y1, fnlo->GetParameter(0), fnlo->GetParameter(1),
                       fnlo->GetParameter(2), 0, 0);
  cout << "par0 "<< fnlos->GetParameter(0) << "y1 "<< y1<<endl<<flush; 

  cout << "Calculate forward smearing and unfold hpt" << endl << flush;

  TGraphErrors *gfold_fwd = new TGraphErrors(0);
  gfold_fwd->SetName(Form("gfold_fwd%s",c));
  TGraphErrors *gcorrpt_fwd = new TGraphErrors(0);
  gcorrpt_fwd->SetName(Form("gcorrpt_fwd%s",c));
  TH1D *hcorrpt_fwd = (TH1D*)hpt->Clone(Form("hcorrpt_fwd%s",c));

  //  for (int i = 0; i != gpt->GetN(); ++i) {
  for (int i = 1; i != gpt->GetN()+1; ++i) {
    double x, y, ex, ey;
    tools::GetPoint(gpt, i, x, y, ex, ey);
    double k = fnlo->Eval(x) / fnlos->Eval(x);

    if (!TMath::IsNaN(k)) {

      tools::SetPoint(gfold_fwd, gfold_fwd->GetN(), x, k, ex, 0.);
      tools::SetPoint(gcorrpt_fwd, gcorrpt_fwd->GetN(), x, k*y, ex, k*ey);

      int j = hpt->FindBin(x);
      hcorrpt_fwd->SetBinContent(j, k*hpt->GetBinContent(j));
      hcorrpt_fwd->SetBinError(j, k*hpt->GetBinError(j));
    }
  }

  // Calculate smearing matrix
  cout << "Generating smearing matrix T..." << endl << flush;

  double tmp_eps = _epsilon;
  _epsilon = 1e-6; // speed up calculations with acceptable loss of precision

  // NB: GetArray only works if custom x binning
  outdir->cd();

  // Deduce range and binning for true and measured spectra
  vector<double> vx; // true, gen
  vector<double> vy; // measured, reco
  for (int i = 1; i != hpt->GetNbinsX()+1; ++i) {

    double x = hpt->GetBinCenter(i);
    double x1 = hpt->GetBinLowEdge(i);
    double x2 = hpt->GetBinLowEdge(i+1);
    double y = hpt->GetBinContent(i); 

    if (x>=jp::unfptmingen && y>0) {
      if (vx.size()==0) vx.push_back(x1);
      vx.push_back(x2);
    }

      if (x>=jp::unfptminreco && y>0) {
       if (vy.size()==0) vy.push_back(x1);
       vy.push_back(x2);
    }
  }

  // copy over relevant part of hpt
  TH1D *hreco = new TH1D(Form("hreco%s",c),";p_{T,reco} (GeV)",  
			 vy.size()-1,&vy[0]);
  for (int i = 1; i != hreco->GetNbinsX()+1; ++i) {  // orig
    int j = hpt->FindBin(hreco->GetBinCenter(i));
    double dpt = hpt->GetBinWidth(j);
    hreco->SetBinContent(i, hpt->GetBinContent(j)*dpt);
    hreco->SetBinError(i, hpt->GetBinError(j)*dpt);
  }

  // copy over relevant part of hnlo
  TH1D *htrue = new TH1D(Form("htrue%s",c),";p_{T,gen} (GeV)",
			 vx.size()-1,&vx[0]);
 for (int i = 1; i != htrue->GetNbinsX()+1; ++i) { 
    int j = hnlo->FindBin(htrue->GetBinCenter(i));
    double dpt = hnlo->GetBinWidth(j);
    htrue->SetBinContent(i, hnlo->GetBinContent(j)*dpt);
    htrue->SetBinError(i, hnlo->GetBinError(j)*dpt);
  }

   TH2D *mt = new TH2D(Form("mt%s",c),"mt;p_{T,reco};p_{T,gen}",           // Construction: x, y
		      vy.size()-1, &vy[0],vx.size()-1, &vx[0]
		       );
  TH1D *mx = new TH1D(Form("mx%s",c),"mx;p_{T,gen};#sigma/dp_{T}",   
		      vx.size()-1, &vx[0]);   // "TRUE"
  TH1D *my = new TH1D(Form("my%s",c),"my;p_{T,reco};#sigma/dp_{T}",
		      vy.size()-1, &vy[0]);   // RECO

  double mtbinsX = mt->GetNbinsX();  double mtbinsY = mt->GetNbinsY();
  
  // From http://hepunx.rl.ac.uk/~adye/software/unfold/RooUnfold.html
  // For 1-dimensional true and measured distribution bins Tj and Mi,
  // the response matrix element Rij gives the fraction of events
  // from bin Tj that end up measured in bin Mi.

  for (int i = 1; i != mt->GetNbinsX()+1; ++i) {
   
    double ptreco1 = mt->GetXaxis()->GetBinLowEdge(i);
    double ptreco2 = mt->GetXaxis()->GetBinLowEdge(i+1);
    double yreco = fnlo->Integral(ptreco1, ptreco2) / (ptreco2 - ptreco1);
    double ptreco = fnlo->GetX(yreco, ptreco1, ptreco2);

    for (int j = 1; j != mt->GetNbinsY()+1; ++j) {

      double ptgen1 = min(jp::emax/cosh(y1), mt->GetYaxis()->GetBinLowEdge(j));
      double ptgen2 = min(jp::emax/cosh(y1), mt->GetYaxis()->GetBinLowEdge(j+1));

      if (ptgen1>=jp::unfptmingen && ptreco>jp::unfptminreco && ptgen1*cosh(y1)<jp::emax) {  // This results in rows and columns of 0 in mt

        fnlos->SetParameter(4, ptgen1);
        fnlos->SetParameter(5, ptgen2);
        // 2D integration over pTreco, pTgen simplified to 1D over pTgen
        mt->SetBinContent(i, j, fnlos->Eval(ptreco) * (ptreco2 - ptreco1));

	fnlos->SetParameter(4, 0);
        fnlos->SetParameter(5, 0);
      }
    } // for j
  } // for i

  for (int j = 1; j != mt->GetNbinsY()+1; ++j) {
  
    double ptgen1 = min(jp::emax/cosh(y1), mt->GetYaxis()->GetBinLowEdge(j));
    double ptgen2 = min(jp::emax/cosh(y1), mt->GetYaxis()->GetBinLowEdge(j+1));
    double ygen = fnlo->Integral(ptgen1, ptgen2);
    mx->SetBinContent(j, ygen);
  }

 for (int i = 1; i != mt->GetNbinsX()+1; ++i) {

    double yrecoprojection(0);
    for (int j = 1; j != mt->GetNbinsY()+1; ++j) {
    
      yrecoprojection += mt->GetBinContent(i, j);
    }
    my->SetBinContent(i, yrecoprojection);
  } // for i

  TH2D *mtu = (TH2D*)mt->Clone(Form("mtu%s",c));
  for (int i = 1; i != mt->GetNbinsX()+1; ++i) {
    for (int j = 1; j != mt->GetNbinsY()+1; ++j) {
        if (mx->GetBinContent(j)!=0) {
	mtu->SetBinContent(i, j, mt->GetBinContent(i,j) / mx->GetBinContent(j));
      }
    } // for j
  } // for i 

  // Check condition number (statcomm recommendation) of the response matrix mt (TDecompSVD)
  // cond(K) = sigma_max/max(0,sigma_min), sigmas are singular values of matrix K

  // TH2D to TMatrixD
  Int_t nbinstotal = mtbinsY*mtbinsX;
  TMatrixD *K = new TMatrixD(mtbinsY,mtbinsX); // K(rows, cols)
  TArrayD mtEntries(nbinstotal);

  for (Int_t j = 1; j <= mtbinsY; j++) {
    for (Int_t i = 1; i <= mtbinsX; i++) {
      mtEntries[(j-1)*mtbinsX+i-1] = mt->GetBinContent(i,j);
    }
  }

  K->SetMatrixArray(mtEntries.GetArray());
  K->Print();

  TDecompSVD *svd = new TDecompSVD(*K);
  svd->Decompose(); 

  // Get singular values
  TVectorD singulars = svd->GetSig();
  singulars.Print();


  // For BinByBin and SVD, need square matrix
  TH2D *mts(0);
  TH1D *mxs(0);

  mts = new TH2D(Form("mts%s",c),"mts;p_{T,reco};p_{T,gen}",
		   vy.size()-1, &vy[0], vy.size()-1, &vy[0]);
  mxs = new TH1D(Form("mxs%s",c),"mxs;p_{T,gen};#sigma/dp_{T}",
		   vy.size()-1, &vy[0]);

  for (int i = 1; i != mts->GetNbinsX()+1; ++i) {
    for (int j = 1; j != mts->GetNbinsY()+1; ++j) {
      double x = ((TAxis*)mts->GetXaxis())->GetBinCenter(i);
      double y = ((TAxis*)mts->GetYaxis())->GetBinCenter(j);
      int i2 = mt->GetXaxis()->FindBin(x);
      int j2 = mt->GetYaxis()->FindBin(y);
      mts->SetBinContent(i, j, mt->GetBinContent(i2, j2));
      mts->SetBinError(i, j, mt->GetBinError(i2, j2));
    }
  }

  for (int i = 1; i != mxs->GetNbinsX()+1; ++i) {
    double x = mxs->GetBinCenter(i);
    int i2 = mx->FindBin(x);
    mxs->SetBinContent(i, mx->GetBinContent(i2));
    mxs->SetBinError(i, mx->GetBinError(i2));
  }
  
  _epsilon = tmp_eps;
  if (jp::debug)
    cout << "done." << endl << flush;

  // Now to actual unfolding business with the d'Agostini method
  if (jp::debug)
    cout << "Unfolding..." << flush;

  // RooUnfoldResponse constructor - create from already-filled histograms
  // "response" (mt) gives the response matrix, measured X truth.
  // "measured" (my) and "truth" (mx) give the projections of "response"
  // onto the X-axis and Y-axis respectively,
  // but with additional entries in "measured" (my) for measurements with
  // no corresponding truth (fakes/background) [not implemented] and
  // in "truth" (mx) for unmeasured events (inefficiency) [is implemented].
  // "measured" and/or "truth" can be specified as 0 (1D case only)
  // or an empty histograms (no entries) as a shortcut
  // to indicate, respectively, no fakes and/or no inefficiency.
  // RooUnfoldResponse(const TH1* measured,
  //                   const TH1* truth, const TH2* response,
  //                   const char* name, const char* title)
  RooUnfoldResponse *uResp = new RooUnfoldResponse(my, mx, mt);

  // RooUnfoldBayes (const RooUnfoldResponse* res, const TH1* meas,
  //                 Int_t niter= 4, Bool_t smoothit= false,
  //                 const char* name= 0, const char* title= 0);
  RooUnfoldBayes *uBayes = new RooUnfoldBayes(uResp, hreco, 4);

// if (jp::debug)
    uBayes->Print();

  TH1D *hTrueBayes = (TH1D*)uBayes->Hreco(RooUnfold::kCovariance);
  assert(hTrueBayes);
  TMatrixD *mCov = new TMatrixD(uBayes->Ereco());
  assert(mCov);

  TH2D *hCov = new TH2D(Form("hCov%s",c), Form("hCov%s;p_{T};p_{T}",c),
			vx.size()-1, &vx[0], vx.size()-1, &vx[0]);
  assert(hCov->GetNbinsX()==mCov->GetNrows());
  assert(hCov->GetNbinsY()==mCov->GetNcols());
  for (int i = 1; i != hCov->GetNbinsX()+1; ++i) {
    for (int j = 1; j != hCov->GetNbinsY()+1; ++j) {
      hCov->SetBinContent(i, j, (*mCov)[i-1][j-1]);
    } // for j
  } // for i

  // copy into original binning for plotting macros to work (e.g. drawClosure)
  TH1D *hcorrpt_dag = (TH1D*)hpt->Clone(Form("hcorrpt_dag%s",c));
  hcorrpt_dag->Reset();
  for (int i = 1; i != hTrueBayes->GetNbinsX()+1; ++i) {

    int j = hpt->FindBin(hTrueBayes->GetBinCenter(i));
    hcorrpt_dag->SetBinContent(j, hTrueBayes->GetBinContent(i));
    hcorrpt_dag->SetBinError(j, hTrueBayes->GetBinError(i));
  }

  // RooUnfoldBinByBin (const RooUnfoldResponse* res, const TH1* meas,
  //                    const char* name=0, const char* title=0);
  //hreco = my;
  // BinByBin and SVD can only handle square matrix

  cout << "Moving to bin by bin" << endl;
  
  TH1D *hcorrpt_bin(0), *hcorrpt_svd(0);
  //  _jet = 0; // added
  if (!_jet) {

    RooUnfoldResponse *uResps = new RooUnfoldResponse(my, mxs, mts);
    RooUnfoldBinByBin *uBin = new RooUnfoldBinByBin(uResps, hreco);
    TH1D *hTrueBin = (TH1D*)uBin->Hreco(RooUnfold::kCovariance);
    assert(hTrueBin);

    hcorrpt_bin = (TH1D*)hpt->Clone(Form("hcorrpt_bin%s",c));
    hcorrpt_bin->Reset();
    for (int i = 1; i != hTrueBin->GetNbinsX()+1; ++i) {

      int j = hpt->FindBin(hTrueBin->GetBinCenter(i));
      hcorrpt_bin->SetBinContent(j, hTrueBin->GetBinContent(i));
      hcorrpt_bin->SetBinError(j, hTrueBin->GetBinError(i));
    }

    cout << "Moving to SVD" << endl; 
    
    int kreg = int(vy.size()/2);
    int ntoys = 0;// default 1000
    RooUnfoldSvd *uSVD = new RooUnfoldSvd(uResps, hreco, kreg, ntoys);
    TH1D *hTrueSVD = (TH1D*)uSVD->Hreco();//RooUnfold::kCovariance);
    assert(hTrueSVD);
    hcorrpt_svd = (TH1D*)hTrueSVD->Clone(Form("hcorrpt_svd%s",c));
    delete uSVD; // ensure static members get destroyed before next instance

    hcorrpt_svd = (TH1D*)hpt->Clone(Form("hcorrpt_svd%s",c));
    hcorrpt_svd->Reset();
    for (int i = 1; i != hTrueSVD->GetNbinsX()+1; ++i) {

      int j = hpt->FindBin(hTrueSVD->GetBinCenter(i));
      hcorrpt_svd->SetBinContent(j, hTrueSVD->GetBinContent(i));
      hcorrpt_svd->SetBinError(j, hTrueSVD->GetBinError(i));
      }
   }

  // TUnfold.

  // cout << "Moving to TUnfold" << endl;



  if (jp::debug)
   cout << "done." << endl << flush;

  // Store "unfolding correction" (unfolded / original) and corrected graph
  //cout << "Store graphs..." << endl << flush;
  outdir->cd();
  TGraphErrors *gfold_dag = new TGraphErrors(0);
  gfold_dag->SetName(Form("gfold_dag%s",c));
  TGraphErrors *gcorrpt_dag = new TGraphErrors(0);
  gcorrpt_dag->SetName(Form("gcorrpt_dag%s",c));
  //
  TGraphErrors *gfold_bin(0), *gcorrpt_bin(0);
  TGraphErrors *gfold_svd(0), *gcorrpt_svd(0);
  if (!_jet) {

    gfold_bin = new TGraphErrors(0);
    gfold_bin->SetName(Form("gfold_bin%s",c));
    gcorrpt_bin = new TGraphErrors(0);
    gcorrpt_bin->SetName(Form("gcorrpt_bin%s",c));
    gfold_svd = new TGraphErrors(0);
    gfold_svd->SetName(Form("gfold_svd%s",c));
    gcorrpt_svd = new TGraphErrors(0);
    gcorrpt_svd->SetName(Form("gcorrpt_svd%s",c));
  }

  // Normalize hcorrpt
  for (int i = 1; i != hcorrpt_dag->GetNbinsX()+1; ++i) {
    double dpt = hcorrpt_dag->GetBinWidth(i);
    hcorrpt_dag->SetBinContent(i, hcorrpt_dag->GetBinContent(i) / dpt);
    hcorrpt_dag->SetBinError(i, hcorrpt_dag->GetBinError(i) / dpt);
  }
  if (!_jet) {

    for (int i = 1; i != hcorrpt_bin->GetNbinsX()+1; ++i) {
      double dpt = hcorrpt_bin->GetBinWidth(i);
      hcorrpt_bin->SetBinContent(i, hcorrpt_bin->GetBinContent(i) / dpt);
      hcorrpt_bin->SetBinError(i, hcorrpt_bin->GetBinError(i) / dpt);
    }
    for (int i = 1; i != hcorrpt_bin->GetNbinsX()+1; ++i) {
      double dpt = hcorrpt_svd->GetBinWidth(i);
      hcorrpt_svd->SetBinContent(i, hcorrpt_svd->GetBinContent(i) / dpt);
      hcorrpt_svd->SetBinError(i, hcorrpt_svd->GetBinError(i) / dpt);
    }
  }

  for (int i = 0; i != gpt->GetN(); ++i) {

    double x, y, ex, ey;
    tools::GetPoint(gpt, i, x, y, ex, ey);

    int j = hcorrpt_dag->FindBin(x);
    double ycorr = hcorrpt_dag->GetBinContent(j);
    double eycorr = hcorrpt_dag->GetBinError(j);

    double k = (ycorr && y ? ycorr / y : 1.);

    if (!TMath::IsNaN(k)) {

      double dk = 0;
      //if (eycorr/ycorr > ey/y)
      dk = eycorr/ycorr*k;
      tools::SetPoint(gfold_dag, gfold_dag->GetN(), x, k, ex, dk);
      tools::SetPoint(gcorrpt_dag, gcorrpt_dag->GetN(), x, ycorr, ex, eycorr);
    }
  } // for i

  if (!_jet) {

    for (int i = 0; i != gpt->GetN(); ++i) {

      double x, y, ex, ey;
      tools::GetPoint(gpt, i, x, y, ex, ey);

      int j = hcorrpt_bin->FindBin(x);
      double ycorr_bin = hcorrpt_bin->GetBinContent(j);
      double eycorr_bin = hcorrpt_bin->GetBinError(j);
      double ycorr_svd = hcorrpt_svd->GetBinContent(j);
      double eycorr_svd = hcorrpt_svd->GetBinError(j);

      double k_bin = (ycorr_bin && y ? ycorr_bin / y : 1.);
      double k_svd = (ycorr_svd && y ? ycorr_svd / y : 1.);

      if (!TMath::IsNaN(k_bin)) {

	double dk = 0;
	//if (eycorr_bin/ycorr_bin > ey/y)
	dk = eycorr_bin/ycorr_bin*k_bin;
	tools::SetPoint(gfold_bin, gfold_bin->GetN(), x, k_bin, ex, dk);
	tools::SetPoint(gcorrpt_bin, gcorrpt_bin->GetN(),
			x, ycorr_bin, ex, eycorr_bin);
      }

      if (!TMath::IsNaN(k_svd)) {

	double dk = 0;
	//if (eycorr_svd/ycorr_svd > ey/y)
	dk = eycorr_svd/ycorr_svd*k_svd;
	tools::SetPoint(gfold_svd, gfold_svd->GetN(), x, k_svd, ex, dk);
	tools::SetPoint(gcorrpt_svd, gcorrpt_svd->GetN(),
			x, ycorr_svd, ex, eycorr_svd);
      }

    } // for i
  }

  outdir->cd();

  // Save resolution function
  TF1 *fres = new TF1(Form("fres%s",c), fPtRes, jp::xmin, jp::xmax, 1);
  fres->SetParameter(0, y1);

  // Store NLO ratio to (unsmeared) fit
  TGraphErrors *grationlo = new TGraphErrors(0);
  grationlo->SetName(Form("grationlo%s",c));
  for (int i = 0; i != gnlo2->GetN(); ++i) {

    double x, y, ex, ey;
    tools::GetPoint(gnlo2, i, x, y, ex, ey);
    double ys = fnlo->Eval(x);
    if (!TMath::IsNaN(ys))
      tools::SetPoint(grationlo, grationlo->GetN(), x, y / ys, ex, ey / ys);
  }

  // Store data ratio to (smeared) fit
  TGraphErrors *gratio = new TGraphErrors(0);
  gratio->SetName(Form("gratio%s",c));
  for (int i = 0; i != gpt->GetN(); ++i) {

    double x, y, ex, ey;
    tools::GetPoint(gpt, i, x, y, ex, ey);
    double ys = fnlos->Eval(x);
    if (!TMath::IsNaN(ys))
      tools::SetPoint(gratio, gratio->GetN(), x, y / ys, ex, ey / ys);
  }

  if (!_jet) {

    // Inputs and central method results
    hpt->Write("hpt");
    hcorrpt_dag->Write("hcorrpt");
    //hpt->Write("hcorrpt"); // TEMP PATCH bypass
    hnlo->Write("hnlo");
    gpt->Write("gpt");
    gcorrpt_dag->Write("gcorrpt");
    //gpt->Write("gcorrpt"); // TEMP PATCH bypass
    gnlo2->Write("gnlo");
    gfold_dag->Write("gfold");

    // Fit functions
    uResp->Write();
    fnlo->SetRange(jp::unfptminnlo, min(jp::xmax, jp::emax/cosh(y1)));
    fnlo->SetNpx(1000); // otherwise ugly on log x-axis after write
    fnlo->Write();
    // Calculating points for fs is taking significant time,
    // and even got stuck at some point when too far out of the range
    fnlos->SetRange(jp::unfptminnlo, min(jp::xmax, jp::emax/cosh(y1)));
    fnlos->SetNpx(1000); // otherwise ugly on log x-axis after write
    fnlos->Write();
    fres->Write();
    grationlo->Write();
    gratio->Write();

    // Unfolding matrices and vectors
    hreco->Write();
    my->Write(); // my: hreco from fit
    mx->Write(); // mx: htrue from fit
    htrue->Write();
    mt->Write();
    mtu->Write();
    
    // Alternative methods
    gfold_dag->Write();
    gfold_fwd->Write();
    gfold_bin->Write();
    gfold_svd->Write();

    gcorrpt_dag->Write();
    gcorrpt_fwd->Write();
    gcorrpt_bin->Write();
    gcorrpt_svd->Write();

    hcorrpt_dag->Write();
    hcorrpt_fwd->Write();
    hcorrpt_bin->Write();
    hcorrpt_svd->Write();

    // Unfolding covariance matrix
    hCov->Write();
  }
  
}