///**
// *  Underlying Event in VB-events analysis
// *
// *  @author Michael Leyton <leyton@cern.ch>
// *  @author Holger Schulz <hschulz@physik.hu-berlin.de>
// * 
// */

#include "UEShapesNoVtx.h"


#include <math.h>
#include <list>
#include <vector>
#include <fstream>
#include <cmath>
#include <dirent.h>
#include <TSystem.h>
#include <TRandom.h>
#include <TString.h>
#include <TMath.h>
#include <TF1.h>
#include <sstream>
#include <algorithm>
#include <iterator>

// histograms generated from par file
#include "UEShapesNoVtx_hist.h"
//#include "TrackD3PDObject_hist.h"

// Some tools
#include "TrackTools.h"
#include "Helpers.h"
#include "NNtree.h"

#include <sys/stat.h>
#include <boost/tokenizer.hpp>
#include <boost/filesystem/path.hpp>
#include <boost/filesystem/operations.hpp>
#include <boost/algorithm/string.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>
#include <boost/archive/binary_oarchive.hpp>
#include <boost/archive/binary_iarchive.hpp>
#include <boost/serialization/base_object.hpp>
#include <boost/serialization/map.hpp>
#include <boost/serialization/vector.hpp>
#pragma link C++ class map< int, vector<vector<vector<double> > > >;

using namespace UEShapesNoVtx_hist;
using boost::lexical_cast;
using namespace boost::filesystem;
//using namespace TrackD3PDObject_hist;

using namespace std;
typedef unsigned int uint;

// 0 = pythia, 1 = sherpa, 2 = alpgen
static int MCGENERATOR = 0;

static double CUT_DRESSINGCONE  = 0.1; // Variations --- 0.1 central, 0.05 down 0.2 up
static double CUT_LEPCONE  = 0.2; // Variations --- 0.1 central, 0.05 down 0.2 up
//static double CUT_LEPCONE  = 0.2;
//static double CUT_LEPCONE  = 0.05;
//static double CUT_LEPPTRES = 0.01; 
//static double CUT_DELTAR   = 0.005; // HS was
//static double CUT_DELTAR   = 0.010; 

const double CUT_PTZ = 10.0;

int ntrk_sig = 0;
int ntrk_add = 0;
int pu_tracks = 0;


#define _USE_RECOIL_
#define _Z_ANA_

inline double fastsin(double x) {
  return x - pow(x,3)/6. + pow(x,5)/120. - pow(x,7)/5040.;
}

double UEShapesNoVtx::Dphi(TVector3 p1, TVector3 p2) { // olegs code
   double dphi, cos_p1_p2, norm_p1, norm_p2, sign;
   //if (p1 && p2) {
       norm_p1 = sqrt(p1.x()*p1.x()+p1.y()*p1.y());
       norm_p2 = sqrt(p2.x()*p2.x()+p2.y()*p2.y());
       if (norm_p1 == 0. || norm_p2 == 0.)
               return 0.;
       sign = ( p1.x()*p2.y() - p1.y()*p2.x() > 0. ? -1. : 1. );
       return sign*acos ( (p1.x()*p2.x()+p1.y()*p2.y()) / norm_p1 / norm_p2 );
   //}
}

//bool UEShapesNoVtx::SelectZ(TVector3 lep1,TVector3 lep2) {
    //bool select=false;
    //double Mee = calculateMeeZ(lep1, lep2, m_muon);
    //if (m_QCD && m_muon) {
        //if ((Mee > 40. && Mee < 66.) || (Mee > 116. && Mee <140.)) select=true;
    //}
    //else if (m_QCD && !m_muon) {
        //if ((Mee > 40. && Mee < 66.) || (Mee > 116. && Mee <140.)) select=true;
    //}
    //else {
        //if ((Mee > 66. && Mee <116.)) select=true;
    //}

    ////cout << "Mee: " << Mee << " select: " << select << endl;

    //return select;

//}
//bool UEShapesNoVtx::SelectZ(double Mee) {
    //bool select=false;
    //if (m_QCD && m_muon) {
        //if ((Mee > 40. && Mee < 66.) || (Mee > 116. && Mee <160.)) select=true;
    //}
    //else if (m_QCD && !m_muon) {
        //if ((Mee > 40. && Mee < 66.) || (Mee > 116. && Mee <160.)) select=true;
    //}
    ////if (m_QCD) {
        ////if ((Mee > 50. && Mee <66.) || (Mee > 116. && Mee <140.)) select=true;
    ////}
    //else {
        //if ((Mee > 66. && Mee <116.)) select=true;
    //}

    ////cout << "Mee: " << Mee << " select: " << select << endl;

    //return select;

//}

TVector3 UEShapesNoVtx::mkLepton(double pt, double eta, double phi, bool ismuon) {
    // if the muon flag is fals, assume we are given the Energy rather than the pt --- needed for electron systematics
    TVector3 lep;
    if (ismuon) lep.SetPtEtaPhi(pt, eta, phi);
    else lep.SetPtEtaPhi(pt/cosh(eta), eta, phi);
    return lep;
}

bool  UEShapesNoVtx::recoCutVB(double VB_m, TVector3 l1, TVector3 l2) {
    if (m_QCD) {
        //if (!( (VB_m > 50. && VB_m < 66.) || (VB_m > 116. && VB_m <160.) )) return false;
        if (!( VB_m > 50.  && VB_m <160. )) return false;
    }
    else {
        if (!(VB_m > m_MLL_MIN && VB_m < m_MLL_MAX)) return false; // CUt on Zmass window
        //if (!(VB_m > 70. && VB_m < 112.)) {
            //cout << "Failed mll cut: " << VB_m << endl;
            //return false; // CUt on Zmass window
        //}
    }
    if (       l1.Pt() <= 20. ||        l2.Pt() <= 20.) return false;
    if (fabs(l1.Eta()) >= 2.4 || fabs(l2.Eta()) >= 2.4) return false;
    if (! m_muon) {
        if ((fabs(l1.Eta()) >1.37 && fabs(l1.Eta()) <1.52) || (fabs(l2.Eta()) >1.37 && fabs(l2.Eta()) <1.52)) return false;
    }
    return true;
}
bool  UEShapesNoVtx::truthCutVB(double VB_m, TVector3 l1, TVector3 l2) { // unused function
    if (!(VB_m > m_MLL_MIN && VB_m < m_MLL_MAX)) return false; // CUt on Zmass window
    if (       l1.Pt() <= 20. ||        l2.Pt() <= 20.) return false;
    if (fabs(l1.Eta()) >= 2.4 || fabs(l2.Eta()) >= 2.4) return false;
    //if (! m_muon) {
        //if ((fabs(l1.Eta()) >1.37 && fabs(l1.Eta()) <1.52) || (fabs(l2.Eta()) >1.37 && fabs(l2.Eta()) <1.52)) return false;
    //}
    //if not m
    //if ((fabs(l1.Eta()) >1.37 && fabs(l1.Eta()) <1.52) || (fabs(l2.Eta()) >1.37 && fabs(l2.Eta()) <1.52)) return false;
    return true;
}

//bool  UEShapesNoVtx::recoCutVB(double VB_m, TVector3 l1, TVector3 l2) {
    //if (!(VB_m > 66. && VB_m < 116.)) return false; // CUt on Zmass window
    //if (       l1.Pt() <= 20. ||        l2.Pt() <= 20.) return false;
    //if (fabs(l1.Eta()) >= 2.4 || fabs(l2.Eta()) >= 2.4) return false;
    //return true;
//}

double UEShapesNoVtx::getRdmDz() {
  double rdmZ = 0.0;
  if (!flag_isMC)                      rdmZ = gRandom->Gaus(-5.263, 57.62);
  else {
      //if      (MCGENERATOR=="pythia8") rdmZ = gRandom->Gaus(-6.244, 57.51);
      if      (MCGENERATOR=="pythia8") rdmZ = gRandom->Gaus(-1.15e-4, 93.02);
      else if (MCGENERATOR=="sherpa")  rdmZ = gRandom->Gaus(-6.304, 57.45);
  }
  // Random deltaZ
  //return fabs(rdmZ - vz);
  return rdmZ;
  //return rdmZ - vz;
}
//______________________________________________________________________

// ctor
UEShapesNoVtx::UEShapesNoVtx(TTree * tree) : 
  nevt(0) { 
  
  cout << "<UEShapesNoVtx::UEShapesNoVtx>" << endl;
  MCGENERATOR="";

  m_tree    = tree;
//  m_tree->SetBranchAddress("mc_vx_z", &mc_vx_z);
  cout << "<UEShapesNoVtx::UEShapesNoVtx> mTree" << endl;
  //m_entries = m_tree->GetEntries();
  //cout << "<UEShapesNoVtx::UEShapesNoVtx> m_entries" << endl;

  //cout << "Number of entries in tree: " << m_entries << endl;

  //Create the D3PDReader objects:
  master = 0; // this counts the events
  mcevt  = new D3PDReader::McEventD3PDObject( master );
  mctruth = new D3PDReader::McTruthD3PDObject ( master );
  mcgen  = new D3PDReader::McGenD3PDObject ( master ); // HS war auskommentiert
  trk    = new D3PDReader::TrackD3PDObject ( master );
  vtx    = new D3PDReader::VertexD3PDObject ( master );
  VB     = new D3PDReader::VectorBosonD3PDObject( master );
  BS     = new D3PDReader::BeamSpotD3PDObject( master );

  // Connect the D3PDReader objects to the TTree:
  mcevt->ReadFrom( m_tree );
  mctruth->ReadFrom( m_tree );
  mcgen->ReadFrom( m_tree ); // HS war auskommentiert
  trk->ReadFrom( m_tree );
  vtx->ReadFrom( m_tree );
  VB->ReadFrom( m_tree );
  BS->ReadFrom( m_tree );

  mcGen = new McGenD3PDObjectAdapter( mcgen );  // HS war auskommentiert
  //mc = new McTruthD3PDObjectAdapter( mcgen ); 
  mc = new McTruthD3PDObjectAdapter( mctruth ); // war von leyton
  
}

void UEShapesNoVtx::LoadHBOMCorrections() {
  l.info("<UEShapes::LoadHBOMCorrections>");
  if (m_HBOMFile.Sizeof()>1) { // Check if HBOM file  name successfully given
    TFile* fHBOM = TFile::Open(m_HBOMFile);
    
    m_HBOM_AvgPt     =(TH1D*)fHBOM->Get("hAvgPtRecSel_correction_w_errs")->Clone("AvgPtcorrectionPU");
    m_HBOM_BeamThrust=(TH1D*)fHBOM->Get("hBeamThrustRecSel_correction_w_errs")->Clone("BeamThrustcorrectionPU");
    m_HBOM_FParameter=(TH1D*)fHBOM->Get("hFParameterRecSel_correction_w_errs")->Clone("FParametercorrectionPU");
    m_HBOM_Minor     =(TH1D*)fHBOM->Get("hMinorRecSel_correction_w_errs")->Clone("MinorcorrectionPU");
    m_HBOM_Ntrk      =(TH1D*)fHBOM->Get("hNtrkRecSel_correction_w_errs")->Clone("NtrkcorrectionPU");
    m_HBOM_Spherocity=(TH1D*)fHBOM->Get("hSpherocityRecSel_correction_w_errs")->Clone("SpherocitycorrectionPU");
    m_HBOM_SumPt     =(TH1D*)fHBOM->Get("hSumPtRecSel_correction_w_errs")->Clone("SumPtcorrectionPU");
    m_HBOM_Thrust    =(TH1D*)fHBOM->Get("hThrustRecSel_correction_w_errs")->Clone("ThrustcorrectionPU");
   
    m_HBOM_AvgPt     ->SetDirectory(0);
    m_HBOM_BeamThrust->SetDirectory(0);
    m_HBOM_FParameter->SetDirectory(0);
    m_HBOM_Minor     ->SetDirectory(0);
    m_HBOM_Ntrk      ->SetDirectory(0);
    m_HBOM_Spherocity->SetDirectory(0);
    m_HBOM_SumPt     ->SetDirectory(0);
    m_HBOM_Thrust    ->SetDirectory(0);

    fHBOM->Close();
  }
}

void UEShapesNoVtx::LoadMeanHBOMCorrections() {
    TFile* fHBOMMean = TFile::Open(m_HBOMMeanFile);
    m_meanCorrNtrk       = (TProfile*)fHBOMMean->Get("hMeanCorrNtrk");
    m_meanCorrSumPt      = (TProfile*)fHBOMMean->Get("hMeanCorrSumPt");
    m_meanCorrAvgPt      = (TProfile*)fHBOMMean->Get("hMeanCorrAvgPt");
    m_meanCorrBeamThrust = (TProfile*)fHBOMMean->Get("hMeanCorrBeamThrust");
    m_meanCorrThrust     = (TProfile*)fHBOMMean->Get("hMeanCorrThrust");
    m_meanCorrMinor      = (TProfile*)fHBOMMean->Get("hMeanCorrMinor");
    m_meanCorrFParameter = (TProfile*)fHBOMMean->Get("hMeanCorrFParameter");
    m_meanCorrSpherocity = (TProfile*)fHBOMMean->Get("hMeanCorrSpherocity");
   
    m_meanCorrNtrk      ->SetDirectory(0);
    m_meanCorrSumPt     ->SetDirectory(0);
    m_meanCorrAvgPt     ->SetDirectory(0);
    m_meanCorrBeamThrust->SetDirectory(0);
    m_meanCorrThrust    ->SetDirectory(0);
    m_meanCorrMinor     ->SetDirectory(0);
    m_meanCorrFParameter->SetDirectory(0);
    m_meanCorrSpherocity->SetDirectory(0);

    fHBOMMean->Close();
}

void UEShapesNoVtx::bookHBOMHistos() {

    double ntrk_bins[21] = {-0.001,2.500,4.500,6.500,8.500,10.500,12.500,15.500,18.500,22.500,26.500,30.500,35.500,40.500,45.500,50.500,60.500,75.500,100.500,125.5,150.5};
    double sumpt_bins[24]= {-0.001,1.,2.,3.,4.,6.,8.,10.,12.,14.,16.,18.,20.,22.,25.,30.,35.,40.,45.,50.,60.,75.,100.,150.};
    double avgpt_bins[19]= {-0.001,0.100,0.200,0.300,0.400,0.500,0.600,0.700,0.800,0.900,1.,1.200,1.400,1.600,1.800,2.500,3.,4.,6.};
    double bt_bins[26]   = {-0.001,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,12.,14.,16.,18.,20.,23.,26.,29.,32.,35.,40.,50.,60.,75.,100.};
    m_meanCorrNtrk       = new TProfile("hMeanCorrNtrk",      "hMeanCorrNtrk",       20, ntrk_bins);
    m_meanCorrSumPt      = new TProfile("hMeanCorrSumPt",     "hMeanCorrSumPt",      23, sumpt_bins); 
    m_meanCorrAvgPt      = new TProfile("hMeanCorrAvgPt",     "hMeanCorrAvgPt",      18, avgpt_bins);
    m_meanCorrBeamThrust = new TProfile("hMeanCorrBeamThrust","hMeanCorrBeamThrust", 25, bt_bins);
    m_meanCorrThrust     = new TProfile("hMeanCorrThrust",    "hMeanCorrThrust",     12, 0.636, 1.001);
    m_meanCorrMinor      = new TProfile("hMeanCorrMinor",     "hMeanCorrMinor",      12, 0.0,   0.7);    
    m_meanCorrFParameter = new TProfile("hMeanCorrFParameter","hMeanCorrFParameter", 12, 0.0,   1.0);
    m_meanCorrSpherocity = new TProfile("hMeanCorrSpherocity","hMeanCorrSpherocity", 12, 0.0,   1.0);

    m_meanCorrNtrk      ->SetDirectory(m_outputFile);
    m_meanCorrSumPt     ->SetDirectory(m_outputFile);
    m_meanCorrAvgPt     ->SetDirectory(m_outputFile);
    m_meanCorrBeamThrust->SetDirectory(m_outputFile);
    m_meanCorrThrust    ->SetDirectory(m_outputFile);
    m_meanCorrMinor     ->SetDirectory(m_outputFile);
    m_meanCorrFParameter->SetDirectory(m_outputFile);
    m_meanCorrSpherocity->SetDirectory(m_outputFile);
}

double UEShapesNoVtx::GetHBOMCorr(TString esname, double esvalue) {
    TH1D* correctionPU;
    if (esname == "AvgPt") correctionPU = m_HBOM_AvgPt;
    if (esname == "BeamThrust") correctionPU = m_HBOM_BeamThrust;
    if (esname == "FParameter") correctionPU = m_HBOM_FParameter;
    if (esname == "Minor") correctionPU = m_HBOM_Minor;
    if (esname == "Ntrk") correctionPU = m_HBOM_Ntrk;
    if (esname == "Spherocity") correctionPU = m_HBOM_Spherocity;
    if (esname == "SumPt") correctionPU = m_HBOM_SumPt;
    if (esname == "Thrust") correctionPU = m_HBOM_Thrust;

    int binid = correctionPU->FindBin(esvalue);
    if (binid > correctionPU->GetNbinsX()) binid = correctionPU->GetNbinsX();
    if (binid <1) binid = 1;

    return correctionPU->GetBinContent(binid);
}

double UEShapesNoVtx::GetHBOMCorrErr(TString esname, double esvalue) {
    TH1D* correctionPU;
    if (esname == "AvgPt") correctionPU = m_HBOM_AvgPt;
    if (esname == "BeamThrust") correctionPU = m_HBOM_BeamThrust;
    if (esname == "FParameter") correctionPU = m_HBOM_FParameter;
    if (esname == "Minor") correctionPU = m_HBOM_Minor;
    if (esname == "Ntrk") correctionPU = m_HBOM_Ntrk;
    if (esname == "Spherocity") correctionPU = m_HBOM_Spherocity;
    if (esname == "SumPt") correctionPU = m_HBOM_SumPt;
    if (esname == "Thrust") correctionPU = m_HBOM_Thrust;

    int binid = correctionPU->FindBin(esvalue);
    if (binid > correctionPU->GetNbinsX()) binid = correctionPU->GetNbinsX();
    if (binid <1) binid = 1;

    return correctionPU->GetBinError(binid);
}

void UEShapesNoVtx::Begin(TTree * /*tree*/) {

  l = Logger();
  //gRandom->SetSeed(0); // internally create a random seed
  gRandom->SetSeed(m_seed); // Set a constant seed
  cout << "Random seed: " << m_seed << endl;
  cout << "HBOM file: " << m_HBOMFile << endl;
  
  // Set up either filling of the mean HBOM correction profiles or reading them in from separate file
  if (m_HBOMMeanFile.Sizeof()>1) {
      LoadMeanHBOMCorrections();
      cout << "Read mean HBOM corrections from " << m_HBOMMeanFile.Sizeof() << endl;
  }
  else {
      cout << "Booking HBOM mean corr histos" << endl;
      if (!m_dumpReco) bookHBOMHistos();
  }
  //exit(1);

  l.setName("UEShapesNoVtx.cxx");

  if (flag_DEBUG) l.debug("<UEShapesNoVtx::Begin>");

  if (m_HBOMFile.Sizeof()>1) LoadHBOMCorrections(); // Check if HBOM file  name successfully given

  m_timer.Reset();
  m_timer.Start();
  m_totaltime = 0;

  dump_vec.clear();
  dump_eventnumber.clear();

  // Initialize global variables
  gRand = new TRandom3();
  gRand->SetSeed(0);
  SetRegions();
  FillValues();
  nevt=0;
  nevt_weighted=0;
  cum_merging_sum = 0;

  m_sidebands_loaded=false;

  // dumpReco
  if (m_dumpReco) {
      cout << "Hellalujah baby" << endl;
      LoadHBOMCorrections();
  }
  //exit(1);

  initTransientVariables();
 
  //readEventAndRunNumbers("/afs/naf.desy.de/user/h/hschulz/UEinVB/run/mcinfo.Merged_Pythia8_nominal_12.txt");
  readEventNumbers("/afs/naf.desy.de/user/h/hschulz/UEinVB/run/mcinfo.intersect.txt");
  //readEventNumbers("/afs/naf.desy.de/user/h/hschulz/UEinVB/run/intersection.txt");

  // Load corrections 
  LoadTrkCorrection("7TeV_ndiff_correctionPlots_version1_noprobcut_rebinned_11June.root"); 

  LoadPUlib();
  bool load_tracklib=true;
  if (!m_writeTrackLib && m_maxAdditionalPileUp>0) LoadTrackLib();

  if (m_writeTrackLib) {
    //m_pu_z = new TH1F("z_pu", "z position of pu vertices", 50, -250, 250);
    m_pu_vtxN = new TH1F("vtxN_pu", "Number of type-3 vertices per event", 30, 0, 30);
    m_pu_vtxDz = new TH2F("vtxDz_pu", "zPrimary - zPU vs zPri", 200, -250, 250, 100, -250, 250);
    m_pu_vtxzPP = new TH2F("vtxzPP", "zPU vs zPri", 700, -350, 350, 700, -350, 350);
    const unsigned int nbins_z = 27;
    double array_z[nbins_z+1] = {-250., -150.,-125., -100.,-80.,-70.,-60.,-50.,-40.,-30.,-20.,-10.,-6.,-2.,2.,6.,10.,20.,30.,40.,50.,60.,70.,80.,100.,125.,150.,250.};
    m_pu_zN = new TH2F("z_puN", "z position of pu vertices vs N tracks", nbins_z, array_z, 120, -0.5, 119.5);
    //m_pu_zN = new TH2F("z_puN", "z position of pu vertices vs N tracks", 15, -250, 250, 120, -0.5, 119.5);
    m_pu_z = new TH1F("z_pu", "z position of pu vertices",nbins_z, array_z);
    m_pu_Dz_z = new TH2F("pu_Dz_z", "DeltaZsintheta vs. z",  30, -1.5, 1.5, 25, -250., 250.);
    m_pu_d_z = new TH2F("pu_d_z", "d vs. z",  30, -1.5, 1.5, 25, -250., 250.);
    m_pu_Dz_N = new TH2F("pu_Dz_N", "DeltaZsintheta vs. Ntracks",  41, -1.5, 1.5, 120, -0.5, 119.5);
  }

  m_outputFile->cd(); 
  
  // booking of autogenerated histograms
  for (int sys=0; sys < MAXSYSTYPES; sys++) {
    for (int h=0; h < 400000; h++) {
      // cout << "sys = " << sys << " and h = " << h << endl;
      harray[sys][h] = 0;
    }
  }
  bookHistograms(m_outputFile, harray);
  //Rebin(harray);

  // imagiro trees and branches
  imagiro_tree_rec = new TTree("imagiroRec", "imagiroRec");
  addImagiroBranches(imagiro_tree_rec);

  imagiro_tree_rec2 = new TTree("imagiroRec2", "imagiroRec2");
  addImagiroBranches2(imagiro_tree_rec2);
  
  imagiro_tree_recTrkEff = new TTree("imagiroRecTrkEff", "imagiroRecTrkEff");
  addImagiroBranches(imagiro_tree_recTrkEff);

  imagiro_tree_rec2TrkEff = new TTree("imagiroRec2TrkEff", "imagiroRec2TrkEff");
  addImagiroBranches2(imagiro_tree_rec2TrkEff);
  
  //imagiro_tree_recTrkSmearUp = new TTree("imagiroRecTrkSmearUp", "imagiroRecTrkSmearUp");
  //addImagiroBranches(imagiro_tree_recTrkSmearUp);

  //imagiro_tree_rec2TrkSmearUp = new TTree("imagiroRec2TrkSmearUp", "imagiroRec2TrkSmearUp");
  //addImagiroBranches2(imagiro_tree_rec2TrkSmearUp);
  
  //imagiro_tree_recTrkSmearDown = new TTree("imagiroRecTrkSmearDown", "imagiroRecTrkSmearDown");
  //addImagiroBranches(imagiro_tree_recTrkSmearDown);

  //imagiro_tree_rec2TrkSmearDown = new TTree("imagiroRec2TrkSmearDown", "imagiroRec2TrkSmearDown");
  //addImagiroBranches2(imagiro_tree_rec2TrkSmearDown);
  // PtVB up
  //imagiro_tree_recPtZUp = new TTree("imagiroRecPtZUp", "imagiroRecPtZUp");
  //addImagiroBranches(imagiro_tree_recPtZUp);

  //imagiro_tree_rec2PtZUp = new TTree("imagiroRec2PtZUp", "imagiroRec2PtZUp");
  //addImagiroBranches2(imagiro_tree_rec2PtZUp);
  // PtVB down
  //imagiro_tree_recPtZDown = new TTree("imagiroRecPtZDown", "imagiroRecPtZDown");
  //addImagiroBranches(imagiro_tree_recPtZDown);

  //imagiro_tree_rec2PtZDown = new TTree("imagiroRec2PtZDown", "imagiroRec2PtZDown");
  //addImagiroBranches2(imagiro_tree_rec2PtZDown);

  if (flag_isMC) {
    imagiro_tree_gen = new TTree("imagiroGen", "imagiroGen");
    addImagiroBranches(imagiro_tree_gen);

    imagiro_tree_gen2 = new TTree("imagiroGen2", "imagiroGen2");
    addImagiroBranches2(imagiro_tree_gen2);
    
  }
  
  // NN tree
  if (flag_dumpNN) {
    nntree = new NNtree("RegTree", "RegTree");  
    l.info("Dumping NN tree");
  }

  
}

void UEShapesNoVtx::addImagiroBranches(TTree* imtree) {

  // all events

  imtree->Branch("EventNumber",    &EventNumber,    "EventNumber/i");
  imtree->Branch("EventWeight",    &EventWeight,    "EventWeight/D");
  imtree->Branch("Ntrk",           &Ntrk,           "Ntrk/D");
  imtree->Branch("SumPt",          &SumPt,          "SumPt/D");
  imtree->Branch("SumPtTrans",     &SumPtTrans,     "SumPtTrans/D");
  imtree->Branch("AvgPt",          &AvgPt,          "AvgPt/D");
  imtree->Branch("EtaVB",          &EtaVB,          "EtaVB/D");
  imtree->Branch("PtVB",           &PtVB,           "PtVB/D");
  imtree->Branch("PtVB_up",        &PtVB_up,        "PtVB_up/D");
  imtree->Branch("PtVB_down",      &PtVB_down,      "PtVB_down/D");
}

void UEShapesNoVtx::addImagiroBranches2(TTree* imtree2) {

  // events with nch or nsel >= 2

  imtree2->Branch("EventNumber",    &EventNumber,    "EventNumber/i");
  imtree2->Branch("EventWeight",    &EventWeight,    "EventWeight/D");
  imtree2->Branch("EtaVB",          &EtaVB,          "EtaVB/D");
  imtree2->Branch("PtVB",           &PtVB,           "PtVB/D");
  imtree2->Branch("PtVB_up",        &PtVB_up,        "PtVB_up/D");
  imtree2->Branch("PtVB_down",      &PtVB_down,      "PtVB_down/D");
  imtree2->Branch("Thrust",         &Thrust,         "Thrust/D");
  imtree2->Branch("ThrustPhi",      &ThrustPhi,      "ThrustPhi/D");
  imtree2->Branch("Minor",          &Minor,          "Minor/D");
  imtree2->Branch("FParameter",     &FParameter,     "FParameter/D");
  imtree2->Branch("Spherocity",     &Spherocity,     "Spherocity/D");
  imtree2->Branch("BeamThrust",     &BeamThrust,     "BeamThrust/D");
 
}

