// File:     match.cpp
// Author:   Zdenek Hubacek, Jiri Kvita
// Modified: $Id: match.cpp,v 1.4 2007/01/27 12:51:48 kvita Exp $ Jiri Kvita, August, September, October 2006 - January 2007
// Modified: Zdenek Hubacek, August 2006

#include <fstream>
#include <sstream>
#include <string>
#include <map>
#include <vector>
#include <TString.h>
#include <iostream>
#include "AnglesUtil.hpp"
#include "TH1D.h"
#include "TH2D.h"
#include "TFile.h"

#include "JetDefs.hpp"
#include "OffsetDeltaFitPars.hpp"
#include "kOkR_FitPars.hpp"
#include "kOkR_FitPars_Utils.hpp"
#include "Binnings.hpp"
#include "TFlags.hpp"

using namespace std;

const Float_t kError_RMPF     = -999;
const Float_t kOnePlusEpsilon = 1.001;
const Float_t kEpsilon        = 1.e-4;
const Float_t kNAverTowers    = 45;
const Float_t kDummyPVz       = 999.;

const int kNFineBins = 2000;
const Double_t kHighestE = 2000.;


// ___________________________________________________

struct zh4vec
{
  float E;
  float px;
  float py;
  float pz;
  float det_eta;
  float det_phi;
};

struct zh4jet
{
  float E;
  float px;
  float py;
  float pz;
  float det_eta;
  float det_phi;
  float offset;
 //  int jet_nitems;
};

struct event_struct_ZB
{
  std::string evtinfo;
  zh4vec photon;
  vector<zh4jet> jets;
  int njets;
  unsigned char pv_size;
  //  float pvz;
  float overlay_instlum;
  //  unsigned char overlay_tick;
  float weight;
  //  float offset;
  //  float met;
  //  int jet_nitems;
  float Response;
  float OrigGamma_calPt;
};

// to save memory:
struct event_struct_noZB
{
  std::string evtinfo;
  zh4vec photon;
  zh4vec jet;
  int njets;
  //  short int pv_size;
  //  float pvz;
  //  float overlay_instlum;
  //  int overlay_tick;
  float weight;
  //  float offset;
  //  float met;
  //  int jet_nitems;
  float Response;
};

// ___________________________________________________

bool PassedPhotonCuts(zh4vec photon, bool doPtPhotonCuts, bool doEtaPhotonCuts)
{

  if (doEtaPhotonCuts)
    if (fabs(photon.det_eta) > 1.)
      return false;

  if (doPtPhotonCuts)
    if ( sqrt(pow(photon.px, 2) + pow(photon.py, 2)) < 7.)
      return false;
  
  return true;
}


// ___________________________________________________

void DumpUsage()
{
  cout << "+---------------------------------------------------------------------+" << endl;
  cout << "|                                                                     |" << endl;
  cout << "| Usage: ./match_x   [ -debug N ]                                     |" << endl;
  cout << "|                                                                     |" << endl;
  cout << "|                  -h -help --help ... prints this help               |" << endl;
  cout << "|                                                                     |" << endl;
  cout << "|                  -doPhotonCuts                                      |" << endl;
  cout << "|                  -doPhotonEtaCuts                                   |" << endl;
  cout << "|                  -doPhotonID [ID=0..3] ]                            |" << endl;
  cout << "|    note: photon pT cut of 7 GeV now default!                        |" << endl;
  cout << "|                                                                     |" << endl;
  cout << "|                  -CorrectOffset                                     |" << endl;
  cout << "|                  -Correct_kOkR                                      |" << endl;
  cout << "|                  -InclusivePVbins ... >=4 PV's as last bin for kO   |" << endl;
  cout << "|                                                                     |" << endl;
  cout << "|                  -doTickCut [tick number within 1st superB)         |" << endl;
  cout << "|                  -doEtaPos  look only at positive jet eta's         |" << endl;
  cout << "|                  -doEtaNeg  look only at negative jet eta's         |" << endl;
  cout << "|                                                                     |" << endl;
  cout << "|                  -AskQuark  ... require a quark as |f> parton       |" << endl;
  cout << "|                  -AskGluon  ... require a gluon as |f> parton       |" << endl;
  cout << "|                                                                     |" << endl;
  cout << "|                  -DeltaRCut   [default 0.025]                       |" << endl;
  cout << "|                  -DeltaPhiCut [default 2.8]                         |" << endl;
  cout << "|                                                                     |" << endl;
  cout << "|                  -zbfile   filename (default: mc_zb.txt)            |" << endl;
  cout << "|                  -nozbfile filename (default: mc_nozb.txt)          |" << endl;
  cout << "|                  -outfile  filename (default: output.root)          |" << endl;
  cout << "|                                                                     |" << endl;
  cout << "|   Possible sensible comparisons:                                    |" << endl;
  cout << "|                                                                     |" << endl;
  cout << "|   -zbfile mc_unsuppzb_JCCA.txt -nozbfile mc_nozb_JCCA.txt           |" << endl;
  cout << "|   -zbfile mc_zb_JCCA.txt       -nozbfile mc_nozb_JCCA.txt           |" << endl;
  cout << "|                                                                     |" << endl;
  cout << "|   -zbfile mc_unsuppzb_JCCB.txt -nozbfile mc_nozb_JCCB.txt           |" << endl;
  cout << "|   -zbfile mc_zb_JCCB.txt       -nozbfile mc_nozb_JCCB.txt           |" << endl;
  cout << "|                                                                     |" << endl;
  cout << "|   -zbfile mc_zb_Loose.txt -nozbfile mc_nozb.txt                     |" << endl;
  cout << "|                                                                     |" << endl;
  cout << "|   (if not explicit, it is JCCA)                                     |" << endl;
  cout << "|                                                                     |" << endl;
  cout << "|                                                                     |" << endl;
  cout << "| Required photon ID must be 0,1,2,3 corresponding to                 |" << endl;
  cout << "|                        Prelim, Loose, Medium, Tight.                |" << endl;
  cout << "|                                                                     |" << endl;
  cout << "|                                                                     |" << endl;
  cout << "|                   Have fun!                                         |" << endl;
  cout << "|                               Jiri and Zdenek (c) 2006              |" << endl;
  cout << "|                                                                     |" << endl;
  cout << "+---------------------------------------------------------------------+" << endl;
}


// ___________________________________________________
// ___________________________________________________
// ___________________________________________________