void UEShapesNoVtx::initTransientVariables() {

  lep1_trk_3vec.Clear();
  lep2_trk_3vec.Clear();

  lep1_cone_sumpt = 0.;
  lep2_cone_sumpt = 0.;
  lep1_cone_n = 0;
  lep2_cone_n = 0;

  trackweights.clear();
}

void UEShapesNoVtx::readEventNumbers(const char* fname) {
  ifstream inFile;
  string line;
  inFile.open(fname);
  using namespace boost;
  while (inFile.good()) {
    getline (inFile,line);
    vector<string> tokens;
    tokens.clear();
    split(tokens, line, is_any_of("\t"), token_compress_on);
    if (tokens.size() == 1 && tokens[0].length() > 0) {
      mcEventNumbers.push_back(lexical_cast<int>(tokens[0]));
    }
  }
  inFile.close();
  cout << "Read " << mcEventNumbers.size() << " event numbers" << endl;
}

void UEShapesNoVtx::openDumpFile() {
    ofstream m_dumpFile;
    cout << "Opening " << m_dumpRecoFile << " for writing" << endl;
    m_dumpFile.open(m_dumpRecoFile);
    cout << "Opened " << m_dumpRecoFile << " for writing" << endl;
}
void UEShapesNoVtx::closeDumpFile() {
    ofstream m_dumpFile;
    cout << "Closing " << m_dumpRecoFile << endl;
    m_dumpFile.close();
    cout << "Closed " << m_dumpRecoFile << endl;
}


bool UEShapesNoVtx::matchEventNumber(int eventnr) {
  bool matched = std::find(mcEventNumbers.begin(), mcEventNumbers.end(), eventnr)!=mcEventNumbers.end();
  return matched;
}

void UEShapesNoVtx::readEventAndRunNumbers(const char* fname) {
  ifstream inFile;
  string line;
  inFile.open(fname);
  using namespace boost;
  while (inFile.good()) {
    getline (inFile,line);
    vector<string> tokens;
    vector<int> temp;
    tokens.clear();
    temp.clear();
    split(tokens, line, is_any_of("\t"), token_compress_on);
    if (tokens.size() == 2) {
      temp.push_back(lexical_cast<int>(tokens[0]));
      temp.push_back(lexical_cast<int>(tokens[1]));
      mcEventAndRunNumbers.push_back(temp);
    }
  }
  inFile.close();
}

bool UEShapesNoVtx::matchEventAndRunNumber(int eventnr, int runnr) {
  vector<int> test;
  test.push_back(eventnr);
  test.push_back(runnr);

  //vector<vector<int> >::iterator it;

  bool matched = std::find(mcEventAndRunNumbers.begin(), mcEventAndRunNumbers.end(), test)!=mcEventAndRunNumbers.end();
  return matched;
}

void UEShapesNoVtx::getSysType(string sysname) {

  if (sysname=="PileUp")           { m_systype = kPileUp; } 
  else if (sysname=="TrkEff")      { m_systype = kTrkEff; } 
  //else if (sysname=="TrkEffUp")    { m_systype = kTrkEffUp; } 
  //else if (sysname=="TrkEffDown")  { m_systype = kTrkEffDown; } 
  else if (sysname=="TrkSmearUp")    { m_systype = kTrkSmearUp; } 
  else if (sysname=="TrkSmearDown")    { m_systype = kTrkSmearDown; } 
  else { m_systype = 0; }
  
  m_outputFile->cd(sysname.c_str());
 
}

void UEShapesNoVtx::LoadTrkCorrection(TString fname) {

  // TRACKING 

  l.info("<UEShapesNoVtx::LoadTrkCorrection>");

  cout << m_corrpath << "(m_corrpath)" << endl;
  TFile* ftrack1 = TFile::Open(m_corrpath+"/"+fname);

  l.info(TString::Format("<LoadTrkCorrection>: Loading file %s", ftrack1->GetName()));

  mTrkEff_pt_eta  = (TH2F*)(ftrack1->Get("TruthCorrection_EfficiencySelectedTrkStableChargedAll_pt_eta"))->Clone();
  mSecondaries_pt = (TH1F*)(ftrack1->Get("TruthCorrection_SecondariesFractionOfSelected_pt"))->Clone();   
  mOutside_pt_eta = (TH2F*)(ftrack1->Get("TruthCorrection_OutsidePhasespaceFractionOfSelected_pt_eta"))->Clone();

  mTrkEff_pt_eta->SetDirectory(0);
  mSecondaries_pt->SetDirectory(0);
  mOutside_pt_eta->SetDirectory(0);

  delete ftrack1;

}

int UEShapesNoVtx::getPtIndex(double pt) {
    int ipt = -1;
    if     (pt > 10 && pt <= 15) ipt  = 0;
    else if(pt > 15 && pt <= 20) ipt  = 1;
    else if(pt > 20 && pt <= 30) ipt  = 2;
    else if(pt > 30)             ipt  = 3;
    else printf("WARNING could not find pt bin for high pT correction\n");
    return ipt;
}

int UEShapesNoVtx::getEtaIndex(double eta, bool highpt=false) {
  // Convenient way to get eta index for Oldrich's code l. 261pp
  //https://svnweb.cern.ch/trac/atlasoff/browser/PhysicsAnalysis/StandardModelPhys/Minbias/InDetTrackCorr/trunk/Root/Corrections09_8TeV.cxx#L210 ?rev=557736  womgl
  int ieta = -1;
  if (!highpt) {
      double aeta = fabs(eta);
      if     (aeta < 1.3) ieta =0;
      else if(aeta < 1.9) ieta =1;
      else if(aeta < 2.1) ieta =2;
      else if(aeta < 2.3) ieta =3;
      else if(aeta < 2.5) ieta =4;
      else printf("WARNING could not find eta bin for material correction\n");
  }
  else {
    if      (eta < -2.25)        ieta = 0;
    else if (fabs(eta) <= 2.25)  ieta = 1;
    else if (eta < 2.5)          ieta = 2;
    else printf("WARNING could not find eta bin for high pT correction\n");
  }
  return ieta;
}

double UEShapesNoVtx::getTrkEffError(double eta, double pt) {
  //https://svnweb.cern.ch/trac/atlasoff/browser/PhysicsAnalysis/StandardModelPhys/Minbias/InDetTrackCorr/trunk/Root/Corrections09_8TeV.cxx#L210
  // Return eta-pT dependent tracking efficiency uncertainty
  double trkEff = getTrkEff(eta, pt, ""); // Nominal tracking efficiency
  // HS2014, neu:
  // NUR wenn pT < 800 |eta|<2.1:
  //  *  0.7% for all tracks with pT>800 MeV  and |eta|<2.1
  //  *  1.5% for all tracks with 500MeV < pT < 800 MeV and |eta|<2.1 
  //
  // TODO check das nochmal mit den relativen unsicherheiten usw.
  double err_highpt = 0.0;
  double err_chi2   = 0.0;
  double err_material = 0.0;
  double err_ptturnon = 0.0;

  bool doNew = true;

  if (doNew) {
    if (fabs(eta)<2.1) {
      if (pt > 0.5 && pt <= 0.8) err_highpt = 0.015*trkEff;
      else if (pt > 0.8) err_highpt = 0.007*trkEff;
      else {
        cout << "wtf?: " << pt << endl;
      }
    }
    else {
      int ieta = getEtaIndex(eta);
      
      // systematics from 10/6/2010
      // material systematics l. 281ff
      double material[5] = { 0.02 , 0.03, 0.04, 0.04 , 0.07  };  //500> MeV
      err_material = trkEff * material[ieta];

      // track selection and pt turn-on l. 295ff
      double ptturnon[5] = { 0.01, 0.01, 0.01, 0.01, 0.01};  //500> MeV
      err_ptturnon = trkEff * ptturnon[ieta];

      // BADLY MEASURED HIGH-PT tracks l. 316ff and 360ff
      double systTrkBad[4][3] =
      {
         {0.200,  0.003,  0.012},
         {0.300,  0.020,  0.250},
         {0.500,  0.040,  0.250},
         {0.800,  0.100,  0.450},
      };
      if (pt > 10.) {
        int ipt2   = getPtIndex(pt);
        int ieta2  = getEtaIndex(eta, true);
        err_highpt = trkEff * systTrkBad[ipt2][ieta2];
        err_chi2   = trkEff * 0.1;
      }
    }
  }
  // Old trk eff shit
  else {
      int ieta = getEtaIndex(eta);
      
      // systematics from 10/6/2010
      // material systematics l. 281ff
      double material[5] = { 0.02 , 0.03, 0.04, 0.04 , 0.07  };  //500> MeV
      err_material = trkEff * material[ieta];

      // track selection and pt turn-on l. 295ff
      double ptturnon[5] = { 0.01, 0.01, 0.01, 0.01, 0.01};  //500> MeV
      err_ptturnon = trkEff * ptturnon[ieta];

      // BADLY MEASURED HIGH-PT tracks l. 316ff and 360ff
      double systTrkBad[4][3] =
      {
         {0.200,  0.003,  0.012},
         {0.300,  0.020,  0.250},
         {0.500,  0.040,  0.250},
         {0.800,  0.100,  0.450},
      };
      if (pt > 10.) {
        int ipt2   = getPtIndex(pt);
        int ieta2  = getEtaIndex(eta, true);
        err_highpt = trkEff * systTrkBad[ipt2][ieta2];
        err_chi2   = trkEff * 0.1;
      }
  }

  // Add all errors in quadrature
  double err = 0.0;
  err += pow(err_material, 2);
  err += pow(err_ptturnon, 2);
  err += pow(err_highpt,   2);
  err += pow(err_chi2  ,   2);

  return sqrt(err);
}

double UEShapesNoVtx::getTrkEff(double eta, double pt, string syst="") {
  // Return eta-pT dependent tracking efficiency from histo mTrkEff read in earlier
  TAxis* pTAxis  = mTrkEff_pt_eta->GetXaxis();
  TAxis* etaAxis = mTrkEff_pt_eta->GetYaxis();

  int i_pt  = pTAxis ->FindBin(pt);
  int i_eta = etaAxis->FindBin(eta);

  // This reads out the central value
  double trkeff = mTrkEff_pt_eta->GetBinContent(i_pt, i_eta);

  // Sytematic variations within the errors
  if (syst == "up") trkeff += mTrkEff_pt_eta->GetBinError(i_pt, i_eta);
  else if (syst == "down") trkeff -= mTrkEff_pt_eta->GetBinError(i_pt, i_eta);
  // Use (1 - Delta TrkEff) as systematic variation (lower efficiency) for use
  // on reco level
  else if (syst == "diffdown") trkeff = 1. - mTrkEff_pt_eta->GetBinError(i_pt, i_eta);
  return trkeff;
}

double UEShapesNoVtx::getTrkEff() {
  // Return a flat Tracking efficiency value --- for testing purposes mainly  
  return 0.8;
}

bool UEShapesNoVtx::randomRemove(double eff) {
  // Return true or false based on Binomial P(n, p) = P(2, eff)
  int test = gRandom->Binomial(1, eff);
  fillh1f(hTrkSysRemoval + hs_TrkEff, 1.-test, 1.);
  fillh1f(hTrkSysEff + hs_TrkEff, eff, 1.);
  //return bool(1 - test);
  return bool(1 - test);
}

// Taken from MuonMomentumCorrections/Root/SmearingClass.cxx
void UEShapesNoVtx::SetRegions(){
  
  m_eta_min.clear();
  m_eta_max.clear();
  m_phi_min.clear();
  m_phi_max.clear();
  
  double tmpval;
  int i=0;
  std::string tmpname;
  std::string regions ="/nfs/dust/atlas/user/hschulz/rjets/2013/MuonMomentumCorrections/share/Regions.dat";
  
  InValues.open(regions.c_str());
  i=0;
  if(!InValues.good()){
    std::cerr<<"SmearingClass ERROR:: File "<<regions<<" NOT FOUND!!"<<std::endl;
    std::cout<<"                      Please provide a correct path!"<<std::endl;
    return;
  }
  else{
    while(!InValues.eof() && i<8){
      tmpval=0;
      if(i==0){
	getline(InValues,tmpname);
      }
      InValues>>tmpval;
      m_eta_min.push_back(tmpval);
      InValues>>tmpval;
      m_eta_max.push_back(tmpval);
      InValues>>tmpval;
      m_phi_min.push_back(tmpval);
      InValues>>tmpval;
      m_phi_max.push_back(tmpval);
      i++;
    }
  }
  InValues.close();
  InValues.clear();
  m_nb_regions=m_eta_min.size();
  
}

// Taken from MuonMomentumCorrections/Root/SmearingClass.cxx
int UEShapesNoVtx::getRegion(double eta,double phi){
  
  for(UInt_t k=0; k<m_nb_regions; k++){
    
    if( eta>=m_eta_min[k] && eta<m_eta_max[k] &&
	phi>=m_phi_min[k] && phi<m_phi_max[k] ){
      return k;
    }
  }
  
  return -1;
}

// Taken from MuonMomentumCorrections/Root/SmearingClass.cxx
void UEShapesNoVtx::FillValues(){
  /* Choose MuID or Staco Algorithm Values */
  
  double tmpval;
  int i=0;
  std::string tmpname;
  std::string data_val  = "/nfs/dust/atlas/user/hschulz/rjets/2013/MuonMomentumCorrections/share/Data11_values_muid_Rel17.dat";

  
  /* filling data vectors*/
  InValues.open(data_val.c_str());
  i=0;
  if(!InValues.good()){
    std::cerr<<"SmearingClass ERROR:: File "<<data_val<<" NOT FOUND!!"<<std::endl;
    std::cout<<"                      Please provide a correct path!"<<std::endl;
    return;
  }
  else{
    while(InValues.good()&& i<m_nb_regions){
      /*main values*/
      tmpval=0;
      if(i==0){
	getline(InValues,tmpname);
      }
      InValues>>tmpval;
      m_p1_ID.push_back(tmpval);
      InValues>>tmpval;
      m_p2_ID.push_back(tmpval);
      InValues>>tmpval;
      m_p2_ID_TAN.push_back(tmpval);
      InValues>>tmpval;
      m_p1_MS.push_back(tmpval);
      InValues>>tmpval;
      m_p2_MS.push_back(tmpval);
      /*errors*/
      InValues>>tmpval;
      m_E_p1_ID.push_back(tmpval);
      InValues>>tmpval;
      m_E_p2_ID.push_back(tmpval);
      InValues>>tmpval;
      m_E_p2_ID_TAN.push_back(tmpval);
      InValues>>tmpval;
      m_E_p1_MS.push_back(tmpval);
      InValues>>tmpval;
      m_E_p2_MS.push_back(tmpval);
      /*systematic*/
      InValues>>tmpval;
      m_S_p1_ID.push_back(tmpval);
      InValues>>tmpval;
      m_S_p2_ID.push_back(tmpval);
      InValues>>tmpval;
      m_S_p2_ID_TAN.push_back(tmpval);
      InValues>>tmpval;
      m_S_p1_MS.push_back(tmpval);
      InValues>>tmpval;
      m_S_p2_MS.push_back(tmpval);
    }
  }
  InValues.close();
  InValues.clear();
}

double UEShapesNoVtx::ErrorMatrix(int detRegion, double pt, double g3, double g4) {
  
  double vid=0.0;
  
  if (detRegion<0 || detRegion>=m_nb_regions)  return vid;

 
  double s1 = pow(m_E_p1_MS[detRegion]*m_E_p1_MS[detRegion]+
      	    m_S_p1_MS[detRegion]*m_S_p1_MS[detRegion],0.5);
  
  double s2 = pow(m_E_p2_MS[detRegion]*m_E_p2_MS[detRegion]+
      	    m_S_p2_MS[detRegion]*m_S_p2_MS[detRegion],0.5);
  
  double s3 = pow(m_E_p1_ID[detRegion]*m_E_p1_ID[detRegion]+
      	    m_S_p1_ID[detRegion]*m_S_p1_ID[detRegion],0.5);
  
  double s4 = pow(m_E_p2_ID[detRegion]*m_E_p2_ID[detRegion]+
      	    m_S_p2_ID[detRegion]*m_S_p2_ID[detRegion],0.5);
  
  
  vid = g3*g3*s3*s3+g4*g4*s4*s4*pt*pt;
  
  vid = pow(fabs(vid),0.5);
 
  return vid;
  
}

void UEShapesNoVtx::smearTrack(TVector3 &vec, int updown) {
  // Smear pT, eta, phi of track vector
  // See p. 304: https://twiki.cern.ch/twiki/pub/Atlas/AtlasTechnicalPaper/Published_version_jinst8_08_s08003.pdf
  // http://cdsweb.cern.ch/record/1379210/files/ATL-COM-PHYS-2011-662__revised.pdf
  double eta = vec.Eta();
  double phi = vec.Phi();
  int detRegion = getRegion(eta, phi);
  double g3 = gRand->Gaus(0,1); 
  double g4 = gRand->Gaus(0,1);

  double error = ErrorMatrix(detRegion, vec.Pt(), g3, g4);

  double dpTneu =  m_p1_ID[detRegion]*g3+m_p2_ID[detRegion]*g4*vec.Pt();
  if (updown <0) error*=-1.;
  double pTsmearneu = (1.+ dpTneu)*vec.Pt();
  //double pTsmearneu = (1.+ error)*vec.Pt();
  //cout << "error: " << error << " muidshit " << dpTneu << " thus " << pTsmearneu << " old : " << vec.Pt() << endl;
  //cout << "Before: " << vec.Pt()  <<  endl;

  //cout << "After neu: " << (1.+ dpTneu)*vec.Pt() << " " <<dpTneu/vec.Pt()*100 << endl;
  //double dpT = gRandom->Gaus(0., vec.Pt()*(vec.Pt()*0.0005 + 0.01));
  //double pTsmear    = vec.Pt()    + dpT;
  //cout << "After old: " << pTsmear << " " << dpT/vec.Pt() *100 << endl;
  //double pTsmear    = vec.Pt()    + gRandom->Gaus(0., vec.Pt()*(vec.Pt()*0.0005 + 0.01));

  // Assume 5 degree angular resolution
  double thetasmear = vec.Theta();// + gRandom->Gaus(0., 5.*M_PI/180.);
  double phismear   = vec.Phi();//   + gRandom->Gaus(0., 5.*M_PI/180.);    
  
  // Modify this tracks properties
  vec.SetPtThetaPhi(pTsmearneu, thetasmear, phismear); // TODO mach ma nur pT
}

void UEShapesNoVtx::smearTrack(vector<double> &trk) {
  // Smear pT, eta, phi of track vector
  // See p. 304: https://twiki.cern.ch/twiki/pub/Atlas/AtlasTechnicalPaper/Published_version_jinst8_08_s08003.pdf
  // http://cdsweb.cern.ch/record/1379210/files/ATL-COM-PHYS-2011-662__revised.pdf
  TVector3 vec = convertDoubletoTVector3(trk);
  double dpT = gRandom->Gaus(0., vec.Pt()*(vec.Pt()*0.0005 + 0.01));
  //cout << dpT << " " << vec.Pt() << " " << 100*dpT/vec.Pt() << endl;
  double pTsmear    = vec.Pt()    + dpT;
  //double pTsmear    = vec.Pt()    + gRandom->Gaus(0., vec.Pt()*(vec.Pt()*0.0005 + 0.01));

  // Assume 5 degree angular resolution
  double thetasmear = vec.Theta();// + gRandom->Gaus(0., 5.*M_PI/180.);
  double phismear   = vec.Phi();//   + gRandom->Gaus(0., 5.*M_PI/180.);    
  
  // Modify this tracks properties
  vec.SetPtThetaPhi(pTsmear, thetasmear, phismear);
  trk[0] = vec.Px();
  trk[1] = vec.Py();
  trk[2] = vec.Pz();
}

void UEShapesNoVtx::LoadPUlib() {
  //TFile* fpu;
  l.info("<UEShapes::LoadPUlib>");
  if (m_genPUlib) {
    //fpu = TFile::Open("test_2d_pulib.root", "recreate");  
    //pileUpLibrary = new TTree("pileUpLibrary","pileUpLibrary");
    //map<int, vector< vector<double> > > temp;
    //pileUpLibrary->Branch("tracks", &temp, 800, 1);
    //map <int, vector< vector<double> > > pumap;
    //for (int i=0; i<200; i++) {
      ////pumap.insert(MyPUMap::value_type(i+1, MyNtrkPULib()));
      ////pumap_vec.push_back(vector < vector<double> >);
      ////pumap[i] = MyNtrkPULib();
      //pumap[i] = vector<vector<double> >;
    //}
    //int num_of_n_bins = 200;
    //int num_of_row = 9;
    //double init_value = 0;

    //vector< vector<double> > matrix;
    ////now we have an empty 2D-matrix of size (0,0). Resizing it with one single command:
    //matrix.resize( num_of col , vector<double>( num_of_row , init_value ) );
    // and we are good to go ...
  }
  //fpu->Write();
  //fpu->Close();
}


void UEShapesNoVtx::WriteTrackLib(int n, vector<vector<vector<double> > > lib) {
  //string s = lexical_cast<string>("tracklibs/tracklib_n_");
  //s+= lexical_cast<string>(n);
  //std::ofstream ofs(s.c_str());
  //boost::archive::text_oarchive oa(ofs);
  //oa << lib;
}

void UEShapesNoVtx::LoadTrackLib() {
    // Code taken from http://www.dreamincode.net/code/snippet2295.htm
    // 
    const char* path = m_tracklibpath.Data();
    cout << "Loading tracklib from " << path << endl;
    // first off, we need to create a pointer to a directory
    DIR *pdir = NULL; // remember, it's good practice to initialise a pointer to NULL!
    pdir = opendir (path); // "." will refer to the current directory
    struct dirent *pent = NULL;
    if (pdir == NULL) printf ("\nERROR! pdir could not be initialised correctly");

    vector<string> contents;
    while ((pent = readdir (pdir)) != NULL) // while there is still something in the directory to list
    {
        if (pent == NULL) printf ("\nERROR! pent could not be initialised correctly");
        string curitem = lexical_cast<string>(path) + lexical_cast<string>("/")  +  lexical_cast<string>(pent->d_name);
        if (! is_directory(curitem)) {
          contents.push_back(lexical_cast<string>(pent->d_name));
          //printf ("%s\n", pent->d_name);
        }
    }
    closedir (pdir);

    // Iterate over files and load contents in PU map
    for (int i=0; i<contents.size();i++) {
      string file_to_read = lexical_cast<string>(path) + lexical_cast<string>("/")  +  contents[i];

      std::ifstream ifs(file_to_read.c_str());
      boost::archive::binary_iarchive ia(ifs);
      vector<vector<vector<double> > > thislib;
      ia >> thislib;

      vector<string> strs;
      boost::split(strs, contents[i], boost::is_any_of("_"));
      int zpos_vtx = lexical_cast<int>(strs[strs.size()-1]);
      //for (int j=0; j<thislib.size(); j++) cout << thislib[j].size() << endl;
      pumap[zpos_vtx] = thislib;
      cout << "\r" << 100.*(i+1)/contents.size() << "\% complete" << std::flush;
    }
    cout << endl << "Done" << endl;
    // Load the corresponding pulib stats
    TString path_pustats = m_tracklibpath + TString("_pulibstats") + TString(".root");
    TFile* fpustats = TFile::Open(path_pustats);
    m_pu_z      = (TH1F*)fpustats->Get("z_pu");
    m_pu_Dz_z   = (TH2F*)fpustats->Get("pu_Dz_z");
    m_pu_Dz_N   = (TH2F*)fpustats->Get("pu_Dz_N");
    m_pu_d_z   = (TH2F*)fpustats->Get("pu_Dz_z");
    m_pu_vtxN   = (TH1F*)fpustats->Get("vtxN_pu");
    //m_pu_vtxDz   = (TH2F*)fpustats->Get("vtxDz_pu");
    m_pu_vtxzPP   = (TH2F*)fpustats->Get("vtxzPP");
    m_pu_zN     = (TH2F*)fpustats->Get("z_puN");
    m_pu_ntrack = (TH1F*)fpustats->Get("h_ntrk_pu");
    m_pu_z     ->SetDirectory(0);
    for (int i=1; i<m_pu_z->GetNbinsX()+1;i++) {
        double newc =m_pu_z->GetBinContent(i)/m_pu_z->GetBinWidth(i);
        m_pu_z->SetBinContent(i, newc);}
    m_pu_z->Scale(1./m_pu_z->Integral()); // Turn it into pdf 
    m_pu_vtxN     ->SetDirectory(0);
    //m_pu_vtxDz     ->SetDirectory(0);
    m_pu_vtxzPP     ->SetDirectory(0);
    m_pu_zN    ->SetDirectory(0);
    m_pu_ntrack->SetDirectory(0);
    m_pu_Dz_z->SetDirectory(0);
    m_pu_Dz_N->SetDirectory(0);
    m_pu_d_z->SetDirectory(0);
    fpustats->Close();

    // The vertex z-resolution histogram #TODO: don't hardcode path
    //TFile* fres = TFile::Open("/scratch/hh/lustre/atlas/users/hschulz/ueinVB/UEinVBTABLR/tables/vertexresolution.root");
    //m_vtx_res = (TH1F*)fres->Get("res_vtx");
    //m_vtx_res->SetDirectory(0);
    //fres->Close();

}    

void UEShapesNoVtx::LoadSideBands() {
  if (m_sidebands_loaded) return;
  else {
    cout << "Loading side bands from " << m_pileuplibpath << endl;
    TFile* fsb = TFile::Open(m_pileuplibpath);
    m_sb_N   = (TH1F*)fsb->Get("hNtrkRecSideBand");
    m_sb_phi = (TH1F*)fsb->Get("hPhiTrkSideBand");
    m_sb_eta = (TH1F*)fsb->Get("hEtaTrkSideBand");
    m_sb_pT  = (TH1F*)fsb->Get("hPtTrkSideBand");
    m_sb_N  ->SetDirectory(0);
    m_sb_phi->SetDirectory(0);
    m_sb_eta->SetDirectory(0);
    m_sb_pT ->SetDirectory(0);
    fsb->Close();
    cout << "Done" << endl;
    m_sidebands_loaded = true;
  }
}

vector<double> UEShapesNoVtx::mkRdmTrack() {
  double phi = m_sb_phi->GetRandom();
  double eta = m_sb_eta->GetRandom();
  double pT =  m_sb_pT ->GetRandom();
  
  double px  = cos(phi) * pT;
  double py  = sin(phi) * pT;
  double theta  = atan(exp(-eta))*2.;
  double pz = sin(theta)*pT;

  vector<double> tv;
  tv.push_back(px);
  tv.push_back(py);
  tv.push_back(pz);

  return tv;
}

vector<double> UEShapesNoVtx::mirrorTrack(vector<double> track) {
  vector<double> mirror;
  mirror.push_back(-1.*track[0]);
  mirror.push_back(-1.*track[1]);
  mirror.push_back(-1.*track[2]);

  return mirror;
}


void UEShapesNoVtx::dumpStuff() {
  ofstream dumpFile;
  dumpFile.open(m_dumpRecoFile);
  cout << dump_vec.size() << " vs. " << dump_eventnumber.size() << endl;
  for (int i=0; i<dump_vec.size(); i++) {
      if (dump_vec[i].size()>0) {
          //dumpFile << lexical_cast<int>(dump_vec[i][0]);
          //dumpFile << dump_eventnumber[i] << " ";
          for (int j=0; j < dump_vec[i].size(); j++) {
              dumpFile << dump_vec[i][j];
              dumpFile << " "; 
          }
          dumpFile << endl;
      }
  }
  dumpFile.close();
  cout << "dumped " << dump_vec.size() << " events" << endl;
}


// CUTS
//___________________________________________________________________

//bool UEShapesNoVtx::isLeptonTrack(TVector3 trkvec, TVector3 lepvec) {
  //// Intended to find one and only one track associated with lepton
  //// Returns true if track is inside deltaR cone around lepton
  //// and has pT within 0.1% of lepton pT

  //double ptdiff = fabs((trkvec.Pt()-lepvec.Pt())/lepvec.Pt());

  //if ( (deltaR(lepvec,trkvec) < CUT_DELTAR) && 
       //(ptdiff < CUT_LEPPTRES) ) return true;

  //else return false;

//}

TVector3 UEShapesNoVtx::getDressedLepton(ITruth& mc, int index) {
    TLorentzVector lep_4vec      = getLepton4Vector(mc, index, m_muon);
    TLorentzVector lep_4vec_temp = getLepton4Vector(mc, index, m_muon);
    //cout << "Found Lepton" << endl;
    //cout << "4: " << lep_4vec.Pt() <<endl;
    //lep_4vec.Print();
    for (uint j = 0; j < mc.n(); j++ ) {
        if (index==j) continue;
        if (!(truthCutPrimary(mc, j)) ) continue;
        if (!(isPhoton(mc, j))) continue;
        //cout << mc[j].status << endl;
        TLorentzVector phot_4vec = getTrack4Vector(mc, j);
        double dR_lep_trk = lep_4vec_temp.DeltaR(phot_4vec);
        if (dR_lep_trk>=CUT_DRESSINGCONE) continue;
        //cout << "Delta R: " << dR_lep_trk << endl;
        //if (phot_4vec.E() > lep_4vec.E()) continue;
        //cout << "Found photon: " <<  endl;
        //phot_4vec.Print();
        else lep_4vec += phot_4vec;
    }
    TVector3 lep_3vec = getTrack3Vector(lep_4vec);
    //cout << "Finished dressing" << endl;
    //lep_4vec.Print();
    //cout << "3: " << lep_3vec.Pt() <<" 4: " << lep_4vec.Pt() << " " << m_muon << endl;
    return lep_3vec;
}

// fill PileUpLibrary
//___________________________________________________________________



//______________________________________________________________________
// Helper function that appends interesting events to a vector defined in main.cxx
void UEShapesNoVtx::rememberEvent(const char* reason, bool silent=false) {
  TString memo = TString::Format("%i %i %s %s", VB->run(), VB->event(), reason, 
          m_tree->GetCurrentFile()->GetName());
  if (!silent && !m_systype && !flag_QUIET) l.warning(memo);
  m_events_for_inspection.push_back(memo);
}

bool UEShapesNoVtx::isIsolatedMuon(TVector3 lep1, TVector3 lep2, vector<TVector3> tracks) {
    bool dec = true;

    double ptsum1 = 0.0;
    double ptsum2 = 0.0;


    //cout << m_CUT_ISO_DR << " " << m_CUT_ISO_PTFRAC << " " << m_CUT_ISO_PT << endl;
    for (uint i=0; i<tracks.size(); ++i) {
        double tpt = tracks[i].Pt();
        if (tpt > m_CUT_ISO_PT) {
            if (lep1.DeltaR(tracks[i]) < m_CUT_ISO_DR) ptsum1 += tpt;
            if (lep2.DeltaR(tracks[i]) < m_CUT_ISO_DR) ptsum2 += tpt;
            //if (lep2.DeltaR(tracks[i]) < 0.2) ptsum2 += tpt;
        }
    }

    if (ptsum1/lep1.Pt() >= m_CUT_ISO_PTFRAC || ptsum2/lep2.Pt() >= m_CUT_ISO_PTFRAC) dec=false;
    //if (ptsum1/lep1.Pt() >= 0.1 || ptsum2/lep2.Pt() >= 0.1) dec=false;

    return dec;
}

// MAIN
//______________________________________________________________________

Bool_t UEShapesNoVtx::Notify()
{
  master = 0; // reset local event counter for d3pdReader objects
  //cout << "Processing file " << this->m_tree->GetDirectory()->GetName() << endl;
  return kTRUE;
}


//______________________________________________________________________