int main( int argc, char *argv[])
{
  
  cout << "+---------------------------------------------------------------------+" << endl;
  cout << "|            Running the executable match_x!                          |" << endl;
  cout << "+---------------------------------------------------------------------+" << endl;

  TFlags Flags(argc, argv);

  if (Flags.Contain("-help") || Flags.Contain("-h") || Flags.Contain("--help") ) {
    DumpUsage();
    return 0;
  }

  // !!!!!!!!! DEBUG DEBUG DEBUG !!!!!!!!
  //  bool _UseAllJets = Flags.Contain("-UseAllJets");
  bool _UseAllJets  = true;
  bool _doPVAverage = true;
  //  bool _doPVAverage = Flags.Contain("-doPVAverage");
  cout << "!!! WARNING: _UseAllJets and _doPVAverage manually set to true !!!" << endl;

  bool _do_inclusive_nPV = Flags.Contain("-InclusivePVbins");


  bool _CorrectOffset  = Flags.Contain("-CorrectOffset");
  if ( _CorrectOffset )
    cout << "OK, will correct Ezb-O with fitted k = 1 + Delta." << endl;

  bool _Correct_kOkR  = Flags.Contain("-Correct_kOkR");
  if ( _Correct_kOkR )
    cout << "OK, will correct for kR/kO" << endl;


  float _DeltaPhiCut = 2.8;
  if ( Flags.GetFloat("-DeltaPhiCut", &_DeltaPhiCut) ) {
    cout << "OK, will use user defined DeltaPhi cut to match jets DeltaPhi=" << _DeltaPhiCut << endl;
  } else {
    cout << "OK, will use default DeltaPhi cut to match jets DeltaPhi=" << _DeltaPhiCut << endl;
  }


  float _InstLumiMin = 0.;
  if ( Flags.GetFloat("-InstLumiMin", &_InstLumiMin) ) {
    cout << "OK, will use user defined cut InstLumiMin=" << _InstLumiMin << endl;
  } else {
    cout << "OK, will use default cut InstLumiMin=" << _InstLumiMin << endl;
  }

  float _InstLumiMax = 999.;
  if ( Flags.GetFloat("-InstLumiMax", &_InstLumiMax) ) {
    cout << "OK, will use user cut defined InstLumiMax=" << _InstLumiMax << endl;
  } else {
    cout << "OK, will use default cut InstLumiMax=" << _InstLumiMax << endl;
  }


  bool _doEtaPos    = Flags.Contain("-doEtaPos");
  bool _doEtaNeg    = Flags.Contain("-doEtaNeg");

  bool _doPhotonCuts    = Flags.Contain("-doPhotonCuts");
  bool _doPtPhotonCuts  = true;
  bool _doEtaPhotonCuts;
  
  if (_doPhotonCuts) {
    _doPtPhotonCuts  = true;
    _doEtaPhotonCuts = true;
    
  } else {
    //    _doPtPhotonCuts  = Flags.Contain("-doPhotonPtCuts");
    _doEtaPhotonCuts = Flags.Contain("-doPhotonEtaCuts");
    if (_doEtaPhotonCuts || _doPtPhotonCuts)
      _doPhotonCuts = true;
  }

  int _doTickCut = -1;
  if ( Flags.GetIntWithinLimits("-doTickCut", 7, 40, &_doTickCut) ) {
    cout << "OK, we will require overlay tick " << _doTickCut << " plus symmetrical ones in next superbunches." << endl;
    if ( (_doTickCut < 7 || _doTickCut > 40 ) || ( (_doTickCut - 7) % 3) ) {
      cerr << "Required tick " << _doTickCut << " must be 7, 10, ... 40, corresponding to a tick in the first superbunch and symmetrically in next two!" << endl;
      return 1;
    }
  } else {
    if ( Flags.Contain("-doTickCut") ) {
      cerr << "ERROR: did not get tick cut within limits! Will not cut on it!" << endl;
      cerr << "Required tick must be 7, 10, ... 40, corresponding to a tick in the first superbunch and symmetrically in next two!" << endl;
      return 1;
    }
  }

  int _doPhotonID    = -1;
  if(  Flags.GetIntWithinLimits("-doPhotonID", 0, 3, &_doPhotonID) )
    cout << "OK, we will require photon ID " << _doPhotonID << endl;
  if (_doPhotonID > 3) {
    cerr << "Required photon ID " << _doPhotonID << " must be 0,1,2,3 corresponding to Prelim, Loose, Medium, Tight!" << endl;
  }
  

  bool _AskQuark = Flags.Contain("-AskQuark");
  bool _AskGluon = Flags.Contain("-AskGluon");


  int debug = 0;
  Flags.GetIntWithinLimits("-debug", 0, 1000, &debug );

  TString mc_nozb = "mc_nozb_JCCA.txt";
  TString mc_zb   = "mc_unsuppzb_JCCA.txt";

  char *mc_zb_ch   = Flags.GetParameter("-zbfile");
  char *mc_nozb_ch = Flags.GetParameter("-nozbfile");


  if ( mc_zb_ch == 0 ) {
    cout << "OK, will use default ZB filename " << mc_zb.Data() << endl;
  } else {
    cout << "OK, will use user specified ZB filename " << mc_zb_ch << endl;
    mc_zb = mc_zb_ch;
  }

  if ( mc_nozb_ch == 0 ) {
    cout << "OK, will use default noZB filename " << mc_nozb.Data() << endl;
  } else {
    cout << "OK, will use user specified NOZB filename " << mc_nozb_ch << endl;
    mc_nozb = mc_nozb_ch;
  }

 if (_UseAllJets) {
   mc_zb.ReplaceAll(TString("JCCA"), "JCCA_alljets");
   mc_zb.ReplaceAll(TString("JCCB"), "JCCB_alljets");
   mc_nozb.ReplaceAll(TString("JCCA"), "JCCA_alljets");
   mc_nozb.ReplaceAll(TString("JCCB"), "JCCB_alljets");

   cout << "WARNING: requested to use dumped all jets, changing input file names to "
	<< mc_zb.Data() << " and " << mc_nozb.Data() << endl;
 }
 

  int jetAlgo  = mc_zb.Contains("JCCA")   ? R2C7 : R2C5;
  int DataType = mc_zb.Contains("unsupp") ? DATA : PLATE;
  float _DeltaRCut = 0.025;
  if ( Flags.GetFloat("-DeltaRCut", &_DeltaRCut) ) {
    cout << "OK, will use user defined DeltaR cut to match jets DeltaR=" << _DeltaRCut << endl;
  } else {
    _DeltaRCut = (mc_zb.Contains("JCCA") ? 0.35 : 0.25);
    cout << "OK, will use default (jetAlgo dependent) DeltaR cut of Rcone/2 to match jets: DeltaR=" << _DeltaRCut << endl;
  }

  TString outputfilename = "output.root";
  char *outputfilename_ch = Flags.GetParameter("-outfile");
  if ( outputfilename_ch == 0 ) {
    cout << "OK, will use default output filename " << outputfilename.Data() << endl;
  } else {
    cout << "OK, will use user specified output filename " << outputfilename_ch << endl;
    outputfilename = outputfilename_ch;
  }

  //Read files
  // default:
  //  TString mc_zb = "mc_zb_Loose.txt";
  // no photon ID here:
  //TString mc_zb = "mc_zb.txt";
  // default:
  //  TString mc_nozb = "mc_nozb.txt";

  // old: (with automatic prelim ID, photon ID in both cases)
  //  TString mc_zb   = "mc_zb_old.txt";
  //  TString mc_nozb = "mc_nozb_old.txt";

  // unsuppressed ZB, new skim:
  //  TString mc_zb   = "mc_unsuppzb.txt";
  // unsuppressed ZB, old skim:
  //  TString mc_zb   = "mc_unsuppzb_old.txt";
  //  TString mc_nozb = "mc_nozb_old.txt";



  std::ifstream f1((const char*) mc_zb,ios::in);
  std::ifstream f2((const char*) mc_nozb, ios::in);

  if (!f1.is_open()) {
    cerr << "ERROR: could not open input ZB text file " << mc_zb.Data() << endl;
    return 1;
  }

  if (!f2.is_open()) {
    cerr << "ERROR: could not open input noZB text file " << mc_nozb.Data() << endl;
    return 1;
  }


  //variables to read
  std::string key;
  zh4vec photon;
  zh4jet zbjet;
  zh4vec nozbjet;
  int njets;
  int pv_size;
  float pvz;
  float overlay_instlum;
  int overlay_tick;
  float weight;
  float offset;
  float dummoffset;
  float met;
  int jet_nitems;
  float _Response;
  float _OrigGamma_calPt;
  std::string dummy;

  //  bool isPrelim, isLoose, isMedium, isTight;
  int PhotonID[4];
  //  bool PhotonID[4];

  // from Mikko's processor:
  bool  hasQuark;
  bool  hasGluon;

  //std structures
  std::vector<string> Keys;
  std::vector<string> Keys2;
  std::map<string, event_struct_ZB> ZB;
  std::map<string, event_struct_noZB> NoZB;


  //Histograms 
  //Fill with weights

  TH1D *InstLumi_h1 = new TH1D("InstLumi","InstLumi",100, 0, 3.);

  TH1D *h_besti = new TH1D("h_besti","Best matched jet index",21, -0.5, 20.5);
  h_besti -> Sumw2();

  TH1D *h_dr_all = new TH1D("dr_all","All ZBJets-noZBJet #DeltaR",1000, 0, 10);
  h_dr_all -> Sumw2();

  TH1D *h_dr = new TH1D("dr_best","Best ZBJet-noZBJet #DeltaR",1000, 0, 10);
  h_dr -> Sumw2();

  TH1D *h_dr_drCut =  new TH1D("dr_best_drCut","Best ZBJet-noZBJet #DeltaR, #DeltaR cut",1000, 0, 10);
  h_dr_drCut -> Sumw2();

  TH1D *h_dr_drCut_dPhiCut =  new TH1D("dr_best_drCut_dPhiCut","Best ZBJet-noZBJet #DeltaR, #DeltaR, #Delta#phi cut",1000, 0, 10);
  h_dr_drCut_dPhiCut -> Sumw2();

  TH1D *h_dphi = new TH1D("dphi","noZB truth #gamma - noZB jet #Delta#phi",100, 0, TMath::Pi());
  h_dphi -> Sumw2();
  TH1D *h_dphi_cut = new TH1D("dphi_cut","noZB truth #gamma - noZB jet #Delta#phi, after cut",100, 0, TMath::Pi());
  h_dphi_cut -> Sumw2();

  TH1D *npv = new TH1D("npv","npv",15, 0, 15);
  npv -> Sumw2();
  TH1D *h_jetpt = new TH1D("jetpt","jetpt",100, 0, 1000);
  h_jetpt -> Sumw2();
  TH1D *filled = new TH1D("filled","filled",300000, 0, 300000);
  filled -> Sumw2();
  TH1D *jetdeta = new TH1D("jetdeta","jedeta",100, -5,5);
  jetdeta -> Sumw2();
  TH1D *ep = new TH1D("eprime","eprime",kNFineBins, 0, kHighestE);
  ep -> Sumw2();


  TH1D *corr[nvertex][neprime][ndeta];
  TH1D *corr2[nvertex][neprime][ndeta];
  TH1D *corr_RMPF[nvertex][neprime][ndeta];
  TH1D *weights[nvertex][neprime][ndeta];
  TH1D *weights2[nvertex][neprime][ndeta];
  TH1D *ehist[nvertex][neprime][ndeta];

  TH1D *corr_pTgamma[nvertex][npTgamma][ndeta];
  TH1D *corr2_pTgamma[nvertex][npTgamma][ndeta];
  TH1D *corr_RMPF_pTgamma[nvertex][npTgamma][ndeta];
  TH1D *weights_pTgamma[nvertex][npTgamma][ndeta];
  TH1D *weights2_pTgamma[nvertex][npTgamma][ndeta];
  TH1D *ehist_pTgamma[nvertex][npTgamma][ndeta];

  TH1D *corr_pTprime[nvertex][npTgamma][ndeta];
  TH1D *corr2_pTprime[nvertex][npTgamma][ndeta];
  TH1D *corr_RMPF_pTprime[nvertex][npTgamma][ndeta];
  TH1D *weights_pTprime[nvertex][npTgamma][ndeta];
  TH1D *weights2_pTprime[nvertex][npTgamma][ndeta];
  TH1D *ehist_pTprime[nvertex][npTgamma][ndeta];

  TH1D *corr_ejet[nvertex][nejet][ndeta];
  TH1D *corr2_ejet[nvertex][nejet][ndeta];
  TH1D *corr_RMPF_ejet[nvertex][nejet][ndeta];
  TH1D *weights_ejet[nvertex][nejet][ndeta];
  TH1D *weights2_ejet[nvertex][nejet][ndeta];
  TH1D *ehist_ejet[nvertex][nejet][ndeta];

  TH1D *corr3[nvertex][neprime][ndeta];
  TH1D *weights3[nvertex][neprime][ndeta];
  TH1D *corr4[nvertex][neprime][ndeta]; //Scale offset by the size of the jet!

  // averages over samples:
  // R_MPF:
  TH2D *response_h2_pTprime_noZB[nvertex][ndeta];
  TH2D *response_h2_pTprime_ZB[nvertex][ndeta];
  TH2D *response_h2_eprime_noZB[nvertex][ndeta];
  TH2D *response_h2_eprime_ZB[nvertex][ndeta];
  // kO numerator and denominator:
  TH2D *kO_numerator_h2_pTprime[nvertex][ndeta];
  TH2D *kO_denominator_h2_pTprime[nvertex][ndeta];
  TH2D *kO_numerator_h2_eprime[nvertex][ndeta];
  TH2D *kO_denominator_h2_eprime[nvertex][ndeta];

  TH1D *corrfac = new TH1D("corrfac_tot","corrfact_tot",kNFineBins, -100, 100);

  TH1D *ntow[nvertex][neprime][ndeta];

  cout << "Preparing histograms..." << endl;

  for (int i = 1; i <= nvertex; i++)
    {
      for (int j = 0; j < ndeta; j++)
	{
	  for (int k = 0; k < neprime; k++)
	    {
	      TString name="corr_fac_";
	      name += i;
	      name += "PV_detetabin";
	      name += j;
	      name += "eprime";
	      name += k;
	      if (debug > 2)
		cout << "Creating histo " << name.Data() << endl;
	      corr[i-1][k][j]     = new TH1D(name, name, kNFineBins, -100, 100);
	      corr[i-1][k][j]     -> Sumw2();

	      name.ReplaceAll("corr_fac","weights");
	      if (debug > 2)
		cout << "Creating histo " << name.Data() << endl;
	      weights[i-1][k][j]  = new TH1D(name, name, 10, 0, 10);
	      weights[i-1][k][j]  -> Sumw2();
	      name.ReplaceAll("weights","sasa");
	      corr2[i-1][k][j]    = new TH1D(name, name, kNFineBins, -1, 1);
	      corr2[i-1][k][j]    -> Sumw2();

	      name.ReplaceAll("sasa","weights2");
	      if (debug > 2)
		cout << "Creating histo " << name.Data() << endl;
	      weights2[i-1][k][j] = new TH1D(name, name, 10, 0, 10);
	      weights2[i-1][k][j] -> Sumw2();
	      name.ReplaceAll("weights2","weights3");
	      //weights3[i-1][k][j] = new TH1D(name, name, 10, 0, 10);
	      //weights3[i-1][k][j] -> Sumw2();

	      name.ReplaceAll("weights3","E");
	      if (debug > 2)
		cout << "Creating histo " << name.Data() << endl;
	      corr3[i-1][k][j]    = new TH1D(name, name, kNFineBins, -1, 1);
	      corr3[i-1][k][j]    -> Sumw2();

	      name.ReplaceAll("E","Scale");
	      if (debug > 2)
		cout << "Creating histo " << name.Data() << endl;
	      corr4[i-1][k][j]    = new TH1D(name, name, kNFineBins, -100, 100);
	      corr4[i-1][k][j]    -> Sumw2();

	      // Response:
	      name.ReplaceAll("Scale", "RMPF");
	      if (debug > 2)
		cout << "Creating histo " << name.Data() << endl;
	      corr_RMPF[i-1][k][j]    = new TH1D(name, name, kNFineBins, -10, 10);
	      corr_RMPF[i-1][k][j]    -> Sumw2();


	      TString ehn="eprime_";
	      ehn += i;
	      ehn += "PV_detetabin";
	      ehn += j;
	      ehn += "eprime";
	      ehn += k;
	      if (debug > 2)
		cout << "Creating histo " << ehn.Data() << endl;
	      ehist[i-1][k][j] = new TH1D(ehn, ehn, kNFineBins, 0, kHighestE);
	      ehist[i-1][k][j] -> Sumw2();

	      TString ntot="Ntow_";
	      ntot += i;
	      ntot += "PV_detetabin";
	      ntot += j;
	      ntot += "eprime";
	      ntot += k;
	      if (debug > 2)
		cout << "Creating histo " << ntot.Data() << endl;
	      ntow[i-1][k][j] = new TH1D(ntot, ntot, kNFineBins, 0, kHighestE);
	      ntow[i-1][k][j] -> Sumw2();

	    } // eprime

	  TString name = "Response_";
	  name += i;
	  name += "PV_detetabin";
	  name += j;
	  TString title = "R_{MPF}, nPV=";
	  title += i;
	  title += ", etabin ";
	  title += j;

	  if (debug > 2)
	    cout << "Creating Response histos for nPV=" << i << endl;

	  response_h2_eprime_noZB[i-1][j] = new TH2D(name + "Eprime_noZB", title + ", noZB;E';R_{MPF}", 
						     neprime-1, eprimebins_forResponse, kNFineBins, -10., 10.);
	  response_h2_eprime_ZB[i-1][j]   = new TH2D(name + "Eprime_ZB", title + ", ZB;E';R_{MPF}", 
						     neprime-1, eprimebins_forResponse, kNFineBins, -10., 10.);
	  response_h2_eprime_noZB[i-1][j] -> Sumw2();
	  response_h2_eprime_ZB[i-1][j]   -> Sumw2();

	  response_h2_pTprime_noZB[i-1][j] = new TH2D(name + "pTprime_noZB", title + ", noZB;p_{T}';R_{MPF}", 
						      npTgamma-1, pTgammabins_forResponse, kNFineBins, -10., 10.);
	  response_h2_pTprime_ZB[i-1][j]   = new TH2D(name + "pTprime_ZB", title + ", ZB;p_{T}';R_{MPF}", 
						      npTgamma-1, pTgammabins_forResponse, kNFineBins, -10., 10.);
	  response_h2_pTprime_noZB[i-1][j] -> Sumw2();
	  response_h2_pTprime_ZB[i-1][j]   -> Sumw2();

	  // kO numerator and denominator, averaged over events:
	  name = "_";
	  name += i;
	  name += "PV_detetabin";
	  name += j;
	  title = ", nPV=";
	  title += i;
	  title += ", etabin ";
	  title += j;

	  if (debug > 2)
	    cout << "Creating kO numerator and denominator histos for nPV=" << i << endl;

	  kO_numerator_h2_pTprime[i-1][j] = new TH2D(TString("kO_numerator") + name + "pTprime", 
						     TString("k_{O} numerator") + title + ";p_{T}';k_{O} numerator", 
						     npTgamma-1, pTgammabins_forResponse, kNFineBins, -100., kHighestE);
	  kO_numerator_h2_pTprime[i-1][j] -> Sumw2();
	  kO_denominator_h2_pTprime[i-1][j] = new TH2D(TString("kO_denominator") + name + "pTprime", 
						       TString("k_{O} denominator") + title + ";p_{T}';k_{O} denominator", 
						       npTgamma-1, pTgammabins_forResponse, kNFineBins, -100., kHighestE);
	  kO_denominator_h2_pTprime[i-1][j] -> Sumw2();
	  kO_numerator_h2_eprime[i-1][j] = new TH2D(TString("kO_numerator") + name + "Eprime", 
						    TString("k_{O} numerator") + title + ";E';k_{O} numerator", 
						    neprime-1, eprimebins_forResponse, kNFineBins, -100., kHighestE);
	  kO_numerator_h2_eprime[i-1][j] -> Sumw2();
	  kO_denominator_h2_eprime[i-1][j] = new TH2D(TString("kO_denominator") + name + "Eprime", 
						      TString("k_{O} denominator") + title + ";E';k_{O} denominator", 
						      neprime-1, eprimebins_forResponse, kNFineBins, -100., kHighestE);
	  kO_denominator_h2_eprime[i-1][j] -> Sumw2();

	  // corrf. histos:

	  for (int k = 0; k < npTgamma; k++)
	    {
	      TString name="corr_fac_";
	      name += i;
	      name += "PV_detetabin";
	      name += j;
	      name += "pTgamma";
	      name += k;
	      if (debug > 2)
		cout << "Creating histo " << name.Data() << endl;
	      corr_pTgamma[i-1][k][j]     = new TH1D(name, name, kNFineBins, -100, 100);
	      corr_pTgamma[i-1][k][j]     -> Sumw2();
	      name.ReplaceAll("corr_fac","weights");
	      if (debug > 2)
		cout << "Creating histo " << name.Data() << endl;
	      weights_pTgamma[i-1][k][j]  = new TH1D(name, name, 10, 0, 10);
	      weights_pTgamma[i-1][k][j]  -> Sumw2();

	      name.ReplaceAll("weights","sasa");
	      if (debug > 2)
		cout << "Creating histo " << name.Data() << endl;
	      corr2_pTgamma[i-1][k][j]    = new TH1D(name, name, kNFineBins, -1, 1);
	      corr2_pTgamma[i-1][k][j]    -> Sumw2();

	      name.ReplaceAll("sasa","weights2");
	      if (debug > 2)
		cout << "Creating histo " << name.Data() << endl;
	      weights2_pTgamma[i-1][k][j] = new TH1D(name, name, 10, 0, 10);
	      weights2_pTgamma[i-1][k][j] -> Sumw2();

	      // Response:
	      name.ReplaceAll("weights2", "RMPF");
	      if (debug > 2)
		cout << "Creating histo " << name.Data() << endl;
	      corr_RMPF_pTgamma[i-1][k][j]    = new TH1D(name, name, kNFineBins, -10, 10);
	      corr_RMPF_pTgamma[i-1][k][j]    -> Sumw2();

	      TString ehn="pTgamma_";
	      ehn += i;
	      ehn += "PV_detetabin";
	      ehn += j;
	      ehn += "pTgamma";
	      ehn += k;
	      if (debug > 2)
		cout << "Creating histo " << ehn.Data() << endl;
	      ehist_pTgamma[i-1][k][j] = new TH1D(ehn, ehn, kNFineBins, 0, kHighestE);
	      ehist_pTgamma[i-1][k][j] -> Sumw2();

	      // pTprime:
	      name="corr_fac_";
	      name += i;
	      name += "PV_detetabin";
	      name += j;
	      name += "pTprime";
	      name += k;
	      if (debug > 2)
		cout << "Creating histo " << name.Data() << endl;
	      corr_pTprime[i-1][k][j]     = new TH1D(name, name, kNFineBins, -100, 100);
	      corr_pTprime[i-1][k][j]     -> Sumw2();
	      name.ReplaceAll("corr_fac","weights");
	      if (debug > 2)
		cout << "Creating histo " << name.Data() << endl;
	      weights_pTprime[i-1][k][j]  = new TH1D(name, name, 10, 0, 10);
	      weights_pTprime[i-1][k][j]  -> Sumw2();

	      name.ReplaceAll("weights","sasa");
	      if (debug > 2)
		cout << "Creating histo " << name.Data() << endl;
	      corr2_pTprime[i-1][k][j]    = new TH1D(name, name, kNFineBins, -1, 1);
	      corr2_pTprime[i-1][k][j]    -> Sumw2();

	      name.ReplaceAll("sasa","weights2");
	      if (debug > 2)
		cout << "Creating histo " << name.Data() << endl;
	      weights2_pTprime[i-1][k][j] = new TH1D(name, name, 10, 0, 10);
	      weights2_pTprime[i-1][k][j] -> Sumw2();

	      // Response:
	      name.ReplaceAll("weights2", "RMPF");
	      if (debug > 2)
		cout << "Creating histo " << name.Data() << endl;
	      corr_RMPF_pTprime[i-1][k][j]    = new TH1D(name, name, kNFineBins, -10, 10);
	      corr_RMPF_pTprime[i-1][k][j]    -> Sumw2();

	      ehn="pTprime_";
	      ehn += i;
	      ehn += "PV_detetabin";
	      ehn += j;
	      ehn += "pTprime";
	      ehn += k;
	      if (debug > 2)
		cout << "Creating histo " << ehn.Data() << endl;
	      ehist_pTprime[i-1][k][j] = new TH1D(ehn, ehn, kNFineBins, 0, kHighestE);
	      ehist_pTprime[i-1][k][j] -> Sumw2();

	    
	    } // pTgamma / pTprime

	  for (int k = 0; k < nejet; k++)
	    {
	      TString name="corr_fac_";
	      name += i;
	      name += "PV_detetabin";
	      name += j;
	      name += "ejet";
	      name += k;
	      if (debug > 2)
		cout << "Creating histo " << name.Data() << endl;
	      corr_ejet[i-1][k][j]     = new TH1D(name, name, kNFineBins, -100, 100);
	      corr_ejet[i-1][k][j]     -> Sumw2();

	      name.ReplaceAll("corr_fac","weights");
	      if (debug > 2)
		cout << "Creating histo " << name.Data() << endl;
	      weights_ejet[i-1][k][j]  = new TH1D(name, name, 10, 0, 10);
	      weights_ejet[i-1][k][j]  -> Sumw2();

	      name.ReplaceAll("weights","sasa");
	      if (debug > 2)
		cout << "Creating histo " << name.Data() << endl;
	      corr2_ejet[i-1][k][j]    = new TH1D(name, name, kNFineBins, -1, 1);
	      corr2_ejet[i-1][k][j]    -> Sumw2();

	      name.ReplaceAll("sasa","weights2");
	      if (debug > 2)
		cout << "Creating histo " << name.Data() << endl;
	      weights2_ejet[i-1][k][j] = new TH1D(name, name, 10, 0, 10);
	      weights2_ejet[i-1][k][j] -> Sumw2();

	      // Response:
	      name.ReplaceAll("weights2", "RMPF");
	      if (debug > 2)
		cout << "Creating histo " << name.Data() << endl;
	      corr_RMPF_ejet[i-1][k][j]    = new TH1D(name, name, kNFineBins, -10, 10);
	      corr_RMPF_ejet[i-1][k][j]    -> Sumw2();


	      TString ehn="ejet_";
	      ehn += i;
	      ehn += "PV_detetabin";
	      ehn += j;
	      ehn += "ejet";
	      ehn += k;
	      if (debug > 2)
		cout << "Creating histo " << ehn.Data() << endl;
	      ehist_ejet[i-1][k][j] = new TH1D(ehn, ehn, kNFineBins, 0, kHighestE);
	      ehist_ejet[i-1][k][j] -> Sumw2();
	    
	    } // ejet


	} // eta
    } // vertex

  cout << "Reading ZB ASCII file..." << endl;

  // +-----------------+
  // |        ZB       |
  // +-----------------+

  //       f1 >> dummy >> dummy
  for (int count; 
       f1 >> key >> photon.E >> photon.px >> photon.py >> photon.pz >> photon.det_eta >> photon.det_phi 
	 >> PhotonID[0] >> PhotonID[1] >> PhotonID[2] >> PhotonID[3]
	 >> _Response >> _OrigGamma_calPt
	 //	 >> hasQuark >> hasGluon
	 >> njets;
	
       )

    {

      // now read the rest:
      event_struct_ZB ev;

      for (int ijet = 0; ijet < njets; ++ijet) {
	f1 >> zbjet.E >> zbjet.px >> zbjet.py >> zbjet.pz >> zbjet.det_eta >> zbjet.det_phi >> zbjet.offset;
	ev.jets.push_back(zbjet);
      }
      f1 >> pv_size >> weight >> overlay_instlum >> dummoffset 
	 >> overlay_tick >> met >> jet_nitems;

      if (overlay_instlum > _InstLumiMax) {
	if (debug > 1) 
	  cout << "Failed maximal InstLumi cut of " << _InstLumiMax << " while lumi=" << overlay_instlum << endl;
	continue;
      }
      if ( overlay_instlum < _InstLumiMin ) {
	if (debug > 1) 
	  cout << "Failed minimal InstLumi cut of " << _InstLumiMin << " while lumi=" << overlay_instlum << endl;
	continue;
      }

      InstLumi_h1 -> Fill(overlay_instlum);

      // we cannot decide heer which jet is the right one...
      /*
	if ( _doEtaPos && zbjet.det_eta < 0.) {
	if (debug > 2) 
	cout << "Failed jet det eta positive cut! Jet eta=" << zbjet.det_eta << endl;
	continue;
	}
	
	if ( _doEtaNeg && zbjet.det_eta > 0.) {
	if (debug > 2) 
	cout << "Failed jet det eta negative cut! Jet eta=" << zbjet.det_eta << endl;
	continue;
	}
      */

      if (_doPhotonCuts)
	if (!PassedPhotonCuts(photon, _doPtPhotonCuts, _doEtaPhotonCuts)) {
	  if (debug > 2) 
	    cout << "Failed photon cut " 
		 << " _doPtPhotonCuts=" << _doPtPhotonCuts
		 << " _doEtaPhotonCuts=" << _doEtaPhotonCuts
		 << endl;
	  continue;
	}

      if (_doPhotonID >= 0) {
	if (! PhotonID[_doPhotonID]) {
	  if (debug > 2) 
	    cout << "Failed required photon ID " << _doPhotonID << " cut!" << endl;
	  continue;
	}
      }

      if (_doTickCut >= 0) {
	if ( overlay_tick != _doTickCut           && 
	     overlay_tick != _doTickCut +  60 - 7 && 
	     overlay_tick != _doTickCut + 113 - 7    ) {
	  if (debug > 2) 
	    cout << "Failed required overlay tick " << _doTickCut << " cut!" << endl;
	  continue;
	}
      }
      
      
      if (_AskQuark && !hasQuark) {
	if (debug > 2) 
	  cout << "Failed AskQuark cut!" << endl;
	continue;
      }

      if (_AskGluon && !hasGluon) {
	if (debug > 2) 
	  cout << "Failed AskGluon cut!" << endl;
	continue;

      }

      pvz = kDummyPVz;
      ev.photon = photon;
      //      ev.jet = jet;
      ev.evtinfo = key;
      ev.pv_size = pv_size;

      if (debug > 2)
	cout << "Read pv_size: " << pv_size << ", stored pv_size: " << (int)(ev.pv_size) << endl;

      //      ev.pvz = pvz;
      ev.overlay_instlum = overlay_instlum;
      //      ev.overlay_tick = overlay_tick;
      ev.weight = weight;
      //      ev.offset = offset;
      //      ev.met = met;
      //      ev.jet_nitems = jet_nitems;
      ev.Response = _Response;
      ev.OrigGamma_calPt = _OrigGamma_calPt;

      Keys.push_back(key);
      std::pair<string, event_struct_ZB> ppair(ev.evtinfo, ev);
      ZB.insert(ppair);
      if (Keys.size() > 0) 
	{
	  //Printout
	  if (debug > 2) {
	    cout << "ZB Event dump" << endl;
	    cout << "Key: " << ev.evtinfo << endl;
	    cout << "Photon: " << ev.photon.E << " " << ev.photon.px << " " << ev.photon.py << " " << ev.photon.pz << " " << ev.photon.det_eta << " " << ev.photon.det_phi << endl;
	    cout << "Response: " << ev.Response << endl;
	    for (int ijet = 0; ijet < ev.jets.size(); ++ijet)
	      cout << "Jet " << ijet << " :" << ev.jets[ijet].E << " " << ev.jets[ijet].px << " " << ev.jets[ijet].py << " " << ev.jets[ijet].pz << " " << ev.jets[ijet].det_eta << " " << ev.jets[ijet].det_phi << " offset: " << ev.jets[ijet].offset << endl;
	    //	    cout << "PV size: " << (int) ev.pv_size << " pv[0].z : " << ev.pvz << endl;
	    cout << "Event weight: " << ev.weight << endl;
	    cout << "Overlay instlum: " << ev.overlay_instlum << endl;
	    //	    cout << "Overlay tick#: " << (int) ev.overlay_tick << endl;
	    //	    cout << "Missing Et: " << ev.met << endl;
	    //	    cout << "Towers in jet: " << ev.jet_nitems << endl;
	  }
	}

    }

  cout << "Reading noZB ASCII file..." << endl;

  // +-----------------+
  // |      No ZB      |
  // +-----------------+

  // >> dummy >> dummy 
  for (int count; 
       f2 >> key >> photon.E >> photon.px >> photon.py >> photon.pz >> photon.det_eta >> photon.det_phi 
	 >> PhotonID[0] >> PhotonID[1] >> PhotonID[2] >> PhotonID[3]
	 >> _Response >> _OrigGamma_calPt
	 >> njets;
	
       )

    {

      // now read the rest:
      
      for (int ijet = 0; ijet < njets; ++ijet) {
	f2 >> nozbjet.E >> nozbjet.px >> nozbjet.py >> nozbjet.pz >> nozbjet.det_eta >> nozbjet.det_phi >> dummoffset;
      }
      f2 >> pv_size >> weight >> overlay_instlum >> dummoffset 
	 >> overlay_tick >> met >> jet_nitems;
      
      if (njets != 1) {
	if (debug > 2) 
	  cout << "Failed 1 jet only cut in ZB sample!" << endl;
      } else {

	//      pvz = kDummyPVz;
	event_struct_noZB ev;
	ev.photon = photon;
	ev.jet = nozbjet;
	ev.evtinfo = key;
	//      ev.pv_size = pv_size;
	//      ev.pvz = pvz;
	//      ev.overlay_instlum = overlay_instlum;
	//      ev.overlay_tick = overlay_tick;
	ev.weight = weight;
	//      ev.offset = offset;
	//      ev.met = met;
	//      ev.jet_nitems = jet_nitems;
	ev.Response = _Response;
	//      ev.OrigGamma_calPt = _OrigGamma_calPt;

	Keys2.push_back(key);
	std::pair<string, event_struct_noZB> ppair(ev.evtinfo, ev);
	NoZB.insert(ppair);
	if (Keys.size() > 0) 
	  {
	    //Printout
	    if (debug > 2) {
	      cout << "noZB Event dump" << endl;
	      cout << "Key: " << ev.evtinfo << endl;
	      cout << "Photon: " << ev.photon.E << " " << ev.photon.px << " " << ev.photon.py << " " << ev.photon.pz << " " << ev.photon.det_eta << " " << ev.photon.det_phi << endl;
	      cout << "Response: " << ev.Response << endl;
	      cout << "Jet: " << ev.jet.E << " " << ev.jet.px << " " << ev.jet.py << " " << ev.jet.pz << " " << ev.jet.det_eta << " " << ev.jet.det_phi << endl;
	      //	    cout << "PV size: " << ev.pv_size << " pv[0].z : " << ev.pvz << endl;
	      cout << "Event weight: " << ev.weight << endl;
	      //	    cout << "Overlay instlum: " << ev.overlay_instlum << endl;
	      //	    cout << "Overlay tick#: " << ev.overlay_tick << endl;
	      //	    cout << "Missing Et: " << ev.met << endl;
	      //	    cout << "Towers in jet: " << ev.jet_nitems << endl;
	    }
	  }
      } // njets == 1
    }

  //Test printout
  cout << "Keys: " << Keys.size() << " " << Keys2.size() << endl;
  cout << "Maps: " << ZB.size() << " " << NoZB.size() << endl;

  //  cout << "Doing the matching between samples..." << endl;

  /*
    //Find 3
    int ff =0;
    for (int i = 0; i< Keys.size(); i++)
    {
    if (ff <3)
    {
    std::map<string, event_struct_noZB>::iterator it = NoZB.find(Keys[i]);
    if (it != NoZB.end() )
    {
    std::map<string, event_struct_ZB>::iterator it2 = ZB.find(Keys[i]);
    
    if (it2 == ZB.end() )
    {
    cerr << "ERROR: Something wierd, did not ind matched events!" << endl;
    }
    if (debug > 2) {
    cout << "Got matched keys!" << endl;
    cout << "Key: " << it -> first << endl;
    cout << "Data: ZB Jet E " << it2 -> second.jet.E << endl;
    cout << "Data: NoZB Jet E " << it -> second.jet.E << endl;
    }
    ff++;	
    }
    }
    }
  */
  
  cout << "Doing the matching between samples, filling histograms..." << endl;

  //Now once more and try matching:


  int Progress = 10000;
  int found = 0;
  int good  = 0;
  int total = 0;
  for (int i = 0; i < Keys.size(); i++)
    {
      std::map<string, event_struct_noZB>::iterator it = NoZB.find(Keys[i]);
      if (it != NoZB.end() )
	{
	  std::map<string, event_struct_ZB>::iterator it2 = ZB.find(Keys[i]);
	  
	  if (it2 == ZB.end() )
	    {
	      cerr << "WARNING: Did not find a matched event in ZB sample!" << endl;
	    }
	  found++;

	  if (! (found % Progress))
	    cout << "Processing " << found << " of " 
		 << (Keys.size() < Keys2.size() ? Keys.size() : Keys2.size() ) << endl;

	  //   /d0dist/dist/packages/kinem_util/devel/kinem_util/AnglesUtil.hpp

	  // returns rapidity:
	  //  double kinem::eta(double x, double y, double z)
	  //  {
	  //    return 0.5*log((sqrt(x*x + y*y + z*z) + z + EPSILON) /
	  //		   (sqrt(x*x + y*y + z*z) - z + EPSILON));
	  //  }


	  // returns pseudorapidity:
	  // double kinem::eta(double th)
	  //{
	  //  if(th == 0) return ETA_LIMIT;
	  //  if(th >= PI-0.0001) return -ETA_LIMIT;
	  //  return -log(tan(th/2.0));
	  //}


	  zh4jet matchedZBjet;
 
	  float theta1 = kinem::theta((double) it  -> second.jet.px, (double) it  -> second.jet.py, (double) it  -> second.jet.pz);
	  float theta2;
	  float eta1 = kinem::eta(theta1);
	  float eta2;
	  float phi1 = kinem::phi((double) it  -> second.jet.px, (double) it  -> second.jet.py);
	  float phi2;
	  float photon_phi2 = kinem::phi((double) it2 -> second.photon.px, (double) it2 -> second.photon.py);
	  float photon_phi1 = kinem::phi((double) it -> second.photon.px, (double) it -> second.photon.py);
	  // will try to find the best dr:
	  float dr = 999.;
	  float bestdr = 9999.;
	  int besti = -1;
	  if (debug > 1) 
	    cout << " --- Finding the best dR match..." << endl;
	  
	  for (int ijet = 0; ijet < it2 -> second.jets.size(); ++ijet) {

	    zh4jet currentZBjet =  it2 -> second.jets[ijet];

	    //Do DR matching	
	    theta2 = kinem::theta((double) currentZBjet.px, (double) currentZBjet.py, (double) currentZBjet.pz);
	    eta2 = kinem::eta(theta2);
	    phi2 = kinem::phi((double) currentZBjet.px, (double) currentZBjet.py);

	  

	    dr = (float) kinem::delta_R((double) eta1, (double) phi1, (double) eta2, (double) phi2);
	    h_dr_all -> Fill(dr);
	    
	    if (debug > 1) {
	      cout << "ijet " << ijet << ", bestdr=" << dr << " dr=" << dr 
		   << " ZBjetpt=" << sqrt(pow(currentZBjet.px, 2) + pow(currentZBjet.py, 2))
		   << " noZBjetpt=" << sqrt(pow(it  -> second.jet.px, 2) + pow(it  -> second.jet.py, 2))
		   << endl;
	      cout << " phi1=" << phi1
		   << " phi2=" << phi2
		   << " eta1=" << eta1
		   << " eta2=" << eta2
		   << endl;
	    }
	    
	    if (dr < bestdr) {
	      matchedZBjet = currentZBjet;
	      bestdr = dr;
	      besti = ijet;
	    }

	    if (debug > 1) 
	      cout << "bestdr=" << dr << " dr=" << dr  
		   << " jetpt=" << sqrt(pow(matchedZBjet.px, 2) + pow(matchedZBjet.py, 2))
		   << " besti=" << besti << endl;

	    //test
	    //	  float rap1 = 0.5* log((double) (it -> second.jet.E+it -> second.jet.pz)/(it -> second.jet.E-it -> second.jet.pz));
	    //	  float rap2 = 0.5* log((double) (currentZBjet.E+currentZBjet.pz)/(currentZBjet.E-currentZBjet.pz));
	    //	  float dr2 = (float) kinem::delta_R((double) rap1, (double) phi1, (double) rap2, (double) phi2);

	  } // loop over ZB jets
	  dr = bestdr;
	  theta2 = kinem::theta((double) matchedZBjet.px, (double) matchedZBjet.py, (double) matchedZBjet.pz);
	  eta2 = kinem::eta(theta2);
	  phi2 = kinem::phi((double) matchedZBjet.px, (double) matchedZBjet.py);


	  // SORRY, N ITEMS KNOCKED OUT AT THE MOMENT!
	  //	  int nit = matchedZBjet_nitems;

	  int nit = 50;


	  //cout << "rap1: " << rap1 << " rap2: " << rap2 << endl;
	  //cout << "dr: " << dr << " dr2: " << dr2 << endl;
	  //cout << "dr ratio: " << dr2/dr << endl;
	  //Ok change to rapidity
	  //	  dr = dr2;

	  h_dr -> Fill(dr);
	  h_besti -> Fill(besti);


	  if (dr < _DeltaRCut )
	    {

	       h_dr_drCut -> Fill(dr);
	      // true ZB Photon - ZB matched jet dphi
	      // float dphi = (float) kinem::delta_phi(phi2, photon_phi2);

	      // true noZB Photon - noZB jet dphi
	      float dphi = (float) kinem::delta_phi(phi1, photon_phi1);
	      h_dphi -> Fill(dphi);
	      if (dphi > _DeltaPhiCut)
		{

		  h_dphi_cut -> Fill(dphi);
		  h_dr_drCut_dPhiCut -> Fill(dr);
		  good++;  

		  float e1 = it -> second.jet.E;  // NoZB
		  float e2 = matchedZBjet.E; // ZB

		  double kO_numerator   = 1.;
		  double kO_denominator = 1.;

		  float response1 = it  -> second.Response;  // NoZB
		  float response2 = it2 -> second.Response;  // ZB

		  float ejet = e1;

		  float offset = matchedZBjet.offset;
		  npv -> Fill((int) it2 -> second.pv_size); // Fill PV
		  float pTgamma = sqrt(pow(it2 -> second.photon.px, 2) + pow(it2 -> second.photon.py, 2));
		  float origGamma_calPt = it2 -> second.OrigGamma_calPt;
		  float jetpt  = sqrt(pow(matchedZBjet.px, 2) + pow(matchedZBjet.py, 2)); 
		  // float eprime  = pTgamma * cosh(kinem::eta(matchedZBjet.px, matchedZBjet.py, matchedZBjet.pz));
		  // new (with photon being corrected to particle level by jes_tools, we still want use
		  // E' defined using the original gamma pT):		  
		  float eprime  = origGamma_calPt * cosh(kinem::eta(matchedZBjet.px, matchedZBjet.py, matchedZBjet.pz));

		  float weight = it2 -> second.weight;
		  h_jetpt -> Fill(jetpt,weight); 
		  ep -> Fill(eprime,weight);
		  //Get the correct bin:
		  int true_npvb = (int) it2 -> second.pv_size;
		  int npvb = true_npvb;
		  int npvb_forMPF = true_npvb;
		  
		  if (_do_inclusive_nPV) {
		    // for arrays binning:
		    if (npvb > merge_nvertex) // see Binnings
		      npvb = merge_nvertex;
		  }

		  // for arrays binning:
		  if (npvb_forMPF > merge_nvertex_forMPF) // see Binnings
		    npvb_forMPF = merge_nvertex_forMPF;
		  
		  // this is cut on 10, co it's safe, we never really look there except for inclusive cases, so it's ok
		  if ( npvb >= nvertex)
		    npvb = nvertex - 1;

		  if (npvb > 0) // this is important
		    //		  if (npvb > 0 && npvb < nvertex)
		    {
		      //jetdeta -> Fill(matchedZBjet.det_eta/10);
		      
		      Float_t feta = matchedZBjet.det_eta/10.;
		      int detabin = Binnings::GetDetaBin(feta);

		      if (debug > 3)
			cout << " feta=" << feta
			    << " detabin=" << detabin
			    << " matchedZBjet.det_eta=" << matchedZBjet.det_eta
			    << endl;

		      //		  float pTprime = eprime / cosh(0.1*(matchedZBjet.det_eta) );
		      // for the new pTprime we need to know in which det eta bin we are
		      // and we divide E' by the center of bin
		      float pTprime = eprime / cosh( eta_center[detabin] ); 
		      
		      if (debug > 3)
			cout << "Eprime: " << eprime << " pTgamma: " << pTgamma 
			     << " pTprime: " << pTprime
			     << " eta_center[detabin=" << detabin << "] " << eta_center[detabin] 
			     << " eta jet: " << kinem::eta(matchedZBjet.px, matchedZBjet.py, matchedZBjet.pz) << endl;
		      


		      if (detabin >= 0 && detabin < ndeta)
			{

			  int eprimebin  = Binnings::GetEprimeBin(eprime);
			  int pTgammabin = Binnings::GetPtGammaBin(pTgamma);
			  int pTprimebin = Binnings::GetPtPrimeBin(pTprime);
			  int ejetbin    = Binnings::GetEjetBin(ejet);
			  int ebin       = Binnings::GetEBin(e2);

			  float corrf = 0.;
			  float corrf2 = 0.;
 
			  if (offset > 0.)
			    {

			      total++;

			      // e1 = Enozb
			      // e2 = Ezb

			      // (Ezb - Enozb) / O 
			      corrf  = ((e2 - e1) / offset);

			      // (O - Ezb + Enozb) / ( Ezb - O) = Enozb / ( Ezb - O) - 1

			      if (! (_CorrectOffset || _Correct_kOkR ) ) {
				corrf2 = (e1) / (e2 - offset) - 1.;
				kO_numerator = e1;
				kO_denominator = e2 - offset;
				if (debug > 1) {
				  cout << " kO_numerator=" << kO_numerator
				       << " kO_denominator=" << kO_denominator
				       << endl;
				}
			      } else
				if (_CorrectOffset || _Correct_kOkR ) {
				  // jiri: OK, we're using the k = 1 + Delta correctin
				  // this applies to E - O_jetcorr
				  // so we have assigned (see above) to e2 = k*(E - O_jetcorr)
				  // as the correction to O_jetcorr is tricky 
				  // (danger of negative energies, fits sensitive to nPV's)
				  corrf2 = (e1 - e2) / e2;
				}
				 
			    } // offset > 0
 
			  float sca    = (float) nit / kNAverTowers;
			  float corrf4 = corrf/sca; //average #towers set to 45...
			  float corrf_RMPF = 1.;
			  if (response2 != 0.)
			    corrf_RMPF = response1/response2; // 1...noZB
			  else
			    corrf_RMPF = kError_RMPF; 

			  // --  Filling for MPF (kR) -- :


			  // fill MPF corrf only for 1 or 2 true nPVs:
			  if (corrf_RMPF > kError_RMPF * kOnePlusEpsilon && true_npvb < nvertex_forMPF_cut ) {

			    // cout << eprimebin << endl;
			    if ((eprimebin >= 0) && (eprimebin < neprime))
			      {
				corr_RMPF[npvb_forMPF-1][eprimebin][detabin]          ->  Fill(corrf_RMPF, weight);
				response_h2_eprime_noZB[npvb_forMPF-1][detabin] -> Fill(eprime, response1, weight);
				response_h2_eprime_ZB[npvb_forMPF-1][detabin]   -> Fill(eprime, response2, weight);
			      }

			    if ((ejetbin >= 0) && (ejetbin < nejet))
			      {
				corr_RMPF_ejet[npvb_forMPF-1][ejetbin][detabin]       ->  Fill(corrf_RMPF, weight);
			      }

			    if ((pTgammabin >= 0) && (pTgammabin < npTgamma))
			      {
				corr_RMPF_pTgamma[npvb_forMPF-1][pTgammabin][detabin] ->  Fill(corrf_RMPF, weight);
			      }
			      
			    if ((pTprimebin >= 0) && (pTprimebin < npTgamma))
			      {
				corr_RMPF_pTprime[npvb_forMPF-1][pTprimebin][detabin] ->  Fill(corrf_RMPF, weight);
				response_h2_pTprime_noZB[npvb_forMPF-1][detabin] -> Fill(pTprime, response1, weight);
				response_h2_pTprime_ZB[npvb_forMPF-1][detabin]   -> Fill(pTprime, response2, weight);
			      }

			  }  // got good response corrf

			  // -- Filling for kO -- :

			  // cout << corrf << " " << corrf4 << " " << nit << endl;
			  // cout << "Correction factor: " << corrf << " difference: " << e2-e1 << " offset: " << offset << endl;
			  //			   cout << "Bin: " << npvb << " " << eprimebin << " " << detabin << endl;

			  if (offset > 0.)
			    {

			      const int maxnFillOptions = 2;
			      int nFillOptions = 2;
			      int pvbin_fill[maxnFillOptions] ={npvb - 1, 0};
			      bool do_pvbin_fill[maxnFillOptions] ={ true, false};
			      // When averaging over PV's,
			      // the idea is to fill vertex bins as follows:
			      // 0 ... 1PV
			      // 1 ... all PV's (average)
			      // 2 ... >= PV's

			      if (_doPVAverage) {
				pvbin_fill[0] = 1;
				if (true_npvb == 1)
				  pvbin_fill[1] = 0;
				else
				  pvbin_fill[1] = 2;
				do_pvbin_fill[1] = true;
			      } else {
				nFillOptions = 1;
			      }

			      for (int ipvbin = 0; ipvbin < nFillOptions; ++ipvbin) {
				//				cout << "ipvbin=" << ipvbin << ", filling bin " << pvbin_fill[ipvbin] << endl;
				if (do_pvbin_fill[ipvbin]) {

				if ((eprimebin >= 0) && (eprimebin < neprime))
				  {
				    kO_numerator_h2_eprime[pvbin_fill[ipvbin]][detabin]   -> Fill(eprime, kO_numerator, weight);
				    kO_denominator_h2_eprime[pvbin_fill[ipvbin]][detabin] -> Fill(eprime, kO_denominator, weight);
				  
				    corr[pvbin_fill[ipvbin]][eprimebin][detabin]      ->  Fill(corrf, weight);
				    corr2[pvbin_fill[ipvbin]][eprimebin][detabin]     ->  Fill(corrf2, weight);
				    weights[pvbin_fill[ipvbin]][eprimebin][detabin]   ->  Fill(1.,(double) weight);
				    weights2[pvbin_fill[ipvbin]][eprimebin][detabin]  ->  Fill(1., (double) (weight*weight));
				    corr4[pvbin_fill[ipvbin]][eprimebin][detabin]     ->  Fill(corrf4, weight);
				    ntow[pvbin_fill[ipvbin]][eprimebin][detabin]      ->  Fill(nit, weight);
			
				    int fill = npvb*10000 + detabin*100 + eprimebin;
				    filled -> Fill(fill);
				    //jetdeta -> Fill(matchedZBjet.det_eta/10);
				    jetdeta -> Fill(matchedZBjet.det_eta/10);
				    corrfac -> Fill(corrf);

				    ehist[pvbin_fill[ipvbin]][eprimebin][detabin]          -> Fill(eprime, weight);
				  } // eprimebin >= 0
			      
				if ((ejetbin >= 0) && (ejetbin < nejet))
				  {
				    corr_ejet[pvbin_fill[ipvbin]][ejetbin][detabin]      ->  Fill(corrf, weight);
				    corr2_ejet[pvbin_fill[ipvbin]][ejetbin][detabin]     ->  Fill(corrf2, weight);
				    weights_ejet[pvbin_fill[ipvbin]][ejetbin][detabin]   ->  Fill(1.,(double) weight);
				    weights2_ejet[pvbin_fill[ipvbin]][ejetbin][detabin]  ->  Fill(1., (double) (weight*weight));
				    ehist_ejet[pvbin_fill[ipvbin]][ejetbin][detabin]       -> Fill(ejet, weight);
				  } // ejetbin >= 0
			      
			      
				if ((pTgammabin >= 0) && (pTgammabin < npTgamma))
				  {
				    corr_pTgamma[pvbin_fill[ipvbin]][pTgammabin][detabin]      ->  Fill(corrf, weight);
				    corr2_pTgamma[pvbin_fill[ipvbin]][pTgammabin][detabin]     ->  Fill(corrf2, weight);
				    weights_pTgamma[pvbin_fill[ipvbin]][pTgammabin][detabin]   ->  Fill(1.,(double) weight);
				    weights2_pTgamma[pvbin_fill[ipvbin]][pTgammabin][detabin]  ->  Fill(1., (double) (weight*weight));
				    ehist_pTgamma[pvbin_fill[ipvbin]][pTgammabin][detabin] -> Fill(pTgamma, weight);
				  } // pTgammabin >= 0

				if ((pTprimebin >= 0) && (pTprimebin < npTgamma))
				  {
				    kO_numerator_h2_pTprime[pvbin_fill[ipvbin]][detabin]   -> Fill(pTprime, kO_numerator, weight);
				    kO_denominator_h2_pTprime[pvbin_fill[ipvbin]][detabin] -> Fill(pTprime, kO_denominator, weight);

				    corr_pTprime[pvbin_fill[ipvbin]][pTprimebin][detabin]      ->  Fill(corrf, weight);
				    corr2_pTprime[pvbin_fill[ipvbin]][pTprimebin][detabin]     ->  Fill(corrf2, weight);
				    weights_pTprime[pvbin_fill[ipvbin]][pTprimebin][detabin]   ->  Fill(1.,(double) weight);
				    weights2_pTprime[pvbin_fill[ipvbin]][pTprimebin][detabin]  ->  Fill(1., (double) (weight*weight));
				    ehist_pTprime[pvbin_fill[ipvbin]][pTprimebin][detabin] -> Fill(pTprime, weight);
				  } // pTprimebin >= 0
			 
				if ((ebin >=0) && (ebin <neprime))
				  {
				    // (O - Ezb + Enozb) / ( Ezb - O) = Enozb / ( Ezb - O) - 1
				    // is this guy the same as corrf2 ?!
				    float corrf3 = (e1) / (e2 - offset) - 1.;
				    corr3[pvbin_fill[ipvbin]][ebin][detabin] -> Fill(corrf3, weight);
				  
				  } // ebin >= 0
			      
				} // if fill

			      } // ipvbin

			    } // offset > 0

			 

			} // detabin
		    } // pvbin

		} // deltaphi cut
	    } // deltaR cut
	  
	} // found match
    } // go through the ZB sample

  cout << "Opening output file..." << endl;
  TFile *output = new TFile(outputfilename.Data(), "RECREATE");

  cout << "Writting histograms..." << endl;

  InstLumi_h1 -> Write();
  h_dr_drCut_dPhiCut -> Write();
  h_dr_drCut -> Write();
  h_dr     -> Write();
  h_dr_all -> Write();
  h_besti -> Write();

  h_dphi_cut -> Write();
  h_dphi -> Write();

  npv     -> Write();
  h_jetpt   -> Write();
  filled  -> Write();
  jetdeta -> Write();
  ep      -> Write();
  corrfac -> Write();

  for (int i = 0; i < nvertex; i++)
    {

      cout << "  ... for vertex " << i + 1 << endl;
      TString pwd;
      pwd += i+1;
      gDirectory -> mkdir(pwd);
      gDirectory -> cd(pwd);
      cout << "       ... for etabin: ";
      for (int k = 0; k < ndeta; k++)
	{
	  cout << "" << k << " " << flush;
	  TString pwd2;
	  pwd2 += k;
	  gDirectory -> mkdir(pwd2);
	  gDirectory -> cd(pwd2);

	  response_h2_eprime_noZB[i][k]   -> Write();
	  response_h2_eprime_ZB[i][k]     -> Write();
	  response_h2_pTprime_noZB[i][k]  -> Write();
	  response_h2_pTprime_ZB[i][k]    -> Write();

	  kO_numerator_h2_pTprime[i][k]   -> Write();
	  kO_denominator_h2_pTprime[i][k] -> Write();
	  kO_numerator_h2_eprime[i][k]    -> Write();
	  kO_denominator_h2_eprime[i][k]  -> Write();


	  for (int j = 0; j < neprime; j++)
	    {
	      corr[i][j][k]      -> Write();
	      corr2[i][j][k]     -> Write();
	      corr_RMPF[i][j][k] -> Write();
	      weights[i][j][k]   -> Write();
	      weights2[i][j][k]  -> Write();
	      ehist[i][j][k]     -> Write();
	      
	      corr3[i][j][k]    -> Write();
	      corr4[i][j][k]    -> Write();
	      ntow[i][j][k]     -> Write();
	      
	      //weights3[i][j][k] -> Write();
	    }
	  
	  for (int j = 0; j < npTgamma; j++) 
	    {
	      corr_pTgamma[i][j][k]      -> Write();
	      corr2_pTgamma[i][j][k]     -> Write();
	      corr_RMPF_pTgamma[i][j][k] -> Write();
	      weights_pTgamma[i][j][k]   -> Write();
	      weights2_pTgamma[i][j][k]  -> Write();
	      ehist_pTgamma[i][j][k]     -> Write();

	      corr_pTprime[i][j][k]      -> Write();
	      corr2_pTprime[i][j][k]     -> Write();
	      corr_RMPF_pTprime[i][j][k] -> Write();
	      weights_pTprime[i][j][k]   -> Write();
	      weights2_pTprime[i][j][k]  -> Write();
	      ehist_pTprime[i][j][k]     -> Write();

	    }
	  
	  for (int j = 0; j < nejet; j++) 
	    {
	      corr_ejet[i][j][k]      -> Write();
	      corr2_ejet[i][j][k]     -> Write();
	      corr_RMPF_ejet[i][j][k] -> Write();
	      weights_ejet[i][j][k]   -> Write();
	      weights2_ejet[i][j][k]  -> Write();
	      ehist_ejet[i][j][k]     -> Write();
	    }
	  gDirectory -> cd("../");
	}
      cout << endl;
      gDirectory -> cd("../");
    }
  

  cout << "Closing the file..." << endl;
  output -> Close();

  cout << "+-----------------------------------------------------------------+" << endl;
  cout << "|  Found same events: "               << found                       << endl;
  cout << "|  Matched, dphi with gamma > " <<  _DeltaPhiCut << ": "  << good    << endl;
  cout << "|  Used for filling: "                << total                       << endl;
  cout << "+-----------------------------------------------------------------+" << endl;
  cout << "Done! (kill me if I hang;-)" << endl;
 return 0;

}