Bool_t UEShapesNoVtx::ProcessCut(Long64_t entry)
{
  // if (flag_DEBUG) cout << "<UEShapesNoVtx::ProcessCut>" << endl;

  extern bool RECVD_KILL_SIGNAL;
  if (RECVD_KILL_SIGNAL) {
      this->Abort("user signal");
      return kFALSE;
  }
  

  // Check if this events number and run number match those required from the list earlier
  //if (!matchEventNumber(VB->event())) {
    //master++;
    //return kTRUE;
  //}
  dump_vec_event.clear();
    //return kTRUE;
  //}
  //if (!matchEventAndRunNumber(VB->event(),VB->run())) {
    //master++;
    //return kTRUE;
  //}

  //// Is this a W or Z event?
//#ifdef _Z_ANA_
  //// totally ignore W events if running Z analysis
  //if (VB->lep2_charge()==0) {
      ////cout << "ignoring W event" << endl;
    //master++;
    //return kTRUE;
  //}
//#else 
  //// totally ignore Z events if running W analysis
  //if (VB->lep2_charge()!=0) {
    //master++;
    //return kTRUE;
  //}
//#endif

  //cout << "Event Number -  " << VB->event() << endl; 
  evt_weight_reco = getEventWeight();
  evt_weight_gen =1.0;
  if (flag_isMC) evt_weight_gen = getEventWeight(true);

  //cout << evt_weight_reco << endl;
  // Compute the progress in data processing
  // ---------------------------------------
  nevt++;
  printEventProcessInfo();
  if (evt_weight_reco < 1.e-6) {
    master++;
    //cout << "W: " << evt_weight_reco << endl;
    return kTRUE;
  }
  //cout << "MASTER: " << master << endl;
  nevt_weighted+=1.*evt_weight_reco;



  // Loop over systematics
  // ---------------------------------------
  // Likely to be pushed deeper in to the Process method

  int begin=0;
  int stop =m_sysNames->size();
  if (m_sysNames->size()>1) begin++;
  for (int s=begin; s<stop; s++) {
  //for (int s=0; s<m_sysNames->size(); s++) {
    getSysType((*m_sysNames)[s]);
//    cout << "sysNames s = " << s << " spass: " << m_systype << " " <<  (*m_sysNames)[s] << endl; 
    
    // RJets reco cut
    bool RJetsSelectZFlag;
    RJetsSelectZFlag =  int(VB->lep2_tightPP())==1;

    bool doReco=true;
    if (RJetsSelectZFlag && doReco) {

    //    getSysType((*m_sysNames)[s]); // Fuck bug 2015-02-06 needed to be above

        initTransientVariables();

        if (!m_systype) {
          trkSelES_nominal.init();
          mcES_nominal.init();
        }

        // set BS displacement in z for pile-up vertex cuts
        double zBS = (flag_isMC) ? -6.11304 : BS->z();


       
        // Annoying declarations
        TVector3 lep1_cl_3vec, lep1_3vec, lep1_3vec_ptUp, lep1_3vec_ptDown;
        TVector3 lep2_cl_3vec, lep2_3vec, lep2_3vec_ptUp, lep2_3vec_ptDown;
        TVector3 lep1_3vec_trigger, lep1_3vec_reco ,lep1_3vec_id, lep2_3vec_trigger, lep2_3vec_reco, lep2_3vec_id;

        lep1_trk_3vec    = mkLepton(VB->lep1_trk_pt(),     VB->lep1_trk_eta(), VB->lep1_trk_phi(), true); //HS changed 20130823 from track pT to cluster pT
        lep2_trk_3vec    = mkLepton(VB->lep2_trk_pt(),     VB->lep2_trk_eta(), VB->lep2_trk_phi(), true);
        /*cout << "lep1_cl_E           : " << VB->lep1_cl_E           () << endl;
        cout << "lep1_f_smear          : " << VB->lep1_f_smear          () << endl;
        cout << "lep1_f_smear_up       : " << VB->lep1_f_smear_up       () << endl;
        cout << "lep1_f_smear_down     : " << VB->lep1_f_smear_down     () << endl;
        cout << "lep2_cl_E           : " << VB->lep2_cl_E           () << endl;
        cout << "lep2_f_smear          : " << VB->lep2_f_smear          () << endl;
        cout << "lep2_f_smear_up       : " << VB->lep2_f_smear_up       () << endl;
        cout << "lep2_f_smear_down     : " << VB->lep2_f_smear_down     () << endl;
        */
        //lep1_trk_3vec    = mkLepton(VB->lep1_cl_pt(),     VB->lep1_trk_eta(), VB->lep1_trk_phi());
        //lep2_trk_3vec    = mkLepton(VB->lep2_cl_pt(),     VB->lep2_trk_eta(), VB->lep2_trk_phi());
        //cout << VB->lep1_trk_eta() << " and " << VB->lep2_trk_eta() << endl;

        //cout << "lep1_E_reco           : " << VB->lep1_E_reco           () << endl;
        //cout << "lep1_cl_E/cosh(eta)   : " << VB->lep1_cl_E()/cosh(VB->lep1_trk_eta()) << endl;
        //cout << "lep1_trk_pt           : " << VB->lep1_trk_pt           () << endl;
        //cout << "lep1_cl_pt           : " << VB->lep1_cl_pt           () << endl;
        //cout << "lep1_E_var_ZeeStatUp  : " << VB->lep1_E_var_ZeeStatUp  () << endl;
        //cout << "lep1_E_var_ZeeMethodUp: " << VB->lep1_E_var_ZeeMethodUp() << endl;
        //cout << "lep1_E_var_ZeeGenUp   : " << VB->lep1_E_var_ZeeGenUp   () << endl;
        //cout << "lep1_E_var_R12StatUp  : " << VB->lep1_E_var_R12StatUp  () << endl;
        //cout << "lep1_E_var_PSUp       : " << VB->lep1_E_var_PSUp       () << endl;
        //cout << "lep2_E_reco           : " << VB->lep2_E_reco           () << endl;
        //cout << "lep2_E_reco           : " << VB->lep2_E_reco           () << endl;
        //cout << "lep2_trk_pt           : " << VB->lep2_trk_pt           () << endl;
        //cout << "lep2_cl_E/cosh(eta)   : " << VB->lep2_cl_E()/cosh(VB->lep2_trk_eta()) << endl;
        //cout << "lep2_cl_pt           : " << VB->lep2_cl_pt           () << endl;
        //cout << "lep2_E_var_ZeeStatUp  : " << VB->lep2_E_var_ZeeStatUp  () << endl;
        //cout << "lep2_E_var_ZeeMethodUp: " << VB->lep2_E_var_ZeeMethodUp() << endl;
        //cout << "lep2_E_var_ZeeGenUp   : " << VB->lep2_E_var_ZeeGenUp   () << endl;
        //cout << "lep2_E_var_PSUp       : " << VB->lep2_E_var_PSUp       () << endl;
        //cout << "lep2_E_var_R12StatUp  : " << VB->lep2_E_var_R12StatUp  () << endl;

        //VB->Print("all");


        // TODO cuts on lepton energies, eta acceptance in systematic variations

        if (m_lepton_sys==0 || m_lepton_sys==501 ||  m_lepton_sys==502 || m_lepton_sys==601 || m_lepton_sys==602
                || m_lepton_sys==503 ||  m_lepton_sys==504 || m_lepton_sys==603 || m_lepton_sys==604 || m_lepton_sys==701 || m_lepton_sys==702|| m_lepton_sys==801 || m_lepton_sys==802|| m_lepton_sys==901 || m_lepton_sys==902) {
            // Careful, in the nominal case we can use the variable VB->lep1_cl_pt directly,
            // in the muon channel this is identical with VB->lep1_trk_pt()
            // and since this is pt already, the bool at the end must be set to true
            // The same is true for the SF systematics
                // First lepton, nominal thingies 
                lep1_3vec        = mkLepton(VB->lep1_cl_pt(),      VB->lep1_trk_eta(), VB->lep1_trk_phi(), true);
                // Second lepton 
                lep2_3vec        = mkLepton(VB->lep2_cl_pt(),      VB->lep2_trk_eta(), VB->lep2_trk_phi(), true);
            //cout << "Using nominal lepton E/pT" << endl;
        }

        // Electron smearing up
        else if (m_lepton_sys ==11) { //
            if (m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
            lep1_3vec        = mkLepton(VB->lep1_cl_E()/VB->lep1_f_smear()*VB->lep1_f_smear_up(),      VB->lep1_trk_eta(), VB->lep1_trk_phi(), m_muon);
//              cout << VB->lep1_cl_E() << " vs: " << VB->lep1_cl_pt() << " vs: " << VB->lep1_cl_E()/cosh(VB->lep1_trk_eta()) << " vs: " << lep1_3vec.Pt() << endl;
            lep2_3vec        = mkLepton(VB->lep2_cl_E()/VB->lep2_f_smear()*VB->lep2_f_smear_up(),      VB->lep2_trk_eta(), VB->lep2_trk_phi(), m_muon);
        }
        // Electron smearing down
        else if (m_lepton_sys ==12) { //
            if (m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
            lep1_3vec        = mkLepton(VB->lep1_cl_E()/VB->lep1_f_smear()*VB->lep1_f_smear_down(),      VB->lep1_trk_eta(), VB->lep1_trk_phi(), m_muon);
            lep2_3vec        = mkLepton(VB->lep2_cl_E()/VB->lep2_f_smear()*VB->lep2_f_smear_down(),      VB->lep2_trk_eta(), VB->lep2_trk_phi(), m_muon);
        }
        // Electron rescaling ZeeStatUp
        else if (m_lepton_sys ==21) { //
            if (m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
            float E_nom1 = VB->lep1_cl_E();
            float E_nom2 = VB->lep2_cl_E();
            float E_syst1 = VB->lep1_f_smear()*VB->lep1_E_var_ZeeStatUp();
            float E_syst2 = VB->lep2_f_smear()*VB->lep2_E_var_ZeeStatUp();
            float DeltaE1 = fabs(E_nom1 - E_syst1);
            float DeltaE2 = fabs(E_nom2 - E_syst2);
            lep1_3vec        = mkLepton(E_nom1 + DeltaE1,      VB->lep1_trk_eta(), VB->lep1_trk_phi(), m_muon);
            lep2_3vec        = mkLepton(E_nom2 + DeltaE2,      VB->lep2_trk_eta(), VB->lep2_trk_phi(), m_muon);
        }
        
        // Electron rescaling ZeeStatUp down
        else if (m_lepton_sys ==22) { //
            if (m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
            float E_nom1 = VB->lep1_cl_E();
            float E_nom2 = VB->lep2_cl_E();
            float E_syst1 = VB->lep1_f_smear()*VB->lep1_E_var_ZeeStatUp();
            float E_syst2 = VB->lep2_f_smear()*VB->lep2_E_var_ZeeStatUp();
            float DeltaE1 = fabs(E_nom1 - E_syst1);
            float DeltaE2 = fabs(E_nom2 - E_syst2);
            lep1_3vec        = mkLepton(E_nom1 - DeltaE1,      VB->lep1_trk_eta(), VB->lep1_trk_phi(), m_muon);
            lep2_3vec        = mkLepton(E_nom2 - DeltaE2,      VB->lep2_trk_eta(), VB->lep2_trk_phi(), m_muon);
        }
        // Electron rescaling ZeeMethodUp
        else if (m_lepton_sys ==31) { //
            if (m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
            float E_nom1 = VB->lep1_cl_E();
            float E_nom2 = VB->lep2_cl_E();
            float E_syst1 = VB->lep1_f_smear()*VB->lep1_E_var_ZeeMethodUp();
            float E_syst2 = VB->lep2_f_smear()*VB->lep2_E_var_ZeeMethodUp();
            float DeltaE1 = fabs(E_nom1 - E_syst1);
            float DeltaE2 = fabs(E_nom2 - E_syst2);
            lep1_3vec        = mkLepton(E_nom1 + DeltaE1,      VB->lep1_trk_eta(), VB->lep1_trk_phi(), m_muon);
            lep2_3vec        = mkLepton(E_nom2 + DeltaE2,      VB->lep2_trk_eta(), VB->lep2_trk_phi(), m_muon);
        }
        
        // Electron rescaling ZeeMethodUp down
        else if (m_lepton_sys ==32) { //
            if (m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
            float E_nom1 = VB->lep1_cl_E();
            float E_nom2 = VB->lep2_cl_E();
            float E_syst1 = VB->lep1_f_smear()*VB->lep1_E_var_ZeeMethodUp();
            float E_syst2 = VB->lep2_f_smear()*VB->lep2_E_var_ZeeMethodUp();
            float DeltaE1 = fabs(E_nom1 - E_syst1);
            float DeltaE2 = fabs(E_nom2 - E_syst2);
            lep1_3vec        = mkLepton(E_nom1 - DeltaE1,      VB->lep1_trk_eta(), VB->lep1_trk_phi(), m_muon);
            lep2_3vec        = mkLepton(E_nom2 - DeltaE2,      VB->lep2_trk_eta(), VB->lep2_trk_phi(), m_muon);
        }
        // Electron rescaling ZeeGenUp
        else if (m_lepton_sys ==41) { //
            if (m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
            float E_nom1 = VB->lep1_cl_E();
            float E_nom2 = VB->lep2_cl_E();
            float E_syst1 = VB->lep1_f_smear()*VB->lep1_E_var_ZeeGenUp();
            float E_syst2 = VB->lep2_f_smear()*VB->lep2_E_var_ZeeGenUp();
            float DeltaE1 = fabs(E_nom1 - E_syst1);
            float DeltaE2 = fabs(E_nom2 - E_syst2);
            lep1_3vec        = mkLepton(E_nom1 + DeltaE1,      VB->lep1_trk_eta(), VB->lep1_trk_phi(), m_muon);
            lep2_3vec        = mkLepton(E_nom2 + DeltaE2,      VB->lep2_trk_eta(), VB->lep2_trk_phi(), m_muon);
        }
        
        // Electron rescaling ZeeGenUp down
        else if (m_lepton_sys ==42) { //
            if (m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
            float E_nom1 = VB->lep1_cl_E();
            float E_nom2 = VB->lep2_cl_E();
            float E_syst1 = VB->lep1_f_smear()*VB->lep1_E_var_ZeeGenUp();
            float E_syst2 = VB->lep2_f_smear()*VB->lep2_E_var_ZeeGenUp();
            float DeltaE1 = fabs(E_nom1 - E_syst1);
            float DeltaE2 = fabs(E_nom2 - E_syst2);
            lep1_3vec        = mkLepton(E_nom1 - DeltaE1,      VB->lep1_trk_eta(), VB->lep1_trk_phi(), m_muon);
            lep2_3vec        = mkLepton(E_nom2 - DeltaE2,      VB->lep2_trk_eta(), VB->lep2_trk_phi(), m_muon);
        }
        // Electron rescaling R12StatUp
        else if (m_lepton_sys ==51) { //
            if (m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
            float E_nom1 = VB->lep1_cl_E();
            float E_nom2 = VB->lep2_cl_E();
            float E_syst1 = VB->lep1_f_smear()*VB->lep1_E_var_R12StatUp();
            float E_syst2 = VB->lep2_f_smear()*VB->lep2_E_var_R12StatUp();
            float DeltaE1 = fabs(E_nom1 - E_syst1);
            float DeltaE2 = fabs(E_nom2 - E_syst2);
            lep1_3vec        = mkLepton(E_nom1 + DeltaE1,      VB->lep1_trk_eta(), VB->lep1_trk_phi(), m_muon);
            lep2_3vec        = mkLepton(E_nom2 + DeltaE2,      VB->lep2_trk_eta(), VB->lep2_trk_phi(), m_muon);
        }
        
        // Electron rescaling R12StatUp down
        else if (m_lepton_sys ==52) { //
            if (m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
            float E_nom1 = VB->lep1_cl_E();
            float E_nom2 = VB->lep2_cl_E();
            float E_syst1 = VB->lep1_f_smear()*VB->lep1_E_var_R12StatUp();
            float E_syst2 = VB->lep2_f_smear()*VB->lep2_E_var_R12StatUp();
            float DeltaE1 = fabs(E_nom1 - E_syst1);
            float DeltaE2 = fabs(E_nom2 - E_syst2);
            lep1_3vec        = mkLepton(E_nom1 - DeltaE1,      VB->lep1_trk_eta(), VB->lep1_trk_phi(), m_muon);
            lep2_3vec        = mkLepton(E_nom2 - DeltaE2,      VB->lep2_trk_eta(), VB->lep2_trk_phi(), m_muon);
        }
        // Electron rescaling PSUp
        else if (m_lepton_sys ==61) { //
            if (m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
            float E_nom1 = VB->lep1_cl_E();
            float E_nom2 = VB->lep2_cl_E();
            float E_syst1 = VB->lep1_f_smear()*VB->lep1_E_var_PSUp();
            float E_syst2 = VB->lep2_f_smear()*VB->lep2_E_var_PSUp();
            float DeltaE1 = fabs(E_nom1 - E_syst1);
            float DeltaE2 = fabs(E_nom2 - E_syst2);
            lep1_3vec        = mkLepton(E_nom1 + DeltaE1,      VB->lep1_trk_eta(), VB->lep1_trk_phi(), m_muon);
            lep2_3vec        = mkLepton(E_nom2 + DeltaE2,      VB->lep2_trk_eta(), VB->lep2_trk_phi(), m_muon);
        }
        
        // Electron rescaling PSUp down
        else if (m_lepton_sys ==62) { //
            if (m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
            float E_nom1 = VB->lep1_cl_E();
            float E_nom2 = VB->lep2_cl_E();
            float E_syst1 = VB->lep1_f_smear()*VB->lep1_E_var_PSUp();
            float E_syst2 = VB->lep2_f_smear()*VB->lep2_E_var_PSUp();
            float DeltaE1 = fabs(E_nom1 - E_syst1);
            float DeltaE2 = fabs(E_nom2 - E_syst2);
            lep1_3vec        = mkLepton(E_nom1 - DeltaE1,      VB->lep1_trk_eta(), VB->lep1_trk_phi(), m_muon);
            lep2_3vec        = mkLepton(E_nom2 - DeltaE2,      VB->lep2_trk_eta(), VB->lep2_trk_phi(), m_muon);
        }

        else if (m_lepton_sys ==921) { //
            if (!m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
            lep1_3vec        = mkLepton(VB->lep1_pt_idup()/1000.,      VB->lep1_trk_eta(), VB->lep1_trk_phi(), m_muon);
            lep2_3vec        = mkLepton(VB->lep2_pt_idup()/1000.,      VB->lep2_trk_eta(), VB->lep2_trk_phi(), m_muon);
            //cout << "Using muon id up syst pTs" << endl;
        }
        else if (m_lepton_sys ==922) { //
            if (!m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
            lep1_3vec        = mkLepton(VB->lep1_pt_iddn()/1000.,      VB->lep1_trk_eta(), VB->lep1_trk_phi(), m_muon);
            lep2_3vec        = mkLepton(VB->lep2_pt_iddn()/1000.,      VB->lep2_trk_eta(), VB->lep2_trk_phi(), m_muon);
            //cout << "Using muon id dn syst pTs" << endl;
        }
        else if (m_lepton_sys ==931) { //
            if (!m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
            lep1_3vec        = mkLepton(VB->lep1_pt_msup()/1000.,      VB->lep1_trk_eta(), VB->lep1_trk_phi(), m_muon);
            lep2_3vec        = mkLepton(VB->lep2_pt_msup()/1000.,      VB->lep2_trk_eta(), VB->lep2_trk_phi(), m_muon);
            //cout << "Using muon ms up syst pTs" << endl;
        }
        else if (m_lepton_sys ==932) { //
            if (!m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
            lep1_3vec        = mkLepton(VB->lep1_pt_msdn()/1000.,      VB->lep1_trk_eta(), VB->lep1_trk_phi(), m_muon);
            lep2_3vec        = mkLepton(VB->lep2_pt_msdn()/1000.,      VB->lep2_trk_eta(), VB->lep2_trk_phi(), m_muon);
            //cout << "Using muon ms dn syst pTs" << endl;
        }
        else if (m_lepton_sys ==941) { //
            if (!m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
            lep1_3vec        = mkLepton(VB->lep1_pt_noscale()/1000.,      VB->lep1_trk_eta(), VB->lep1_trk_phi(), m_muon);
            lep2_3vec        = mkLepton(VB->lep2_pt_noscale()/1000.,      VB->lep2_trk_eta(), VB->lep2_trk_phi(), m_muon);
        }
        // SCALEKUP
        else if (m_lepton_sys ==951) { //
            if (!m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
            lep1_3vec        = mkLepton(VB->lep1_f_smear_up()/1000.,      VB->lep1_trk_eta(), VB->lep1_trk_phi(), m_muon);
            lep2_3vec        = mkLepton(VB->lep2_f_smear_up()/1000.,      VB->lep2_trk_eta(), VB->lep2_trk_phi(), m_muon);
        }
        // SCALEKDOWN
        else if (m_lepton_sys ==952) { //
            if (!m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
            lep1_3vec        = mkLepton(VB->lep1_f_smear_down()/1000.,      VB->lep1_trk_eta(), VB->lep1_trk_phi(), m_muon);
            lep2_3vec        = mkLepton(VB->lep2_f_smear_down()/1000.,      VB->lep2_trk_eta(), VB->lep2_trk_phi(), m_muon);
        }
        // SCAECUP
        else if (m_lepton_sys ==961) { //
            if (!m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
            lep1_3vec        = mkLepton(VB->lep1_E_var_ZeeStatUp()/1000.,      VB->lep1_trk_eta(), VB->lep1_trk_phi(), m_muon);
            lep2_3vec        = mkLepton(VB->lep2_E_var_ZeeStatUp()/1000.,      VB->lep2_trk_eta(), VB->lep2_trk_phi(), m_muon);
        }
        // SCAECUP
        else if (m_lepton_sys ==962) { //
            if (!m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
            lep1_3vec        = mkLepton(VB->lep1_E_var_ZeeMethodUp()/1000.,      VB->lep1_trk_eta(), VB->lep1_trk_phi(), m_muon);
            lep2_3vec        = mkLepton(VB->lep2_E_var_ZeeMethodUp()/1000.,      VB->lep2_trk_eta(), VB->lep2_trk_phi(), m_muon);
        }

        // A posteriori cut on lepton phase-space
        
        //// Systematic bullshit
        //if (flag_isMC && !m_muon) { // electron channel
            //lep1_3vec_ptUp   = mkLepton(VB->lep1_cl_pt_up(),   VB->lep1_trk_eta(), VB->lep1_trk_phi());
            //lep1_3vec_ptDown = mkLepton(VB->lep1_cl_pt_down(), VB->lep1_trk_eta(), VB->lep1_trk_phi()); 
            //lep2_3vec_ptUp   = mkLepton(VB->lep2_cl_pt_up(),   VB->lep2_trk_eta(), VB->lep2_trk_phi());
            //lep2_3vec_ptDown = mkLepton(VB->lep2_cl_pt_down(), VB->lep2_trk_eta(), VB->lep2_trk_phi());
        //}
        //if (flag_isMC && m_muon) { // Muon channel
            //lep1_3vec_trigger = mkLepton(VB->corr_trigger_lep1(), VB->lep1_trk_eta(), VB->lep1_trk_phi());
            //lep1_3vec_reco    = mkLepton(VB->corr_reco_lep1(),    VB->lep1_trk_eta(), VB->lep1_trk_phi()); 
            //lep1_3vec_id      = mkLepton(VB->corr_id_lep1(),      VB->lep1_trk_eta(), VB->lep1_trk_phi()); 
            //lep2_3vec_trigger = mkLepton(VB->corr_trigger_lep2(), VB->lep2_trk_eta(), VB->lep2_trk_phi());
            //lep2_3vec_reco    = mkLepton(VB->corr_reco_lep2(),    VB->lep2_trk_eta(), VB->lep2_trk_phi());
            //lep2_3vec_id      = mkLepton(VB->corr_id_lep2(),      VB->lep2_trk_eta(), VB->lep2_trk_phi());
        //}

     //corr_trigger_lep1
     //corr_reco_lep1;
     //corr_id_lep1;
     //corr_trigger_lep2
     //corr_reco_lep2;
     //corr_id_lep2;
        // MET with and without HR
        //TVector3 met_3vec, met_hr_3vec;
        //met_3vec.SetPtEtaPhi(VB->met_et(),0.0,VB->met_phi());
        //met_hr_3vec.SetPtEtaPhi(VB->hr_met_et(),0.0,VB->hr_met_phi());

        double VB_mit_rec, VB_pt_rec, VB_pt_rec_hr, VB_pt_rec_up, VB_pt_rec_down, VB_rap_rec, VB_rap_rec_up, VB_rap_rec_down, VB_phi_rec;
        vector<double> VB_vec;
#ifdef _Z_ANA_
        // Nominal reco values
        //VB_mit_rec           = VB->mee();
        VB_pt_rec            = calculatePtVB(lep1_3vec,lep2_3vec,   m_muon);
        VB_rap_rec           = calculateYZ(lep1_3vec,lep2_3vec,     m_muon);
        VB_phi_rec           = calculatePhiZ(lep1_3vec,lep2_3vec,   m_muon);
        VB_vec               = getZThreeVector(lep1_3vec,lep2_3vec, m_muon);

        //cout << "Mu: " << VB_mit_rec << " " << calculateMeeZ(lep1_3vec, lep2_3vec, m_muon) << endl;
        //cout << "e : " << VB_mit_rec << " " << calculateMeeZ(lep1_3vec, lep2_3vec) << endl;
      
        //if (!flag_isMC)
        // Lepton cuts a posteriori
        
        VB_mit_rec = calculateMeeZ(lep1_3vec, lep2_3vec, m_muon);

        //if (m_muon) {
            //RJetsSelectZFlag = true;
        //}
        //else {
        //}

        if ((recoCutVB(VB_mit_rec, lep1_3vec, lep2_3vec) && RJetsSelectZFlag)) {
            // Now messing with reco weight
            //
            // 01 ... up up
            // 02 ... dn dn
            // 03 ... up dn
            // 04 ... dn up
            //
            if (m_lepton_sys==501) {
            if (m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
                double w_temp = evt_weight_reco/(VB->lep1_SFRECO()*VB->lep2_SFRECO());
                double f_syst = (VB->lep1_SFRECO() + VB->lep1_SFRECO_err())*(VB->lep2_SFRECO() + VB->lep2_SFRECO_err());
//                cout << "w_temp: " << w_temp << " f_syst: " << f_syst << endl;
                evt_weight_reco = w_temp*f_syst;
            }
            else if (m_lepton_sys==502) {
            if (m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
                double w_temp = evt_weight_reco/(VB->lep1_SFRECO()*VB->lep2_SFRECO());
                double f_syst = (VB->lep1_SFRECO() - VB->lep1_SFRECO_err())*(VB->lep2_SFRECO() - VB->lep2_SFRECO_err());
                evt_weight_reco = w_temp*f_syst;
            }
            else if (m_lepton_sys==503) {
            if (m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
                double w_temp = evt_weight_reco/(VB->lep1_SFRECO()*VB->lep2_SFRECO());
                double f_syst = (VB->lep1_SFRECO() + VB->lep1_SFRECO_err())*(VB->lep2_SFRECO() - VB->lep2_SFRECO_err());
                evt_weight_reco = w_temp*f_syst;
            }
            else if (m_lepton_sys==504) {
            if (m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
                double w_temp = evt_weight_reco/(VB->lep1_SFRECO()*VB->lep2_SFRECO());
                double f_syst = (VB->lep1_SFRECO() - VB->lep1_SFRECO_err())*(VB->lep2_SFRECO() + VB->lep2_SFRECO_err());
                evt_weight_reco = w_temp*f_syst;
            }
            if (m_lepton_sys==601) {
            if (m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
                double w_temp = evt_weight_reco/(VB->lep1_SFID()*VB->lep2_SFID());
                double f_syst = (VB->lep1_SFID() + VB->lep1_SFID_err())*(VB->lep2_SFID() + VB->lep2_SFID_err());
                evt_weight_reco = w_temp*f_syst;
            }
            else if (m_lepton_sys==602) {
            if (m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
                double w_temp = evt_weight_reco/(VB->lep1_SFID()*VB->lep2_SFID());
                double f_syst = (VB->lep1_SFID() - VB->lep1_SFID_err())*(VB->lep2_SFID() - VB->lep2_SFID_err());
                evt_weight_reco = w_temp*f_syst;
            }
            else if (m_lepton_sys==603) {
            if (m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
                double w_temp = evt_weight_reco/(VB->lep1_SFID()*VB->lep2_SFID());
                double f_syst = (VB->lep1_SFID() + VB->lep1_SFID_err())*(VB->lep2_SFID() - VB->lep2_SFID_err());
                evt_weight_reco = w_temp*f_syst;
            }
            else if (m_lepton_sys==604) {
            if (m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
                double w_temp = evt_weight_reco/(VB->lep1_SFID()*VB->lep2_SFID());
                double f_syst = (VB->lep1_SFID() - VB->lep1_SFID_err())*(VB->lep2_SFID() + VB->lep2_SFID_err());
                evt_weight_reco = w_temp*f_syst;
            }
            // Dielectron trigger SF up
            if (m_lepton_sys==901) {
              if (m_muon) { // Protection against srambling of muon and electron systematics
                  cout << "Error, do not use electron systematics with muons!!!" << endl;
                  exit(1);
              }
                evt_weight_reco *= VB->corr_trigger_lep1();
            }
            // Dielectron trigger SF dn
            else if (m_lepton_sys==902) {
            if (m_muon) { // Protection against srambling of muon and electron systematics
                cout << "Error, do not use electron systematics with muons!!!" << endl;
                exit(1);
            }
                evt_weight_reco *= VB->corr_trigger_lep2();
            }

            // MC muon scale factors 
            if (m_muon && flag_isMC) { // Protection against srambling of muon and electron systematics
                //cout <<"before: " << evt_weight_reco << endl;
                if (m_lepton_sys==0) {
                    evt_weight_reco *= VB->lep1_SFRECO(); // this is actually combined for both muons the recoSF
                    evt_weight_reco *= VB->lep1_SFID();   // this is actually combined for both muons the triggerSF
                }
                else if (m_lepton_sys==701) {
                    evt_weight_reco *= VB->lep1_SFRECO_err(); // this is actually combined for both muons the recoSF up variation
                    evt_weight_reco *= VB->lep1_SFID();   // this is actually combined for both muons the triggerSF
                }
                else if (m_lepton_sys==702) {
                    evt_weight_reco *= VB->lep2_SFRECO_err(); // this is actually combined for both muons the recoSF dn variation
                    evt_weight_reco *= VB->lep1_SFID();   // this is actually combined for both muons the triggerSF
                }
                else if (m_lepton_sys==801) {
                    evt_weight_reco *= VB->lep1_SFID_err(); // this is actually combined for both muons the recoSF up variation
                    evt_weight_reco *= VB->lep1_SFRECO();   // this is actually combined for both muons the triggerSF
                }
                else if (m_lepton_sys==802) {
                    evt_weight_reco *= VB->lep2_SFID_err(); // this is actually combined for both muons the recoSF dn variation
                    evt_weight_reco *= VB->lep1_SFRECO();   // this is actually combined for both muons the triggerSF
                }
                //cout <<"after:  " << evt_weight_reco << endl;
            }

            //cout << "passed reco cuts" << endl;
//            double pTl1 = VB->lep1_cl_pt();
//            double pTl2 = VB->lep2_cl_pt();
            double pTl1 = lep1_3vec.Pt();
            double pTl2 = lep2_3vec.Pt();
            fillh1f(hPtLept1 + hs_VBRecAll, pTl1,  evt_weight_reco);
            fillh1f(hPtLept2 + hs_VBRecAll, pTl2,  evt_weight_reco);
            fillh1f(hPtLept  + hs_VBRecAll, pTl1,  evt_weight_reco);
            fillh1f(hPtLept  + hs_VBRecAll, pTl2,  evt_weight_reco);
//            double etal1 = VB->lep1_trk_eta();
//            double etal2 = VB->lep2_trk_eta();
            double etal1 = lep1_3vec.Eta();
            double etal2 = lep2_3vec.Eta();
            fillh1f(hEtaLept1 + hs_VBRecAll, etal1,  evt_weight_reco);
            fillh1f(hEtaLept2 + hs_VBRecAll, etal2,  evt_weight_reco);
            fillh1f(hEtaLept  + hs_VBRecAll, etal1,  evt_weight_reco);
            fillh1f(hEtaLept  + hs_VBRecAll, etal2,  evt_weight_reco);
         


            // Back to back lepton pairs 
            double l_dphi = Dphi(lep1_3vec, lep2_3vec);
            fillh1f( hLeptonDPhi + hs_VBRecAll, l_dphi,  evt_weight_reco);
        
            if (fabs(l_dphi) >3.0) {
                fillh1f(hB2BLeptonEta  + hs_VBRecAll, lep1_3vec.Eta(),  evt_weight_reco);
                fillh1f(hB2BLeptonEta  + hs_VBRecAll, lep2_3vec.Eta(),  evt_weight_reco);
                fillh1f(hB2BLeptonPhi  + hs_VBRecAll, lep1_3vec.Phi(),  evt_weight_reco);
                fillh1f(hB2BLeptonPhi  + hs_VBRecAll, lep2_3vec.Phi(),  evt_weight_reco);
            }


            //if (VB_pt_rec > CUT_PTZ) {
                //continue; // HS: hier pTZ cut
            //}
            //else {
                //cout << "Keeping " << VB_pt_rec << endl;
            //}

            // Electron systematics
            //if (flag_isMC &! m_muon) {
                //VB_pt_rec_up   = calculatePtVB(lep1_3vec_ptUp,lep2_3vec_ptUp);
                //VB_pt_rec_down = calculatePtVB(lep1_3vec_ptDown,lep2_3vec_ptDown);
                
                //VB_rap_rec_up   = calculateYZ(lep1_3vec_ptUp,lep2_3vec_ptUp);
                //VB_rap_rec_down = calculateYZ(lep1_3vec_ptDown,lep2_3vec_ptDown);
            //}
            //
            
            VB_pt_rec_hr   = VB_pt_rec; // Leave in!!!!
#else 
            //VB_mit_rec     = VB->mt(); 
            //VB_pt_rec      = calculatePtVB(lep1_3vec,met_3vec);
            //VB_pt_rec_up   = calculatePtVB(lep1_3vec_ptUp,met_3vec);
            //VB_pt_rec_down = calculatePtVB(lep1_3vec_ptDown,met_3vec);
            //VB_pt_rec_hr   = calculatePtVB(lep1_3vec,met_hr_3vec);
#endif


            // fill trk-cl correlation plots
            fillh2f(hTrkClCorrPt  + hs_lep1, VB->lep1_trk_pt(),  VB->lep1_cl_pt(),  evt_weight_reco);
            fillh2f(hTrkClCorrEta + hs_lep1, VB->lep1_trk_eta(), VB->lep1_cl_eta(), evt_weight_reco);
            fillh2f(hTrkClCorrPhi + hs_lep1, VB->lep1_trk_phi(), VB->lep1_cl_phi(), evt_weight_reco);
            fillh2f(hTrkClCorrPt  + hs_lep2, VB->lep2_trk_pt(),  VB->lep2_cl_pt(),  evt_weight_reco);
            fillh2f(hTrkClCorrEta + hs_lep2, VB->lep2_trk_eta(), VB->lep2_cl_eta(), evt_weight_reco);
            fillh2f(hTrkClCorrPhi + hs_lep2, VB->lep2_trk_phi(), VB->lep2_cl_phi(), evt_weight_reco);

#ifdef _USE_RECOIL_
            VB_pt_rec  =  VB_pt_rec_hr;
            //met_3vec   =  met_hr_3vec;
#endif


            vector<TVector3>  track3vectors_all;
            vector<TVector3>  track3vectors_sel;
            vector<TVector3>  track3vectors_selall;
            vector<TVector3>  track3vectors_sideband;
            vector<TVector3>  track3vectors_PU;

            // Main loop over all vertices 
            // ---------------------------------------------------
            int nvtx = vtx->n();
            int mynvcs = 0;
            int mysvcs = 0;

            vector<int> PV_track_idcs = (*vtx)[0].trk_index();
            //cout << "tracks: " << PV_track_idcs.size() << endl;
            //for (uint t=0; t<PV_track_idcs.size(); ++t) {
                //TVector3 trk3vec = getTrack3Vector(*trk, t);
                //cout << t <<": " << trk3vec.Pt() << endl;
            //}
    //vtx->trk_index()->at(id).at(0);





          bool doTrackAtBeamLine=true;

          if (m_writeTrackLib) {
          
            double mindist = 60.;
            double DZCUT = 3.0;
            // A vertex loop over type-3 vertices only to find isolated vertices
            vector<int> PU_vtx_ids;
            vector<double> PU_vtx_zs;
            int nType3 = 0;
            for (int i=0; i<nvtx; i++) {
              int type      = (*vtx)[i].type();
              int vtx_ntrack = (*vtx)[i].nTracks();

              // Exclude vertices without tracks and non-type3-vertices
              if (type!=3 || vtx_ntrack ==0) continue;
              nType3++;
              //if (type!=3) continue;
              double vz = 0.0;
              //if (doTrackAtBeamLine) {
                  vz = vtx->z()->at(i) -zBS; // Assume this to be the z position of the electron pair
              //}
              //else {
                  //vz = vtx->z()->at(i); // Assume this to be the z position of the electron pair
              //}
              double deltaz = vtx->z()->at(i) - vtx->z()->at(0);
              // Fill deltaZ histo for merging probability estimation
              //m_pu_vtxDz->Fill(vtx->z()->at(0) - vtx->z()->at(i), vtx->z()->at(0));
              // Require the PU vertex to be far away from the Primary vertex
              if(fabs(deltaz) < mindist) continue;
              PU_vtx_ids.push_back(i);
              PU_vtx_zs.push_back(vz);
            } // nvtx
            m_pu_vtxN->Fill(nType3);

            assert(PU_vtx_zs.size() == PU_vtx_ids.size());

            // Iterate over type3 vertices to find isolated vertices
            vector<int> isolated_vtx_ids;
            isolated_vtx_ids.clear();

            for(int i=0;i<PU_vtx_zs.size();i++) {
              bool isol = true;
              for(int j=0;j<PU_vtx_zs.size();j++) {
                // Prevent double counting and comparison with itself
                if (i!=j && j>i) {
                  double distance = fabs(PU_vtx_zs[PU_vtx_ids[i]] - PU_vtx_zs[PU_vtx_ids[j]]);
                  if (distance < mindist) {
                    isol=false;
                    continue;
                  }
                }
              }
              if (isol) isolated_vtx_ids.push_back(PU_vtx_ids[i]);
            } // vertices

            // The vector 'isolated_vtx_ids' contains the vtx indices of isolated type-3 vertices
            for (int i=0;i<isolated_vtx_ids.size();i++) {
              int id = isolated_vtx_ids[i];
              int first_index = vtx->trk_index()->at(id).at(0);
              //int count = (*vtx)[id].trk_n();
              double zpos = 0.0;
              //if (doTrackAtBeamLine) {
                  zpos = (*vtx)[id].z() -zBS;
                  m_pu_vtxzPP->Fill(zpos, vtx->z()->at(0)-zBS);
              //}
              //else {
                  //zpos = (*vtx)[id].z();
                  //m_pu_vtxzPP->Fill(zpos, vtx->z()->at(0));
              //}
              m_pu_z->Fill(zpos);
              
              vector<TVector3> temp;
              vector<double> temp_z;
              vector<double> temp_theta;
              // This loops over the  tracks associated with the pu vertex
              //for(int n = first_index; n < first_index+count; n++) {
                //// Apply InDet, PhaseSpace and d0 cuts  
                //if (trkCutSel(*trk, n, zpos)) {
                  //TVector3 trk3vec = getTrack3Vector(*trk, n);
                  //temp.push_back(trk3vec);
                //}
              //}
              // This loops over all tracks
              for (int n=0; n<trk->n(); n++) {
                // Apply InDet, PhaseSpace and d0 cuts  
                
                if (doTrackAtBeamLine) {
                    if (trkCutPUVertexBL(*trk, n, zpos, DZCUT)) {
                      TVector3 trk3vec = getTrack3Vector(*trk, n);
                      temp.push_back(trk3vec);
                      m_pu_Dz_z->Fill(trkDZPUVertexBL(*trk, n, zpos), zpos);
                      m_pu_d_z->Fill((*trk)[n].d0_wrtBL(), zpos);
                      //cout << (*trk)[n].d0_wrtBL() << " " << (*trk)[n].d0_wrtPV() << " " << (*trk)[n].d0_wrtBL() - (*trk)[n].d0_wrtPV() << endl;
                      //temp_z.push_back(trkDZPUVertexBL(*trk, n, zpos));
                      //temp_z.push_back(zpos - (*trk)[n].z0_wrtBL());
                      temp_z.push_back((*trk)[n].z0_wrtBL()-zpos);
                      temp_theta.push_back((*trk)[n].theta());
                    }
                }
                else {
                    if (trkCutSel(*trk, n, zpos, DZCUT)) {
                      TVector3 trk3vec = getTrack3Vector(*trk, n);
                      temp.push_back(trk3vec);
                      m_pu_Dz_z->Fill(trkDZPUVertexBL(*trk, n, zpos), zpos);
                      m_pu_d_z->Fill((*trk)[n].d0_wrtBL(), zpos);
                      //cout << (*trk)[n].d0_wrtBL() << " " << (*trk)[n].d0_wrtPV() << " " << (*trk)[n].d0_wrtBL() - (*trk)[n].d0_wrtPV() << endl;
                      //temp_z.push_back(trkDZPUVertexBL(*trk, n, zpos));
                      //temp_z.push_back(zpos - (*trk)[n].z0_wrtBL());
                      temp_z.push_back((*trk)[n].z0_wrtPV()-zpos);
                      temp_theta.push_back((*trk)[n].theta());
                    }
                } // not track at beam line
              } // track loop

              // Fill resolution histo
              for (int n=0; n<temp.size(); n++) {
                m_pu_Dz_N->Fill(temp_z[n]*sin(temp_theta[n]), temp.size());
              }

              //m_pu_z->Fill(zpos);
              m_pu_zN->Fill(zpos, temp.size());
              vector< vector<double> > vtx_trackvectors = convertTVector3toDouble(&temp);

              int zpos_idx = int(floor(zpos));
              pumap[zpos_idx].push_back(zipVectors(vtx_trackvectors, temp_z, temp_theta));
              //pumap[temp.size()].push_back(vtx_trackvectors);
              //cout << temp.size() << " " << pumap[temp.size()].size() << " " << pumap[temp.size()][pumap[temp.size()V].size()].size()
              //assert(temp.size() == pumap[temp.size()][pumap[temp.size()].size()].size());
              pu_tracks += temp.size();
            } // isolated vertices
          } // write Tracklib

          // The actual analysis
          else {
            vector<int> trkindices;
            trkindices.clear(); 
            double z_vtx = vtx->z()->at(0) -zBS;
            for (int i=0; i<trk->n(); i++) {
          
              // Loop over the tracks associated to the vertices
              // -----------------------------------------------------
              //bool trackAtBeamLine=true;
              bool trackAtBeamLine=false; // AUA!!!
              if (trackAtBeamLine) { 
                  double z_trk = (*trk)[i].z0_wrtBL();
                  //double deltaz_trk = vtx->z()->at(i) - vtx->z()->at(0);
                  double theta = (*trk)[i].theta();
                  double deltaz_trk = (z_vtx - z_trk)*sin(theta);
                  //double deltaz_trk = (z_vtx - z_trk -zBS);

                  /*
                  double z_trk = (*trk)[i].z0_wrtPV();
                  double theta = (*trk)[i].theta();
                  double deltaz_trk = z_trk*sin(theta); // The sin theta thing
                  */
                  //double cut_trk = z_trk*sin(theta);
                  //double deltaz_trk = z_vtx - z_trk;
                  //double deltaz_trk = z_trk;              // the wrtPV thing
                  //double deltaz_trk = z_vtx - z_trk;
                  //double deltaz_trk = (z_vtx - z_trk)*sin(theta);
               
                  //doTrackAnalysis(hsTrkSel, i, track3vectors_sel, deltaz_trk);
                  //doTrackAnalysis(hsTrkSideBand, i, track3vectors_sideband, deltaz_trk);
                  //
                  // HS: hier ist der z0 cut!!!
                  //
                  if (fabs(deltaz_trk) < 0.5 ) doTrackAnalysis(hsTrkSel, i, track3vectors_sel, deltaz_trk); // This cuts on delta z
                  if (fabs(deltaz_trk) < m_sbcut) doTrackAnalysis(hsTrkSideBand, i, track3vectors_sideband, deltaz_trk);
              }
              else {
                  if (doTrackAnalysisLegacy(hsTrkSel, i, track3vectors_sel)) {; // This cuts on delta z
                      trkindices.push_back(i); // this is ony for the print function
                  }
              }
            } // end loop over tracks
            
             ////Cut on muon isolation
            //bool cutMuonIso=false; // false per default
            bool cutMuonIso=true; // false per default --- set to true for lepton cone stuff

//            if (cutMuonIso &! m_muon) {
//              cout << "ERROR, asking for muon cut in electron analysis" << endl;
//              exit(1);
//            }

            if (cutMuonIso &! isIsolatedMuon(lep1_3vec, lep2_3vec, track3vectors_sel) ) {
                ////cout << "Event failed muon isolation cut" << endl;
                break;
            }


            //printEvent(VB_mit_rec, VB_pt_rec, VB_phi_rec, trkindices);

            fillh1f(hVtxN+hs_RecVtxAll,  mynvcs,  evt_weight_reco);
            fillh1f(hVtxN+hs_RecVtxPU,   mysvcs,  evt_weight_reco);

            // fill trk-cl correlation plots for dressed leptons
            fillh2f(hTrkClCorrSumPt + hs_lep1, lep1_cone_sumpt,  VB->lep1_cl_pt(),  evt_weight_reco);
            fillh2f(hTrkClCorrSumPt + hs_lep2, lep2_cone_sumpt,  VB->lep2_cl_pt(),  evt_weight_reco);

            // fill VB mass and pT histos
            fillh1f(hMitVB +hs_VBRecAll, VB_mit_rec, evt_weight_reco);
            fillh2f(hMitVBPtVB +hs_VBRecAll, VB_mit_rec, VB_pt_rec, evt_weight_reco);
            fillh1f(hPtVB  +hs_VBRecAll, VB_pt_rec,  evt_weight_reco);

            if (track3vectors_sel.size() >= 2) {
              fillh1f(hMitVB +hs_VBRecNsel2, VB_mit_rec, evt_weight_reco);
              fillh1f(hPtVB  +hs_VBRecNsel2, VB_pt_rec,  evt_weight_reco);
            }



            vector< vector<double> > stl_vec_pileUpLibrary;
            if (m_systype == kPileUp) {
              ntrk_sig+=track3vectors_sel.size();
              //LoadSideBands();
              // Add PU tracks  N times 
              for(int i = 0; i < m_maxAdditionalPileUp; i++) {
                  //cout << " i " << i << endl;

                  double z_vtx = vtx->z()->at(0) -zBS;

                  int n_addVertices = int(floor(m_pu_vtxN->GetRandom()+0.5));
                  fillh1f(hNaddVtx + hsTrkAfter, n_addVertices, 1.0);
                  
                  map<int, vector<vector<vector<double> > > >::const_iterator itr_pumap;
                  for (int i=0;i<n_addVertices; i++) {
                    TH1D* projn = m_pu_zN->ProjectionY();

                    //double rdmZ=m_pu_vtxzPP->ProjectionX()->GetRandom();
                    double rdmZ=m_pu_vtxzPP->ProjectionY()->GetRandom();
                    while ( pumap.find(rdmZ) == pumap.end() ) {
                      rdmZ=m_pu_vtxzPP->ProjectionY()->GetRandom(); // AUa!
                    }

                    bool ok = false;
                    int n_zidx = int(floor(rdmZ));
                    // Pick a random vertex from that vertex position
                    int rdmEntry = int(floor(gRandom->Uniform(pumap[n_zidx].size())));
                    vector<vector<double> > vtx_add;
                    //cout << "n_zidx" <<  n_zidx << endl;
                    vector<vector<double> > temp = pumap[n_zidx][rdmEntry];

                    // Loop ovetr its tracks and cut on |(ztrk - zPri)*sin(theta_trk)| <1.5
                    for (int t=0;t<temp.size();t++) {
                      double trk_deltaZ = temp[t][3];
                      double trk_theta = temp[t][4]; 
                      //double zstwrtPV =  (z_vtx + rdmZ + trk_deltaZ)*fastsin(trk_theta);
                      double zstwrtPV =  (z_vtx + rdmZ + trk_deltaZ)*sin(trk_theta);
                      //cout << z_vtx << " " << rdmZ << " " << trk_deltaZ << " " << zstwrtPV << endl;
                      if (fabs(zstwrtPV) < 1.5) {
                        vector<double> mytrack;
                        for (int k=0; k<3;k++) mytrack.push_back(temp[t][k]);
                        vtx_add.push_back(mytrack);
                        fillh1f(hDeltaZwrtPVPeak + hsTrkAfter, zstwrtPV,1.0);
                      }
                    } // end track loop
                    fillh1f(hZPri + hsTrkAfter, z_vtx, 1.0);
                    fillh2f(hNaddZ + hsTrkAfter,  z_vtx, vtx_add.size(), 1.0);
                    stl_vec_pileUpLibrary.insert(stl_vec_pileUpLibrary.end(), vtx_add.begin(), vtx_add.end());
                  } // end vertices loop

                  ntrk_add+= stl_vec_pileUpLibrary.size();
                  ntrk_sig+= stl_vec_pileUpLibrary.size();
              } // end loop over additional HBOM steps
            } // kPileUp

            // add additional track weights!
            for (int pui=0; pui<stl_vec_pileUpLibrary.size(); pui++) {
              TVector3 tmp = TVector3(stl_vec_pileUpLibrary.at(pui)[0],
                                      stl_vec_pileUpLibrary.at(pui)[1],
                                      stl_vec_pileUpLibrary.at(pui)[2]);

              trackweights.push_back(getTrackWeight(TVector3(stl_vec_pileUpLibrary.at(pui)[0],
                                                             stl_vec_pileUpLibrary.at(pui)[1],
                                                             stl_vec_pileUpLibrary.at(pui)[2])));
            }
            
            //EventShapes trkRdmES(stl_vec_pileUpLibrary);
            //fillEventShapes(hsRecRdm,  &trkRdmES,  VB_pt_rec, vtx->z()->at(0)-zBS ); 
              
            // Reco EventShape calculations
            // ---------------------------------------
            vector< vector<double> > stl_vec = convertTVector3toDouble(&track3vectors_all);
            //vector< vector<double> > stl_vec = convertTVector3toDouble(&track3vectors_all);
            stl_vec.insert(stl_vec.end(), stl_vec_pileUpLibrary.begin(), stl_vec_pileUpLibrary.end());

            EventShapes trkAllES(stl_vec);
            stl_vec.clear();
            
            // This part does the event shapes for the sideband region
            stl_vec = convertTVector3toDouble(&track3vectors_sideband);
            //stl_vec.insert(stl_vec.end(), stl_vec_pileUpLibrary.begin(), stl_vec_pileUpLibrary.end());
            EventShapes trkSideBandES(stl_vec);
            stl_vec.clear();
            

            // Reco Tracking efficiency application
            if (m_systype == kTrkEff) {
                int n_before = track3vectors_sel.size();
                for (unsigned int pos=0; pos < track3vectors_sel.size(); pos++) {
                    double test_eta = track3vectors_sel[pos].Eta();
                    double test_pt  = track3vectors_sel[pos].Pt();
                    double prob = getTrkEffError(test_eta, test_pt);
                    if (randomRemove(1. - prob)) {
                        //cout << prob << " " << test_eta << " " << test_pt << endl;
                        track3vectors_sel.erase(track3vectors_sel.begin() + pos);
                        trackweights.erase(trackweights.begin() + pos);
                    }
                }
                int n_after = track3vectors_sel.size();
                  //cout << "Throwing away track" << endl;
                //cout << " Before: " << n_before << " After: " << n_after << " removed " << n_before - n_after << " tracks" << endl;
            } // kTrkEff
            else if (m_systype == kTrkSmearUp) {
              vector<TVector3>  temp;
              int n_before = track3vectors_sel.size();
              for (int it_ntrk=0; it_ntrk<track3vectors_sel.size(); it_ntrk++) {
                //cout << "Before: " << track3vectors_sel[it_ntrk].Pt() << endl;
                smearTrack(track3vectors_sel[it_ntrk], 1);
                //cout << "After: " << track3vectors_sel[it_ntrk].Pt() << endl;
                double smeared = track3vectors_sel[it_ntrk].Pt();
                if (smeared > 0.5) temp.push_back(track3vectors_sel[it_ntrk]);
              }
              track3vectors_sel = temp;
            }
            else if (m_systype == kTrkSmearDown) {
              vector<TVector3>  temp;
              for (int it_ntrk=0; it_ntrk<track3vectors_sel.size(); it_ntrk++) {
                smearTrack(track3vectors_sel[it_ntrk], -1);
                double smeared = track3vectors_sel[it_ntrk].Pt();
                if (smeared > 0.5) temp.push_back(track3vectors_sel[it_ntrk]);
              }
              track3vectors_sel = temp;
            }


            /*// HS 2014 
            bool smearTracks=true;

            if (smearTracks) {
              vector<TVector3>  temp;
              for (int it_ntrk=0; it_ntrk<track3vectors_sel.size(); it_ntrk++) {
                smearTrack(track3vectors_sel[it_ntrk]);
                double smeared = track3vectors_sel[it_ntrk].Pt();
                if (smeared > 0.5) temp.push_back(track3vectors_sel[it_ntrk]);
              }
              stl_vec =convertTVector3toDouble(&temp);
            }
            else {
              stl_vec = convertTVector3toDouble(&track3vectors_sel);
            }
            */
            // HS
            for (int i_holger=0; i_holger<track3vectors_sel.size(); i_holger++) {
              fillTrackHistos(hsTrkSel,  track3vectors_sel[i_holger], evt_weight_reco);
            }


            stl_vec = convertTVector3toDouble(&track3vectors_sel);
            stl_vec.insert(stl_vec.end(), stl_vec_pileUpLibrary.begin(), stl_vec_pileUpLibrary.end());
            //// Forward closure cut on PV tracks
            //if (stl_vec.size() < 10 || stl_vec.size() > 30) continue;

            //EventShapes trkSelES(stl_vec);
            //EventShapes trkSelES(stl_vec, VB_phi_rec);
            EventShapes trkSelES(stl_vec, VB_vec, true);
            //trkSelES.print();
            EventShapes trkCorrES(stl_vec, trackweights);

            //trkSelES.print();

            // HS 2014 --- pT distribution of low multiplicity events
            if (trkSelES.Ntrk()>0 && trkSelES.Ntrk()<=5) {
              //cout << trkSelES.Ntrk() << " vs: " << stl_vec.size() << endl;
              for (int i=0; i<stl_vec.size();i++) {
                //cout << stl_vec[i][0] << " " << stl_vec[i][1] << stl_vec[i][2] << endl;
                double pt = sqrt(pow(stl_vec[i][0],2) + pow(stl_vec[i][1],2));
                //cout << "pt: " << pt << " " << endl;
                fillh1f(hPtLowNtrk+hsRecSel, pt, evt_weight_reco );
              }
            }

            // HS2014 --- Oleg request, n,sumpt cone removal stuff
            fillh1f(hConeN+hsRecSel, lep1_cone_n+lep2_cone_n, evt_weight_reco );
            fillh1f(hConeSumPt+hsRecSel, lep1_cone_sumpt+lep2_cone_sumpt, evt_weight_reco );

            if ( (lep1_cone_n+lep2_cone_n) <2) {
              cout << "############ UPSI = " << lep1_cone_n+lep2_cone_n <<endl;
              cout << "lep1 trk pt = " << lep1_trk_3vec.Pt() << endl;
              cout << "lep1 trk eta = " << lep1_trk_3vec.Eta() << endl;
              cout << "lep1 trk phi = " << lep1_trk_3vec.Phi() << endl;

              cout << "lep1  pt = " << lep1_3vec.Pt() << endl;
              cout << "lep1  eta = " << lep1_3vec.Eta() << endl;
              cout << "lep1  phi = " << lep1_3vec.Phi() << endl;

              cout << "lep2 trk pt = " << lep2_trk_3vec.Pt() << endl;
              cout << "lep2 trk eta = " << lep2_trk_3vec.Eta() << endl;
              cout << "lep2 trk phi = " << lep2_trk_3vec.Phi() << endl;

              cout << "lep2  pt = " << lep2_3vec.Pt() << endl;
              cout << "lep2  eta = " << lep2_3vec.Eta() << endl;
              cout << "lep2  phi = " << lep2_3vec.Phi() << endl;
            }

            if (!m_systype) trkSelES_nominal = trkSelES; // store event shapes for systematics


            //cout << "sys: " << m_systype<< " Ntrk: " <<trkSelES.Ntrk() << endl;


            fillEventShapes(hsRecSel,  &trkSelES,  VB_pt_rec, evt_weight_reco, 0.0, true); 
            fillEventShapes3D(hsRecSel3D,  &trkSelES,  VB_pt_rec, VB_mit_rec, evt_weight_reco, 0.0, true); 
        //////    
            //fillEventShapes(hsRecAll,  &trkAllES,  VB_pt_rec); 
            //fillEventShapes(hsRecCorr, &trkCorrES, VB_pt_rec); 
            //fillEventShapes(hsRecSideBand,  &trkSideBandES,  VB_pt_rec, vtx->z()->at(0)-zBS ); 
        //////// HS switched off April 16 2013
         
            // Imagiro Trees
            // ---------------------------------------
            //fillImagiroTree(imagiro_tree_rec, imagiro_tree_rec2, &trkSelES, VB_pt_rec, 
                            //make_pair(VB_pt_rec_up,VB_pt_rec_down));
            if (m_systype == kTrkEff) {
                fillImagiroTree(imagiro_tree_recTrkEff, imagiro_tree_rec2TrkEff,  &trkSelES, VB_rap_rec,  VB_pt_rec, evt_weight_reco, 
                                make_pair(VB_pt_rec_up,VB_pt_rec_down));
            }
           /* else if (m_systype == kTrkSmearUp) {
                fillImagiroTree(imagiro_tree_recTrkSmearUp, imagiro_tree_rec2TrkSmearUp,  &trkSelES, VB_rap_rec,  VB_pt_rec, evt_weight_reco, 
                                make_pair(VB_pt_rec_up,VB_pt_rec_down));
            }
            else if (m_systype == kTrkSmearDown) {
                fillImagiroTree(imagiro_tree_recTrkSmearDown, imagiro_tree_rec2TrkSmearDown,  &trkSelES, VB_rap_rec,  VB_pt_rec, evt_weight_reco, 
                                make_pair(VB_pt_rec_up,VB_pt_rec_down));
            }
            */
            else {
                fillImagiroTree(imagiro_tree_rec,        imagiro_tree_rec2,        &trkSelES, VB_rap_rec, VB_pt_rec, evt_weight_reco, 
                                make_pair(VB_pt_rec_up,VB_pt_rec_down));
                // PtZ up
                //fillImagiroTree(imagiro_tree_recPtZUp,   imagiro_tree_rec2PtZUp,   &trkSelES, VB_rap_rec, VB_pt_rec_up, evt_weight_reco, 
                    //            make_pair(VB_pt_rec_up,VB_pt_rec_down));
                // PtZ down
                //fillImagiroTree(imagiro_tree_recPtZDown, imagiro_tree_rec2PtZDown, &trkSelES, VB_rap_rec, VB_pt_rec_down, evt_weight_reco, 
                  //              make_pair(VB_pt_rec_up,VB_pt_rec_down));
            }

            // Contamination Correction
            // ---------------------------------------
            fillCTCorrection(&trkSelES, evt_weight_reco);

            // Smearing of event shapes due to systematics
            // -------------------------------------------
            fillEventShapeSmear(hs_RecSyst, &trkSelES_nominal, &trkSelES, VB_pt_rec, VB_pt_rec, evt_weight_reco);
            fillEventShapeDiff (hs_RecSyst, &trkSelES_nominal, &trkSelES, VB_pt_rec, VB_pt_rec, evt_weight_reco);
        } // actual analysis
      } // VB reco cut
    } // RJetSelectZFlag

    // Loop over the mc particles
    // ---------------------------------------
    if (flag_isMC) {

      TVector3 lep_pos_3vec, lep_neg_3vec; 


      //m_tree->GetEntry(master);
      //cout << mc_vx_z[1] << endl;
      const int nmcevt = (mcevt->pileUpType())->size();

      vector<TVector3> mcparticle3vectors_temp; // for primary vertex ... these contain also the tracks close to leptons
      vector<TVector3> mcparticle3vectors; // for primary vertex
      vector< vector<TVector3> > mc3vecs_pileup(nmcevt); // for pile-up vertices
      mcparticle3vectors_temp.clear();
      mcparticle3vectors.clear();
      vector<int> barcodes;

      bool doPU = false;
      //bool doPU = true;

      int n_VB_leptons =0;
      int n_VB_leptons_neg =0;
      int n_VB_leptons_pos =0;
      
      //cout << "Next EVENT, " << evt_weight_gen <<endl;
      std::pair<vector<TLorentzVector>, std::pair<vector<int>, vector<int> > > opfer = getDressedLeptons(*mc, MCGENERATOR, m_muon, true, true);
      std::pair<vector<TLorentzVector>, std::pair<vector<int>, vector<int> > > opfer_nohadrons = getDressedLeptons(*mc, MCGENERATOR, m_muon, false, true);
      std::pair<vector<TLorentzVector>, std::pair<vector<int>, vector<int> > > opfer_bare = getDressedLeptons(*mc, MCGENERATOR, m_muon, true, false);
      vector<TLorentzVector> dressed_leptons = opfer.first ;// getDressedLeptons(*mc, MCGENERATOR, m_muon);
      vector<TLorentzVector> dressed_leptons_nohadrons = opfer_nohadrons.first ;// getDressedLeptons(*mc, MCGENERATOR, m_muon);
      vector<TLorentzVector> dressed_leptons_bare = opfer_bare.first ;// getDressedLeptons(*mc, MCGENERATOR, m_muon);
      vector<int> charges        = opfer.second.first;
      vector<int> lepton_indices = opfer.second.second;
      vector<int> accepted;
      
      // other lepton definitions
      double VB_pt_gen_nohadrons, VB_pt_gen_bare, VB_mit_gen_nohadrons, VB_mit_gen_bare;
              
      // Leptons dressed omitting photns off meson decays
      if (dressed_leptons_nohadrons.size() == 2) {
          TLorentzVector z_4vec_nohadrons = opfer_nohadrons.first[0] + opfer_nohadrons.first[1];
          VB_mit_gen_nohadrons  = z_4vec_nohadrons.M();
          VB_pt_gen_nohadrons   = z_4vec_nohadrons.Pt();
          //if (70. < VB_mit_gen_nohadrons && 112. > VB_mit_gen_nohadrons) {
          if (m_MLL_MIN < VB_mit_gen_nohadrons && m_MLL_MAX > VB_mit_gen_nohadrons) {
              //cout << "No hadrons: " << VB_pt_gen_nohadrons << endl;
              fillh1f(hPtVB+hs_VBGenNoHadrons, VB_pt_gen_nohadrons, evt_weight_gen);
          }
      }
      // Bare leptons aka Born level
      if (dressed_leptons_bare.size() == 2) {
          TLorentzVector z_4vec_bare = opfer_bare.first[0] + opfer_bare.first[1];
          VB_mit_gen_bare  = z_4vec_bare.M();
          VB_pt_gen_bare   = z_4vec_bare.Pt();
          //if (70. < VB_mit_gen_bare && 112. > VB_mit_gen_bare) {
          if (m_MLL_MIN < VB_mit_gen_bare && m_MLL_MAX > VB_mit_gen_bare) {
              //cout << "Bare: " << VB_pt_gen_bare << endl;
              fillh1f(hPtVB+hs_VBGenBare, VB_pt_gen_bare, evt_weight_gen);
          }
      }
      

      // The stuff that is really used as truth definition
      //if (dressed_leptons.size() == 2) {
          //if (charges[0] > 0) {
              //lep_pos_3vec = dressed_leptons[0].Vect();
              //lep_neg_3vec = dressed_leptons[1].Vect();
          //}
          //else {
              //lep_pos_3vec = dressed_leptons[1].Vect();
              //lep_neg_3vec = dressed_leptons[0].Vect();
          //}

          //TLorentzVector z_4vec = dressed_leptons[0] + dressed_leptons[1];
      if (dressed_leptons_nohadrons.size() == 2) {
          if (charges[0] > 0) {
              lep_pos_3vec = dressed_leptons_nohadrons[0].Vect();
              lep_neg_3vec = dressed_leptons_nohadrons[1].Vect();
          }
          else {
              lep_pos_3vec = dressed_leptons_nohadrons[1].Vect();
              lep_neg_3vec = dressed_leptons_nohadrons[0].Vect();
          }

          TLorentzVector z_4vec = dressed_leptons_nohadrons[0] + dressed_leptons_nohadrons[1];
          double VB_mit_gen, VB_pt_gen, VB_rap_gen, VB_phi_gen;
          VB_mit_gen  = z_4vec.M();
          VB_pt_gen   = z_4vec.Pt();
          VB_phi_gen  = z_4vec.Phi();
          VB_rap_gen  = z_4vec.Rapidity();

          if (m_MLL_MIN < VB_mit_gen && m_MLL_MAX > VB_mit_gen) {
          //if (70. < VB_mit_gen && 112. > VB_mit_gen) {
              //cout << "All photons: " << VB_pt_gen << endl;
              //cout << "event accepted" << endl;
     
              // Loop over mc particles to find charged
              for (uint i = 0; i < mc->n(); i++ ) {
                
                // is PU information available?
                //if (i==0) doPU = (mc->mcevt_index(i) != -999); 


                //statuscodes[(int)(*mc).status(i)]++;
                // cut on primary final-state particles
                if (!truthCutPrimary(*mc,i)) continue;

                //cout <<" i: " << i << endl;
                // Do not reuse leptons
                if (i == lepton_indices[0] || i == lepton_indices[1]) {
                    //cout << "Not using ";
                    //printMCinfo(*mc, i);//(*mc).eta.at(i) << endl;// << " " << *mc.pt(i) << " " << *mc.pdgId(i) << endl;
                    continue;
                }

                TVector3 mc3vec = getTrack3Vector(*mc,i);
                
                // Track smearing
                //if (m_systype == kTrkSmear) {
                  //smearTrack(mc3vec); // TODO upsi? 
                //}
                
                //// Tracking efficiency systematics ---randomly remove non-lepton track if
                //// binomial B(1,eff) evaluates to 1 where eff is the eta-pT bin entry in the
                //// tracking efficiency histo
                //if (!isVBLepton(*mc,i, MCGENERATOR)) { // Make sure we don't throw away the VB leptons
                  //if (m_systype == kTrkEffUp) {
                    //if (randomRemove(getTrkEff(mc3vec.Eta(), mc3vec.Pt(), "up"))) continue;
                  //}
                  //else if (m_systype == kTrkEffDown) {
                    //if (randomRemove(getTrkEff(mc3vec.Eta(), mc3vec.Pt(), "down"))) continue;
                  //}
                  //else if (m_systype == kTrkEff) {
                    //if (randomRemove(getTrkEff(mc3vec.Eta(), mc3vec.Pt()))) {
                        ////cout << "Throwing away track" << endl;
                        //continue;
                    //}
                  //}
                //}

                // cut on pile-up particles
                //if (doPU) { // skip if PU info not available (TruthParticle)
                  //if (!truthCutPileUp(*mc,i,mcevt)) {
                    //if (truthCutCharged(*mc,i) && truthCutPhaseSpace(*mc,i)) {
                      //fillTrackHistos(hsGenPartPileUp, mc3vec, evt_weight); 
                      //// int currindex = mc->mcevt_index(i); 
                      //mc3vecs_pileup[mc->mcevt_index(i)].push_back(mc3vec);
                    //}
                    //continue;
                  //}
                //}
                
                // MAL TODO: DRESS LEPTONS !!! Adding four vectors of all photons in cone of DR=0.1 around leptons
                if (truthCutCharged(*mc,i)  && truthCutPhaseSpace(*mc,i)) {
                  mcparticle3vectors_temp.push_back(mc3vec);
                  accepted.push_back(i);
                }
/*
                int pid = (int)(*mc).pdgId(i); 
                if (isChargedLepton(pid, m_muon) && isVBLepton(*mc,i, MCGENERATOR, m_muon)) {
                  //mcparticle3vectors_temp.push_back(mc3vec);
                  if (pid >0) {
                    lep_pos_3vec = getDressedLepton(*mc, i);
                    n_VB_leptons++;
                    n_VB_leptons_pos++;
                    //lep_pos_3vec = mc3vec;
                  }
                  else if (pid < 0) {
                    //	    cout << "saving neg lepton" << endl;
                    lep_neg_3vec = getDressedLepton(*mc, i);
                    n_VB_leptons++;
                    n_VB_leptons_neg++;
                    //lep_neg_3vec = mc3vec;
                  }
                  else cout << "ERROR: VB lepton is neither positively nor negatively charged!" << endl;
                }
                //else if (isNeutrino(pid) && isVBLepton(*mc,i, MCGENERATOR))  nu3vec = mc3vec;
                else if (truthCutCharged(*mc,i) && truthCutPhaseSpace(*mc,i) &! isVBLepton(*mc,i, MCGENERATOR, m_muon) &! (fabs(pid) == 11 || fabs(pid) ==13 )) {
                  mcparticle3vectors_temp.push_back(mc3vec);
                  barcodes.push_back((int)(*mc).barcode(i));
                }
                //else {
                    //continue;
                    ////cout << "shut the fuck up" << endl;
                    ////if(truthCutCharged(*mc,i)) printMCinfo(*mc,i);
                //}
*/
              } // end loop over mc particles
              //if (!(n_VB_leptons ==2)) {
                  //return kTRUE;
              //}
                  //cout << n_VB_leptons << "VB Leptons" << endl;
              // remove tracks close to leptons!!! as in reco!!!
              //for (unsigned int i=0; i<mcparticle3vectors_temp.size(); i++) {
                  //if ( (deltaR(mcparticle3vectors_temp[i], lep_pos_3vec) < CUT_LEPCONE) ||  (deltaR(mcparticle3vectors_temp[i], lep_neg_3vec) < CUT_LEPCONE)   ) {
                      ////cout << "Skip mc track" << endl;
                  //}
                  //else {
                  //fillTrackHistos(hsGenPart, mcparticle3vectors_temp[i], evt_weight);
                      //mcparticle3vectors.push_back(mcparticle3vectors_temp[i]);
                  //}
              //}
              //cout << m_systype << " systemactic" << endl;
              //cout << mcparticle3vectors.size() << " mc tracks" << endl;

              // Calculate truth mT, mee and pT
              // ------------------------------

              //cout << "11: " << lep_pos_3vec.Pt() << endl;
              //cout << "12: " << lep_neg_3vec.Pt() << endl;

        //#ifdef _Z_ANA_
              if (lep_pos_3vec.Pt()<1.e-6 || lep_neg_3vec.Pt()<1.e-6) {
                cout << "ERROR: Did not find at least one of the leptons!" << endl;
              }
              //VB_mit_gen  = calculateMeeZ(lep_pos_3vec,lep_neg_3vec, m_muon);
              //VB_pt_gen   = calculatePtVB(lep_pos_3vec,lep_neg_3vec, m_muon);
              ////cout << VB_pt_gen << " " << VB_mit_gen << " " << m_muon << endl;
              //VB_phi_gen  = calculatePhiZ(lep_pos_3vec,lep_neg_3vec, m_muon);
              //VB_rap_gen  = calculateYZ(lep_pos_3vec,lep_neg_3vec, m_muon);
              //cout << "uff: " << VB_pt_gen << endl;
              //if (fabs(z_4vec.Pt() -VB_pt_gen) > .0001)  cout << "True that - uff: " << z_4vec.Pt() -VB_pt_gen << endl;
        //#else 
              //if (lep_pos_3vec.Mag()==0) lep_pos_3vec = lep_neg_3vec;
              //VB_mit_gen = calculateMtW(lep_pos_3vec,nu3vec);
              //VB_pt_gen  = calculatePtVB(lep_pos_3vec,nu3vec):
        //#endif
         
              //if (truthCutVB(VB_mit_gen, lep_pos_3vec, lep_neg_3vec) && (n_VB_leptons ==2) && (n_VB_leptons_pos ==1) && (n_VB_leptons_neg==1)) {

              double pTl1 = lep_pos_3vec.Pt();
              double pTl2 = lep_neg_3vec.Pt();
              fillh1f(hPtLept1 + hs_VBGenAll, pTl1,  evt_weight_gen);
              fillh1f(hPtLept2 + hs_VBGenAll, pTl2,  evt_weight_gen);
              fillh1f(hPtLept  + hs_VBGenAll, pTl1,  evt_weight_gen);
              fillh1f(hPtLept  + hs_VBGenAll, pTl2,  evt_weight_gen);
              double etal1 = lep_pos_3vec.Eta() ;
              double etal2 = lep_neg_3vec.Eta() ;
              fillh1f(hEtaLept1 + hs_VBGenAll, etal1,  evt_weight_gen);
              fillh1f(hEtaLept2 + hs_VBGenAll, etal2,  evt_weight_gen);
              fillh1f(hEtaLept  + hs_VBGenAll, etal1,  evt_weight_gen);
              fillh1f(hEtaLept  + hs_VBGenAll, etal2,  evt_weight_gen);
              double pxl1 = lep_pos_3vec.Px();
              double pxl2 = lep_neg_3vec.Px();
              fillh1f(hPxLept  + hs_VBGenAll, pxl1,  evt_weight_gen);
              fillh1f(hPxLept  + hs_VBGenAll, pxl2,  evt_weight_gen);
              double pyl1 = lep_pos_3vec.Py();
              double pyl2 = lep_neg_3vec.Py();
              fillh1f(hPyLept  + hs_VBGenAll, pyl1,  evt_weight_gen);
              fillh1f(hPyLept  + hs_VBGenAll, pyl2,  evt_weight_gen);
              double pzl1 = lep_pos_3vec.Pz();
              double pzl2 = lep_neg_3vec.Pz();
              fillh1f(hPzLept  + hs_VBGenAll, pzl1,  evt_weight_gen);
              fillh1f(hPzLept  + hs_VBGenAll, pzl2,  evt_weight_gen);
              fillh1f(hLeptDeltaR  + hs_VBGenAll, lep_pos_3vec.DeltaR(lep_neg_3vec),  evt_weight_gen);
                
              for (unsigned int i=0; i<mcparticle3vectors_temp.size(); i++) {
                  //if ( (deltaR(mcparticle3vectors_temp[i], lep_pos_3vec) < CUT_LEPCONE) ||  (deltaR(mcparticle3vectors_temp[i], lep_neg_3vec) < CUT_LEPCONE)   ) {
                      ////cout << "Skip mc track" << endl;
                      //continue;
                  //}
                  //else {
                      fillTrackHistos(hsGenPart, mcparticle3vectors_temp[i], 1.0);
                      //fillh1f(hBarcode + hsGenPart, barcodes[i], 1.0);
                      mcparticle3vectors.push_back(mcparticle3vectors_temp[i]);
                  //}
              }
              //if (mcparticle3vectors_temp.size() <2) {
                  ////cout << mcparticle3vectors_temp.size() << " Particles" << endl;
                  //for (uint a =0; a< mcparticle3vectors_temp.size(); a++) {
                    ////printMCinfo(*mc,accepted[a]);
                    ////nch2particles[pdgMCinfo(*mc, accepted[a]).first].push_back(pdgMCinfo(*mc, accepted[a]).second);
                    //int pdgid = pdgMCinfo(*mc, accepted[a]).first;
                    //int binid =10;
                    //switch (pdgid) {
                        //case 11: binid=0; break;
                        //case 13: binid=1;break;
                        //case 211: binid=2;break;
                        //case 321: binid=3;break;
                        //case 2212: binid=4;break;
                        //case 3112: binid=5;break;
                        //case 3222: binid=6;break;
                        //case 3312: binid=7;break;
                        //case 3334: binid=8;break;
                    //}
                    //fillh1f(hPDGIDNch2 + hs_VBGenAll, binid,  evt_weight_gen);
                    //// pt histogram for those
                    //double Etanch2 = EtaMCinfo(*mc, accepted[a]);
                    //fillh1f(hEtaPartNch2 + hs_VBGenAll, Etanch2,  evt_weight_gen);
                    //double pTnch2 = pTMCinfo(*mc, accepted[a]);
                    //fillh1f(hPtPartNch2 + hs_VBGenAll, pTnch2,  evt_weight_gen);
                  //}
                  //fillh1f(hPtLept1Nch2 + hs_VBGenAll, pTl1,  evt_weight_gen);
                  //fillh1f(hPtLept2Nch2 + hs_VBGenAll, pTl2,  evt_weight_gen);
                  //fillh1f(hPtLeptNch2  + hs_VBGenAll, pTl1,  evt_weight_gen);
                  //fillh1f(hPtLeptNch2  + hs_VBGenAll, pTl2,  evt_weight_gen);
                  //fillh1f(hEtaLept1Nch2 + hs_VBGenAll, etal1,  evt_weight_gen);
                  //fillh1f(hEtaLept2Nch2 + hs_VBGenAll, etal2,  evt_weight_gen);
                  //fillh1f(hEtaLeptNch2  + hs_VBGenAll, etal1,  evt_weight_gen);
                  //fillh1f(hEtaLeptNch2  + hs_VBGenAll, etal2,  evt_weight_gen);
                  //fillh1f(hPxLeptNch2  + hs_VBGenAll, pxl1,  evt_weight_gen);
                  //fillh1f(hPxLeptNch2  + hs_VBGenAll, pxl2,  evt_weight_gen);
                  //fillh1f(hPyLeptNch2  + hs_VBGenAll, pyl1,  evt_weight_gen);
                  //fillh1f(hPyLeptNch2  + hs_VBGenAll, pyl2,  evt_weight_gen);
                  //fillh1f(hPzLeptNch2  + hs_VBGenAll, pzl1,  evt_weight_gen);
                  //fillh1f(hPzLeptNch2  + hs_VBGenAll, pzl2,  evt_weight_gen);
                  ////Delta R
                  //fillh1f(hLeptDeltaRNch2  + hs_VBGenAll, lep_pos_3vec.DeltaR(lep_neg_3vec),  evt_weight_gen);

              //} // if n particles <2
              
              fillh1f(hNGen + hsGenPart, mcparticle3vectors.size(), evt_weight_gen);

              fillh1f(hMitVB+hs_VBGenAll, VB_mit_gen, evt_weight_gen);
              fillh1f(hPtVB+hs_VBGenAll, VB_pt_gen, evt_weight_gen);

              if (mcparticle3vectors.size() >= 2) {
                fillh1f(hMitVB+hs_VBGenNch2, VB_mit_gen, evt_weight_gen);
                fillh1f(hPtVB+hs_VBGenNch2, VB_pt_gen, evt_weight_gen);
              }

                  //fillMcMitVBHistos(VB_mit_gen, VB_mit_rec);
                  //fillMcPtVBHistos(VB_pt_gen, VB_pt_rec);
                  // MC event shape calculations
                  // ---------------------------------------

                  //int nmc_PUvtx = 0;
                  // First the pile-up vertices
                  //for (int i=0; i < nmcevt; i++) {

                    //short PUT = (mcevt->pileUpType())->at(i);

                    //if (!(PUT==1 || PUT==2)) continue; // keep only in-time and out-of-time pile-up

                    //nmc_PUvtx++;

                    //vector< vector<double> > stl_vec_mc_pu = convertTVector3toDouble(&mc3vecs_pileup[i]);
                    //EventShapes mcES_pileup = EventShapes(stl_vec_mc_pu);
                    //fillEventShapes(hsGenPileUp, &mcES_pileup, VB_pt_gen);

                    //dynPileUpLibrary_MCtruth.push_back(stl_vec_mc_pu);
                  //}

                  //vector< vector<double> > stl_vec_pileUpLibrary_MCtruth;
                  //if ((m_systype == kPileUp) && dynPileUpLibrary_MCtruth.size()>0 && addPileUp)
                    //stl_vec_pileUpLibrary_MCtruth = dynPileUpLibrary_MCtruth.at((unsigned int)gRandom->Uniform(0.,dynPileUpLibrary_MCtruth.size()));

                  // Now the primary vertex
                  vector< vector<double> > stl_vec_mc = convertTVector3toDouble(&mcparticle3vectors);
                  //stl_vec_mc.insert(stl_vec_mc.end(),stl_vec_pileUpLibrary_MCtruth.begin(),stl_vec_pileUpLibrary_MCtruth.end()); //HS commented out 20130716
                  //EventShapes mcES = EventShapes(stl_vec_mc, VB_phi_gen, true);
                  vector<double> VB_vec_mc;
                  VB_vec_mc   = getZThreeVector(lep_pos_3vec, lep_neg_3vec, m_muon);
                  EventShapes mcES = EventShapes(stl_vec_mc, VB_vec_mc, true);
                  //EventShapes mcES = EventShapes(stl_vec_mc, VB_phi_gen, true);
                  //EventShapes trkSelES(stl_vec, VB_vec, true);

                  if (!m_systype) mcES_nominal = mcES; // store event shapes for systematics

                  // Imagiro Trees
                  // ---------------------------------------
                  //if (m_systype == kTrkEff) {
                      ////fillImagiroTree(imagiro_tree_genTrkEff, imagiro_tree_gen2TrkEff, &mcES, VB_pt_gen);
                      //fillImagiroTree(imagiro_tree_genTrkEff, imagiro_tree_gen2TrkEff, &mcES, VB_Eta_gen, VB_pt_gen);
                  //}
                  //else {
                  if (!m_systype || m_systype==kTrkEff) { // HS careful when reediting 2015
                      //if (truthCutVB(VB_mit_gen, lep_pos_3vec, lep_neg_3vec) && (n_VB_leptons ==2)) {
                  //if (!(n_VB_leptons ==2)) {
                      //return kTRUE;
                  //}
                          fillImagiroTree(imagiro_tree_gen, imagiro_tree_gen2, &mcES, VB_rap_gen, VB_pt_gen, evt_weight_gen);
                          fillEventShapes(hsGen, &mcES, VB_pt_gen, evt_weight_gen, 0.0, true);
                          fillEventShapesDump(hsGen, &mcES, VB_pt_gen, stl_vec_mc); 
                          //}
                  }
                      //}
                  
                  // Out-of-Acceptance Correction
                  // ---------------------------------------
                  //fillOACorrection(&trkSelES, &mcES, VB_pt_gen);
            ////////////
                  // Smearing of event shapes (Truth vs. Reco)
                  // -------------------------------------------
                  //fillEventShapeSmear(hs_GenRec, &mcES, &trkSelES, VB_pt_gen, VB_pt_rec); // TODO wieder einbauen
                  //fillEventShapeDiff (hs_GenRec, &mcES, &trkSelES, VB_pt_gen, VB_pt_rec);
                  
                  // Smearing of event shapes due to systematics
                  // -------------------------------------------
                  //fillEventShapeSmear(hs_GenSyst, &mcES_nominal, &mcES, VB_pt_rec, VB_pt_rec);
                  //fillEventShapeDiff (hs_GenSyst, &mcES_nominal, &mcES, VB_pt_rec, VB_pt_rec);
                  
                  // That other thing, dumping  
                  //fillEventShapesDump(hsRecSel, &mcES, &trkSelES, VB_pt_gen); 
                  //fillEventShapeSmear(hsRecSel, &trkSelES, VB_pt_rec); 
            //// HS switched off April 16 2013
                  
                  // NNTree 
                  //if (flag_dumpNN && !m_systype) {
                    //UEShapesNoVtx::fillNNtreeRecoAndGen(&trkSelES, VB_pt_rec, &mcES, VB_pt_gen, nntree, evt_weight);
                  //}

            } // Z cuts // lept0n, VB acceptance
                  //else {
                      //cout << "failed truth cut mvb :" << VB_mit_gen << endl;
                  //}
    } // n dressed leptons == 2
    } //isMC
    //else {
      //// NNTree 
      //if (flag_dumpNN && !m_systype) {
        //UEShapesNoVtx::fillNNtreeReco(&trkSelES, VB_pt_rec, nntree, evt_weight);
      //}
    //}

  } // systematics loop

  //}
  master++;

  
  return kTRUE;
}

//______________________________________________________________________ 

// returns 0 if track was not selected; 1 if track was selected

bool UEShapesNoVtx::doTrackAnalysis(int suffix, int itrk, vector<TVector3>& trkvec, double deltaZ, double z0) {

  // returns 0

  if (!(suffix==hsTrkAll || suffix==hsTrkSel || suffix==hsTrkSelPileUp || hsTrkSideBand )) {
    cout << "ERROR: Unknown histo suffix!" << endl;
    return 0;
  }
  TVector3 trk3vec = getTrack3Vector(*trk, itrk);

  if (suffix==hsTrkAll) {
    trkvec.push_back(trk3vec);
    fillTrackHistos(hsTrkAll, trk3vec, evt_weight_reco); // fill tracks
    if (flag_doTrackControlPlots ) fillTrackD3PDObjectControlPlots(hs_All, itrk, evt_weight_reco);
    return 1;
  }

  // can now safely assume that (suffix==hsTrkSelPileUp || suffix==hsTrkSel)

  // Track selection

  if (!trkCutPhaseSpace(*trk, itrk)) return 0;
  if (!trkCutInnDet(*trk, itrk))     return 0;
  if (!trkCutD0(*trk, itrk))     return 0;

  // remove (dressed) lep1 and lep2 tracks

  // HS 2015 --- not removing leptons to see if there is a bug in the cone algo --- off
  //
  
  if (deltaR(trk3vec,lep1_trk_3vec) < CUT_LEPCONE) {
    if (suffix==hsTrkSel) lep1_cone_sumpt += trk3vec.Pt();
    return 0;
  }
#ifdef _Z_ANA_
  if (deltaR(trk3vec,lep2_trk_3vec) < CUT_LEPCONE) {
    if (suffix==hsTrkSel) lep2_cone_sumpt += trk3vec.Pt();
    return 0;
  }
#endif
 

  fillh1f(hZ0wrtPUV + suffix, (*trk)[itrk].z0_wrtBL(), evt_weight_reco);

  // The new pu extraction plots
  if (fabs(deltaZ) <1.5) fillh1f(hDeltaZwrtPVPeak + suffix, deltaZ, evt_weight_reco);
  fillh1f(hDeltaZwrtPVFull + suffix, deltaZ, evt_weight_reco);
  fillh2f(hDeltaZwrtPVFullZ0 + suffix, deltaZ, z0, evt_weight_reco);

  //if (suffix==hsTrkSel && (!trkCutVertex(*trk, itrk))) return 0;
  //else if (suffix==hsTrkSelPileUp && (!trkCutPUVertex(*trk, itrk, deltaZ))) return 0;
  if (suffix==hsTrkSelPileUp && (!trkCutPUVertex(*trk, itrk, deltaZ))) return 0;

  trkvec.push_back(trk3vec);
  fillTrackHistos(suffix,  trk3vec, evt_weight_reco);

  if (suffix==hsTrkSel) {
    double trk_weight = getTrackWeight(trk3vec);
    trackweights.push_back(trk_weight);

    fillTrackHistos(hsTrkCorr, trk3vec, trk_weight*evt_weight_reco);

    if (flag_doTrackControlPlots) {
      fillTrackD3PDObjectControlPlots(hsTrkSel,  itrk, evt_weight_reco);
      fillTrackD3PDObjectControlPlots(hsTrkCorr, itrk, trk_weight*evt_weight_reco);
    }

  }

  return 1;

}

bool UEShapesNoVtx::doTrackAnalysisLegacy(int suffix, int itrk, vector<TVector3>& trkvec) {

  // returns 0

  if (!(suffix==hsTrkAll || suffix==hsTrkSel || suffix==hsTrkSelPileUp || hsTrkSideBand )) {
    cout << "ERROR: Unknown histo suffix!" << endl;
    return 0;
  }
  TVector3 trk3vec = getTrack3Vector(*trk, itrk);

  if (suffix==hsTrkAll) {
    trkvec.push_back(trk3vec);
    fillTrackHistos(hsTrkAll, trk3vec, evt_weight_reco); // fill tracks
    if (flag_doTrackControlPlots ) fillTrackD3PDObjectControlPlots(hs_All, itrk, evt_weight_reco);
    return 1;
  }
  // remove (dressed) lep1 and lep2 tracks
  if (deltaR(trk3vec,lep1_trk_3vec) < CUT_LEPCONE) {
    if (suffix==hsTrkSel) {
      lep1_cone_sumpt += trk3vec.Pt();
      lep1_cone_n++;
      if ( (deltaR(trk3vec,lep1_trk_3vec) <0.01) && trk->pt()->at(itrk) > 20000.) {
        //cout << trk->z0_wrtPV()->at(itrk) << endl;
        fillh1f(hLeptonDZ + hsRecSel, trk->z0_wrtPV()->at(itrk), evt_weight_reco);
      }
    }
    return 0;
  }
#ifdef _Z_ANA_
  if (deltaR(trk3vec,lep2_trk_3vec) < CUT_LEPCONE) {
    if (suffix==hsTrkSel) {
      lep2_cone_sumpt += trk3vec.Pt();
      lep2_cone_n++;
      if ( (deltaR(trk3vec,lep1_trk_3vec) <0.01) && trk->pt()->at(itrk) > 20000.) {
        //cout << trk->z0_wrtPV()->at(itrk) << endl;
        fillh1f(hLeptonDZ + hsRecSel, trk->z0_wrtPV()->at(itrk), evt_weight_reco);
      }
    }
    return 0;
  }
#endif

  // can now safely assume that (suffix==hsTrkSelPileUp || suffix==hsTrkSel)

  // Track selection

  if (!trkCutPhaseSpace(*trk, itrk)) return 0;
  if (!trkCutInnDet(*trk, itrk))     return 0;
  if (!trkCutVertex(*trk, itrk))     return 0;


  //fillh1f(hZ0wrtPUV + suffix, (*trk)[itrk].z0_wrtBL(), evt_weight_reco);

  // The new pu extraction plots
  //if (fabs(deltaZ) <1.5) fillh1f(hDeltaZwrtPVPeak + suffix, deltaZ, evt_weight_reco);
  //fillh1f(hDeltaZwrtPVFull + suffix, deltaZ, evt_weight_reco);
  //fillh2f(hDeltaZwrtPVFullZ0 + suffix, deltaZ, z0, evt_weight_reco);

  //if (suffix==hsTrkSelPileUp && (!trkCutPUVertex(*trk, itrk, deltaZ))) return 0;

  trkvec.push_back(trk3vec);
  //fillTrackHistos(suffix,  trk3vec, evt_weight_reco);

  if (suffix==hsTrkSel) {
    double trk_weight = getTrackWeight(trk3vec);
    trackweights.push_back(trk_weight);

    fillTrackHistos(hsTrkCorr, trk3vec, trk_weight*evt_weight_reco);

    if (flag_doTrackControlPlots) {
      fillTrackD3PDObjectControlPlots(hsTrkSel,  itrk, evt_weight_reco);
      fillTrackD3PDObjectControlPlots(hsTrkCorr, itrk, trk_weight*evt_weight_reco);
    }

  }

  return 1;

}
//______________________________________________________________________ 
void UEShapesNoVtx::bumpIdentityInMigrationMatrix(double percentage_merged, TH2F& inputHisto){
    /* assuming percentage_merged is the number we are allowed to merge
    but we merged all events.
    the rest (meaning 1-percentage_merged) is now filled as identity in the histo. */
    double processed = 0;
    double factor = (1./percentage_merged-1.);
    double binCenter = 0;
    for(int binNr = 1; binNr <= inputHisto.GetNbinsX(); binNr++){
        processed = inputHisto.Integral(binNr,binNr,0,-1);
        binCenter = inputHisto.GetBinCenter(binNr);
        inputHisto.Fill(binCenter, binCenter, processed*factor);
    }
}


//______________________________________________________________________ 

 void UEShapesNoVtx::Terminate()
{

  
  if (m_dumpReco) {
      dumpStuff();
  }
  map<int, vector<int> >::const_iterator itr_nch2;
  cout << "\n" << endl;
  //for(itr_nch2 = nch2particles.begin(); itr_nch2 != nch2particles.end(); ++itr_nch2){
      //cout << "PDGID " << (*itr_nch2).first << "  Found " << (*itr_nch2).second.size() << endl;
    ////for (uint l=0;l < (*itr_nch2).second.size();l++) {
        ////cout << "\tBARCODE " << (*itr_nch2).second[l] << endl;
      ////}
  //}
  //map<int, int >::const_iterator itr_statuscodes;
  //cout << "\n" << endl;
  //for(itr_statuscodes = statuscodes.begin(); itr_statuscodes != statuscodes.end(); ++itr_statuscodes){
      //cout << "SC " << (*itr_statuscodes).first << "  Found " << (*itr_statuscodes).second << endl;
    ////for (uint l=0;l < (*itr_statuscodes).second.size();l++) {
        ////cout << "\tBARCODE " << (*itr_statuscodes).second[l] << endl;
      ////}
  //}
  

  if (m_writeTrackLib) {
    cout << "Found a total of " << pu_tracks << " tracks associated with isolated type3 vertices" << endl;  
    
    if (! exists(m_tracklibpath.Data())) {
        cout << "Creating new directory for track library: " << m_tracklibpath << endl;
        mkdir(m_tracklibpath.Data(), S_IRWXU | S_IRWXG | S_IRWXO);
    }
    
    map<int, vector<vector<vector<double> > > >::const_iterator itr;
    vector<vector<vector<double> > >::const_iterator itr_vec;
    
    int n_vertices = 0;
    int n_tracks = 0;
    for(itr = pumap.begin(); itr != pumap.end(); ++itr){
      cout << "Vertices with " << (*itr).first << " tracks: Found " << (*itr).second.size() << endl;
      //for(itr_vec = pumap[(*itr).first].begin(); itr_vec != pumap[(*itr).first].end(); ++itr_vec){
          ////cout << (*itr).first << " vs. " << (*itr_vec).size() << endl;
          //assert((*itr).first == (*itr_vec).size());
      //}
      n_vertices += (*itr).second.size();
      string s = lexical_cast<string>(m_tracklibpath) + lexical_cast<string>("/tracklib_n_");
      s+= lexical_cast<string>((*itr).first);
      //cout << "Wrote tracklib to " << s << endl;
      std::ofstream ofs(s.c_str());
      boost::archive::binary_oarchive oa(ofs);
      oa << (*itr).second;
      for (int i=0;i <(*itr).second.size();i++) n_tracks+=(*itr).second[i].size();
    }
    cout << "Found " << n_tracks << " tracks in " << n_vertices << " vertices" << endl;

    // Write some histograms with information on tracks in tracklib
    TString histo_out = m_tracklibpath + TString("_pulibstats") + TString(".root");
    TFile* fawesome = TFile::Open(histo_out, "recreate");
    fawesome->cd();
    TH1F* pu_n = new TH1F("h_ntrk_pu", "Ntracks per vertex", 130, 0, 130);
    TH2F* pu_eta = new TH2F("h_ntrk_eta", "eta", 150, -0.5, 149.5, 21, -2.5, 2.5);
    TH2F* pu_phi = new TH2F("h_ntrk_phi", "phi", 150, -0.5, 149.5, 21, -0.5, 3.5);
    TH2F* pu_pT = new TH2F("h_ntrk_pT", "pT", 150, -0.5, 149.5, 81, 0., 10.);
    // Loop over vertex classes
    for(itr = pumap.begin(); itr != pumap.end(); ++itr){
      // Loop over vertices
      for (int j=0; j<pumap[(*itr).first].size(); j++) {
        int ntrkatv =  (*itr).first;
        pu_n->Fill(ntrkatv);
        // Loop over tracks at that vertex
        for (int k=0;k<pumap[ntrkatv][j].size(); k++) {
          TVector3 vec = convertDoubletoTVector3(pumap[ntrkatv][j][k]);
          pu_eta->Fill(ntrkatv, vec.Eta());
          pu_phi->Fill(ntrkatv, vec.Phi());
          pu_pT->Fill(ntrkatv, vec.Pt());
        }
 
      }
    }
    pu_n->SetDirectory(fawesome);
    pu_n->Write();
    pu_eta->SetDirectory(fawesome);
    pu_eta->Write();
    pu_phi->SetDirectory(fawesome);
    pu_phi->Write();
    pu_pT->SetDirectory(fawesome);
    pu_pT->Write();
    m_pu_z->SetDirectory(fawesome);
    m_pu_z->Write();  
    m_pu_zN->SetDirectory(fawesome);
    m_pu_zN->Write();  
    m_pu_vtxN->SetDirectory(fawesome);
    m_pu_vtxN->Write();  
    //m_pu_vtxDz->SetDirectory(fawesome);
    //m_pu_vtxDz->Write();  
    m_pu_vtxzPP->SetDirectory(fawesome);
    m_pu_vtxzPP->Write();  
    m_pu_Dz_z->SetDirectory(fawesome);
    m_pu_Dz_z->Write();
    m_pu_Dz_N->SetDirectory(fawesome);
    m_pu_Dz_N->Write();
    m_pu_d_z->SetDirectory(fawesome);
    m_pu_d_z->Write();
    fawesome->Close();

    cout << "Wrote track library to " << m_tracklibpath << endl;
    cout << "Wrote track library summary to " << histo_out << endl;
  }




  cout << "Total tracks is signal region: " << ntrk_sig << endl;
  cout << "Total tracks added to signal region: " << ntrk_add << endl;
  cout << "Fraction of added tracks: " << 100.*ntrk_add/ntrk_sig << endl;
  
  if (m_writeSideBands) {
    TFile* fsb = TFile::Open(m_pileuplibpath, "recreate");
    harray[0][hNtrk+hsRecSideBand]->SetDirectory(fsb);
    harray[0][hNtrkZ0+hsRecSideBand]->SetDirectory(fsb);

    harray[0][15020]->SetDirectory(fsb);  // pt
    harray[0][15023]->SetDirectory(fsb); // eta
    harray[0][15025]->SetDirectory(fsb); //phi
    harray[0][15029]->SetDirectory(fsb); // dz vs z
      //hNtrkRecSideBand->SetDirectory(fsb);
//    hPtTrkSideBand->SetDirectory(fsb);
 //   hPhiTrkSideBand->SetDirectory(fsb);
   // hEtaTrkSideBand->SetDirectory(fsb);
    fsb->Write();
    cout << endl << "Wrote side-band plots to " << m_pileuplibpath << endl;
    //cout << "Exiting" << endl;
    //exit(1);
  }


  m_outputFile->cd();

  m_outputFile->WriteObject(&pumap, "PUMAP");

  //imagiro_tree_rec        ->Write("imagiroRec",         TObject::kOverwrite);
  //imagiro_tree_rec2       ->Write("imagiroRec2",        TObject::kOverwrite);
  //imagiro_tree_recTrkEff  ->Write("imagiroRecTrkEff",   TObject::kOverwrite);
  //imagiro_tree_rec2TrkEff ->Write("imagiroRec2TrkEff",  TObject::kOverwrite);
  //imagiro_tree_recTrkSmearUp  ->Write("imagiroRecTrkSmearUp",   TObject::kOverwrite);
  //imagiro_tree_rec2TrkSmearUp ->Write("imagiroRec2TrkSmearUp",  TObject::kOverwrite);
  //imagiro_tree_recTrkSmearDown  ->Write("imagiroRecTrkSmearDown",   TObject::kOverwrite);
  //imagiro_tree_rec2TrkSmearDown ->Write("imagiroRec2TrkSmearDown",  TObject::kOverwrite);
  //imagiro_tree_recPtZUp   ->Write("imagiroRecPtZUp",    TObject::kOverwrite);
  //imagiro_tree_rec2PtZUp  ->Write("imagiroRec2PtZUp",   TObject::kOverwrite);
  //imagiro_tree_recPtZDown ->Write("imagiroRecPtZDown",  TObject::kOverwrite);
  //imagiro_tree_rec2PtZDown->Write("imagiroRec2PtZDown", TObject::kOverwrite);

  if (flag_isMC) {
   // imagiro_tree_gen ->Write("imagiroGen",  TObject::kOverwrite);
   // imagiro_tree_gen2->Write("imagiroGen2", TObject::kOverwrite);
    //imagiro_tree_genTrkEff ->Write("imagiroGenTrkEff",  TObject::kOverwrite);
    //imagiro_tree_gen2TrkEff->Write("imagiroGen2TrkEff", TObject::kOverwrite);
  }
  
  if (flag_DEBUG) l.debug("<UEShapesNoVtx::Terminate()>");
  cout << "" << endl;
  cout << "Events Processed:         " << nevt << endl;
  cout << "No. Weighted Events:      " << nevt_weighted << endl;


  if (m_sysNames) delete m_sysNames;

}

//______________________________________________________________________ 

double UEShapesNoVtx::getEventWeight(bool gen) {

  double w_evt;  
  //double w_evt = VB->evt_weight(); 
  if (gen)  w_evt = VB->weight_gen(); // Correcte
  else      w_evt = VB->evt_weight(); // Correcte
  //double w_evt =1.0; 
  //cout << w_evt << endl;

  // MAL: temporary hack!
  //if (flag_isMC) return w_evt;
  //else return 1.0;

  return w_evt;
}

//______________________________________________________________________ 

double UEShapesNoVtx::getTrackWeight(TVector3 trkvec) {

  return 1.0;

  /*
  double pt   = trkvec.Pt();
  double eta  = trkvec.Eta();
  double aeta = fabs(eta);

  double ptcorrmax = 50.; // for pT > ptcorrmax, apply correction from bin (ptcorrmax-1.)
  
  // TRACK EFFICIENCY /////////////////////////////////////////

  double etrk = -1;
  double err  = 0.;

  if ( aeta < CUT_ETA && pt > CUT_PT) {

    Int_t etrkbin;

    if (pt < ptcorrmax) etrkbin = mTrkEff_pt_eta->FindBin(pt,eta);
    else etrkbin = mTrkEff_pt_eta->FindBin(ptcorrmax-1.,eta);

    etrk = mTrkEff_pt_eta->GetBinContent(etrkbin);

    // SYSTEMATICS

    // MC STATISTICS
    if (m_systype == kStatUp   ) err =  mTrkEff_pt_eta->GetBinError(etrkbin); 
    if (m_systype == kStatDown ) err = -mTrkEff_pt_eta->GetBinError(etrkbin);

    // MATERIAL 
    // (based on SCT extension efficiency + 10% sample)
    double deta=0.;
    if (m_systype == kMatUp || m_systype == kMatDown) {

      int ieta = -1;
      if(aeta<1.3) ieta = 0;
      else if(aeta < 1.9) ieta =1;
      else if(aeta < 2.1) ieta =2;
      else if(aeta < 2.3) ieta =3;
      else if(aeta < 2.5) ieta =4;
      //else printf("WARNING could not find eta bin for material correction\n");
      l.warning("Could not find eta bin for material correcion");

      // systematics from 10/6/2010
      // material systematics
      double static syst[5] = {0.02 , 0.03, 0.04, 0.04 , 0.07  }; // > 500 MeV 
      // track selection and pt turn-on
      double static systTrk = 0.01;

      deta = sqrt( pow(syst[ieta], 2)  + pow(systTrk, 2) );

      if (m_systype == kMatUp   ) err = +deta*etrk; 
      if (m_systype == kMatDown ) err = -deta*etrk; 
    }


    // FLAT SYSTEMATICS
    if (m_systype == kFlatUp || m_systype == kFlatDown) {
      double flatErr = sqrt( 
			   pow(0.002,2)   // Vertex+Trigger  /// conservative estimate 
			   +pow(0.002 ,2));  // Particle Composition

      if(m_systype == kFlatUp)   err =  flatErr;
      if(m_systype == kFlatDown) err = -flatErr;
    }
  }

  // PT RESOLUTION LOSS ////////////////////////////////////////////////////
  double ptloss=0;   
  if ( aeta < CUT_ETA ) {

    int ptlbin = mOutside_pt_eta->FindBin(pt,eta);
    if (pt > CUT_PT && pt < ptcorrmax) ptloss = mOutside_pt_eta->GetBinContent(ptlbin);
  }

  // CONTAMINATIONS (SECONDARY TRACKS) ////////////////////////////

  double secScaleFactor = 0.;
  double sec    = 0.;
  double secErr = 0.;

  if ( fabs(eta) < CUT_ETA ) {

    Int_t cbin = mSecondaries_pt->FindBin(pt);

    if ( flag_isMC) secScaleFactor = 1.; // no scaling factor in MC
    else secScaleFactor = 1.21;

    if ( pt > CUT_PT && pt < 20. ) sec  = secScaleFactor * mSecondaries_pt->GetBinContent(cbin);
    else sec = secScaleFactor * mSecondaries_pt->GetBinContent(mSecondaries_pt->FindBin(19.)); 

    secErr = sec; 
    if (m_systype == kSecUp   ) { secErr = sec*1.25; } 
    if (m_systype == kSecDown ) { secErr = sec*0.75; } 
  }

  // add errors
  etrk += err;
  double deltaSec = secErr - sec;
  sec = secErr;
  double w_trk_noevent = (1.-ptloss) * (1.-sec) / etrk;

  if( deltaSec != 0 && err != 0 ) l.error("Error: Trying to vary two systematics at the same time");

  return w_trk_noevent;
  */
}

//______________________________________________________________________ 

void UEShapesNoVtx::printEvent(double mee, double pTZ, double PhiZ, vector<int> trkindices) {
    ofstream m_printFile;
    bool doMuon=true;
    m_printFile.open("awesome.txt", ofstream::app);
    m_printFile << TString::Format("%i\t%i\t%f\t%f\t%f\n",VB->run(), VB->event(), vtx->z()->at(0), vtx->x()->at(0), vtx->y()->at(0));
    m_printFile << TString::Format("\t%f\t%f\t%f\n#\n", mee, pTZ, PhiZ);
    m_printFile << TString::Format("\t%f\t%f\t%f\n", VB->lep1_cl_pt(), VB->lep1_cl_eta(), VB->lep1_cl_phi());
    m_printFile << TString::Format("\t%f\t%f\t%f\n", VB->lep2_cl_pt(), VB->lep2_cl_eta(), VB->lep2_cl_phi());
    m_printFile << TString::Format("\t%f\t%f\t%f\n", VB->lep1_trk_pt(), VB->lep1_trk_eta(), VB->lep1_trk_phi());
    m_printFile << TString::Format("\t%f\t%f\t%f\n", VB->lep2_trk_pt(), VB->lep2_trk_eta(), VB->lep2_trk_phi());
    if (doMuon) {
        m_printFile << TString::Format("\t%f\t%f\n"   , VB->lep1_mu_staco_nBLHits(), VB->lep1_mu_staco_ptcone20());
        m_printFile << TString::Format("\t%f\t%f\n#\n", VB->lep2_mu_staco_nBLHits(), VB->lep2_mu_staco_ptcone20());
    }
    else{
        m_printFile << TString::Format("\t%i\t%i\t%i\n", VB->lep1_loosePP(), VB->lep1_mediumPP(), VB->lep1_tightPP());
        m_printFile << TString::Format("\t%i\t%i\t%i\n#\n", VB->lep2_loosePP(), VB->lep2_mediumPP(), VB->lep2_tightPP());
    }
    m_printFile << TString::Format("\t%i tracks\n", trkindices.size());
    vector<int> transTracks;
    vector<int> otherTracks;
    for (unsigned int i=0;i<trkindices.size();i++) {
        double phiTrack = trk->phi_wrtPV()->at(trkindices[i]);
      double dPhi = fabs(PhiZ - phiTrack);
      if (dPhi > M_PI) {
          dPhi*= -1.;
          dPhi+=2.*M_PI;
      }
      if (dPhi < 2.*M_PI/3.0 && dPhi > M_PI/3.) transTracks.push_back(trkindices[i]);
      else otherTracks.push_back(trkindices[i]);
    }
    m_printFile << TString::Format("\t%i transverse\n", transTracks.size());
    for (unsigned int i=0;i<transTracks.size();i++) {
        int idx=transTracks[i];
        m_printFile << TString::Format("**\t%f\t%f\t%f\t%f\t%f\n", trk->d0_wrtPV()->at(idx),
                trk->z0_wrtPV()->at(idx)*sin(trk->theta()->at(idx)),
                trk->pt()->at(idx)/1.e3, trk->eta()->at(idx), trk->phi()->at(idx)
                );
        int nbl = (trk->expectBLayerHit()->at(idx)>=1) ? trk->nBLHits()->at(idx) : -42;
        m_printFile << TString::Format("\t%i\t%i\t%i\n", trk->nPixHits()->at(idx), nbl, trk->nSCTHits()->at(idx));

  //if (!(trk[index].nPixHits() >= 1)) return 0;
  //if ( (trk[index].nBLHits() < 1) && (trk[index].expectBLayerHit() >= 1)) return 0;
  //if (!(trk[index].nSCTHits() >= 6)) return 0;
    }
    for (unsigned int i=0;i<otherTracks.size();i++) {
        int idx=otherTracks[i];
        m_printFile << TString::Format("*\t%f\t%f\t%f\t%f\t%f\n",  trk->d0_wrtPV()->at(idx),
                trk->z0_wrtPV()->at(idx)*sin(trk->theta()->at(idx)),
                trk->pt()->at(idx)/1.e3, trk->eta()->at(idx), trk->phi()->at(idx)
                );
        int nbl = (trk->expectBLayerHit()->at(idx)>=1) ? trk->nBLHits()->at(idx) : -42;
        m_printFile << TString::Format("\t%i\t%i\t%i\n", trk->nPixHits()->at(idx), nbl, trk->nSCTHits()->at(idx));
    }
    m_printFile.close();
}

void UEShapesNoVtx::printEventProcessInfo() {

    double elapsed, togo;
    if ((int)nevt%m_nprint==0) {  
      elapsed = m_timer.RealTime();
      togo = elapsed * (double(evtmax)/double(nevt) -1.);
      l.info(TString::Format("Process(): %.0f/%.0f %.2f %% processed. Elapsed: %is, (Done in: %is)", double(nevt), double(evtmax), (100.*nevt)/evtmax, int(elapsed), int(togo)), true);
      m_timer.Continue();
    
      if (flag_DEBUG) {
          ProcInfo_t p;
          gSystem->GetProcInfo(&p);
          l.debug(TString::Format("Processing event %10d Rmem=%8.3fMb Vmem=%8.f3Mb Read %8.3fMb (Calls %6d) CPU %6.1fs SYS %6.1fs TOT %6.1fs %f Hz",
                      (int)nevt, p.fMemResident*1e-3, p.fMemVirtual*1e-3,TFile::GetFileBytesRead()*1.e-6,TFile::GetFileReadCalls(), p.fCpuUser,p.fCpuSys,m_timer.RealTime(),(nevt/m_timer.RealTime())));
          m_timer.Continue();
          fflush( stdout );
      }
    }
}


// HISTOGRAM FILLING METHODS
//______________________________________________________________________ 

void UEShapesNoVtx::fillTrackHistos(int suffix, TVector3 t, double w_trk) {

  double pT  = t.Pt(); 
  double phi = t.Phi();
  double eta = t.Eta();
 
  // Translate pT and phi into px, py pairs
  double px  = cos(phi) * pT;
  double py  = sin(phi) * pT;

  fillh1f(hPt    + suffix,  pT,   w_trk);
  fillh1f(hPx    + suffix,  px,   w_trk);
  fillh1f(hPy    + suffix,  py,   w_trk);
  fillh1f(hEta   + suffix,  eta,  w_trk);
  fillh1f(hPhi   + suffix,  phi,  w_trk);
  
  fillh2f(hEtaPt + suffix,  eta, pT, w_trk);

}

//void UEShapesNoVtx::fillMcMitVBHistos(double genMit, double recMit) {

  //int suffix = hs_VB;

  //// fill resolution histos
  //double VB_mit_res = (recMit - genMit)/genMit;
  //fillh1f(hMitRes+suffix,                 VB_mit_res, evt_weight);
  //fillh2f(hMitGenMitRes+suffix,  genMit,  VB_mit_res, evt_weight);
  //fillh2f(hMitRecMitRes+suffix,  recMit,  VB_mit_res, evt_weight);

  //// fill migration matrices
  //fillh2f(hMitRecMitGen+suffix,  recMit,  genMit,    evt_weight);
//}

//void UEShapesNoVtx::fillMcPtVBHistos(double genPt, double recPt) {

  //int suffix = hs_VB;

  //// fill resolution histos
  //double VB_pt_diff = recPt - genPt;
  //double VB_pt_res  = VB_pt_diff/genPt;
  //fillh1f(hPtRes+suffix,               VB_pt_res,   evt_weight);
  //fillh2f(hPtGenPtRes+suffix,  genPt,  VB_pt_res,   evt_weight);
  //fillh2f(hPtRecPtRes+suffix,  recPt,  VB_pt_res,   evt_weight);
  //fillh2f(hPtGenPtDiff+suffix, genPt,  VB_pt_diff,  evt_weight);
  //fillh2f(hPtRecPtDiff+suffix, recPt,  VB_pt_diff,  evt_weight); 

  //// fill migration matrices
  //fillh2f(hPtRecPtGen+suffix,  recPt,  genPt,      evt_weight);
//}

//___________________________________________________________________ 
void UEShapesNoVtx::fillNNtreeReco(EventShapes* ES, double pt_VB, NNtree* nntree, double weight) {
  //nntree->VB_run   = (UInt_t)VB->event();
  nntree->VB_evt_weight = weight;;

  if (!(ES->hasFParam() && ES->hasSpheroThrust())) return;
  nntree->VB_ptReco      = pt_VB;
  nntree->ThrustReco     = (double)ES->ThrustV();
  nntree->ThrustPhiReco  = (double)ES->ThrustPhiV();
  nntree->MinorReco      = (double)ES->MinorV();
  nntree->FParameterReco = (double)ES->FParameterV();
  nntree->SpherocityReco = (double)ES->SpherocityV();
  nntree->BeamThrustReco = (double)ES->BeamThrustV();

  nntree->VB_ptGen      = 0.0;
  nntree->ThrustGen     = 0.0;
  nntree->ThrustPhiGen  = 0.0;
  nntree->MinorGen      = 0.0;
  nntree->NtrkGen       = 0.0;
  nntree->SumPtGen      = 0.0;
  nntree->AvgPtGen      = 0.0;
  nntree->FParameterGen = 0.0;
  nntree->SpherocityGen = 0.0;
  nntree->BeamThrustGen = 0.0;

  //
  nntree->Fill();
}

//___________________________________________________________________ 
void UEShapesNoVtx::fillNNtreeRecoAndGen(EventShapes* ES, double pt_VB, EventShapes* ESgen, double pt_VBgen, NNtree* nntree, double weight) {

  if (!(ES->hasFParam()    && ES->hasSpheroThrust())) return;
  if (!(ESgen->hasFParam() && ESgen->hasSpheroThrust())) return;
  
  nntree->VB_evt_weight = weight;;

  nntree->VB_ptReco      = pt_VB;
  nntree->ThrustReco     = (double)ES->ThrustV();
  nntree->ThrustPhiReco  = (double)ES->ThrustPhiV();
  nntree->MinorReco      = (double)ES->MinorV();
  nntree->NtrkReco       = (double)ES->Ntrk() ;
  nntree->SumPtReco      = (double)ES->SumPt();
  nntree->AvgPtReco      = (double)ES->AvgPt();
  nntree->FParameterReco = (double)ES->FParameterV();
  nntree->SpherocityReco = (double)ES->SpherocityV();
  nntree->BeamThrustReco = (double)ES->BeamThrustV();
  
  nntree->VB_ptGen      = pt_VBgen;
  nntree->ThrustGen     = (double)ESgen->ThrustV();
  nntree->ThrustPhiGen  = (double)ESgen->ThrustPhiV();
  nntree->MinorGen      = (double)ESgen->MinorV();
  nntree->NtrkGen       = (double)ESgen->Ntrk() ;
  nntree->SumPtGen      = (double)ESgen->SumPt();
  nntree->AvgPtGen      = (double)ESgen->AvgPt();
  nntree->FParameterGen = (double)ESgen->FParameterV();
  nntree->SpherocityGen = (double)ESgen->SpherocityV();
  nntree->BeamThrustGen = (double)ESgen->BeamThrustV();

  nntree->Fill();
}

void UEShapesNoVtx::fillImagiroTree(TTree* imtree, TTree* imtree2, EventShapes* ES, double pt_Eta,
			       double pt_VB, double evt_weight, pair<double,double> pt_VB_updown) {


  //if (m_systype) return; // HS this was in as of April 21st 2013 
  if (!(m_systype == 0 || m_systype == kTrkEff || m_systype == kTrkSmearUp || m_systype == kTrkSmearDown)) return; 

  // Prevent multiple filling of reco Imagiro tree
  if (m_systype == kTrkEff && TString(imtree->GetName()) == TString("imagiroRec")) {
      return;
  }
  if (m_systype == kTrkSmearUp && TString(imtree->GetName()) == TString("imagiroRec")) {
      return;
  }
  if (m_systype == kTrkSmearDown && TString(imtree->GetName()) == TString("imagiroRec")) {
      return;
  }
  /*if (m_systype == kTrkEff && TString(imtree->GetName()) == TString("imagiroGen")) { // HS careful when reediting 2015
    cout << "No way" << endl;
      return;
  }
  */
  if (m_systype == kTrkSmearUp && TString(imtree->GetName()) == TString("imagiroGen")) {
      return;
  }
  if (m_systype == kTrkSmearDown && TString(imtree->GetName()) == TString("imagiroGen")) {
      return;
  }
  
//  cout << "Fill method 1: " << m_systype << imtree->GetName()<<endl;
//  cout << "Fill method 1: " << m_systype << imtree2->GetName()<<endl;

  // set event-level imagiro variables
  
  // all events
  EventNumber = (UInt_t)VB->event();
  EventWeight = (double)evt_weight;
  Ntrk        = (double)ES->Ntrk();
  SumPt       = (double)ES->SumPt();
  SumPtTrans  = (double)ES->SumPtTrans();
  AvgPt       = (double)ES->AvgPt();
  PtVB        = (double)pt_VB;
  EtaVB       = fabs((double)pt_Eta);
  if (pt_VB_updown.first > 0. || pt_VB_updown.second > 0.) {
    PtVB_up     = (double)pt_VB_updown.first; 
    PtVB_down   = (double)pt_VB_updown.second;
  }

  imtree->Fill();

  if (!(ES->hasFParam() && ES->hasSpheroThrust())) return;

  // events with nch or nsel >= 2
  FParameter  = (double)ES->FParameterV();
  Thrust      = (double)ES->ThrustV();
  ThrustPhi   = (double)ES->ThrustPhiV();
  Minor       = (double)ES->MinorV();
  Spherocity  = (double)ES->SpherocityV();
  BeamThrust  = (double)ES->BeamThrustV();

  imtree2->Fill();

}

// not used anymore
/*void UEShapesNoVtx::fillImagiroTree(TTree* imtree, TTree* imtree2, EventShapes* ES, 
			       double pt_VB, double evt_weight, pair<double,double> pt_VB_updown) {

  //if (m_systype) return; // HS this was in as of April 21st 2013 

  // set event-level imagiro variables
  cout << "Fill method 2: " << m_systype << imtree->GetName()<<endl;
  cout << "Fill method 2: " << m_systype << imtree2->GetName()<<endl;
  
  // all events
  EventNumber = (UInt_t)VB->event();
  EventWeight = (double)evt_weight;
  Ntrk        = (double)ES->Ntrk();
  SumPt       = (double)ES->SumPt();
  AvgPt       = (double)ES->AvgPt();
  PtVB        = (double)pt_VB;
  if (pt_VB_updown.first > 0. || pt_VB_updown.second > 0.) {
    PtVB_up     = (double)pt_VB_updown.first; 
    PtVB_down   = (double)pt_VB_updown.second;
  }

  imtree->Fill();

  if (!(ES->hasFParam() && ES->hasSpheroThrust())) return;

  // events with nch or nsel >= 2
  FParameter  = (double)ES->FParameterV();
  Thrust      = (double)ES->ThrustV();
  ThrustPhi   = (double)ES->ThrustPhiV();
  Minor       = (double)ES->MinorV();
  Spherocity  = (double)ES->SpherocityV();
  BeamThrust  = (double)ES->BeamThrustV();

  imtree2->Fill();

}
*/

void UEShapesNoVtx::fillEventShapes(int suffix, EventShapes* ES, double pt_VB, double weight, double z0, bool doSumPtTrans) {

  int suffix2d = ERRORVAL;  
  
  if      (suffix==hsGen)          suffix2d = hsGen2d;
  else if (suffix==hsGenPileUp)    suffix2d = hsGenPileUp2d;
  else if (suffix==hsRecAll)       suffix2d = hsRecAll2d;
  else if (suffix==hsRecSel)       suffix2d = hsRecSel2d;
  else if (suffix==hsRecCorr)      suffix2d = hsRecCorr2d;
  else if (suffix==hsRecSelPileUp) suffix2d = hsRecSelPileUp2d;
  else if (suffix==hsRecSideBand) suffix2d = hsRecSideBand2d;
  else return;

  //if ( ES->hasFParam() ) {
  if (ES->Ntrk()>=0) {
      fillh1f(hNtrk+suffix,                ES->Ntrk(),          weight);
      fillh2f(hNtrkZ0+suffix,              ES->Ntrk(),      z0, weight);
      fillh2f(hNtrkPtVB+suffix2d,          ES->Ntrk(),   pt_VB, weight);
      fillh2f(hNtrkZpt+suffix,             ES->Ntrk(),   pt_VB, weight);
      fillh1f(hSumPt+suffix,               ES->SumPt(),         weight);
      fillh2f(hSumPtPtVB+suffix2d,         ES->SumPt(),  pt_VB, weight);
      fillh2f(hSumPtZpt+suffix,            ES->SumPt(),  pt_VB, weight);
  }
  if (ES->Ntrk()>0) {
      fillh1f(hAvgPt+suffix,               ES->AvgPt(),         weight);
      fillh2f(hAvgPtPtVB+suffix2d,         ES->AvgPt(),  pt_VB, weight);
      fillh2f(hAvgPtZpt+suffix,         ES->AvgPt(),  pt_VB, weight);
  }

  if (doSumPtTrans && ES->Ntrk()>0) {
      fillh2f(hSumPtTransZpt+suffix,         ES->SumPtTrans(),  pt_VB, weight);
  }

  if ( ES->hasFParam() ) {
    fillh1f(hFParameter+suffix,        ES->FParameterV(),   weight); 
    fillh2f(hFParameterPtVB+suffix2d,   ES->FParameterV(),   pt_VB,   weight);
    // 2D for HBOM
    fillh2f(hFParameterZpt+suffix,   ES->FParameterV(),   pt_VB,   weight);
    
    // Write out some memorable events
    if (fuzzyEquals(ES->FParameterV(), 0.0, 1e-4)) rememberEvent(Form("F = %f", ES->FParameterV()));
    if (fuzzyEquals(ES->FParameterV(), 1.0, 1e-4)) rememberEvent(Form("F = %f", ES->FParameterV()));

  }
  if ( ES->hasSpheroThrust() ) {
    fillh1f(hSpherocity+suffix,        ES->SpherocityV(),   weight);
    fillh1f(hThrust+suffix,            ES->ThrustV(),       weight);
    fillh1f(hThrustPhi+suffix,         ES->ThrustPhiV(),    weight);
    fillh1f(hBeamThrust+suffix,        ES->BeamThrustV(),   weight);
    fillh1f(hMinor+suffix,             ES->MinorV(),        weight); 

    fillh2f(hSpherocityPtVB+suffix2d,   ES->SpherocityV(),   pt_VB,   weight);
    fillh2f(hThrustPtVB+suffix2d,       ES->ThrustV(),       pt_VB,   weight);
    fillh2f(hThrustPhiPtVB+suffix2d,    ES->ThrustPhiV(),    pt_VB,   weight);
    fillh2f(hBeamThrustPtVB+suffix2d,   ES->BeamThrustV(),   pt_VB,   weight);
    fillh2f(hMinorPtVB+suffix2d,        ES->MinorV(),        pt_VB,   weight);
    // 2D for HBOM
    fillh2f(hSpherocityZpt+suffix,   ES->SpherocityV(),   pt_VB,   weight);
    fillh2f(hThrustZpt+suffix,       ES->ThrustV(),       pt_VB,   weight);
    fillh2f(hThrustPhiZpt+suffix,    ES->ThrustPhiV(),    pt_VB,   weight);
    fillh2f(hBeamThrustZpt+suffix,   ES->BeamThrustV(),   pt_VB,   weight);
    fillh2f(hMinorZpt+suffix,        ES->MinorV(),        pt_VB,   weight);
    
    // Correlations of ES
    //
    //hThrustFParameter    
    //hThrustBeamThrust    
    //hThrustMinor         
    //hThrustSpherocity    
    //hFParameterBeamThrust
    //hFParameterMinor     
    //hFParameterSpherocity
    //hBeamThrustMinor     
    //hBeamThrustSpherocity
    //hMinorSpherocity     
    //
    //
    fillh2f(hThrustFParameter     +suffix,  ES->ThrustV()    , ES->FParameterV(),   weight);
    fillh2f(hThrustBeamThrust     +suffix,  ES->ThrustV()    , ES->BeamThrustV(),   weight);
    fillh2f(hThrustMinor          +suffix,  ES->ThrustV()    , ES->MinorV()     ,   weight);
    fillh2f(hThrustSpherocity     +suffix,  ES->ThrustV()    , ES->SpherocityV(),   weight);
    fillh2f(hFParameterBeamThrust +suffix,  ES->FParameterV(), ES->BeamThrustV(),   weight);
    fillh2f(hFParameterMinor      +suffix,  ES->FParameterV(), ES->MinorV()     ,   weight);
    fillh2f(hFParameterSpherocity +suffix,  ES->FParameterV(), ES->SpherocityV(),   weight);
    fillh2f(hBeamThrustMinor      +suffix,  ES->BeamThrustV(), ES->MinorV()     ,   weight);
    fillh2f(hBeamThrustSpherocity +suffix,  ES->BeamThrustV(), ES->SpherocityV(),   weight);
    fillh2f(hMinorSpherocity      +suffix,  ES->MinorV()     , ES->SpherocityV(),   weight);
    
    
    // Write out some memorable events
    if (fuzzyEquals(ES->SpherocityV(), 0.0, 1e-4)) rememberEvent(Form("S = %f", ES->SpherocityV()));
    if (fuzzyEquals(ES->SpherocityV(), 1.0, 1e-4)) rememberEvent(Form("S = %f", ES->SpherocityV()));
    if (fuzzyEquals(ES->ThrustV(), 2.0/M_PI, 1e-3)) rememberEvent(Form("T = %f", ES->ThrustV()));
    if (fuzzyEquals(ES->ThrustV(), 1.0, 1e-3)) rememberEvent("T = %f", ES->ThrustV());

  }                                    
                                       
}                                      

void UEShapesNoVtx::fillEventShapes3D(int suffix, EventShapes* ES, double pt_VB, double mll,double weight, double z0, bool doSumPtTrans) {

  if (suffix != hsRecSel3D) return;
  
  if (ES->Ntrk()>=0) {
    fillh3f(h3DNtrkZptMll       + suffix, ES->Ntrk(),        pt_VB, mll, weight);
    fillh3f(h3DSumPtZptMll      + suffix, ES->SumPt(),       pt_VB, mll, weight);
  }                                                          
  //if (ES->Ntrk()>0) {                                        
    //fillh3f(h3DAvgPtZptMll      + suffix, ES->AvgPt(),       pt_VB, mll, weight);
  //}
  if (doSumPtTrans && ES->Ntrk()>0) {
      fillh3f(hSumPtTransZpt+suffix,         ES->SumPtTrans(),  pt_VB, mll, weight);
  }

  if ( ES->hasFParam() ) {
    fillh3f(h3DFParameterZptMll + suffix, ES->FParameterV(), pt_VB, mll, weight);
  }
  if ( ES->hasSpheroThrust() ) {
    fillh3f(h3DSpherocityZptMll + suffix, ES->SpherocityV(), pt_VB, mll, weight);
    fillh3f(h3DThrustZptMll     + suffix, ES->ThrustV(),     pt_VB, mll, weight);
    //fillh3f(h3DThrustPhiZptMll  + suffix, ES->ThrustPhiV(),  pt_VB, mll, weight);
    fillh3f(h3DBeamThrustZptMll + suffix, ES->BeamThrustV(), pt_VB, mll, weight);
    fillh3f(h3DMinorZptMll      + suffix, ES->MinorV(),      pt_VB, mll, weight);
  }                                    
}                                      
                                       
//void UEShapesNoVtx::fillEventShapesDump(int suffix, EventShapes* ESgen, EventShapes* ESreco, double pt_vb) {
void UEShapesNoVtx::fillEventShapesDump(int suffix, EventShapes* ESreco,  double pt_vb, vector < vector<double> > stl_vec) {
    // Dump values to vector
    if ( m_dumpReco) {
        if ( ESreco->hasSpheroThrust() && pt_vb < 6) {
            double dThrust = (1- 2./M_PI)/3.;
            //if (ESreco->SpherocityV() <0.01 ||  ESreco->SpherocityV() >0.9) {
            //if (  ESreco->ThrustV() < 2./M_PI + dThrust) {
            //if (  ESreco->ThrustV() > 1 - dThrust) {
                //cout << 2./M_PI + dThrust<< " "  << 1. -dThrust << " " << dThrust << "("<<ESreco->ThrustV()<<")" <<endl;
            if (  ESreco->ThrustV() > 2./M_PI + dThrust && ESreco->ThrustV() < 1. -dThrust ) {
                dump_vec_event.push_back(VB->event()); // MC event number
                dump_vec_event.push_back(ESreco->ThrustV()    ); // 
                dump_vec_event.push_back(ESreco->MinorV()     ); //
                dump_vec_event.push_back(ESreco->SpherocityV()); // 
                dump_vec_event.push_back(ESreco->FParameterV()); // 
                dump_vec_event.push_back(ESreco->ThrustPhiV()    ); // 
                dump_vec_event.push_back(ESreco->SpherocityPhiV()); // 
                dump_vec_event.push_back(pt_vb); // MC event number
                dump_vec_event.push_back(ESreco->Ntrk() ); // Ntrk
                dump_vec_event.push_back(ESreco->SumPt()); // Ntrk
                dump_vec_event.push_back(ESreco->AvgPt()); // Ntrk
            
                dump_vec_event.push_back(ESreco->BeamThrustV()); // 

    dump_eventnumber.push_back(VB->event());    
    dump_vec.push_back(dump_vec_event);
    for (int i=0;i<stl_vec.size();i++){
        dump_vec.push_back(stl_vec[i]);
    }

    dump_vec_event.clear();            
            }
        }
    }
}

void UEShapesNoVtx::fillEventShapeSmear(int suffix, EventShapes* ES1, double pt_VB_1, double evt_weight) {

  if (m_HBOMFile.Sizeof()>1 &! m_dumpReco) {
     //cout << "Fill da shit " << suffix << endl; 
    // Fill the actual resolutions histos RecSel
      double ntrk  = ES1->Ntrk();
      double sumpt = ES1->SumPt();
      double avgpt = ES1->AvgPt();
      double c_ntrk  = GetHBOMCorr("Ntrk" , ntrk );
      double c_sumpt = GetHBOMCorr("SumPt", sumpt);
      double c_avgpt = GetHBOMCorr("AvgPt", avgpt);
      double e_ntrk  = GetHBOMCorrErr("Ntrk" , ntrk );
      double e_sumpt = GetHBOMCorrErr("SumPt", sumpt);
      double e_avgpt = GetHBOMCorrErr("AvgPt", avgpt);
      fillh2f(hNtrkHBOMSmear   + suffix, c_ntrk *ntrk , ntrk *(gRandom->Gaus(0., e_ntrk ) + c_ntrk ), evt_weight);
      fillh2f(hSumPtHBOMSmear  + suffix, c_sumpt*sumpt, sumpt*(gRandom->Gaus(0., e_sumpt) + c_sumpt), evt_weight);
      fillh2f(hAvgPtHBOMSmear  + suffix, c_avgpt*avgpt, avgpt*(gRandom->Gaus(0., e_avgpt) + c_avgpt), evt_weight);

      if (ES1->hasFParam()) {
        double fp = ES1->FParameterV();
        double c_HBOM = GetHBOMCorr("FParameter",fp);
        double e_HBOM = GetHBOMCorrErr("FParameter",fp);
        fillh2f(hFParameterHBOMSmear  + suffix, c_HBOM*fp, fp*(gRandom->Gaus(0.,e_HBOM) + c_HBOM), evt_weight);
      }
      
      if ( ES1->hasSpheroThrust() ) {
        double sp = ES1->SpherocityV();
        double th = ES1->ThrustV()    ;
        double bt = ES1->BeamThrustV();
        double mn = ES1->MinorV()     ;
        double c_sp = GetHBOMCorr("Spherocity", sp);
        double c_th = GetHBOMCorr("Thrust"    , th);
        double c_bt = GetHBOMCorr("BeamThrust", bt);
        double c_mn = GetHBOMCorr("Minor"     , mn);
        double e_sp = GetHBOMCorrErr("Spherocity", sp);
        double e_th = GetHBOMCorrErr("Thrust"    , th);
        double e_bt = GetHBOMCorrErr("BeamThrust", bt);
        double e_mn = GetHBOMCorrErr("Minor"     , mn);
        fillh2f(hSpherocityHBOMSmear + suffix, c_sp*sp, sp*(gRandom->Gaus(0.,e_sp) + c_sp), evt_weight);
        fillh2f(hThrustHBOMSmear     + suffix, c_th*th, th*(gRandom->Gaus(0.,e_th) + c_th), evt_weight);
        fillh2f(hBeamThrustHBOMSmear + suffix, c_bt*bt, bt*(gRandom->Gaus(0.,e_bt) + c_bt), evt_weight);
        fillh2f(hMinorHBOMSmear      + suffix, c_mn*mn, mn*(gRandom->Gaus(0.,e_mn) + c_mn), evt_weight);
      }
  }
}

void UEShapesNoVtx::fillEventShapeSmear(int suffix, EventShapes* ES1, EventShapes* ES2,
				   double pt_VB_1, double pt_VB_2, double evt_weight) {

  fillh2f(hPtVBSmear  + suffix, pt_VB_1,      pt_VB_2,      evt_weight);
  fillh2f(hNtrkSmear  + suffix, ES1->Ntrk(),  ES2->Ntrk(),  evt_weight);
  fillh2f(hSumPtSmear + suffix, ES1->SumPt(), ES2->SumPt(), evt_weight);
  fillh2f(hAvgPtSmear + suffix, ES1->AvgPt(), ES2->AvgPt(), evt_weight);

  if ( ES1->hasFParam() && ES2->hasFParam() )
    fillh2f(hFParameterSmear + suffix, ES1->FParameterV(), ES2->FParameterV(), evt_weight);
 
  if ( ES1->hasSpheroThrust() && ES2->hasSpheroThrust() ) {
    
    fillh2f(hSpherocitySmear + suffix, ES1->SpherocityV(), ES2->SpherocityV(), evt_weight);
    fillh2f(hThrustSmear     + suffix, ES1->ThrustV(),     ES2->ThrustV(),     evt_weight);
    fillh2f(hThrustPhiSmear  + suffix, ES1->ThrustPhiV(),  ES2->ThrustPhiV(),  evt_weight);
    fillh2f(hBeamThrustSmear + suffix, ES1->BeamThrustV(), ES2->BeamThrustV(), evt_weight);
    fillh2f(hMinorSmear      + suffix, ES1->MinorV(),      ES2->MinorV(),      evt_weight);
  }


  // Also use this function to fill the HBOM correction smearing histos
  // Keep in mind htat we store relative differences here
  // hMeanCorrFP->GetBinContent(hMeanCorrFP->FindBin(.5)) 
  // ES1 is MC truth
  // ES2 is reco
  //
  if (m_HBOMFile.Sizeof()>1 &! m_dumpReco) {
      // On demand filling of mean HBOM correction histo
      double ntrk  = ES2->Ntrk();
      double sumpt = ES2->SumPt();
      double avgpt = ES2->AvgPt();
      double c_ntrk  = GetHBOMCorr("Ntrk" , ntrk );
      double c_sumpt = GetHBOMCorr("SumPt", sumpt);
      double c_avgpt = GetHBOMCorr("AvgPt", avgpt);
      double e_ntrk  = GetHBOMCorrErr("Ntrk" , ntrk );
      double e_sumpt = GetHBOMCorrErr("SumPt", sumpt);
      double e_avgpt = GetHBOMCorrErr("AvgPt", avgpt);
      //fillh2f(hNtrkHBOMSmear   + suffix, c_ntrk *ntrk , ntrk *(gRandom->Gaus(0., e_ntrk ) + c_ntrk ), evt_weight);
      //fillh2f(hSumPtHBOMSmear  + suffix, c_sumpt*sumpt, sumpt*(gRandom->Gaus(0., e_sumpt) + c_sumpt), evt_weight);
      //fillh2f(hAvgPtHBOMSmear  + suffix, c_avgpt*avgpt, avgpt*(gRandom->Gaus(0., e_avgpt) + c_avgpt), evt_weight);
      fillh2f(hNtrkHBOMSmear   + suffix, ES1->Ntrk() , ntrk *(gRandom->Gaus(0., e_ntrk ) + c_ntrk ), evt_weight);
      fillh2f(hSumPtHBOMSmear  + suffix, ES1->SumPt(), sumpt*(gRandom->Gaus(0., e_sumpt) + c_sumpt), evt_weight);
      fillh2f(hAvgPtHBOMSmear  + suffix, ES1->AvgPt(), avgpt*(gRandom->Gaus(0., e_avgpt) + c_avgpt), evt_weight);

      if (ES2->hasFParam()) {
        double fp = ES2->FParameterV();
        double c_HBOM = GetHBOMCorr("FParameter",fp);
        double e_HBOM = GetHBOMCorrErr("FParameter",fp);
        //fillh2f(hFParameterHBOMSmear  + suffix, c_HBOM*fp, fp*(gRandom->Gaus(0.,e_HBOM) + c_HBOM), evt_weight);
        fillh2f(hFParameterHBOMSmear  + suffix,  ES1->FParameterV(), fp*(gRandom->Gaus(0.,e_HBOM) + c_HBOM), evt_weight);
      }
      if ( ES2->hasSpheroThrust() ) {
        double sp = ES2->SpherocityV();
        double th = ES2->ThrustV()    ;
        double bt = ES2->BeamThrustV();
        double mn = ES2->MinorV()     ;
        double c_sp = GetHBOMCorr("Spherocity", sp);
        double c_th = GetHBOMCorr("Thrust"    , th);
        double c_bt = GetHBOMCorr("BeamThrust", bt);
        double c_mn = GetHBOMCorr("Minor"     , mn);
        double e_sp = GetHBOMCorrErr("Spherocity", sp);
        double e_th = GetHBOMCorrErr("Thrust"    , th);
        double e_bt = GetHBOMCorrErr("BeamThrust", bt);
        double e_mn = GetHBOMCorrErr("Minor"     , mn);
        //fillh2f(hSpherocityHBOMSmear + suffix, c_sp*sp, sp*(gRandom->Gaus(0.,e_sp) + c_sp), evt_weight);
        //fillh2f(hThrustHBOMSmear     + suffix, c_th*th, th*(gRandom->Gaus(0.,e_th) + c_th), evt_weight);
        //fillh2f(hBeamThrustHBOMSmear + suffix, c_bt*bt, bt*(gRandom->Gaus(0.,e_bt) + c_bt), evt_weight);
        //fillh2f(hMinorHBOMSmear      + suffix, c_mn*mn, mn*(gRandom->Gaus(0.,e_mn) + c_mn), evt_weight);
        fillh2f(hSpherocityHBOMSmear + suffix, ES1->SpherocityV(), sp*(gRandom->Gaus(0.,e_sp) + c_sp), evt_weight);
        fillh2f(hThrustHBOMSmear     + suffix, ES1->ThrustV()    , th*(gRandom->Gaus(0.,e_th) + c_th), evt_weight);
        fillh2f(hBeamThrustHBOMSmear + suffix, ES1->BeamThrustV(), bt*(gRandom->Gaus(0.,e_bt) + c_bt), evt_weight);
        fillh2f(hMinorHBOMSmear      + suffix, ES1->MinorV()     , mn*(gRandom->Gaus(0.,e_mn) + c_mn), evt_weight);
      }
  }
}

void UEShapesNoVtx::fillEventShapeDiff(int suffix, EventShapes* ES1, EventShapes* ES2,
				  double pt_VB_1, double pt_VB_2, double evt_weight) {

  fillh2f(hPtVBDiff  +suffix,   pt_VB_1,      pt_VB_2-pt_VB_1,            evt_weight);
  fillh2f(hNtrkDiff  +suffix,   ES1->Ntrk(),  ES2->Ntrk()-ES1->Ntrk(),    evt_weight);
  fillh2f(hSumPtDiff +suffix,   ES1->SumPt(), ES2->SumPt()-ES1->SumPt(),  evt_weight);
  fillh2f(hAvgPtDiff +suffix,   ES1->AvgPt(), ES2->AvgPt()-ES1->AvgPt(),  evt_weight);

  if ( ES1->hasFParam() && ES2->hasFParam() )
    fillh2f(hFParameterDiff+suffix, ES1->FParameterV(), ES2->FParameterV()-ES1->FParameterV(), evt_weight); 

  if ( ES1->hasSpheroThrust() && ES2->hasSpheroThrust() ) {

    fillh2f(hSpherocityDiff+suffix, ES1->SpherocityV(), ES2->SpherocityV()-ES1->SpherocityV(), evt_weight);
    fillh2f(hThrustDiff    +suffix, ES1->ThrustV(),     ES2->ThrustV()-ES1->ThrustV(),         evt_weight);
    fillh2f(hThrustPhiDiff +suffix, ES1->ThrustPhiV(),  ES2->ThrustPhiV()-ES1->ThrustPhiV(),   evt_weight);
    fillh2f(hBeamThrustDiff+suffix, ES1->BeamThrustV(), ES2->BeamThrustV()-ES1->BeamThrustV(), evt_weight);
    fillh2f(hMinorDiff     +suffix, ES1->MinorV(),      ES2->MinorV()-ES1->MinorV(),           evt_weight); 
  }
}

void UEShapesNoVtx::fillCTCorrection(EventShapes* recES, double evt_weight) {

  double Nsel = recES->Ntrk();

  fillh2f(hNtrkN+hssel,              recES->Ntrk(),          Nsel,   evt_weight);
  fillh2f(hSumPtN+hssel,             recES->SumPt(),         Nsel,   evt_weight);
  fillh2f(hAvgPtN+hssel,             recES->AvgPt(),         Nsel,   evt_weight);

  if (recES->hasFParam()) {
    fillh2f(hFParameterN+hssel,      recES->FParameterV(),   Nsel,   evt_weight);
  }
  if (recES->hasSpheroThrust()) {
    fillh2f(hSpherocityN+hssel,      recES->SpherocityV(),   Nsel,   evt_weight);
    fillh2f(hThrustN+hssel,          recES->ThrustV(),       Nsel,   evt_weight);
    fillh2f(hThrustPhiN+hssel,       recES->ThrustPhiV(),    Nsel,   evt_weight);
    fillh2f(hBeamThrustN+hssel,      recES->BeamThrustV(),   Nsel,   evt_weight);
    fillh2f(hMinorN+hssel,           recES->MinorV(),        Nsel,   evt_weight);
  }
}

void UEShapesNoVtx::fillOACorrection(EventShapes* recES, EventShapes* genES, double pt_VB, double evt_weight) {
  
  // Out-of-Acceptance and Contamination corrections
  
  int suffOA = ERRORVAL;
  int suffCT = ERRORVAL;

  // should always be an integer for mcparticles or selected tracks with no weights
  int Nch  = (int)genES->Ntrk(); 

  fillh2f(hNtrkN+hsch,              genES->Ntrk(),          Nch,   evt_weight);
  fillh2f(hSumPtN+hsch,             genES->SumPt(),         Nch,   evt_weight);
  fillh2f(hAvgPtN+hsch,             genES->AvgPt(),         Nch,   evt_weight);

  if (genES->hasFParam()) {
    fillh2f(hFParameterN+hsch,      genES->FParameterV(),   Nch,   evt_weight);
  }
  if (genES->hasSpheroThrust()) {
    fillh2f(hSpherocityN+hsch,      genES->SpherocityV(),   Nch,   evt_weight);
    fillh2f(hThrustN+hsch,          genES->ThrustV(),       Nch,   evt_weight);
    fillh2f(hThrustPhiN+hsch,       genES->ThrustPhiV(),    Nch,   evt_weight);
    fillh2f(hBeamThrustN+hsch,      genES->BeamThrustV(),   Nch,   evt_weight);
    fillh2f(hMinorN+hsch,           genES->MinorV(),        Nch,   evt_weight);
  }

  switch((int)recES->Ntrk()) {
  case 0: 
    suffOA = hsGenNsel0;
    break;
  case 1: 
    suffOA = hsGenNsel1;
    break;
  case 2: 
    suffOA = hsGenNsel2;
    break;
  }

  switch((int)genES->Ntrk()) {
  case 0: 
    suffCT = hsRecNch0;
    break;
  case 1: 
    suffCT = hsRecNch1;
    break;
  case 2: 
    suffCT = hsRecNch2;
    break;
  }

  if (!recES->hasFParam() && genES->hasFParam() ) { // nch>=2 && nsel<2
    fillh2f(hFParameterPtVB+suffOA,    genES->FParameterV(),   pt_VB,   evt_weight);
  }
  else if (recES->hasFParam() && !genES->hasFParam() ) { // nch<2 && nsel>=2
    fillh2f(hFParameterPtVB+suffCT,    recES->FParameterV(),   pt_VB,   evt_weight);
  }

  if (!recES->hasSpheroThrust() && genES->hasSpheroThrust() ) {
    fillh2f(hSpherocityPtVB+suffOA,    genES->SpherocityV(),   pt_VB,   evt_weight);
    fillh2f(hThrustPtVB+suffOA,        genES->ThrustV(),       pt_VB,   evt_weight);
    fillh2f(hThrustPhiPtVB+suffOA,     genES->ThrustPhiV(),    pt_VB,   evt_weight);
    fillh2f(hBeamThrustPtVB+suffOA,    genES->BeamThrustV(),   pt_VB,   evt_weight);
    fillh2f(hMinorPtVB+suffOA,         genES->MinorV(),        pt_VB,   evt_weight); 
  }
  else if (recES->hasSpheroThrust() && !genES->hasSpheroThrust() ) {
    fillh2f(hSpherocityPtVB+suffCT,    recES->SpherocityV(),   pt_VB,   evt_weight);
    fillh2f(hThrustPtVB+suffCT,        recES->ThrustV(),       pt_VB,   evt_weight);
    fillh2f(hThrustPhiPtVB+suffCT,     recES->ThrustPhiV(),    pt_VB,   evt_weight);
    fillh2f(hBeamThrustPtVB+suffCT,    recES->BeamThrustV(),   pt_VB,   evt_weight);
    fillh2f(hMinorPtVB+suffCT,         recES->MinorV(),        pt_VB,   evt_weight); 
  }
}


// void UEShapesNoVtx::fillMigMatrix(int hES, int suff, double valRec, double valGen, double VBpt) {

//   // note that VBpt here is the reconstructed value!

//   if (suff==hsRecSel) {

//     fillh2f(hES+hsMat_Sel, valRec, valGen, evt_weight);

//     /*
//     if (VBpt < 4.0)  fillh2f(hES+hsSel_PtVB4,  valRec, valGen, evt_weight);
//     if (VBpt < 8.0)  fillh2f(hES+hsSel_PtVB8,  valRec, valGen, evt_weight);
//     if (VBpt < 12.0) fillh2f(hES+hsSel_PtVB12, valRec, valGen, evt_weight);
//     if (VBpt < 16.0) fillh2f(hES+hsSel_PtVB16, valRec, valGen, evt_weight);
//     if (VBpt < 20.0) fillh2f(hES+hsSel_PtVB20, valRec, valGen, evt_weight);
//     */
//   }
//   else if (suff==hsRecCorr) {

//     fillh2f(hES+hsMat_Corr, valRec, valGen, evt_weight);

//     /*
//     if (VBpt < 4.0)  fillh2f(hES+hsCorr_PtVB4,  valRec, valGen, evt_weight);
//     if (VBpt < 8.0)  fillh2f(hES+hsCorr_PtVB8,  valRec, valGen, evt_weight);
//     if (VBpt < 12.0) fillh2f(hES+hsCorr_PtVB12, valRec, valGen, evt_weight);
//     if (VBpt < 16.0) fillh2f(hES+hsCorr_PtVB16, valRec, valGen, evt_weight);
//     if (VBpt < 20.0) fillh2f(hES+hsCorr_PtVB20, valRec, valGen, evt_weight);
//     */
//   }
// }


// This little friend is automatically generated with d3pdreader/makepar -w
void UEShapesNoVtx::fillTrackD3PDObjectControlPlots(int suffix, int index, double w_trk) {

  /*  
  fillh1f(htrk_d0 +suffix, (float) (*trk)[index].d0(), w_trk);
  fillh1f(htrk_z0 +suffix, (float) (*trk)[index].z0(), w_trk);
  fillh1f(htrk_phi +suffix, (float) (*trk)[index].phi(), w_trk);
  fillh1f(htrk_theta +suffix, (float) (*trk)[index].theta(), w_trk);
  fillh1f(htrk_qoverp +suffix, (float) (*trk)[index].qoverp(), w_trk);
  fillh1f(htrk_pt +suffix, (float) (*trk)[index].pt(), w_trk);
  fillh1f(htrk_eta +suffix, (float) (*trk)[index].eta(), w_trk);
  fillh1f(htrk_d0_wrtPV +suffix, (float) (*trk)[index].d0_wrtPV(), w_trk);
  fillh1f(htrk_z0_wrtPV +suffix, (float) (*trk)[index].z0_wrtPV(), w_trk);
  fillh1f(htrk_phi_wrtPV +suffix, (float) (*trk)[index].phi_wrtPV(), w_trk);
  fillh1f(htrk_chi2 +suffix, (float) (*trk)[index].chi2(), w_trk);
  fillh1f(htrk_ndof +suffix, (float) (*trk)[index].ndof(), w_trk);
  fillh1f(htrk_nBLHits +suffix, (float) (*trk)[index].nBLHits(), w_trk);
  fillh1f(htrk_nPixHits +suffix, (float) (*trk)[index].nPixHits(), w_trk);
  fillh1f(htrk_nSCTHits +suffix, (float) (*trk)[index].nSCTHits(), w_trk);
  fillh1f(htrk_nTRTHits +suffix, (float) (*trk)[index].nTRTHits(), w_trk);
  fillh1f(htrk_nTRTHighTHits +suffix, (float) (*trk)[index].nTRTHighTHits(), w_trk);
  fillh1f(htrk_nPixHoles +suffix, (float) (*trk)[index].nPixHoles(), w_trk);
  fillh1f(htrk_nSCTHoles +suffix, (float) (*trk)[index].nSCTHoles(), w_trk);
  fillh1f(htrk_nTRTHoles +suffix, (float) (*trk)[index].nTRTHoles(), w_trk);
  fillh1f(htrk_expectBLayerHit +suffix, (float) (*trk)[index].expectBLayerHit(), w_trk);
  //fillh1f(htrk_nMDTHits +suffix, (float) (*trk)[index].nMDTHits(), w_trk);
  //fillh1f(htrk_nCSCEtaHits +suffix, (float) (*trk)[index].nCSCEtaHits(), w_trk);
  //fillh1f(htrk_nCSCPhiHits +suffix, (float) (*trk)[index].nCSCPhiHits(), w_trk);
  //fillh1f(htrk_nRPCEtaHits +suffix, (float) (*trk)[index].nRPCEtaHits(), w_trk);
  //fillh1f(htrk_nRPCPhiHits +suffix, (float) (*trk)[index].nRPCPhiHits(), w_trk);
  //fillh1f(htrk_nTGCEtaHits +suffix, (float) (*trk)[index].nTGCEtaHits(), w_trk);
  //fillh1f(htrk_nTGCPhiHits +suffix, (float) (*trk)[index].nTGCPhiHits(), w_trk);
  fillh1f(htrk_nHits +suffix, (float) (*trk)[index].nHits(), w_trk);
  fillh1f(htrk_nHoles +suffix, (float) (*trk)[index].nHoles(), w_trk);
  fillh1f(htrk_hitPattern +suffix, (float) (*trk)[index].hitPattern(), w_trk);
  //fillh1f(htrk_TRTHighTHitsRatio +suffix, (float) (*trk)[index].TRTHighTHitsRatio(), w_trk);
  //fillh1f(htrk_TRTHighTOutliersRatio +suffix, (float) (*trk)[index].TRTHighTOutliersRatio(), w_trk);
  fillh1f(htrk_fitter +suffix, (float) (*trk)[index].fitter(), w_trk);
  fillh1f(htrk_patternReco1 +suffix, (float) (*trk)[index].patternReco1(), w_trk);
  fillh1f(htrk_patternReco2 +suffix, (float) (*trk)[index].patternReco2(), w_trk);
  //fillh1f(htrk_seedFinder +suffix, (float) (*trk)[index].seedFinder(), w_trk);
  //fillh1f(htrk_mc_probability +suffix, (float) (*trk)[index].mc_probability(), w_trk);
  
  */
}


     /*
      bool foundW = false;

      for (uint i = 0; i < mc->n(); i++ ) {

	// cut on primary final-state particles
	if (!truthCutPrimary(*mc,i)) continue;

	int pid = (int)fabs((*mc)[i].pdgId());

	if (isChargedLepton(pid)) {

	  for (uint j = 0; j < mc->n(); j++ ) {

	    if (i==j) continue;
	    if (!truthCutPrimary(*mc,j)) continue;

	    int pid2 = (int)fabs((*mc)[j].pdgId());

	    if (isNeutrino(pid2)) {
	      bool isW = false;
	      isW = isWLeptonNeutrinoPair(*mc,i,j);

	      if (isW) {     
		if (foundW) l.error(TString::Format("More than one W found in event (Run #%i, Evt #%i)", VB->run(), VB->event()));
		foundW = true;
		lep3vec = getTrack3Vector(*mc,i);
		nu3vec  = getTrack3Vector(*mc,j);
	      }
	    }
	  }
	}
      }

      if (!foundW) l.error(TString::Format("Did not find W in event (Run #%i, Evt #%i)", VB->run(), VB->event()));
      */


/*
int l_tracks = 0; // number of removed 'lepton' tracks

    int index_ptdiff1, index_deltaR1, index_ptdiff2, index_deltaR2;
    float smallest_ptdiff1 = 999.;
    float deltaR_smallest_ptdiff1 = 999.;
    float smallest_deltaR1 = 999.;
    float ptdiff_smallest_deltaR1 = 999.;
    float smallest_ptdiff2 = 999.;
    float deltaR_smallest_ptdiff2 = 999.;
    float smallest_deltaR2 = 999.;
    float ptdiff_smallest_deltaR2 = 999.;


    
    // Find lep1 and lep2 tracks
    // -------------------------
    for (uint i = 0; i < trk->n(); i++ ) {

      TVector3 trk3vec = getTrack3Vector(*trk, i);

      double ptdiff1 = fabs((trk3vec.Pt()-lep1_trk_3vec.Pt())/lep1_trk_3vec.Pt());
      double dR1     = deltaR(lep1_trk_3vec,trk3vec);

      double ptdiff2 = fabs((trk3vec.Pt()-lep2_trk_3vec.Pt())/lep2_trk_3vec.Pt());
      double dR2     = deltaR(lep2_trk_3vec,trk3vec);

      if (ptdiff1 < smallest_ptdiff1) {
	smallest_ptdiff1 = ptdiff1;
	deltaR_smallest_ptdiff1 = dR1;
	index_ptdiff1 = i;
      }
      if (dR1 < smallest_deltaR1) {
	smallest_deltaR1 = dR1;
	ptdiff_smallest_deltaR1 = ptdiff1;
	index_deltaR1 = i;
      }
      if (ptdiff2 < smallest_ptdiff2) {
	smallest_ptdiff2 = ptdiff2;
	deltaR_smallest_ptdiff2 = dR2;
	index_ptdiff2 = i;
      }
      if (dR2 < smallest_deltaR2) {
	smallest_deltaR2 = dR2;
	ptdiff_smallest_deltaR2 = ptdiff2;
	index_deltaR2 = i;
      }
    }

    if (index_ptdiff1 != index_deltaR1) {
      rememberEvent("ptdiff1 and deltaR1 indices don't match!");
      TVector3 trk3vecA = getTrack3Vector(*trk,index_ptdiff1);
      TVector3 trk3vecB = getTrack3Vector(*trk,index_deltaR1);

      cout << "lep1 trk pt = " << lep1_trk_3vec.Pt() << endl;
      cout << "lep1 trk eta = " << lep1_trk_3vec.Eta() << endl;
      cout << "lep1 trk phi = " << lep1_trk_3vec.Phi() << endl;

      cout << "lep1 cl pt = " << lep1_cl_3vec.Pt() << endl;
      cout << "lep1 cl eta = " << lep1_cl_3vec.Eta() << endl;
      cout << "lep1 cl phi = " << lep1_cl_3vec.Phi() << endl;

      cout << "trkA pt = " << trk3vecA.Pt() << endl;
      cout << "trkA eta = " << trk3vecA.Eta() << endl;
      cout << "trkA phi = " << trk3vecA.Phi() << endl;

      cout << "trkB pt = " << trk3vecB.Pt() << endl;
      cout << "trkB eta = " << trk3vecB.Eta() << endl;
      cout << "trkB phi = " << trk3vecB.Phi() << endl;

      cout << "smallest_ptdiff1 = " << smallest_ptdiff1 << endl;
      cout << "deltaR for smallest_ptdiff1 = " << deltaR_smallest_ptdiff1 << endl;
      cout << "smallest_deltaR1 = " << smallest_deltaR1 << endl;
      cout << "ptdiff_smallest_deltaR1 = " << ptdiff_smallest_deltaR1 << endl;
    }

    if (index_ptdiff2 != index_deltaR2) {
      rememberEvent("ptdiff2 and deltaR2 indices don't match!");
      TVector3 trk3vecA = getTrack3Vector(*trk,index_ptdiff2);
      TVector3 trk3vecB = getTrack3Vector(*trk,index_deltaR2);

      cout << "lep2 trk pt = " << lep2_trk_3vec.Pt() << endl;
      cout << "lep2 trk eta = " << lep2_trk_3vec.Eta() << endl;
      cout << "lep2 trk phi = " << lep2_trk_3vec.Phi() << endl;

      cout << "lep2 cl pt = " << lep2_cl_3vec.Pt() << endl;
      cout << "lep2 cl eta = " << lep2_cl_3vec.Eta() << endl;
      cout << "lep2 cl phi = " << lep2_cl_3vec.Phi() << endl;

      cout << "trkA pt = " << trk3vecA.Pt() << endl;
      cout << "trkA eta = " << trk3vecA.Eta() << endl;
      cout << "trkA phi = " << trk3vecA.Phi() << endl;

      cout << "trkB pt = " << trk3vecB.Pt() << endl;
      cout << "trkB eta = " << trk3vecB.Eta() << endl;
      cout << "trkB phi = " << trk3vecB.Phi() << endl;

      cout << "smallest_ptdiff2 = " << smallest_ptdiff2 << endl;
      cout << "deltaR for smallest_ptdiff2 = " << deltaR_smallest_ptdiff2 << endl;
      cout << "smallest_deltaR2 = " << smallest_deltaR2 << endl;
      cout << "ptdiff_smallest_deltaR2 = " << ptdiff_smallest_deltaR2 << endl;
    }
*/