#define ElecTauAnalysis_cxx
#include "ElecTauAnalysis.h"
#include "pdfUtils.h"
#include "efficiencies.h"
#include <TStyle.h>
#include <TCanvas.h>
#include <stdlib.h>
#include <fstream>
#include <iomanip>

using namespace std;

void ElecTauAnalysis::Loop(){
  // Book histos and initialize out tree
  bookHistos();
  
  if (fChain == 0) return;
  //Long64_t nentries = _nEvents;
  if(_nEvents < 0 || _nEvents > fChain->GetEntries()) _nEvents = fChain->GetEntries();
  double eventsFraction = (double)_nEvents/(double)fChain->GetEntries();
  //double eventsFraction = 1.;
  //cross section in pico barns
  if (_source == "Signal"){
   _theXSection = make_pair<float, float>(_signalXSection, 0); // 0% error on xsection
   skimmedEvents = make_pair<int, int>((int)_nEvents, (int)(50400.*eventsFraction));
  }
  else if (_source == "ZTauTau"){
   _theXSection = make_pair<float, float>(1915.08, 0); // NNLO Z -> ll (M>20)				         
   skimmedEvents = make_pair<int, int>(51046, (int)(3295238.*eventsFraction));  		         
   _doFactorization = false; 		         
  }													         
  else if (_source == "ZEE"){										         
   _theXSection = make_pair<float, float>(1915.08, 0); // NNLO Z -> ll (M>20)				         
   skimmedEvents = make_pair<int, int>(4716085, (int)(18162170.*eventsFraction));
   _doFactorization = false; 		         
  }													         
  else if (_source == "TTJets"){									         
   _theXSection = make_pair<float, float>(245.8, 0); //NNLO					         
   skimmedEvents = make_pair<int, int>(1687519, (int)(6900671.*eventsFraction));  		         
   _doFactorization = false; 		         
  }													         
  else if(_source == "WJets"){  									         
   _theXSection = make_pair<float, float>(37509.0*0.9595, 0);  //from 0-50 GeV						         
   //_theXSection = make_pair<float, float>(37509.0, 0);  //inclusive						         
   skimmedEvents = make_pair<double, double>(247809, (int)(17166890.*eventsFraction));
   if(_doWJetsCR > 0) _doFactorization = false;
   else _doFactorization = true;	         
  }													         
  else if(_source == "WJetsPtW50-70"){  									         
   _theXSection = make_pair<float, float>(811.2, 0);  //						         
   skimmedEvents = make_pair<double, double>(4362176, (int)(47739087.*eventsFraction));		         
   if(_doWJetsCR > 0) _doFactorization = false;
   else _doFactorization = true;	         
  }													         
  else if(_source == "WJetsPtW70-100"){  									         
   _theXSection = make_pair<float, float>(428.9, 0);  //						         
   skimmedEvents = make_pair<double, double>(2880289, (int)(21762635.*eventsFraction));		         
   if(_doWJetsCR > 0) _doFactorization = false;
   else _doFactorization = true;	         
  }													         
  else if(_source == "WJetsPtW100"){  									         
   _theXSection = make_pair<float, float>(228.9, 0);  //						         
   skimmedEvents = make_pair<double, double>(2379929, (int)(12742382.*eventsFraction));		         
   if(_doWJetsCR > 0) _doFactorization = false;
   else _doFactorization = true;	         
  }													         
  else if(_source == "WW"){										         
   _theXSection = make_pair<float, float>(54.838, 0); // MCFM NLO inclusive				         
   skimmedEvents = make_pair<int, int>(870228, (int)(9980431.*eventsFraction)); 		         
   _doFactorization = false; 		         
  }													         
  else if(_source == "WZ"){										         
   _theXSection = make_pair<float, float>(32.3161, 0);  //MCFM NLO M_ll > 40 inclusive  		         
   //_theXSection = make_pair<float, float>(0.7192, 0.0277);  //WZ->3leptons+nu with 3.8% error 	         
   skimmedEvents = make_pair<double, double>(748174, (int)(9933615.*eventsFraction));		         
   _doFactorization = false; 		         
  }													         
  else if(_source == "ZZ"){										         
   _theXSection = make_pair<float, float>(17.654, 0);  // NLO inclusive 				         
   skimmedEvents = make_pair<double, double>(758397, (int)(9799908.*eventsFraction));		         
   _doFactorization = false; 		         
  }							         
  else if(_source == "QCD20to30"){  									               
   _theXSection = make_pair<float, float>(2.886e8*0.0101, 0.);  // no error for qcd 		         
   skimmedEvents = make_pair<double, double>(25310, (int)(35040695.*eventsFraction)); 	     	         
   _doFactorization = true; 		         
  }													         
  else if(_source == "QCD30to80"){  									               
   _theXSection = make_pair<float, float>(7.433e7*0.0621, 0.);  // no error for qcd 		         
   skimmedEvents = make_pair<double, double>(411319, (int)(32956534.*eventsFraction)); 	     	         
   _doFactorization = true; 		         
  }													         
  else if(_source == "QCD80to170"){  									               
   _theXSection = make_pair<float, float>(1191000.0*0.1539, 0.);  // no error for qcd 		         
   skimmedEvents = make_pair<double, double>(2529357, (int)(34404593.*eventsFraction)); 	     	         
   _doFactorization = true; 		         
  }													         
  else if(_source == "QCD170to250"){  									               
   _theXSection = make_pair<float, float>(30990.0*0.148, 0.);  // no error for qcd 		         
   skimmedEvents = make_pair<double, double>(3349014, (int)(31506884.*eventsFraction)); 	     	         
   _doFactorization = true; 		         
  }													         
  
  else if(_source == "QCD250to350"){  									               
   _theXSection = make_pair<float, float>(4250.0*0.131, 0.);  // no error for qcd 		         
   skimmedEvents = make_pair<double, double>(3709410, (int)(33965250.*eventsFraction)); 	     	         
   _doFactorization = true; 		         
  }													         
  else if(_source == "QCD350"){  									               
   _theXSection = make_pair<float, float>(810.0*0.11, 0.);  // no error for qcd 		         
   skimmedEvents = make_pair<double, double>(3527663, (int)(33898802.*eventsFraction)); 	     	         
   _doFactorization = true; 		         
  }													         
 
  													         
  //get events normilized to lumi									         
  eventsAtLumi = make_pair<double, double>(_theXSection.first * _theLumi, _theXSection.second * _theLumi);       
  skimEff = GetEfficiency(skimmedEvents.first, skimmedEvents.second);
  nSkimmedAtLumi = GetSurvivingEvents(eventsAtLumi, skimEff);						         
  
  GetPUWeights();
  
  if(_doWJetsCR > 0) _doFactorization = false;
  
  Long64_t nbytes = 0, nb = 0;
  for (Long64_t jentry=0; jentry<_nEvents;jentry++) {
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    										
    resetCandCounters();
        
    resetHistosDefault();

    //get PileUp Weights
    _eventPUWeight = 1.;										       
    if(_source != "Data"){										       
      if(InTimePU < (int)_thePUWeights.size()) _eventPUWeight = _thePUWeights.at(InTimePU);		       
      if(InTimePU >= (int)_thePUWeights.size()) _eventPUWeight = _thePUWeights.at(_thePUWeights.size() -1);    
    }													       
    _puWeightedEvents += _eventPUWeight;

    if(_source != "Signal" && !passedTrigger()) continue;
    _triggerCounter++;
    
    if(_source == "WJets" && wCandGenPt > 50) continue;
    
    getCandCounters();
    
    getEventCounters();
    fillHistos();
  }
  getReport();
  writeOutFile();
}

void ElecTauAnalysis::getCandCounters(){
  unsigned int vItBegin = 0;
  unsigned int vItEnd = elecPt->size();

  for(unsigned int vIt = vItBegin; vIt < vItEnd; ++vIt){ 
    //get N-1 counters first
    getNMinus1CandCounters(vIt);
    
    if(_source == "Signal"){
      if(!passedElecMatched(vIt)) continue;
      _elecMatchCounter++;
      if(!passedTauMatched(vIt)) continue;
      _tauMatchCounter++;
    }
    
    bool passedZTTCRSel = passedElecAccCuts(vIt) && passedElecId(vIt) 
    && passedHadTauAcc(vIt) && passedHadTauId(vIt) && passedElecTauCharge(vIt) && passedElecTauCosDPhi(vIt)
    && MEt > 50. && passed2DZetaCut(vIt) && passedBTagVeto(vIt)
    && elecHadTauMass->at(vIt) <= 100. && elecMEtMt->at(vIt) < 40.; 
    if(passedZTTCRSel)_zttCRCounter++;
    if(passedZTTCRSel && passedVBFSel(vIt)){
      _zttCRVBFCounter++;
      if(diJetMass->at(vIt) > _maxZTTCRDiJetMass) _maxZTTCRDiJetMass = diJetMass->at(vIt);
    }
    
    // Set central cuts control region and vbf cuts 
    bool passedDYEECR1Sel = passedElecAccCuts(vIt) && passedElecId(vIt) && passedHadTauAcc(vIt)
    && passedHadTauDecayModeFinding(vIt) && passedHadTauIsoDisc(vIt) && passedHadTauMuonDisc(vIt)
    && hadTauIdByLooseElectronRejection->at(vIt) < 0.5
    && passedElecTauCosDPhi(vIt) && passed2DZetaCut(vIt) && passedBTagVeto(vIt) 
    && passedWJetsVeto(vIt) && passedElecTauMass(vIt);
    if(passedDYEECR1Sel) _zeeCR1CentralCutsCounter++;
    if(passedDYEECR1Sel && elecCharge->at(vIt)*hadTauCharge->at(vIt) < 0) _zeeCR1OSCounter++;
    if(passedDYEECR1Sel && elecCharge->at(vIt)*hadTauCharge->at(vIt) > 0) _zeeCR1LSCounter++;
    if(passedDYEECR1Sel && passedMET()) _zeeCR1METCounter++;
    if(passedDYEECR1Sel && passedVBFSel(vIt)){
      _zeeCR1VBFCounter++;
      if(diJetMass->at(vIt) > _maxZEECRDiJetMass) _maxZEECRDiJetMass = diJetMass->at(vIt);
    }
    
    //double zWidth = 2.4952;
    bool passedDYEECR2Sel = passedElecAccCuts(vIt) && passedElecId(vIt) && passedHadTauAcc(vIt)
    && passedHadTauDecayModeFinding(vIt) && passedTauNProngs(vIt) && passedHadTauIsoDisc(vIt) && passedHadTauMuonDisc(vIt)
    //&& (fabs(zMassDiff) < _zMassWidthCut * zWidth || fabs(zPtAsymm) < _zPtAsymmCut)
    && passedElecTauCosDPhi(vIt) && passed2DZetaCut(vIt) && passedBTagVeto(vIt) && passedWJetsVeto(vIt)
    && MEt < 30. && elecHadTauMass->at(vIt) <= 120. && elecHadTauMass->at(vIt) >= 80.;
    
    if(passedDYEECR2Sel) _zeeCR2Counter++;
    if(passedDYEECR2Sel && hadTauIdByLooseElectronRejection->at(vIt) < 0.5) _zeeCR2AntiElecDiscCounter++;
    if(passedDYEECR2Sel && passedHadTauElecDisc(vIt)) _zeeCR2ElecDiscCounter++;
    
    
    //define CR1 using MT cut, no cosDPhi and no PZeta and VLooseIso
    bool passedWJetsCR1Sel = passedElecAccCuts(vIt) && passedElecId(vIt) 
    && passedHadTauAcc(vIt) && passedHadTauDecayModeFinding(vIt) && passedTauNProngs(vIt) 
    //&& hadTauIdByVLooseCombinedIsolationDBSumPtCorr->at(vIt) > 0.5
    && hadTauIdByLooseCombinedIsolationDBSumPtCorr3Hits->at(vIt) > 0.5
    && passedHadTauElecDisc(vIt) && passedHadTauMuonDisc(vIt)
    && MEt > 30 && passedBTagVeto(vIt) && passedElecTauMass(vIt)
    && elecMEtMt->at(vIt) < _elecMETMaxMassCut && elecMEtMt->at(vIt) > _elecMETMinMassCut;
    if(passedWJetsCR1Sel) _wjetsCR1Counter++;
    
    if(passedWJetsCR1Sel && elecCharge->at(vIt)*hadTauCharge->at(vIt) > 0) _wjetsCR1LSCounter++;
    if(passedWJetsCR1Sel && elecCharge->at(vIt)*hadTauCharge->at(vIt) < 0) _wjetsCR1OSCounter++;
    //if(passedWJetsCR1Sel && passedHadTauIsoDisc(vIt)) _wjetsCR1TauIsoCounter++;
    if(passedWJetsCR1Sel && passedVBFSel(vIt)){
      _wjetsCR1VBFCounter++;
      if(diJetMass->at(vIt) > _maxWJetsCRDiJetMass) _maxWJetsCRDiJetMass = diJetMass->at(vIt);
    }

    //get efficiency of cosDPhi and Pzeta together due to correlations 
    if(passedWJetsCR1Sel && passedElecTauCosDPhi(vIt) && passed2DZetaCut(vIt)) _wjetsCR1ETauCosDPhiPZetaCounter++;
    
    if(passedWJetsCR1Sel && passedElecTauCosDPhi(vIt) && passed2DZetaCut(vIt)) _wjetsCRCentralCutsCounter++;
    if(passedWJetsCR1Sel && passedElecTauCosDPhi(vIt) && passed2DZetaCut(vIt) && passedVBFSel(vIt)) _wjetsCRVBFCounter++;
    
    bool passedWJetsCR2Sel =  passedElecAccCuts(vIt) && passedElecId(vIt) 
    && passedHadTauAcc(vIt) && passedHadTauId(vIt)
    && passedElecTauCharge(vIt) 
    && passedZEEVeto() && passedMET() && passedBTagVeto(vIt) && passedElecTauMass(vIt)
    && elecHadTauCosDPhi->at(vIt) >= _elecTauCosDPhiCut
    && elecHadTauPZeta->at(vIt) - _zeta2DSlope * elecHadTauPZetaVis->at(vIt) <= _zeta2DCut;
    if(passedWJetsCR2Sel) _wjetsCR2Counter++;
    if(passedWJetsCR2Sel && elecMEtMt->at(vIt) < _elecMETMaxMassCut && elecMEtMt->at(vIt) > _elecMETMinMassCut) _wjetsCR2MtVetoCounter++;
    
    bool passedTTbarCR1Sel = passedElecAccCuts(vIt) && passedElecId(vIt) 
    && passedHadTauAcc(vIt) && passedHadTauId(vIt) && passedElecTauCharge(vIt)
    && passedMET() && passedWJetsVeto(vIt) && passedElecTauMass(vIt) && nBTagsCombSecVtxLoose->at(vIt) >= 2;
    if(passedTTbarCR1Sel){
      _ttbarCR1Counter++;
       if(hadTauPt->at(vIt) > _maxTauPtTTbarCR1) _maxTauPtTTbarCR1 = hadTauPt->at(vIt);
    }  
    if(passedTTbarCR1Sel && passedElecTauCosDPhi(vIt) && passed2DZetaCut(vIt)) _ttbarCR1ETauCosDPhiPZetaCounter++;
    if(passedTTbarCR1Sel && elecCharge->at(vIt)*hadTauCharge->at(vIt) > 0) _ttbarCR1LSCounter++;
    if(passedTTbarCR1Sel && elecCharge->at(vIt)*hadTauCharge->at(vIt) < 0) _ttbarCR1OSCounter++;
    if(passedTTbarCR1Sel && passedVBFSel(vIt)){
      _ttbarCRVBFCounter++;
      if(diJetMass->at(vIt) > _maxTTJetsCRDiJetMass) _maxTTJetsCRDiJetMass = diJetMass->at(vIt);
    }
   
    //make pt histo before and after tau id for scale factor
    bool passedTTbarCR1NoTauId = passedElecAccCuts(vIt) && passedElecId(vIt) 
    && passedHadTauAcc(vIt) && passedElecTauCharge(vIt)
    && passedZEEVeto() && passedMET()
    && passedWJetsVeto(vIt) && passedElecTauMass(vIt) && nBTagsCombSecVtxMedium->at(vIt) >= 2;
    if(passedTTbarCR1NoTauId && hadTauPt->at(vIt) > _maxTauPtTTbarCRNoTauId) _maxTauPtTTbarCRNoTauId = hadTauPt->at(vIt);
    
    bool passedTTbarCR2Sel = passedElecAccCuts(vIt) && passedElecId(vIt) 
    && passedHadTauAcc(vIt) && passedHadTauId(vIt) && passedElecTauCharge(vIt)
    && passedZEEVeto() && MEt > 50.
    && passedWJetsVeto(vIt) && passedElecTauMass(vIt) && jetEtSum->at(vIt) >= 250.;
    if(passedTTbarCR2Sel){
    //make btag plot
      _ttbarCR2Counter++;
      if(nBtagsHiEffTrkCnt->at(vIt) < _minBtagsTCHEMTTJetsCR)_minBtagsTCHEMTTJetsCR = nBtagsHiEffTrkCnt->at(vIt);
      if(nBtagsHiPurityTrkCnt->at(vIt) < _minBtagsTCHPTTTJetsCR) _minBtagsTCHPTTTJetsCR = nBtagsHiPurityTrkCnt->at(vIt);
      if(nBTagsCombSecVtxLoose->at(vIt) < _minBtagsCSVLTTJetsCR) _minBtagsCSVLTTJetsCR = nBTagsCombSecVtxLoose->at(vIt);
      if(nBTagsCombSecVtxMedium->at(vIt) < _minBtagsCSVMTTJetsCR) _minBtagsCSVMTTJetsCR = nBTagsCombSecVtxMedium->at(vIt);
      if(nBTagsCombSecVtxTight->at(vIt) < _minBtagsCSVTTTJetsCR) _minBtagsCSVTTTJetsCR = nBTagsCombSecVtxTight->at(vIt);
    }
    if(passedTTbarCR2Sel && passedBTagVeto(vIt)) _ttbarCR2BtagVetoCounter++;
    if(passedTTbarCR2Sel && nBTagsCombSecVtxMedium->at(vIt) >= 2) _ttbarCR2TwoBtagCounter++;
    
    bool passedQCDLSSel = passedElecAccCuts(vIt) && passedElecId(vIt)
    && passedHadTauAcc(vIt) && passedHadTauId(vIt)
    && passedElecTauCosDPhi(vIt) && passedZEEVeto() && passedMET() && passed2DZetaCut(vIt)
    && passedBTagVeto(vIt) && passedWJetsVeto(vIt) && passedElecTauMass(vIt)
    && elecCharge->at(vIt)*hadTauCharge->at(vIt) > 0;
    if(passedQCDLSSel){
    //make mass plots for ls events
      _qcdLSCounter++;
      if(elecHadTauMass->at(vIt) > _maxElecTauLSMass) _maxElecTauLSMass = elecHadTauMass->at(vIt);
      if(elecHadTauMetMass->at(vIt) > _maxElecTauMetLSMass) _maxElecTauMetLSMass = elecHadTauMetMass->at(vIt);
    }
    //set qcd region for OS/LS ratio by requiring taus failing  loose iso and MT < 50
    bool passedRevTauIso = passedElecAccCuts(vIt) && passedElecId(vIt)
    && passedHadTauAcc(vIt) &&  passedHadTauDecayModeFinding(vIt) && passedTauNProngs(vIt)
    && passedHadTauMuonDisc(vIt) && passedHadTauElecDisc(vIt)
    && hadTauIdByLooseCombinedIsolationDBSumPtCorr3Hits->at(vIt) < 1
    && passedElecTauCosDPhi(vIt) && passedZEEVeto() && passedMET() && passed2DZetaCut(vIt)
    && passedBTagVeto(vIt) && passedElecTauMass(vIt) && elecMEtMt->at(vIt) < 50.;
    if(passedRevTauIso) _qcdRevTauIsoCounter++;
    if(passedRevTauIso && elecCharge->at(vIt)*hadTauCharge->at(vIt) > 0) _qcdRevTauIsoLSCounter++;
    if(passedRevTauIso && elecCharge->at(vIt)*hadTauCharge->at(vIt) < 0) _qcdRevTauIsoOSCounter++;
    //make mass plots for LS and OS for failing iso events
    
    bool _passedDiTauChargeCut = passedElecTauCharge(vIt);
    bool _passedTauIso = passedHadTauIsoDisc(vIt);
    if(_doFactorization){
      _passedDiTauChargeCut = true;
      //_passedMETCut = true; //met and MT are correlated -- can't mix them
      _passedTauIso = hadTauIdByVLooseIsolationMVA3newDMwLT->at(vIt) > 0;
    }
    
    if(elecPt->at(vIt) > _maxElecPt) _maxElecPt = elecPt->at(vIt);
    if(!passedElecPtCut(vIt)) continue;
    _elecPtCounter++;
    
    if(!elecInGaps(vIt) && fabs(elecEta->at(vIt)) < _minElecEta) _minElecEta = elecEta->at(vIt);
    if(!passedElecEtaCut(vIt)) continue;
    _elecEtaCounter++;
    if(passedElecAccCuts(vIt)) _elecAccCounter++;
    
    if(elecInBarrel(vIt) && fabs(elecDPhiIn->at(vIt)) < _minElecDPhiEB) 
    _minElecDPhiEB = elecDPhiIn->at(vIt);
    if(elecInEndCap(vIt) && fabs(elecDPhiIn->at(vIt)) < _minElecDPhiEE) 
    _minElecDPhiEE = elecDPhiIn->at(vIt);    
    if(!passedElecDPhiIn(vIt)) continue;
    _elecDPhiInCounter++;
   
    if(elecInBarrel(vIt) && fabs(elecDEtaIn->at(vIt)) < _minElecDEtaEB) 
    _minElecDEtaEB = elecDEtaIn->at(vIt);
    if(elecInEndCap(vIt) && fabs(elecDEtaIn->at(vIt)) < _minElecDEtaEE) 
    _minElecDEtaEE = elecDEtaIn->at(vIt);
    if(!passedElecDEtaIn(vIt)) continue;
    _elecDEtaInCounter++;
  
    if(elecInBarrel(vIt) && elecSigmaIEtaIEta->at(vIt) < _minElecSigmaIEtaIEtaEB) 
    _minElecSigmaIEtaIEtaEB = elecSigmaIEtaIEta->at(vIt); 
    if(elecInEndCap(vIt) && elecSigmaIEtaIEta->at(vIt) < _minElecSigmaIEtaIEtaEE) 
    _minElecSigmaIEtaIEtaEE = elecSigmaIEtaIEta->at(vIt); 
    if(!passedElecSigmaIEtaIEta(vIt)) continue;
    _elecSigmaIEtaIEtaCounter++;
  
    if(elecInBarrel(vIt) && elecHadOverE->at(vIt) < _minElecHOverEEB) 
    _minElecHOverEEB = elecHadOverE->at(vIt);
    if(elecInEndCap(vIt) && elecHadOverE->at(vIt) < _minElecHOverEEE) 
    _minElecHOverEEE = elecHadOverE->at(vIt);
    if(!passedElecHOverE(vIt)) continue;
    _elecHadFracCounter++;
  
    if(elecInBarrel(vIt) && fabs(elecOneOverEMinusOneOverPIn->at(vIt)) < _minElecOneOverEMinusOneOverPEB)
    _minElecOneOverEMinusOneOverPEB = elecOneOverEMinusOneOverPIn->at(vIt);
    if(elecInEndCap(vIt) && fabs(elecOneOverEMinusOneOverPIn->at(vIt)) < _minElecOneOverEMinusOneOverPEE)
    _minElecOneOverEMinusOneOverPEE = elecOneOverEMinusOneOverPIn->at(vIt);
    if(!passedElecOneOverEMinusOneOverP(vIt)) continue;
    _elecEOverPCounter++;
  
    if(fabs(elecDxyVtx->at(vIt)) < _minElecDxy) _minElecDxy = elecDxyVtx->at(vIt);
    if(!passedElecDxyVtx(vIt)) continue;
    _elecDxyCounter++;
  
    if(fabs(elecDzVtx->at(vIt)) < _minElecDz) _minElecDz = elecDzVtx->at(vIt);
    if(!passedElecDzVtx(vIt)) continue;
    _elecDzCounter++;
  
    if(elecMissingHits->at(vIt) < _minElecMissHits) _minElecMissHits = elecMissingHits->at(vIt);
    if(!passedElecMissingHits(vIt)) continue;
    _elecMissHitsCounter++;
    
    double effArea = 0.;
    if(_source == "Data") effArea = getEffAreaData(vIt);
    else effArea = getEffAreaMC(vIt);
 
    if((elecPFIsoDR03SumChargedHadronPt->at(vIt) 
    + std::max(0., elecPFIsoDR03SumNeutralHadronPt->at(vIt)
    + elecPFIsoDR03SumPhotonPt->at(vIt) - elecRhoIso * effArea))/elecPt->at(vIt) < _minElecPFRelIsoDR03)
    _minElecPFRelIsoDR03 = (elecPFIsoDR03SumChargedHadronPt->at(vIt) 
    + std::max(0., elecPFIsoDR03SumNeutralHadronPt->at(vIt)
    + elecPFIsoDR03SumPhotonPt->at(vIt) - elecRhoIso * effArea))/elecPt->at(vIt);
  
    if((elecPFIsoDR04SumChargedHadronPt->at(vIt) 
    + std::max(0., elecPFIsoDR04SumNeutralHadronPt->at(vIt)
    + elecPFIsoDR04SumPhotonPt->at(vIt) - elecRhoIso * effArea))/elecPt->at(vIt) < _minElecPFRelIsoDR04)
    _minElecPFRelIsoDR04 = (elecPFIsoDR04SumChargedHadronPt->at(vIt) 
    + std::max(0., elecPFIsoDR04SumNeutralHadronPt->at(vIt)
    + elecPFIsoDR04SumPhotonPt->at(vIt) - elecRhoIso * effArea))/elecPt->at(vIt);
    
    if(!passedElecPFRelIso(vIt)) continue;
    _elecPFRelIsoCounter++;
  
    if(passedElecId(vIt)) _elecIdCounter++;
    
    //if(hadTauIsGenElec->at(vIt) < 1) continue;
    
    if(elecHadTauDelatR->at(vIt) > _maxElecTauDR) _maxElecTauDR = elecHadTauDelatR->at(vIt);
    if(!passedElecTauDR(vIt)) continue;
    _elecTauDRCounter++;
    
    if(hadTauPt->at(vIt) > _maxTauPt) _maxTauPt = hadTauPt->at(vIt);
    if(hadTauAltLVPt->at(vIt) > _maxTauAltLVPt) _maxTauAltLVPt = hadTauAltLVPt->at(vIt);
    if(!passedHadTauMinPt(vIt)) continue;
    _hadTauPtCounter++;  
    
    if(fabs(hadTauEta->at(vIt)) < _minTauEta) _minTauEta = hadTauEta->at(vIt);
    if(!passedHadTauMaxEta(vIt)) continue;
    _hadTauEtaCounter++;
    
    if(passedHadTauAcc(vIt)) _hadTauAccCounter++;
    
    if(hadTauGenPt->at(vIt) > _maxTauGenPtBeforeDecayMode) _maxTauGenPtBeforeDecayMode = hadTauGenPt->at(vIt);
    if(hadTauPt->at(vIt) > _maxTauPtBeforeDecayMode) _maxTauPtBeforeDecayMode = hadTauPt->at(vIt);
    
    if((int)hadTauIdByDecayModeFinding->at(vIt) > _maxTauDecayMode) 
    _maxTauDecayMode = hadTauIdByDecayModeFinding->at(vIt);
    if(!passedHadTauDecayModeFinding(vIt)) continue;
    _hadTauDecayModeCounter++;
    if(hadTauGenPt->at(vIt) > _maxTauGenPtAfterDecayMode) _maxTauGenPtAfterDecayMode = hadTauGenPt->at(vIt);
    if(hadTauPt->at(vIt) > _maxTauPtAfterDecayMode) _maxTauPtAfterDecayMode = hadTauPt->at(vIt);
    
    if((int) hadTauNProngs->at(vIt) > _maxTauNProngs) _maxTauNProngs = hadTauNProngs->at(vIt);
    if(!passedTauNProngs(vIt)) continue;
    _hadTauNProngsCounter++;    
    
    if(hadTauPt->at(vIt) > _maxTauPtBeforeIso) _maxTauPtBeforeIso = hadTauPt->at(vIt);
    if(hadTauGenPt->at(vIt) > _maxTauGenPtBeforeIso) _maxTauGenPtBeforeIso = hadTauGenPt->at(vIt);
    
    if(hadTauIdByLooseCombinedIsolationDBSumPtCorr3Hits->at(vIt) < 0.5) _tauIsoDB = 0;
    if(hadTauIdByLooseCombinedIsolationDBSumPtCorr3Hits->at(vIt) > 0) _tauIsoDB = 1;
    if(hadTauIdByMediumCombinedIsolationDBSumPtCorr3Hits->at(vIt) > 0) _tauIsoDB = 2;
    if(hadTauIdByTightCombinedIsolationDBSumPtCorr3Hits->at(vIt) > 0) _tauIsoDB = 3;
    
    if(hadTauIdByVLooseIsolationMVA3newDMwLT->at(vIt) < 0.5) _tauMVAIso = 0;
    if(hadTauIdByVLooseIsolationMVA3newDMwLT->at(vIt) > 0) _tauMVAIso = 1;
    if(hadTauIdByLooseIsolationMVA3newDMwLT->at(vIt) > 0) _tauMVAIso = 2;
    if(hadTauIdByMediumIsolationMVA3newDMwLT->at(vIt) > 0) _tauMVAIso = 3;
    if(hadTauIdByTightIsolationMVA3newDMwLT->at(vIt) > 0) _tauMVAIso = 4;
    if(hadTauIdByVTightIsolationMVA3newDMwLT->at(vIt) > 0) _tauMVAIso = 5;
    if(passedHadTauIsoDisc(vIt)) _hadTauTightIsoCounter++;
    if(!_passedTauIso) continue;
    _hadTauIsoDiscCounter++;
    
    if(hadTauPt->at(vIt) > _maxTauPtAfterIso) _maxTauPtAfterIso = hadTauPt->at(vIt);
    if(hadTauGenPt->at(vIt) > _maxTauGenPtAfterIso) _maxTauGenPtAfterIso = hadTauGenPt->at(vIt);
    
    if(hadTauIdByLooseElectronRejection->at(vIt) < 0.5) _tauElecDisc = 0;
    if(hadTauIdByLooseElectronRejection->at(vIt) > 0) _tauElecDisc = 1;
    if(hadTauIdByMediumElectronRejection->at(vIt) > 0) _tauElecDisc = 2;
    if(hadTauIdByTightElectronRejection->at(vIt) > 0) _tauElecDisc = 3;
    if(hadTauIdByMVA5VLooseElectronRejection->at(vIt) < 0.5) _tauMVA5ElecDisc = 0;
    if(hadTauIdByMVA5VLooseElectronRejection->at(vIt) > 0)  _tauMVA5ElecDisc = 1;
    if(hadTauIdByMVA5LooseElectronRejection->at(vIt) > 0)  _tauMVA5ElecDisc = 2;
    if(hadTauIdByMVA5MediumElectronRejection->at(vIt) > 0)  _tauMVA5ElecDisc = 3;
    if(hadTauIdByMVA5TightElectronRejection->at(vIt) > 0)  _tauMVA5ElecDisc = 4;
    if(hadTauIdByMVA5VTightElectronRejection->at(vIt) > 0)  _tauMVA5ElecDisc = 5;
    if(!passedHadTauElecDisc(vIt)) continue;
    _hadTauElecDiscCounter++;
    if(hadTauGenPt->at(vIt) > _maxTauGenPtAfterElecDisc) _maxTauGenPtAfterElecDisc = hadTauGenPt->at(vIt);
    if(hadTauPt->at(vIt) > _maxTauPtAfterElecDisc) _maxTauPtAfterElecDisc = hadTauPt->at(vIt);
    
    if(hadTauIdByLooseMuonRejection3->at(vIt) < 0.5) _tauMuonDisc = 0;
    if(hadTauIdByLooseMuonRejection3->at(vIt) > 0) _tauMuonDisc = 1;
    if(hadTauIdByTightMuonRejection3->at(vIt) > 0) _tauMuonDisc = 2;
    if(!passedHadTauMuonDisc(vIt)) continue;
    _hadTauMuonDiscCounter++;
    if(hadTauGenPt->at(vIt) > _maxTauGenPtAfterMuonDisc) _maxTauGenPtAfterMuonDisc = hadTauGenPt->at(vIt);
    if(hadTauPt->at(vIt) > _maxTauPtAfterMuonDisc) _maxTauPtAfterMuonDisc = hadTauPt->at(vIt);
    if(passedHadTauId(vIt)) _hadTauIdCounter++;
    
    //if(!passedZEEVeto()) continue;
    //_zeeVetoCounter++;
    
    if((int)(elecCharge->at(vIt)*hadTauCharge->at(vIt)) < _minElecTauCharge) _minElecTauCharge =
    (int)(elecCharge->at(vIt)*hadTauCharge->at(vIt));
    if(passedElecTauCharge(vIt)) _tightDiTauChargeCounter++;
    if(!_passedDiTauChargeCut) continue;
    _elecHadTauChargeCounter++;
   
    if(elecHadTauCosDPhi->at(vIt) < _minElecTauCosDPhi) _minElecTauCosDPhi = elecHadTauCosDPhi->at(vIt);
    if(!passedElecTauCosDPhi(vIt)) continue;
    _elecHadTauCosDPhiCounter++;
    
    _maxMET = MEt;
    if(!passedMET()) continue;
    _metCounter++;
    
    if(elecHadTauPZeta->at(vIt) > _maxPZeta) _maxPZeta = elecHadTauPZeta->at(vIt);
    if(elecHadTauPZetaVis->at(vIt) < _minPZetaVis) _minPZetaVis = elecHadTauPZetaVis->at(vIt);
    if(!passed2DZetaCut(vIt)) continue;
    _elecHadTauPZetaCounter++;
    
    if(elecMEtMt->at(vIt) > _maxElecMETMass) _maxElecMETMass = elecMEtMt->at(vIt);
    if(elecMetCosDPhi->at(vIt) < _minElecMETCosDPhi) _minElecMETCosDPhi = elecMetCosDPhi->at(vIt);
    if(!passedWJetsVeto(vIt)) continue;
    _wjetsVetoCounter++;
    
    if(nBtagsHiEffTrkCnt->at(vIt) < _minBtagsTCHEM)_minBtagsTCHEM = nBtagsHiEffTrkCnt->at(vIt);
    if(nBtagsHiPurityTrkCnt->at(vIt) < _minBtagsTCHPT) _minBtagsTCHPT = nBtagsHiPurityTrkCnt->at(vIt);
    if(nBTagsCombSecVtxLoose->at(vIt) < _minBtagsCSVL) _minBtagsCSVL = nBTagsCombSecVtxLoose->at(vIt);
    if(nBTagsCombSecVtxMedium->at(vIt) < _minBtagsCSVM) _minBtagsCSVM = nBTagsCombSecVtxMedium->at(vIt);
    if(nBTagsCombSecVtxTight->at(vIt) < _minBtagsCSVT) _minBtagsCSVT = nBTagsCombSecVtxTight->at(vIt);
    if(nJets->at(vIt) < _minNJets) _minNJets = nJets->at(vIt);
    if(jetEtSum->at(vIt) > _maxJetEtSum) _maxJetEtSum = jetEtSum->at(vIt);
    if(!passedBTagVeto(vIt)) continue;
    _bTagsCounter++;
    
    if(hadTauPt->at(vIt) > _maxTauPtAfterAllCuts) _maxTauPtAfterAllCuts = hadTauPt->at(vIt);
    if(elecPt->at(vIt) > _maxElecPtAfterAllCuts) _maxElecPtAfterAllCuts = elecPt->at(vIt);
    _maxMETAfterAllCuts = MEt;
    
    if(elecHadTauMass->at(vIt) > _maxElecTauMass) _maxElecTauMass = elecHadTauMass->at(vIt);
    if(elecHadTauMetMass->at(vIt) > _maxElecTauMetMass) _maxElecTauMetMass = elecHadTauMetMass->at(vIt);
    if(!passedElecTauMass(vIt)) continue;
    _elecTauMassCounter++;
    
    if((int)nVBFJets->at(vIt) > _maxVBFJets) _maxVBFJets = (int)nVBFJets->at(vIt);
    if(passedTopologyCuts(vIt)) _topologyCounter++;
    
    if(!passedNVBFJets(vIt)) continue;
    _nVBFJetsCounter++;
    
    if(!passedVBFJetsPt(vIt)) continue;
    _vbfJetsPtCounter++;
    
    if(!passedDiJetDEta(vIt)) continue;
    _diJetDEtaCounter++;
    
    if(diJetMass->at(vIt) > _maxDiJetMass) _maxDiJetMass = diJetMass->at(vIt);
    if(!passedDiJetMass(vIt)) continue;
    _diJetMassCounter++;
    
    if(passedVBFSel(vIt)) _vbfCounter++;
  }
}

void ElecTauAnalysis::getNMinus1CandCounters(unsigned int theIndex){
  if(passedElecEtaCut(theIndex) && passedElecId(theIndex) && passedHadTauAcc(theIndex)
  && passedHadTauId(theIndex) && passedTopologyCuts(theIndex) && passedVBFSel(theIndex)){
    _nMinus1ElecPtCounter++;
    if(elecPt->at(theIndex) > _maxElecPtNMinus1) _maxElecPtNMinus1 = elecPt->at(theIndex);
  }
  
  if(passedElecPtCut(theIndex) && passedElecId(theIndex) && passedHadTauAcc(theIndex)
  && passedHadTauId(theIndex) && passedTopologyCuts(theIndex) && passedVBFSel(theIndex)){
    _nMinus1ElecEtaCounter++;
    if(!elecInGaps(theIndex) && fabs(elecEta->at(theIndex)) < _minElecEtaNMinus1) _minElecEtaNMinus1 = elecEta->at(theIndex);
  }
  
  if(passedElecAccCuts(theIndex) && passedElecDEtaIn(theIndex) && passedElecSigmaIEtaIEta(theIndex) && passedElecHOverE(theIndex)
  && passedElecOneOverEMinusOneOverP(theIndex) && passedElecDxyVtx(theIndex) && passedElecDzVtx(theIndex)
  && passedElecMissingHits(theIndex) && passedElecPFRelIso(theIndex) && passedHadTauAcc(theIndex) && passedHadTauId(theIndex) 
  && passedTopologyCuts(theIndex) && passedVBFSel(theIndex)){
    _nMinus1ElecDPhiInCounter++;
    if(elecInBarrel(theIndex) && fabs(elecDPhiIn->at(theIndex)) < _minElecDPhiEBNMinus1) 
    _minElecDPhiEBNMinus1 = elecDPhiIn->at(theIndex);
    if(elecInEndCap(theIndex) && fabs(elecDPhiIn->at(theIndex)) < _minElecDPhiEENMinus1) 
    _minElecDPhiEENMinus1 = elecDPhiIn->at(theIndex);
  }
  
  if(passedElecAccCuts(theIndex) && passedElecDPhiIn(theIndex) && passedElecSigmaIEtaIEta(theIndex) && passedElecHOverE(theIndex)
  && passedElecOneOverEMinusOneOverP(theIndex) && passedElecDxyVtx(theIndex) && passedElecDzVtx(theIndex)
  && passedElecMissingHits(theIndex) && passedElecPFRelIso(theIndex) && passedHadTauAcc(theIndex) && passedHadTauId(theIndex) 
  && passedTopologyCuts(theIndex) && passedVBFSel(theIndex)){
    _nMinus1ElecDEtaInCounter++;
    if(elecInBarrel(theIndex) && fabs(elecDEtaIn->at(theIndex)) < _minElecDEtaEBNMinus1) 
    _minElecDEtaEBNMinus1 = elecDEtaIn->at(theIndex);
    if(elecInEndCap(theIndex) && fabs(elecDEtaIn->at(theIndex)) < _minElecDEtaEENMinus1) 
    _minElecDEtaEENMinus1 = elecDEtaIn->at(theIndex);
  }
  
  if(passedElecAccCuts(theIndex) && passedElecDPhiIn(theIndex) && passedElecDEtaIn(theIndex) && passedElecHOverE(theIndex)
  && passedElecOneOverEMinusOneOverP(theIndex) && passedElecDxyVtx(theIndex) && passedElecDzVtx(theIndex)
  && passedElecMissingHits(theIndex) && passedElecPFRelIso(theIndex) && passedHadTauAcc(theIndex) && passedHadTauId(theIndex) 
  && passedTopologyCuts(theIndex) && passedVBFSel(theIndex)){
    _nMinus1ElecSigmaIEtaIEtaCounter++;
    if(elecInBarrel(theIndex) && elecSigmaIEtaIEta->at(theIndex) < _minElecSigmaIEtaIEtaEBNMinus1) 
    _minElecSigmaIEtaIEtaEBNMinus1 = elecSigmaIEtaIEta->at(theIndex); 
    if(elecInEndCap(theIndex) && elecSigmaIEtaIEta->at(theIndex) < _minElecSigmaIEtaIEtaEENMinus1) 
    _minElecSigmaIEtaIEtaEENMinus1 = elecSigmaIEtaIEta->at(theIndex); 
  }
  
  if(passedElecAccCuts(theIndex) && passedElecDPhiIn(theIndex) && passedElecDEtaIn(theIndex) && passedElecSigmaIEtaIEta(theIndex)
  && passedElecOneOverEMinusOneOverP(theIndex) && passedElecDxyVtx(theIndex) && passedElecDzVtx(theIndex)
  && passedElecMissingHits(theIndex) && passedElecPFRelIso(theIndex) && passedHadTauAcc(theIndex) && passedHadTauId(theIndex) 
  && passedTopologyCuts(theIndex) && passedVBFSel(theIndex)){
    _nMinus1ElecHOverECounter++;
    if(elecInBarrel(theIndex) && elecHadOverE->at(theIndex) < _minElecHOverEEBNMinus1) 
    _minElecHOverEEBNMinus1 = elecHadOverE->at(theIndex);
    if(elecInEndCap(theIndex) && elecHadOverE->at(theIndex) < _minElecHOverEEENMinus1) 
    _minElecHOverEEENMinus1 = elecHadOverE->at(theIndex);
  }
  
  if(passedElecAccCuts(theIndex) && passedElecDPhiIn(theIndex) && passedElecDEtaIn(theIndex) && passedElecSigmaIEtaIEta(theIndex)
  && passedElecHOverE(theIndex) && passedElecDxyVtx(theIndex) && passedElecDzVtx(theIndex)
  && passedElecMissingHits(theIndex) && passedElecPFRelIso(theIndex) && passedHadTauAcc(theIndex) && passedHadTauId(theIndex) 
  && passedTopologyCuts(theIndex) && passedVBFSel(theIndex)){
    _nMinus1OneOverEMinusOneOverPCounter++;
    if(elecInBarrel(theIndex) && fabs(elecOneOverEMinusOneOverPIn->at(theIndex)) < _minElecOneOverEMinusOneOverPEBNMinus1)
    _minElecOneOverEMinusOneOverPEBNMinus1 = elecOneOverEMinusOneOverPIn->at(theIndex);
    if(elecInEndCap(theIndex) && fabs(elecOneOverEMinusOneOverPIn->at(theIndex)) < _minElecOneOverEMinusOneOverPEENMinus1)
    _minElecOneOverEMinusOneOverPEENMinus1 = elecOneOverEMinusOneOverPIn->at(theIndex);
  }
  
  if(passedElecAccCuts(theIndex) && passedElecDPhiIn(theIndex) && passedElecDEtaIn(theIndex) && passedElecSigmaIEtaIEta(theIndex)
  && passedElecHOverE(theIndex) && passedElecOneOverEMinusOneOverP(theIndex) && passedElecDzVtx(theIndex)
  && passedElecMissingHits(theIndex) && passedElecPFRelIso(theIndex) && passedHadTauAcc(theIndex) && passedHadTauId(theIndex) 
  && passedTopologyCuts(theIndex) && passedVBFSel(theIndex)){
    _nMinus1ElecDxyVtxCounter++;
    if(fabs(elecDxyVtx->at(theIndex)) < _minElecDxyNMinus1) _minElecDxyNMinus1 = elecDxyVtx->at(theIndex);
  } 
  
  if(passedElecAccCuts(theIndex) && passedElecDPhiIn(theIndex) && passedElecDEtaIn(theIndex) && passedElecSigmaIEtaIEta(theIndex)
  && passedElecHOverE(theIndex) && passedElecOneOverEMinusOneOverP(theIndex) && passedElecDxyVtx(theIndex)
  && passedElecMissingHits(theIndex) && passedElecPFRelIso(theIndex) && passedHadTauAcc(theIndex) && passedHadTauId(theIndex) 
  && passedTopologyCuts(theIndex) && passedVBFSel(theIndex)){
    _nMinus1ElecDzVtxCounter++;
    if(fabs(elecDzVtx->at(theIndex)) < _minElecDzNMinus1) _minElecDzNMinus1 = elecDzVtx->at(theIndex);
  }
  
  if(passedElecAccCuts(theIndex) && passedElecDPhiIn(theIndex) && passedElecDEtaIn(theIndex) && passedElecSigmaIEtaIEta(theIndex)
  && passedElecHOverE(theIndex) && passedElecOneOverEMinusOneOverP(theIndex) && passedElecDxyVtx(theIndex)
  && passedElecDzVtx(theIndex) && passedElecPFRelIso(theIndex) && passedHadTauAcc(theIndex) && passedHadTauId(theIndex) 
  && passedTopologyCuts(theIndex) && passedVBFSel(theIndex)){
    _nMinus1ElecMissHitsCounter++;
    if(elecMissingHits->at(theIndex) < _minElecMissHitsNMinus1) _minElecMissHitsNMinus1 = elecMissingHits->at(theIndex);
  }
  
  if(passedElecAccCuts(theIndex) && passedElecDPhiIn(theIndex) && passedElecDEtaIn(theIndex) && passedElecSigmaIEtaIEta(theIndex)
  && passedElecHOverE(theIndex) && passedElecOneOverEMinusOneOverP(theIndex) && passedElecDxyVtx(theIndex)
  && passedElecDzVtx(theIndex) && passedElecMissingHits(theIndex) && passedHadTauAcc(theIndex) && passedHadTauId(theIndex) 
  && passedTopologyCuts(theIndex) && passedVBFSel(theIndex)){
    double effArea = 0.;
    if(_source == "Data") effArea = getEffAreaData(theIndex);
    else effArea = getEffAreaMC(theIndex);
    _nMinus1ElecPFRelIsoCounter++;
    
    if((elecPFIsoDR03SumChargedHadronPt->at(theIndex) 
    + std::max(0., elecPFIsoDR03SumNeutralHadronPt->at(theIndex)
    + elecPFIsoDR03SumPhotonPt->at(theIndex) - elecRhoIso * effArea))/elecPt->at(theIndex) < _minElecPFRelIsoDR03NMinus1)
    _minElecPFRelIsoDR03NMinus1 = (elecPFIsoDR03SumChargedHadronPt->at(theIndex) 
    + std::max(0., elecPFIsoDR03SumNeutralHadronPt->at(theIndex)
    + elecPFIsoDR03SumPhotonPt->at(theIndex) - elecRhoIso * effArea))/elecPt->at(theIndex);
  
    if((elecPFIsoDR04SumChargedHadronPt->at(theIndex) 
    + std::max(0., elecPFIsoDR04SumNeutralHadronPt->at(theIndex)
    + elecPFIsoDR04SumPhotonPt->at(theIndex) - elecRhoIso * effArea))/elecPt->at(theIndex) < _minElecPFRelIsoDR04NMinus1)
    _minElecPFRelIsoDR04NMinus1 = (elecPFIsoDR04SumChargedHadronPt->at(theIndex) 
    + std::max(0., elecPFIsoDR04SumNeutralHadronPt->at(theIndex)
    + elecPFIsoDR04SumPhotonPt->at(theIndex) - elecRhoIso * effArea))/elecPt->at(theIndex);

  }
  
  if(passedElecAccCuts(theIndex) && passedElecId(theIndex) && passedHadTauMinPt(theIndex) && passedHadTauMaxEta(theIndex) 
  && passedHadTauId(theIndex) && passedTopologyCuts(theIndex) && passedVBFSel(theIndex)){
    _nMinus1ElecHadTauDRCounter++;
    if(elecHadTauDelatR->at(theIndex) > _maxElecTauDRNMinus1) _maxElecTauDRNMinus1 = elecHadTauDelatR->at(theIndex);
  }
  
  if(passedElecAccCuts(theIndex) && passedElecId(theIndex) && passedElecTauDR(theIndex) && passedHadTauMaxEta(theIndex) 
  && passedHadTauId(theIndex) && passedTopologyCuts(theIndex) && passedVBFSel(theIndex)){
    _nMinus1TauPtCounter++;
    if(hadTauPt->at(theIndex) > _maxTauPtNMinus1) _maxTauPtNMinus1 = hadTauPt->at(theIndex);
    if(hadTauAltLVPt->at(theIndex) > _maxTauAltLVPtNMinus1) _maxTauAltLVPtNMinus1 = hadTauAltLVPt->at(theIndex);
  }
  
  if(passedElecAccCuts(theIndex) && passedElecId(theIndex) && passedElecTauDR(theIndex) && passedHadTauMinPt(theIndex) 
  && passedHadTauId(theIndex) && passedTopologyCuts(theIndex) && passedVBFSel(theIndex)){
    _nMinus1TauEtaCounter++;
    if(fabs(hadTauEta->at(theIndex)) < _minTauEtaNMinus1) _minTauEtaNMinus1 = hadTauEta->at(theIndex);
  }
  
  if(passedElecAccCuts(theIndex) && passedElecId(theIndex) && passedHadTauAcc(theIndex) && passedTauNProngs(theIndex)
  && passedHadTauIsoDisc(theIndex) && passedHadTauMuonDisc(theIndex) && passedHadTauElecDisc(theIndex) 
  && passedTopologyCuts(theIndex) && passedVBFSel(theIndex)){
    _nMinus1TauDecayModeCounter++;
    if((int)hadTauIdByDecayModeFinding->at(theIndex) > _maxTauDecayModeNMinus1) 
    _maxTauDecayModeNMinus1 = hadTauIdByDecayModeFinding->at(theIndex);
  }
  
  if(passedElecAccCuts(theIndex) && passedElecId(theIndex) && passedHadTauAcc(theIndex) 
  && passedHadTauDecayModeFinding(theIndex)
  && passedHadTauIsoDisc(theIndex) && passedHadTauMuonDisc(theIndex) && passedHadTauElecDisc(theIndex) 
  && passedTopologyCuts(theIndex) && passedVBFSel(theIndex)){
    _nMinus1TauNProngsCounter++;
    if((int)hadTauNProngs->at(theIndex) > _maxTauNProngsNMinus1) _maxTauNProngsNMinus1 = hadTauNProngs->at(theIndex);
  }
  
  if(passedElecAccCuts(theIndex) && passedElecId(theIndex) && passedHadTauAcc(theIndex)
  && passedHadTauDecayModeFinding(theIndex) && passedTauNProngs(theIndex)
  && passedHadTauIsoDisc(theIndex) && passedHadTauElecDisc(theIndex) 
  && passedTopologyCuts(theIndex) && passedVBFSel(theIndex)){
    _nMinus1TauMuonDiscCounter++;
    if(hadTauIdByLooseMuonRejection3->at(theIndex) < 0.5) _tauMuonDiscNMinus1 = 0;
    if(hadTauIdByLooseMuonRejection3->at(theIndex) > 0) _tauMuonDiscNMinus1 = 1;
    if(hadTauIdByTightMuonRejection3->at(theIndex) > 0) _tauMuonDiscNMinus1 = 2;
  }
  
  if(passedElecAccCuts(theIndex) && passedElecId(theIndex) && passedHadTauAcc(theIndex)
  && passedHadTauDecayModeFinding(theIndex) && passedTauNProngs(theIndex) 
  && passedHadTauIsoDisc(theIndex) && passedHadTauMuonDisc(theIndex) 
  && passedTopologyCuts(theIndex) && passedVBFSel(theIndex)){
    _nMinus1TauElecDiscCounter++;
    if(hadTauIdByLooseElectronRejection->at(theIndex) < 0.5) _tauElecDiscNMinus1 = 0;
    if(hadTauIdByLooseElectronRejection->at(theIndex) > 0) _tauElecDiscNMinus1 = 1;
    if(hadTauIdByMediumElectronRejection->at(theIndex) > 0) _tauElecDiscNMinus1 = 2;
    if(hadTauIdByTightElectronRejection->at(theIndex) > 0) _tauElecDiscNMinus1 = 3;
    if(hadTauIdByMVA5VLooseElectronRejection->at(theIndex) < 0.5) _tauMVA5ElecDiscNMinus1 = 0;
    if(hadTauIdByMVA5VLooseElectronRejection->at(theIndex) > 0)  _tauMVA5ElecDiscNMinus1 = 1;
    if(hadTauIdByMVA5LooseElectronRejection->at(theIndex) > 0)  _tauMVA5ElecDiscNMinus1 = 2;
    if(hadTauIdByMVA5MediumElectronRejection->at(theIndex) > 0)  _tauMVA5ElecDiscNMinus1 = 3;
    if(hadTauIdByMVA5TightElectronRejection->at(theIndex) > 0)  _tauMVA5ElecDiscNMinus1 = 4;
    if(hadTauIdByMVA5VTightElectronRejection->at(theIndex) > 0)  _tauMVA5ElecDiscNMinus1 = 5;
  }
  
  if(passedElecAccCuts(theIndex) && passedElecId(theIndex) && passedHadTauAcc(theIndex)
  && passedHadTauDecayModeFinding(theIndex) && passedTauNProngs(theIndex) 
  && passedHadTauMuonDisc(theIndex) && passedHadTauElecDisc(theIndex) 
  && passedTopologyCuts(theIndex) && passedVBFSel(theIndex)){
    _nMinus1TauIsoCounter++;
    
    if(hadTauIdByVLooseIsolationMVA3newDMwLT->at(theIndex) < 0.5) _tauMVAIsoNMinus1 = 0;
    if(hadTauIdByVLooseIsolationMVA3newDMwLT->at(theIndex) > 0) _tauMVAIsoNMinus1 = 1;
    if(hadTauIdByLooseIsolationMVA3newDMwLT->at(theIndex) > 0) _tauMVAIsoNMinus1 = 2;
    if(hadTauIdByMediumIsolationMVA3newDMwLT->at(theIndex) > 0) _tauMVAIsoNMinus1 = 3;
    if(hadTauIdByTightIsolationMVA3newDMwLT->at(theIndex) > 0) _tauMVAIsoNMinus1 = 4;
    if(hadTauIdByVTightIsolationMVA3newDMwLT->at(theIndex) > 0) _tauMVAIsoNMinus1 = 5;
    
    if(hadTauIdByLooseCombinedIsolationDBSumPtCorr3Hits->at(theIndex) < 0.5) _tauIsoDB3HitsNMinus1 = 0;
    if(hadTauIdByLooseCombinedIsolationDBSumPtCorr3Hits->at(theIndex) > 0.5) _tauIsoDB3HitsNMinus1 = 1;
    if(hadTauIdByMediumCombinedIsolationDBSumPtCorr3Hits->at(theIndex) > 0.5) _tauIsoDB3HitsNMinus1 = 2;
    if(hadTauIdByTightCombinedIsolationDBSumPtCorr3Hits->at(theIndex) > 0.5) _tauIsoDB3HitsNMinus1 = 3;
  }
  
  if(passedElecAccCuts(theIndex) && passedElecId(theIndex) 
  && passedHadTauAcc(theIndex) && passedHadTauId(theIndex)
  && passedElecTauCosDPhi(theIndex) && passedMET()
  && passed2DZetaCut(theIndex) && passedBTagVeto(theIndex) && passedWJetsVeto(theIndex) && passedVBFSel(theIndex)){
    _nMinus1ElecTauChargeCounter++;
    if(elecCharge->at(theIndex)*hadTauCharge->at(theIndex) < _minElecTauChargeNMinus1) _minElecTauChargeNMinus1 =
    elecCharge->at(theIndex)*hadTauCharge->at(theIndex);
  }
  
  if(passedElecAccCuts(theIndex) && passedElecId(theIndex) 
  && passedHadTauAcc(theIndex) && passedHadTauId(theIndex)
  && passedElecTauCharge(theIndex) && passedMET()
  && passed2DZetaCut(theIndex) && passedBTagVeto(theIndex) && passedWJetsVeto(theIndex) && passedVBFSel(theIndex)){
    _nMinus1ElecTauCosDPhiCounter++;
    if(elecHadTauCosDPhi->at(theIndex) < _minElecTauCosDPhiNMinus1) _minElecTauCosDPhiNMinus1 =
    elecHadTauCosDPhi->at(theIndex);
  }
  
  if(passedElecAccCuts(theIndex) && passedElecId(theIndex) 
  && passedHadTauAcc(theIndex) && passedHadTauId(theIndex)
  && passedElecTauCharge(theIndex) && passedElecTauCosDPhi(theIndex)
  && passed2DZetaCut(theIndex) && passedBTagVeto(theIndex) && passedWJetsVeto(theIndex) && passedVBFSel(theIndex)){
    _nMinus1METCounter++;
    if(MEt > _maxMETNMinus1) _maxMETNMinus1 = MEt;
  }
  
  if(passedElecAccCuts(theIndex) && passedElecId(theIndex) 
  && passedHadTauAcc(theIndex) && passedHadTauId(theIndex)
  && passedElecTauCharge(theIndex) && passedElecTauCosDPhi(theIndex)
  && passedMET() && passedBTagVeto(theIndex) && passedWJetsVeto(theIndex) && passedVBFSel(theIndex)){
    _nMinus1PZeta2DCounter++;
    if(elecHadTauPZeta->at(theIndex) > _maxPZetaNMinus1) _maxPZetaNMinus1 = elecHadTauPZeta->at(theIndex);
    if(elecHadTauPZetaVis->at(theIndex) < _minPZetaVisNMinus1) _minPZetaVisNMinus1 = elecHadTauPZetaVis->at(theIndex);
  }
  
  bool passedBtagging = true;
  if(_btaggingType == "ReversedJets" || _btaggingType == "JetSumEt") passedBtagging = passedBTagVeto(theIndex);
  if(passedElecAccCuts(theIndex) && passedElecId(theIndex) 
  && passedHadTauAcc(theIndex) && passedHadTauId(theIndex)
  && passedElecTauCharge(theIndex) && passedElecTauCosDPhi(theIndex) && passedBtagging
  && passedMET() && passed2DZetaCut(theIndex) && passedWJetsVeto(theIndex) && passedVBFSel(theIndex)){
    _nMinus1NBtagsCounter++;
    if((int)nBtagsHiEffTrkCnt->at(theIndex) < _minBtagsTCHEMNMinus1) _minBtagsTCHEMNMinus1 = nBtagsHiEffTrkCnt->at(theIndex);
    if((int)nBtagsHiPurityTrkCnt->at(theIndex) < _minBtagsTCHPTNMinus1) _minBtagsTCHPTNMinus1 = nBtagsHiPurityTrkCnt->at(theIndex);
    if((int)nBTagsCombSecVtxLoose->at(theIndex) < _minBtagsCSVLNMinus1) _minBtagsCSVLNMinus1 = nBTagsCombSecVtxLoose->at(theIndex);
    if((int)nBTagsCombSecVtxMedium->at(theIndex) < _minBtagsCSVMNMinus1) _minBtagsCSVMNMinus1 = nBTagsCombSecVtxMedium->at(theIndex);
    if((int)nBTagsCombSecVtxTight->at(theIndex) < _minBtagsCSVTNMinus1) _minBtagsCSVTNMinus1 = nBTagsCombSecVtxTight->at(theIndex);
    if((int)nJets->at(theIndex) < _minNJetsNminus1) _minNJetsNminus1 = nJets->at(theIndex);
    if(jetEtSum->at(theIndex) > _maxJetEtSumNMinus1) _maxJetEtSumNMinus1 = jetEtSum->at(theIndex);
  }
  if(passedElecAccCuts(theIndex) && passedElecId(theIndex) 
  && passedHadTauAcc(theIndex) && passedHadTauId(theIndex)
  && passedElecTauCharge(theIndex) && passedElecTauCosDPhi(theIndex)
  && passedMET() && passed2DZetaCut(theIndex) && passedWJetsVeto(theIndex) && passedVBFSel(theIndex)){
    if((int)nJets->at(theIndex) < _minNJetsNminus1) _minNJetsNminus1 = nJets->at(theIndex);
    if(jetEtSum->at(theIndex) > _maxJetEtSumNMinus1) _maxJetEtSumNMinus1 = jetEtSum->at(theIndex);
  }
  
  if(passedElecAccCuts(theIndex) && passedElecId(theIndex) 
  && passedHadTauAcc(theIndex) && passedHadTauId(theIndex) && passedTopologyCuts(theIndex) 
  && passedVBFJetsPt(theIndex) && passedDiJetDEta(theIndex) && passedDiJetMass(theIndex)){
    _nMinus1NVBFJetsCounter++;
    if((int)nVBFJets->at(theIndex) > _maxNVBFJetsNMinus1)_maxNVBFJetsNMinus1 = nVBFJets->at(theIndex);
  } 
  
  if(passedElecAccCuts(theIndex) && passedElecId(theIndex) 
  && passedHadTauAcc(theIndex) && passedHadTauId(theIndex) && passedTopologyCuts(theIndex) 
  && passedNVBFJets(theIndex) && passedDiJetDEta(theIndex) && passedDiJetMass(theIndex)){
    _nMinus1JetsPtCounter++;
    if(jet1Pt->at(theIndex) > _maxJet1PtNMinus1) _maxJet1PtNMinus1 = jet1Pt->at(theIndex);
    if(jet2Pt->at(theIndex) > _maxJet2PtNMinus1) _maxJet2PtNMinus1 = jet2Pt->at(theIndex);
  } 
  
  if(passedElecAccCuts(theIndex) && passedElecId(theIndex) 
  && passedHadTauAcc(theIndex) && passedHadTauId(theIndex) && passedTopologyCuts(theIndex) 
  && passedNVBFJets(theIndex) && passedVBFJetsPt(theIndex) && passedDiJetMass(theIndex)){
    _nMinus1DiJetDEtaCounter++;
    if(diJetDEta->at(theIndex) > _maxDiJetDEtaNMinus1) _maxDiJetDEtaNMinus1 = diJetDEta->at(theIndex);
  } 
  
  if(passedElecAccCuts(theIndex) && passedElecId(theIndex) 
  && passedHadTauAcc(theIndex) && passedHadTauId(theIndex) && passedTopologyCuts(theIndex) 
  && passedNVBFJets(theIndex) && passedVBFJetsPt(theIndex) && passedDiJetDEta(theIndex)){
    _nMinus1DiJetMassCounter++;
    if(diJetMass->at(theIndex) > _maxDiJetMassNMinus1) _maxDiJetMassNMinus1 = diJetMass->at(theIndex);
  } 
  
}

void ElecTauAnalysis::getReport(){
  //start xml file with root element
  
  std::ofstream outLogFile;
  outLogFile.open(_outLogFile.c_str());    
  outLogFile.setf(ios::fixed,ios::floatfield);
  
  outLogFile
  << "<?xml version=\"1.0\" encoding=\"ISO-8859-1\"?>\n"
  << "<?xml-stylesheet type=\"text/xsl\" href=\"elecTauLog.xsl\"?>\n"
  << "<etauAnalysis>\n";
  
  double elecPtDeno = _puWeightedEvents;
  double nMatchedElec = _nEvents;
  double nMatchedTau = _nEvents;
  string matchingEff = PrintEfficiency(_puWeightedEvents,_puWeightedEvents);
  if(_source == "Signal"){
    nMatchedElec = _nElecMatched;
    nMatchedTau = _nTauMatched;
    elecPtDeno = nMatchedTau;
    _nTriggers = elecPtDeno; //Don't apply trigger to signal => get this eff to 100%
    outLogFile
    << PrintXMLFormat("cutEff", "Elec Match", nMatchedElec, _puWeightedEvents, true)
    << PrintXMLFormat("cutEff", "Tau Match", nMatchedTau, nMatchedElec, true)
    << "\n";
    /*
    outLogFile << "\t Matching Efficiencies \n"
    << " Cut \t \t \t \t"		<< "NEvents \t" 			<< "Efficiency "
    << "\n -------------------------------------------------------------"
    << "\n All \t \t \t \t"	   	<< setprecision(0) << _puWeightedEvents	   
    << "\n Elec Matching \t \t \t" 	<< setprecision(0) << nMatchedElec << "\t \t"	   << PrintEfficiency(nMatchedElec, _puWeightedEvents)
    << "\n Tau Matching \t \t \t" 	<< setprecision(0) << nMatchedTau << "\t \t"	   << PrintEfficiency(nMatchedTau, nMatchedElec)
    << "\n"
    << "\n";
    */
    matchingEff = PrintEfficiency(nMatchedTau, _puWeightedEvents);
  }
  
  // correction to raw numbers due to factorization
  std::pair<double, std::pair<double, double> > isoLooseToTightCorr = GetEfficiency(_nTightTauIso,_nHadTauIsoDisc);
  double isoCorr = isoLooseToTightCorr.first;
  
  std::pair<double, std::pair<double, double> > chargeLooseToTightCorr = GetEfficiency(_nTightDiTauCharge, _nElecHadTauCharge);
  double chargeCorr = chargeLooseToTightCorr.first; 
  
  std::pair<double, double> fullCorrEff = GetEfficiencyProd(
  	std::make_pair<double, double>(isoLooseToTightCorr.first, std::max(isoLooseToTightCorr.second.first,isoLooseToTightCorr.second.second)), 
  	std::make_pair<double, double>(chargeLooseToTightCorr.first, std::max(chargeLooseToTightCorr.second.first,chargeLooseToTightCorr.second.second))
  );
  
  pair<double,pair<double,double> > efficiency = GetEfficiency(_nVBFSel, _puWeightedEvents);
  std::pair<double, double> correctedEff = GetEfficiencyProd(fullCorrEff, 
  	std::make_pair<double,double>(efficiency.first, std::max(efficiency.second.first, efficiency.second.second)));
  
  pair<double, double> survivingEventsAtLumi = GetSurvivingEvents(nSkimmedAtLumi, 
  	std::pair<double,std::pair<double,double> > (correctedEff.first, std::make_pair<double,double>(correctedEff.second, correctedEff.second)));

  float finalNEvents = survivingEventsAtLumi.first;
  float finalNEventError = survivingEventsAtLumi.second;
  if(_source == "Data"){
    finalNEvents = _nVBFSel;
    finalNEventError = sqrt(_nVBFSel);
  }
  
  outLogFile
  << "\t <sample> \n"
  << "\t \t <name> " << _source << " </name> \n"
  << "\t \t <xsection> " << _theXSection.first << " </xsection> \n"
  << "\t \t <lumi> " << PrintWeightedEvent(_theLumi) << " </lumi> \n"
  << "\t \t <allEvents> " << PrintWeightedEvent(_nEvents) << " </allEvents> \n"
  << "\t \t <reweightedEvents> " << PrintWeightedEvent(_puWeightedEvents) << "</reweightedEvents>\n"
  << "\t \t <passingEvents> " << PrintWeightedEvent(_nVBFSel * isoCorr * chargeCorr) << " </passingEvents>\n"
  << "\t \t <expectedEvents> " << PrintWeightedEvent(finalNEvents) << " </expectedEvents>\n"
  << "\t \t <expectedEventsError> " << PrintWeightedEvent(finalNEventError) << " </expectedEventsError> \n"
  << "\t </sample> \n"
  << "\n"
  
  << PrintXMLFormat("cutEff","All", elecPtDeno, elecPtDeno, true)
  << PrintXMLFormat("cutEff","Trigger",_nTriggers, elecPtDeno, true)
  << PrintXMLFormat("cutEff","Elec Pt", _nElecPt, _nTriggers, true)
  << PrintXMLFormat("cutEff","Elec Eta", _nElecEta, _nElecPt, true)
  << PrintXMLFormat("cutEff","Elec DPhiIn", _nElecDPhiIn, _nElecEta, true)
  << PrintXMLFormat("cutEff","Elec DEtaIn", _nElecDEtaIn,_nElecDPhiIn, true)
  << PrintXMLFormat("cutEff","Elec SigmaIEtaIEta", _nElecSigmaIEtaIEta,_nElecDEtaIn, true)
  << PrintXMLFormat("cutEff","Elec H/E", _nElecHadFrac,_nElecSigmaIEtaIEta, true)
  << PrintXMLFormat("cutEff","Elec (1/E - 1/P)", _nElecEOverP,_nElecHadFrac, true)
  << PrintXMLFormat("cutEff","Elec Dxy", _nElecDxy,_nElecEOverP, true)
  << PrintXMLFormat("cutEff","Elec Dz",_nElecDz,_nElecDxy, true)
  << PrintXMLFormat("cutEff","Elec Miss Hits", _nElecMissHits,_nElecDz, true)
  << PrintXMLFormat("cutEff","Elec PFRelIso", _nElecPFRelIso,_nElecMissHits, true)
  << PrintXMLFormat("cutEff","DR(e,tau)", _nElecTauDR,_nElecPFRelIso, true)
  << PrintXMLFormat("cutEff","Tau Pt", _nHadTauPt,_nElecTauDR, true)
  << PrintXMLFormat("cutEff","Tau Eta", _nHadTauEta,_nHadTauPt, true)
  << PrintXMLFormat("cutEff","Tau DecayMode", _nHadTauDecayMode,_nHadTauEta, true)
  << PrintXMLFormat("cutEff","Tau NProngs", _nHadTauNProngs,_nHadTauDecayMode, true)
  << PrintXMLFormat("cutEff","Tau IsoDisc", _nTightTauIso, _nHadTauNProngs, true)
  << PrintXMLFormat("cutEff","Tau ElecDisc", _nHadTauElecDisc,isoCorr,_nHadTauIsoDisc,isoCorr, true)
  << PrintXMLFormat("cutEff","Tau MuonDisc", _nHadTauMuonDisc,isoCorr,_nHadTauElecDisc,isoCorr, true)
  << PrintXMLFormat("cutEff","Charge(e,tau)",  _nTightDiTauCharge,isoCorr,_nHadTauMuonDisc,isoCorr, true)
  << PrintXMLFormat("cutEff","CosDPhi(e,tau)", _nElecHadTauCosDPhi,isoCorr*chargeCorr,_nElecHadTauCharge,isoCorr*chargeCorr, true)
  << PrintXMLFormat("cutEff","MET", _nMEt,isoCorr*chargeCorr, _nElecHadTauCosDPhi,isoCorr*chargeCorr, true)
  << PrintXMLFormat("cutEff","2D PZeta",  _nElecHadTauPZeta,isoCorr*chargeCorr,_nMEt,isoCorr*chargeCorr, true)
  << PrintXMLFormat("cutEff","WJets Veto", _nWJetsVeto,isoCorr*chargeCorr,_nElecHadTauPZeta,isoCorr*chargeCorr, true)
  << PrintXMLFormat("cutEff","BTag Veto",  _nBTags,isoCorr*chargeCorr, _nWJetsVeto,isoCorr*chargeCorr, true)
  << PrintXMLFormat("cutEff","ETau Mass", _nElecTauMass,isoCorr*chargeCorr, _nBTags, isoCorr*chargeCorr, true)
  << "\n"
  << PrintXMLFormat("cutEffVBF","Two VBF Jets", _nVBFJets,isoCorr*chargeCorr,_nElecTauMass,isoCorr*chargeCorr, true)
  << PrintXMLFormat("cutEffVBF","VBF Jets Pt", _nVBFJetsPt,isoCorr*chargeCorr,_nVBFJets,isoCorr*chargeCorr, true)
  << PrintXMLFormat("cutEffVBF","DiJet Delta Eta",  _nDiJetDEta,isoCorr*chargeCorr,_nVBFJetsPt,isoCorr*chargeCorr, true)
  << PrintXMLFormat("cutEffVBF","DiJet Mass", _nDiJetsMass,isoCorr*chargeCorr, _nDiJetDEta,isoCorr*chargeCorr, true)
  << "\n"
  << PrintXMLFormat("summary", "Skim", skimmedEvents.first, skimmedEvents.second, true)
  << PrintXMLFormat("summary", "Elec Acc", _nElecAcc, elecPtDeno, true)
  << PrintXMLFormat("summary", "Elec Id", _nElecId, _nElecAcc, true)
  << PrintXMLFormat("summary", "Tau Acc", _nHadTauAcc, _nElecId, true)
  << PrintXMLFormat("summary", "Tau Id", _nHadTauId, isoCorr, _nHadTauAcc, 1., true)
  << PrintXMLFormat("summary", "Topology", _nElecTauMass,isoCorr*chargeCorr,_nHadTauId, isoCorr, true)
  << PrintXMLFormat("summary", "VBFSel", _nVBFSel , isoCorr * chargeCorr, _nElecTauMass,isoCorr*chargeCorr, true)
  << PrintXMLFormat("summary", "Total", _nVBFSel, isoCorr * chargeCorr, _puWeightedEvents, 1., true)
  << "\n"
  << "\n";
  
  if(_source == "Data") outLogFile << PrintBgFromDataInfo();
  outLogFile << "</etauAnalysis>\n";
  outLogFile.close();
}

std::string ElecTauAnalysis::PrintBgFromDataInfo(){
  //get efficiencies
  
  //zee data driven estimations
  double zeeMETEff = GetEfficiency(_nZEECR1MET, _nZEECR1CentralCuts).first;
  double zeeMETEffError = max(GetEfficiency(_nZEECR1MET, _nZEECR1CentralCuts).second.first, 
  			GetEfficiency(_nZEECR1MET, _nZEECR1CentralCuts).second.second);
  pair<double, double> zeeMETEffWithError = make_pair<double,double>(zeeMETEff, zeeMETEffError);
  
  double zeeOSEff = GetEfficiency(_nZEECR1OS, _nZEECR1CentralCuts).first;
  double zeeOSEFFError = max(GetEfficiency(_nZEECR1OS, _nZEECR1CentralCuts).second.first, 
  			GetEfficiency(_nZEECR1OS, _nZEECR1CentralCuts).second.second);
  pair<double, double> zeeOSEffWithError = make_pair<double,double>(zeeOSEff, zeeOSEFFError);
  
  double zeeLSEff = GetEfficiency(_nZEECR1LS, _nZEECR1CentralCuts).first;
  double zeeLSEFFError = max(GetEfficiency(_nZEECR1LS, _nZEECR1CentralCuts).second.first, 
  			GetEfficiency(_nZEECR1LS, _nZEECR1CentralCuts).second.second);
  pair<double, double> zeeLSEffWithError = make_pair<double,double>(zeeLSEff, zeeLSEFFError);
  
  double zeeVBFEff = GetEfficiency(_nZEECR1VBFSel, _nZEECR1CentralCuts).first;
  double zeeVBFEFFError = max(GetEfficiency(_nZEECR1VBFSel, _nZEECR1CentralCuts).second.first, 
  			GetEfficiency(_nZEECR1VBFSel, _nZEECR1CentralCuts).second.second);
  pair<double, double> zeeVBFEffWithError = make_pair<double,double>(zeeVBFEff, zeeVBFEFFError);
  
  double zeeElecDiscEff = GetEfficiency(_nZEECR2ElecDisc, _nZEECR2).first;
  double zeeElecDiscEFFError = max(GetEfficiency(_nZEECR2ElecDisc, _nZEECR2).second.first, 
  			GetEfficiency(_nZEECR2ElecDisc, _nZEECR2).second.second);
  pair<double, double> zeeElecDiscEffWithError = make_pair<double,double>(zeeElecDiscEff, zeeElecDiscEFFError);
  
  double zeeAntiElecDiscEff = GetEfficiency(_nZEECR2AntiElecDisc, _nZEECR2).first;
  double zeeAntiElecDiscEFFError = max(GetEfficiency(_nZEECR2AntiElecDisc, _nZEECR2).second.first, 
  			GetEfficiency(_nZEECR2AntiElecDisc, _nZEECR2).second.second);
  pair<double, double> zeeAntiElecDiscEffWithError = make_pair<double,double>(zeeAntiElecDiscEff, zeeAntiElecDiscEFFError);
  
  pair<double, double> zeeElecAntiElecEffRatio = GetEfficiencyRatio(zeeElecDiscEffWithError, zeeAntiElecDiscEffWithError);
  
  pair<double, double> zeeLSCentralEff = GetEfficiencyProd(GetEfficiencyProd(zeeMETEffWithError, zeeLSEffWithError),zeeElecAntiElecEffRatio);
  pair<double, double> zeeOSCentralEff = GetEfficiencyProd(GetEfficiencyProd(zeeMETEffWithError, zeeOSEffWithError),zeeElecAntiElecEffRatio);
  
  pair<double,double> zeeLSBG =  GetSurvivingEvents(_nZEECR1CentralCuts, zeeLSCentralEff);
  pair<double,double> zeeOSBG =  GetSurvivingEvents(_nZEECR1CentralCuts, zeeOSCentralEff);

  pair<double,double> zeeLSCentralVBFEff = GetEfficiencyProd(zeeLSCentralEff,zeeVBFEffWithError);
  pair<double,double> zeeOSCentralVBFEff = GetEfficiencyProd(zeeOSCentralEff,zeeVBFEffWithError);
  
  pair<double,double> zeeLSCentralVBFBG = GetSurvivingEvents(_nZEECR1CentralCuts,zeeLSCentralVBFEff);
  pair<double,double> zeeOSCentralVBFBG = GetSurvivingEvents(_nZEECR1CentralCuts,zeeOSCentralVBFEff);

  //wjets data driven estimation
  double wjetsCosDPhiEff = GetEfficiency(_nWJetsCR1ElecTauCosDPhiPZeta,_nWJetsCR1).first;
  double wjetsCosDPhiEffError = max(GetEfficiency(_nWJetsCR1ElecTauCosDPhiPZeta,_nWJetsCR1).second.first, 
                                    GetEfficiency(_nWJetsCR1ElecTauCosDPhiPZeta,_nWJetsCR1).second.second);
  pair<double,double> wjetsCosDPhiEffWithError = make_pair<double,double>(wjetsCosDPhiEff,wjetsCosDPhiEffError);
  
  double wjetsMTEff = GetEfficiency(_nWJetsCR2MtVeto,_nWJetsCR2).first;
  double wjetsMTEffError = max(GetEfficiency(_nWJetsCR2MtVeto,_nWJetsCR2).second.first,
                               GetEfficiency(_nWJetsCR2MtVeto,_nWJetsCR2).second.second);
  pair<double,double> wjetsMTEffWithError = make_pair<double,double>(wjetsMTEff, wjetsMTEffError);
  
  //double wjetsMTAntiEff = 1. - wjetsMTEff;
  //double wjetsMTRatio = wjetsMTEff/wjetsMTAntiEff;
  //double wjetsMTRatioError = ((1./wjetsMTAntiEff + (wjetsMTEff)/(wjetsMTAntiEff*wjetsMTAntiEff))*wjetsMTEffError);
  //pair<double,double> wjetsMTRatioEffWithError = make_pair<double,double>(wjetsMTRatio,wjetsMTRatioError);
  
  pair<double,double> wjetsFinalEff = GetEfficiencyRatio(wjetsCosDPhiEffWithError,wjetsMTEffWithError);
  
  pair<double,double> wjetsBG =  GetSurvivingEvents(_nWJetsCR1,wjetsFinalEff);
  
  //ttjets data driven estimation
  pair<double,pair<double,double> > efficiencyTTbarCR2 = GetEfficiency(_nTTbarCR2, _puWeightedEvents);
  pair<double, double> survivingEventsAtLumiTTbarCR2 = GetSurvivingEvents(nSkimmedAtLumi, efficiencyTTbarCR2);
  
  double ttjetsCosDPhiEff = GetEfficiency(_nTTbarCR1ElecTauCosDPhiPZeta,_nTTbarCR1).first;
  double ttjetsCosDPhiEffError = max(GetEfficiency(_nTTbarCR1ElecTauCosDPhiPZeta,_nTTbarCR1).second.first, 
                                    GetEfficiency(_nTTbarCR1ElecTauCosDPhiPZeta,_nTTbarCR1).second.second);
  double ttjetsBtagEff = GetEfficiency(_nTTbarCR2BtagVeto,_nTTbarCR2).first;
  double ttjetsBtagEffError = max(GetEfficiency(_nTTbarCR2BtagVeto,_nTTbarCR2).second.first,
                               GetEfficiency(_nTTbarCR2BtagVeto,_nTTbarCR2).second.second);
  pair<double,double> ttjetsBtagEffWithError = make_pair<double,double>(ttjetsBtagEff,ttjetsBtagEffError);
  
  double ttjetsTwoBtagEff = GetEfficiency(_nTTbarCR2TwoBtag,_nTTbarCR2).first;
  double ttjetsTwoBtagEffError = max(GetEfficiency(_nTTbarCR2TwoBtag,_nTTbarCR2).second.first,
                               GetEfficiency(_nTTbarCR2TwoBtag,_nTTbarCR2).second.second);
  pair<double,double> ttjetsTwoBtagEffWithError = make_pair<double,double>(ttjetsTwoBtagEff, ttjetsTwoBtagEffError);
  
  pair<double,double> ttjetsBtagRatioEffWithError = GetEfficiencyRatio(ttjetsBtagEffWithError,ttjetsTwoBtagEffWithError);
    
  pair<double,double> ttjetsCosDPhiEffWithError = make_pair<double,double>(ttjetsCosDPhiEff,ttjetsCosDPhiEffError);
  pair<double,double> ttjetsFinalEff = GetEfficiencyProd(ttjetsCosDPhiEffWithError,ttjetsBtagRatioEffWithError);
  
  pair<double,double> ttjetsBG =  GetSurvivingEvents(_nTTbarCR1,ttjetsFinalEff);
  
  pair<double,pair<double,double> > qcdLSSelEff = GetEfficiency(_nQCDLSSel, _puWeightedEvents);
  pair<double,double> qcdLSOSRatio = GetRatio(_nQCDTauRevIsoOS, _nQCDTauRevIsoLS);
  pair<double, double> lsEventsAtLumi = GetSurvivingEvents(nSkimmedAtLumi, qcdLSSelEff);
  if(_source == "Data"){
    lsEventsAtLumi.first = _nQCDLSSel;
    lsEventsAtLumi.second = sqrt(_nQCDLSSel);
    survivingEventsAtLumiTTbarCR2.first = _nTTbarCR2;
    survivingEventsAtLumiTTbarCR2.second = sqrt(_nTTbarCR2);
  }
  
  pair<double, double> qcdBGEst = GetSurvivingEvents(lsEventsAtLumi.first, qcdLSOSRatio);
 
  
  std::stringstream theOutputStream;
  theOutputStream.setf(ios::fixed);
  theOutputStream.setf(ios::showpoint);
  theOutputStream.precision(4);
  theOutputStream << std::scientific
  //theOutputStream 
    << PrintXMLFormat("zeeCR1", "ZEECR1 All", _nZEECR1CentralCuts, _nZEECR1CentralCuts, false)
    << PrintXMLFormat("zeeCR1", "ZEECR1 LS", _nZEECR1LS, _nZEECR1CentralCuts, false)
    << PrintXMLFormat("zeeCR1", "ZEECR1 OS", _nZEECR1OS, _nZEECR1CentralCuts, false)
    << PrintXMLFormat("zeeCR1", "ZEECR1 MET", _nZEECR1MET, _nZEECR1CentralCuts, false)
    << PrintXMLFormat("zeeCR1", "ZEECR1 VBF", _nZEECR1VBFSel, _nZEECR1CentralCuts, false)
    << PrintXMLFormat("zeeCR2", "ZEECR2 All",  _nZEECR2, _nZEECR2, false)
    << PrintXMLFormat("zeeCR2", "ZEECR2 AntiElecDisc", _nZEECR2AntiElecDisc, _nZEECR2, false)
    << PrintXMLFormat("zeeCR2", "ZEECR2 ElecDisc", _nZEECR2ElecDisc, _nZEECR2, false)
    << PrintXMLFormat("wjetsCR1", "WJetsCR1 All", _nWJetsCR1, _nWJetsCR1, false)
    << PrintXMLFormat("wjetsCR1", "WJetsCR1 LS", _nWJetsCR1LS, _nWJetsCR1, false)
    << PrintXMLFormat("wjetsCR1", "WJetsCR1 OS", _nWJetsCR1OS, _nWJetsCR1, false)
    << PrintXMLFormat("wjetsCR1", "WJets CR1 TauIso", _nWJetsCR1TauIso, _nWJetsCR1, false)
    << PrintXMLFormat("wjetsCR1", "WJets CR1 VBF", _nWJetsCR1VBF, _nWJetsCR1, false)
    //<< PrintXMLFormat("wjetsCR1", "WJetsCR1 CosDPhi &amp; PZeta", _nWJetsCR1ElecTauCosDPhiPZeta,_nWJetsCR1, false)
    //<< PrintXMLFormat("wjetsCR2", "WJetsCR2 All", _nWJetsCR2,_nWJetsCR2, false)
    //<< PrintXMLFormat("wjetsCR2", "WJetsCR2 Mt(e,MET)", _nWJetsCR2MtVeto,_nWJetsCR2, false)
    << PrintXMLFormat("ttjetsCR1", "TTJetsCR1 All", _nTTbarCR1,_nTTbarCR1, false)
    << PrintXMLFormat("ttjetsCR1", "TTJetsCR1 OS", _nTTbarCR1OS,_nTTbarCR1, false)
    << PrintXMLFormat("ttjetsCR1", "TTJetsCR1 LS", _nTTbarCR1LS,_nTTbarCR1, false)
    << PrintXMLFormat("ttjetsCR1", "TTJetsCR1 VBF", _nTTbarVBFSel,_nTTbarCR1, false)
    //<< PrintXMLFormat("ttjetsCR1", "TTJetsCR1 CosDPhi &amp; PZeta Eff", _nTTbarCR1ElecTauCosDPhiPZeta,_nTTbarCR1, false)
    //<< PrintXMLFormat("ttjetsCR2", "TTJetsCR2 All", _nTTbarCR2, _nTTbarCR2, false)
    //<< PrintXMLFormat("ttjetsCR2", "TTJetsCR2 0 BTags", _nTTbarCR2BtagVeto,_nTTbarCR2, false)
    //<< PrintXMLFormat("ttjetsCR2", "TTJetsCR2 &gt; 1 BTag", _nTTbarCR2TwoBtag,_nTTbarCR2, false)
    << PrintXMLFormat("zttCR", "ZTTCR All", _nZTTCR, _nZTTCR, false)
    << PrintXMLFormat("zttCR", "ZTTCR VBF", _nZTTCRVBF, _nZTTCR, false)
    << "\n"
    << "\t <bgEst> \n"
    << "\t \t <zee> \n"
    << "\t \t \t <zeeOSBGEst> " << zeeOSBG.first << "</zeeOSBGEst> \n"
    << "\t \t \t <zeeOSBGEstError> " << zeeOSBG.second << "</zeeOSBGEstError> \n"
    << "\t \t \t <zeeLSBGEst> " << zeeLSBG.first << "</zeeLSBGEst> \n"
    << "\t \t \t <zeeLSBGEstError> " << zeeLSBG.second << "</zeeLSBGEstError> \n"
    << "\t \t \t <zeeOSVBFBGEst> " << zeeOSCentralVBFBG.first << "</zeeOSVBFBGEst> \n"
    << "\t \t \t <zeeOSVBFBGEstError> " << zeeOSCentralVBFBG.second << "</zeeOSVBFBGEstError> \n"
    << "\t \t \t <zeeLSVBFBGEst> " << zeeLSCentralVBFBG.first << "</zeeLSVBFBGEst> \n"
    << "\t \t \t <zeeLSVBFBGEstError> " << zeeLSCentralVBFBG.second << "</zeeLSVBFBGEstError> \n"
    << "\t \t </zee> \n"
    << "\t \t <wjets> \n"
    << "\t \t \t  <wjetsBGEst> " << wjetsBG.first << " </wjetsBGEst> \n"
    << "\t \t \t  <wjetsBGEstError> " << wjetsBG.second << " </wjetsBGEstError> \n"
    << "\t \t </wjets> \n"
    << "\t \t <ttjets> \n"
    << "\t \t \t <normTTJetsCR2> " << survivingEventsAtLumiTTbarCR2.first << " </normTTJetsCR2> \n"
    << "\t \t \t <normTTJetsCR2Error> " << survivingEventsAtLumiTTbarCR2.second << " </normTTJetsCR2Error> \n"
    << "\t \t \t <ttjetsBtagRatio> " << ttjetsBtagRatioEffWithError.first << " </ttjetsBtagRatio> \n"
    << "\t \t \t <ttjetsBtagRatioError> " << ttjetsBtagRatioEffWithError.second << " </ttjetsBtagRatioError> \n"
    << "\t \t \t <ttjetsBGEst> " << ttjetsBG.first << " </ttjetsBGEst> \n"
    << "\t \t \t <ttjetsBGEstError> " << ttjetsBG.second << " </ttjetsBGEstError> \n"
    << "\t \t </ttjets> \n"
    << "\t \t <qcd> \n"
    << "\t \t \t <qcdLS> " << lsEventsAtLumi.first << " </qcdLS> \n"
    << "\t \t \t <qcdRevTauIso> " << PrintWeightedEvent(_nQCDTauRevIso) << " </qcdRevTauIso> \n"
    << "\t \t \t <qcdRevIsoLS> " << PrintWeightedEvent(_nQCDTauRevIsoLS) << " </qcdRevIsoLS> \n"
    << "\t \t \t <qcdRevIsoOS> " << PrintWeightedEvent(_nQCDTauRevIsoOS) << " </qcdRevIsoOS> \n"
    << "\t \t \t <qcdBGEst> " << PrintWeightedEvent(qcdBGEst.first) << " </qcdBGEst> \n"
    << "\t \t \t <qcdBGEstError> " << PrintWeightedEvent(qcdBGEst.second) << " </qcdBGEstError> \n"
    << "\t \t </qcd> \n"
    << "\t </bgEst> \n"  
    << "\n";
    return theOutputStream.str();
}

std::string ElecTauAnalysis::PrintXMLFormat(std::string tag, std::string name, double nume, double corrNume, 
					   double deno, double corrDeno, bool printAbsEff = false){
  std::pair<double, std::pair<double, double> > theEff = GetEfficiency(nume, deno);
  
  std::stringstream theOutputStream;
  theOutputStream.setf(ios::showpoint);
  theOutputStream.setf(std::ios::fixed);
  theOutputStream.precision(2);
  double theAbsEff = -100.;
  if(printAbsEff) theAbsEff = GetEfficiency(nume*corrNume, _puWeightedEvents).first*100.;
  
  theOutputStream
  << "\t <" << tag << "> \n"
  << "\t \t <name> " << name << " </name> \n"
  << "\t \t <nume> " << PrintWeightedEvent(nume*corrNume) << " </nume> \n"
  << "\t \t <deno> " << PrintWeightedEvent(deno*corrDeno) << " </deno> \n"
  << "\t \t <eff> " << theEff.first*100. << " </eff> \n";	
  if(theEff.first < 0.01){
    theOutputStream 
    << "\t \t <effErrorLow> " << std::scientific << theEff.second.first*100. << " </effErrorLow> \n"
    << "\t \t <effErrorHi> " << theEff.second.second*100. << " </effErrorHi> \n";	
  }
  else {
    theOutputStream 
    << "\t \t <effErrorLow> " << theEff.second.first*100. << " </effErrorLow> \n"
    << "\t \t <effErrorHi> " << theEff.second.second*100. << " </effErrorHi> \n";
  }
  
  if(theAbsEff < 0.01) theOutputStream << "\t \t <absEff> " << std::scientific << theAbsEff << " </absEff>\n"; 
  else theOutputStream << "\t \t <absEff> " << theAbsEff << " </absEff>\n"; 
  theOutputStream << "\t </" << tag <<"> \n";
  return theOutputStream.str();
}

std::string ElecTauAnalysis::PrintXMLFormat(std::string tag, std::string name, double nume, double deno, bool printAbsEff){
  return PrintXMLFormat(tag, name, nume, 1., deno, 1., printAbsEff);
}

void ElecTauAnalysis::writeOutFile(){
  
  outRootFile->Write();
  outRootFile->Close();
} 

void ElecTauAnalysis::fillHistos(){
  _pileUPDistro->Fill(InTimePU);									       
  _pileUPReweighted->Fill(InTimePU,_eventPUWeight);
  
  _elecPtHisto->Fill(_maxElecPt,_eventPUWeight);
  _elecEtaHisto->Fill(_minElecEta,_eventPUWeight);
  _elecDPhiInEBHisto->Fill(_minElecDPhiEB,_eventPUWeight);
  _elecDPhiInEEHisto->Fill(_minElecDPhiEE,_eventPUWeight);
  _elecDEtaInEBHisto->Fill(_minElecDEtaEB,_eventPUWeight);
  _elecDEtaInEEHisto->Fill(_minElecDEtaEE,_eventPUWeight);
  _elecSigmaIEtaIEtaEBHisto->Fill(_minElecSigmaIEtaIEtaEB,_eventPUWeight);
  _elecSigmaIEtaIEtaEEHisto->Fill(_minElecSigmaIEtaIEtaEE,_eventPUWeight);
  _elecHOverEEBHisto->Fill(_minElecHOverEEB,_eventPUWeight);
  _elecHOverEEEHisto->Fill(_minElecHOverEEE,_eventPUWeight);
  _elecOneOverEMinusOneOverPEBHisto->Fill(_minElecOneOverEMinusOneOverPEB,_eventPUWeight);
  _elecOneOverEMinusOneOverPEEHisto->Fill(_minElecOneOverEMinusOneOverPEE,_eventPUWeight);
  _elecDxyVtxHisto->Fill(_minElecDxy,_eventPUWeight);
  _elecDzVtxHisto->Fill(_minElecDz,_eventPUWeight);
  _elecMissHisto->Fill(_minElecMissHits,_eventPUWeight);
  _elecPFRelIsoDR03Histo->Fill(_minElecPFRelIsoDR03,_eventPUWeight);
  _elecPFRelIsoDR04Histo->Fill(_minElecPFRelIsoDR04,_eventPUWeight);
  
  _elecTauDRHisto->Fill(_maxElecTauDR,_eventPUWeight);
  _tauPtHisto->Fill(_maxTauPt,_eventPUWeight);
  _tauAltLVPtHisto->Fill(_maxTauAltLVPt,_eventPUWeight);
  _tauEtaHisto->Fill(_minTauEta,_eventPUWeight);
  _tauDecayModeHisto->Fill(_maxTauDecayMode,_eventPUWeight);
  _tauNProngsHisto->Fill(_maxTauNProngs,_eventPUWeight);
  _tauMuonDiscHisto->Fill(_tauMuonDisc,_eventPUWeight);  							           
  _tauElecDiscHisto->Fill(_tauElecDisc,_eventPUWeight);
  _tauElecMVA5DiscHisto->Fill(_tauMVA5ElecDisc,_eventPUWeight);
  _tauIsoDB3HitsHisto->Fill(_tauIsoDB,_eventPUWeight);
  _tauIsoMVAHisto->Fill(_tauMVAIso,_eventPUWeight);
  _tauGenPtBeforeDecayModeHisto->Fill(_maxTauGenPtBeforeDecayMode,_eventPUWeight);
  _tauGenPtAfterDecayModeHisto->Fill(_maxTauGenPtAfterDecayMode,_eventPUWeight);
  _tauGenPtAfterMuonDiscHisto->Fill(_maxTauGenPtAfterMuonDisc,_eventPUWeight);
  _tauGenPtAfterElecDiscHisto->Fill(_maxTauGenPtAfterElecDisc,_eventPUWeight);
  _tauGenPtBeforeIsoHisto->Fill(_maxTauGenPtBeforeIso,_eventPUWeight);
  _tauGenPtAfterIsoHisto->Fill(_maxTauGenPtAfterIso,_eventPUWeight);
  _tauPtBeforeDecayModeHisto->Fill(_maxTauPtBeforeDecayMode,_eventPUWeight);
  _tauPtAfterDecayModeHisto->Fill(_maxTauPtAfterDecayMode,_eventPUWeight);
  _tauPtAfterMuonDiscHisto->Fill(_maxTauPtAfterMuonDisc,_eventPUWeight);
  _tauPtAfterElecDiscHisto->Fill(_maxTauPtAfterElecDisc,_eventPUWeight);
  _tauPtBeforeIsoHisto->Fill(_maxTauPtBeforeIso,_eventPUWeight);
  _tauPtAfterIsoHisto->Fill(_maxTauPtAfterIso,_eventPUWeight);
  _ptResAfterIsoHisto->Fill((_maxTauPtAfterIso-_maxTauGenPtAfterIso)/_maxTauGenPtAfterIso,_eventPUWeight);
  _ptResBeforeTauIdHisto->Fill((_maxTauPtBeforeDecayMode-_maxTauGenPtBeforeDecayMode)/_maxTauGenPtBeforeDecayMode,_eventPUWeight);
  _elecTauChargeHisto->Fill(_minElecTauCharge,_eventPUWeight);
  _elecTauCosDPhiHisto->Fill(_minElecTauCosDPhi,_eventPUWeight);
  _metHisto->Fill(_maxMET,_eventPUWeight);
  _diElecAsymmVsdiElecMassDiff->Fill(fabs(zPtAsymm),fabs(zMassDiff),_eventPUWeight);
  _pZetaVsPZetaVis2DHisto->Fill(_minPZetaVis,_maxPZeta,_eventPUWeight);
  _pZetaHisto->Fill(_maxPZeta,_eventPUWeight);
  _pZetaVisHisto->Fill(_minPZetaVis,_eventPUWeight);
  _nBTagsTCHEMHisto->Fill(_minBtagsTCHEM,_eventPUWeight);
  _nBTagsTCHPTHisto->Fill(_minBtagsTCHPT,_eventPUWeight);
  _nBTagsCSVLHisto->Fill(_minBtagsCSVL,_eventPUWeight);
  _nBTagsCSVMHisto->Fill(_minBtagsCSVM,_eventPUWeight);
  _nBTagsCSVTHisto->Fill(_minBtagsCSVT,_eventPUWeight);
  _nJetsHisto->Fill(_minNJets, _eventPUWeight);
  _jetEtSumHisto->Fill(_maxJetEtSum,_eventPUWeight);
  _vbfDiJetMassHisto->Fill(_maxDiJetMass,_eventPUWeight);
  _elecMETMassHisto->Fill(_maxElecMETMass,_eventPUWeight);
  _elecMETCosDPhi->Fill(_minElecMETCosDPhi,_eventPUWeight);
  _elecMETMassVsElecMETCosDPhi->Fill(_minElecMETCosDPhi,_maxElecMETMass,_eventPUWeight);
  _elecMETMassVsElecTauCosDPhi->Fill(_minElecTauCosDPhi,_maxElecMETMass,_eventPUWeight);
  _nVBFJetsHisto->Fill(_maxVBFJets,_eventPUWeight);
  _elecTauMassHisto->Fill(_maxElecTauMass,_eventPUWeight);
  _elecTauLSMassHisto->Fill(_maxElecTauLSMass,_eventPUWeight);
  _elecTauMetMassHisto->Fill(_maxElecTauMetMass,_eventPUWeight);
  _elecTauMetLSMassHisto->Fill(_maxElecTauMetLSMass,_eventPUWeight);
  _nBTagsTCHEMTTJetsCRHisto->Fill(_minBtagsTCHEMTTJetsCR,_eventPUWeight);
  _nBTagsTCHPTTTJetsCRHisto->Fill(_minBtagsTCHPTTTJetsCR,_eventPUWeight);
  _nBTagsCSVLTTJetsCRHisto->Fill(_minBtagsCSVLTTJetsCR,_eventPUWeight);
  _nBTagsCSVMTTJetsCRHisto->Fill(_minBtagsCSVMTTJetsCR,_eventPUWeight);
  _nBTagsCSVTTTJetsCRHisto->Fill(_minBtagsCSVTTTJetsCR,_eventPUWeight);
  _tauPtTTbarCR1Histo->Fill(_maxTauPtTTbarCR1,_eventPUWeight);
  _tauPtTTbarCR1HistoNoTauId->Fill(_maxTauPtTTbarCRNoTauId,_eventPUWeight);
  _tauPtAfterAllCutsHisto->Fill(_maxTauPtAfterAllCuts,_eventPUWeight);
  _elecPtAfterAllCutsHisto->Fill(_maxElecPtAfterAllCuts,_eventPUWeight);
  _metAfterAllCutsHisto->Fill(_maxMETAfterAllCuts,_eventPUWeight);
  
  _elecPtNMinus1Histo->Fill(_maxElecPtNMinus1,_eventPUWeight);
  _elecEtaNMinus1Histo->Fill(_minElecEtaNMinus1,_eventPUWeight);
  _elecDPhiInEBNMinus1Histo->Fill(_minElecDPhiEBNMinus1,_eventPUWeight);
  _elecDPhiInEENMinus1Histo->Fill(_minElecDPhiEENMinus1,_eventPUWeight);
  _elecDEtaInEBNMinus1Histo->Fill(_minElecDEtaEBNMinus1,_eventPUWeight);
  _elecDEtaInEENMinus1Histo->Fill(_minElecDEtaEENMinus1,_eventPUWeight);
  _elecSigmaIEtaIEtaEBNMinus1Histo->Fill(_minElecSigmaIEtaIEtaEBNMinus1,_eventPUWeight);
  _elecSigmaIEtaIEtaEENMinus1Histo->Fill(_minElecSigmaIEtaIEtaEENMinus1,_eventPUWeight);
  _elecHOverEEBNMinus1Histo->Fill(_minElecHOverEEBNMinus1,_eventPUWeight);
  _elecHOverEEENMinus1Histo->Fill(_minElecHOverEEENMinus1,_eventPUWeight);
  _elecOneOverEMinusOneOverPEBNMinus1Histo->Fill(_minElecOneOverEMinusOneOverPEBNMinus1,_eventPUWeight);
  _elecOneOverEMinusOneOverPEENMinus1Histo->Fill(_minElecOneOverEMinusOneOverPEENMinus1,_eventPUWeight);
  _elecDxyVtxNMinus1Histo->Fill(_minElecDxyNMinus1,_eventPUWeight);
  _elecDzVtxNMinus1Histo->Fill(_minElecDzNMinus1,_eventPUWeight);
  _elecMissNMinus1Histo->Fill(_minElecMissHitsNMinus1,_eventPUWeight);
  _elecPFRelIsoDR03NMinus1Histo->Fill(_minElecPFRelIsoDR03NMinus1,_eventPUWeight);
  _elecPFRelIsoDR04NMinus1Histo->Fill(_minElecPFRelIsoDR04NMinus1,_eventPUWeight);
  
  _elecTauDRNMinus1Histo->Fill(_maxElecTauDRNMinus1,_eventPUWeight);
  _tauPtNMinus1Histo->Fill(_maxTauPtNMinus1,_eventPUWeight);
  _tauAltLVPtNMinus1Histo->Fill(_maxTauAltLVPtNMinus1,_eventPUWeight);
  _tauEtaNMinus1Histo->Fill(_minTauEtaNMinus1,_eventPUWeight);
  _tauDecayModeNMinus1Histo->Fill(_maxTauDecayModeNMinus1,_eventPUWeight);
  _tauNProngsNMinus1Histo->Fill(_maxTauNProngsNMinus1,_eventPUWeight);
  _tauMuonDiscNMinus1Histo->Fill(_tauMuonDiscNMinus1,_eventPUWeight);  							           
  _tauElecDiscNMinus1Histo->Fill(_tauElecDiscNMinus1,_eventPUWeight);
  _tauElecMVA5DiscNMinus1Histo->Fill(_tauMVA5ElecDiscNMinus1,_eventPUWeight);
  _tauIsoDB3HitsNMinus1Histo->Fill(_tauIsoDB3HitsNMinus1,_eventPUWeight);
  _tauIsoMVANMinus1Histo->Fill(_tauMVAIsoNMinus1,_eventPUWeight);
  _elecTauChargeNMinus1Histo->Fill(_minElecTauChargeNMinus1,_eventPUWeight);
  _elecTauCosDPhiNMinus1Histo->Fill(_minElecTauCosDPhiNMinus1,_eventPUWeight);
  _metNMinus1Histo->Fill(_maxMETNMinus1,_eventPUWeight);
  _pZetaVsPZetaVis2DNMinus1Histo->Fill(_minPZetaVisNMinus1,_maxPZetaNMinus1,_eventPUWeight);
  _pZetaNMinus1Histo->Fill(_maxPZetaNMinus1,_eventPUWeight);
  _pZetaVisNMinus1Histo->Fill(_minPZetaVisNMinus1,_eventPUWeight);
  _pZetaMinusPZetaVisNMinus1Histo->Fill(_maxPZetaNMinus1 - _zeta2DSlope * _minPZetaVisNMinus1,_eventPUWeight);
  _nBTagsTCHEMNMinus1Histo->Fill(_minBtagsTCHEMNMinus1,_eventPUWeight);
  _nBTagsTCHPTNMinus1Histo->Fill(_minBtagsTCHPTNMinus1,_eventPUWeight);
  _nBTagsCSVLNMinus1Histo->Fill(_minBtagsCSVLNMinus1,_eventPUWeight);
  _nBTagsCSVMNMinus1Histo->Fill(_minBtagsCSVMNMinus1,_eventPUWeight);
  _nBTagsCSVTNMinus1Histo->Fill(_minBtagsCSVTNMinus1,_eventPUWeight);
  _nJetsNMinus1Histo->Fill(_minNJetsNminus1,_eventPUWeight);
  _jetEtSumNMinus1Histo->Fill(_maxJetEtSumNMinus1,_eventPUWeight);
  _nJetsVsNBtagsCSVM2DHisto->Fill(_minBtagsCSVMNMinus1, _minNJetsNminus1, _eventPUWeight);
  _jetEtSumVsNBtagsCSVM2DHisto->Fill(_minBtagsCSVMNMinus1, _maxJetEtSumNMinus1, _eventPUWeight);
  _nVBFJetsNMinus1Histo->Fill(_maxNVBFJetsNMinus1,_eventPUWeight);
  _jet1PtNMinus1Histo->Fill(_maxJet1PtNMinus1,_eventPUWeight);
  _jet2PtNMinus1Histo->Fill(_maxJet2PtNMinus1,_eventPUWeight);
  _vbfDiJetDEtaNMinus1Histo->Fill(_maxDiJetDEtaNMinus1,_eventPUWeight);
  _vbfDiJetMassNMinus1Histo->Fill(_maxDiJetMassNMinus1,_eventPUWeight);
  
  _zeeCRVBFDiJetMassHisto->Fill(_maxZEECRDiJetMass,_eventPUWeight);
  _ztautauCRVBFDiJetMassHisto->Fill(_maxZTTCRDiJetMass,_eventPUWeight);
  _wjetsCRVBFDiJetMassHisto->Fill(_maxWJetsCRDiJetMass,_eventPUWeight);
  _ttjetsCRVBFDiJetMassHisto->Fill(_maxTTJetsCRDiJetMass,_eventPUWeight);
}

void ElecTauAnalysis::bookHistos(){
  outRootFile = new TFile(_outRootFileName.c_str(), "Recreate");
  TDirectory* theDir; 
  std::string histosDir = _source + "HistosDirectory"; 
  if(outRootFile->FindObjectAny(histosDir.c_str()) == 0) theDir = outRootFile->mkdir(histosDir.c_str());
  else theDir = (TDirectory*)outRootFile->FindObjectAny(histosDir.c_str());
  theDir->cd();
  
  float ptAndMetBinning[25] = {0,10,20,30,40,50,60,70,80,95,110,125,140,160,180,200,220,250,280,310,340,390,440,500,600};
  float massBinning[27] = {0,20,40,60,80,100,120,140,160,180,200,225,250,275,300,333,366,400,450,500,600,700,800,1000,1250,1500,2000};
  float vbfMassBinning[35] = {0,5,10,50,100,150,200,250,300,350,400,450,500,550,600,650,700,750,800,850,900,950,1000,1050,1100,1150,1200,1250,1300,1400,1500,1750,2000,2500,3000};
  
  _pileUPDistro = new TH1F("_pileUPDistro", "PileUp", 50, 0., 50.);
  _pileUPReweighted = new TH1F("_pileUPReweighted", "PileUp Reweighted", 50, 0., 50.);
    
  _elecPtHisto = new TH1F("_elecPtHisto", "Elec P_{T}", sizeof(ptAndMetBinning)/sizeof(float)-1, ptAndMetBinning);
  _elecEtaHisto = new TH1F("_elecEtaHisto", "Elec #eta", 100, -5., 5.);
  _elecDPhiInEBHisto = new TH1F("_elecDPhiInEBHisto", "Elec #Delta#phi_{in} EB", 100, -0.5, 0.5);
  _elecDPhiInEEHisto = new TH1F("_elecDPhiInEEHisto", "Elec #Delta#phi_{in} EE", 100, -0.5, 0.5);
  _elecDEtaInEBHisto = new TH1F("_elecDEtaInEBHisto", "Elec #Delta#eta_{in} EB", 200, -0.1, 0.1);
  _elecDEtaInEEHisto = new TH1F("_elecDEtaInEEHisto", "Elec #Delta#eta_{in} EE", 200, -0.1, 0.1);
  _elecSigmaIEtaIEtaEBHisto = new TH1F("_elecSigmaIEtaIEtaEBHisto", "Elec #sigma_{i#eta i#eta} EB", 100, 0., 0.1); 
  _elecSigmaIEtaIEtaEEHisto = new TH1F("_elecSigmaIEtaIEtaEEHisto", "Elec #sigma_{i#eta i#eta} EE", 100, 0., 0.1);
  _elecHOverEEBHisto = new TH1F("_elecHOverEEBHisto", "Elec H/E EB", 200, 0, 2);
  _elecHOverEEEHisto = new TH1F("_elecHOverEEEHisto", "Elec H/E EE", 200, 0, 2);
  _elecOneOverEMinusOneOverPEBHisto = new TH1F("_elecOneOverEMinusOneOverPEBHisto", "Elec 1/E - 1/p EB", 100, -0.5, 0.5);
  _elecOneOverEMinusOneOverPEEHisto = new TH1F("_elecOneOverEMinusOneOverPEEHisto", "Elec 1/E - 1/p EE", 100, -0.5, 0.5);
  _elecDxyVtxHisto = new TH1F("_elecDxyVtxHisto", "Elec d_{xy}", 200, -2, 2);
  _elecDzVtxHisto = new TH1F("_elecDzVtxHisto", "Elec d_{z}", 200, -2, 2);
  _elecMissHisto = new TH1F("_elecMissHisto", "Elec Missing Hits", 10, 0, 10);
  _elecPFRelIsoDR03Histo = new TH1F("_elecPFRelIsoDR03Histo", "Elec PFRelIso #Delta#beta #Delta R = 0.3", 250, 0, 10);
  _elecPFRelIsoDR04Histo = new TH1F("_elecPFRelIsoDR04Histo", "Elec PFRelIso #Delta#beta #Delta R = 0.4", 250, 0, 10);
  _elecTauDRHisto = new TH1F("_elecTauDRHisto", "#Delta R(e,#tau)", 100, 0, 6.28);
  _tauPtHisto = new TH1F("_tauPtHisto", "Tau p_{T}", sizeof(ptAndMetBinning)/sizeof(float)-1, ptAndMetBinning);
  _tauAltLVPtHisto = new TH1F("_tauAltLVPtHisto", "Tau p_{T}", sizeof(ptAndMetBinning)/sizeof(float)-1, ptAndMetBinning);
  _tauEtaHisto = new TH1F("_tauEtaHisto", "Tau #eta",  100, -5., 5.);
  _tauDecayModeHisto = new TH1F("_tauDecayModeHisto", "Tau Decay Mode Finding", 2, 0, 2);
  _tauNProngsHisto = new TH1F("_tauNProngsHisto", "Tau NProngs", 5, 0, 5);
  _tauMuonDiscHisto = new TH1F("_tauMuonDiscHisto", "Tau Muon Disc", 3, 0, 3);
  _tauElecDiscHisto = new TH1F("_tauElecDiscHisto", "Tau Elec Disc", 5, 0, 5);
  _tauElecMVA5DiscHisto = new TH1F("_tauElecMVA5DiscHisto", "Tau Elec Disc MVA", 6, 0, 6);
  _tauIsoDB3HitsHisto = new TH1F("_tauIsoDB3HitsHisto", "Tau Iso DB", 4, 0, 4);
  _tauIsoMVAHisto = new TH1F("_tauIsoMVAHisto", "Tau Iso MVA", 6, 0, 6);
  _tauGenPtBeforeIsoHisto = new TH1F("_tauGenPtBeforeIsoHisto", "Tau GenPt", 30, 0, 1500);
  _tauGenPtAfterIsoHisto = new TH1F("_tauGenPtAfterIsoHisto", "Tau GenPt", 30, 0, 1500);
  _tauPtBeforeIsoHisto = new TH1F("_tauPtBeforeIsoHisto", "Tau p_{T}", 30, 0, 1500);
  _tauPtAfterIsoHisto = new TH1F("_tauPtAfterIsoHisto", "Tau p_{T}", 30, 0, 1500);
  _tauGenPtBeforeDecayModeHisto = new TH1F("_tauGenPtBeforeDecayModeHisto", "Tau GenPt", 30, 0, 1500);
  _tauGenPtAfterDecayModeHisto = new TH1F("_tauGenPtAfterDecayModeHisto", "Tau GenPt", 30, 0, 1500);
  _tauPtBeforeDecayModeHisto = new TH1F("_tauPtBeforeDecayModeHisto", "Tau Pt", 30, 0, 1500);
  _tauPtAfterDecayModeHisto = new TH1F("_tauPtAfterDecayModeHisto", "Tau Pt", 30, 0, 1500);
  _tauGenPtAfterMuonDiscHisto = new TH1F("_tauGenPtAfterMuonDiscHisto", "Tau GenPt", 30, 0, 1500);
  _tauGenPtAfterElecDiscHisto = new TH1F("_tauGenPtAfterElecDiscHisto", "Tau GenPt", 30, 0, 1500);
  _tauPtAfterMuonDiscHisto = new TH1F("_tauPtAfterMuonDiscHisto", "Tau Pt", 30, 0, 1500);
  _tauPtAfterElecDiscHisto = new TH1F("_tauPtAfterElecDiscHisto", "Tau Pt", 30, 0, 1500);
  _ptResAfterIsoHisto = new TH1F("_ptResAfterIsoHisto", "(Reco Pt - GenPt)/GenPt", 100, -1, 1);
  _ptResBeforeTauIdHisto = new TH1F("_ptResBeforeTauIdHisto", "(Reco Pt - GenPt)/GenPt", 100, -1, 1);
  _elecTauChargeHisto = new TH1F("_elecTauChargeHisto", "Q_{e} * Q_{#tau}", 3, -1, 2);
  _elecTauCosDPhiHisto = new TH1F("_elecTauCosDPhiHisto", "cos#Delta#phi(e,#tau)", 200, -1, 1);
  _metHisto = new TH1F("_metHisto","#slash{E}_{T}", sizeof(ptAndMetBinning)/sizeof(float)-1,ptAndMetBinning);
  //_metHisto = new TH1F("_metHisto","#slash{E}_{T}", 100, 0, 500);
  _diElecAsymmVsdiElecMassDiff = new TH2F("_diElecAsymmVsdiElecMassDiff", "diElec Asymm vs diElec Mass Diff", 100, 0, 1, 100, 0,100);
  _pZetaVsPZetaVis2DHisto = new TH2F("_pZetaVsPZetaVis2DHisto", "p_{#zeta} Vs p_{#zeta}^{Vis}", 100, 0, 100, 200, -100, 100);
  _pZetaHisto = new TH1F("_pZetaHisto", "p_{#zeta}", 50, -50, 50);
  _pZetaVisHisto = new TH1F("_pZetaVisHisto", "p_{#zeta}^{Vis}", 100, 0, 100);
  _nBTagsTCHEMHisto = new TH1F("_nBTagsTCHEHisto", "N Btags TCHEM", 5, 0, 5);
  _nBTagsTCHPTHisto = new TH1F("_nBTagsTCHPTHisto", "N Btags TCHPT", 5, 0, 5);
  _nBTagsCSVLHisto = new TH1F("_nBTagsCSVLHisto", "N Btags CSVL", 5, 0, 5);
  _nBTagsCSVMHisto = new TH1F("_nBTagsCSVMHisto", "N Btags CSVM", 5, 0, 5);
  _nBTagsCSVTHisto = new TH1F("_nBTagsCSVTHisto", "N Btags CSVT", 5, 0, 5);
  _nJetsHisto = new TH1F("_nJetsHisto", "NJets", 10, 0, 10);
  _jetEtSumHisto = new TH1F("_jetEtSumHisto", "Jet Sum E_{T}", 100, 0, 1000);
  _vbfDiJetMassHisto = new TH1F("_vbfDiJetMassHisto", "VBF DiJet Mass", sizeof(vbfMassBinning)/sizeof(float)-1,vbfMassBinning);
  _elecMETCosDPhi = new TH1F("_elecMETCosDPhi", "cos#Delta#phi(e,#slash{E}_{T})", 200, -1, 1);
  _elecMETMassVsElecMETCosDPhi = new TH2F("_elecMETMassVsElecMETCosDPhi", "M(e,#slash{E}_{T}) vs cos#Delta#phi(e,#slash{E}_{T})", 200, -1, 1, 500, 0, 500);
  _elecMETMassVsElecTauCosDPhi = new TH2F("_elecMETMassVsElecTauCosDPhi", "M(e,#slash{E}_{T}) vs cos#Delta#phi(e,#tau)", 200, -1, 1, 500, 0, 500);
  _nVBFJetsHisto = new TH1F("_nVBFJetsHisto", "N VBFJets", 10, 0, 10);
  _elecMETMassHisto = new TH1F("_elecMETMassHisto", "M(e,#slash{E}_{T})", 20, massBinning);
  if(_maxElecTauMassCut <= 500){
    _elecTauMassHisto = new TH1F("_elecTauMassHisto", "M(e,#tau)", (int)(_maxElecTauMassCut - _minElecTauMassCut)/10.,_minElecTauMassCut, _maxElecTauMassCut);
    _elecTauLSMassHisto = new TH1F("_elecTauLSMassHisto", "M(e,#tau) LS", (int)(_maxElecTauMassCut - _minElecTauMassCut)/10.,_minElecTauMassCut, _maxElecTauMassCut);
    _elecTauMetMassHisto = new TH1F("_elecTauMetMassHisto", "M(e,#tau,#slash{E}_{T})", (int)(_maxElecTauMassCut - _minElecTauMassCut)/10.,_minElecTauMassCut, _maxElecTauMassCut);
    _elecTauMetLSMassHisto = new TH1F("_elecTauMetLSMassHisto", "M(e,#tau,#slash{E}_{T}) LS", (int)(_maxElecTauMassCut - _minElecTauMassCut)/10.,_minElecTauMassCut, _maxElecTauMassCut);
  }
  else{
    _elecTauMassHisto = new TH1F("_elecTauMassHisto", "M(e,#tau)", sizeof(massBinning)/sizeof(float)-1,massBinning);
    _elecTauLSMassHisto = new TH1F("_elecTauLSMassHisto", "M(e,#tau) LS", sizeof(massBinning)/sizeof(float)-1,massBinning);
    _elecTauMetMassHisto = new TH1F("_elecTauMetMassHisto", "M(e,#tau,#slash{E}_{T})", sizeof(massBinning)/sizeof(float)-1,massBinning);
    _elecTauMetLSMassHisto = new TH1F("_elecTauMetLSMassHisto", "M(e,#tau,#slash{E}_{T}) LS", sizeof(massBinning)/sizeof(float)-1,massBinning);
  }
  _nBTagsTCHEMTTJetsCRHisto = new TH1F("_nBTagsTCHETTJetsCRHisto", "N Btags TCHEM (TTJets CR)", 5, 0, 5);
  _nBTagsTCHPTTTJetsCRHisto = new TH1F("_nBTagsTCHPTTTJetsCRHisto", "N Btags TCHPT (TTJets CR)", 5, 0, 5);
  _nBTagsCSVLTTJetsCRHisto = new TH1F("_nBTagsCSVLTTJetsCRHisto", "N Btags CSVL (TTJets CR)", 5, 0, 5);
  _nBTagsCSVMTTJetsCRHisto = new TH1F("_nBTagsCSVMTTJetsCRHisto", "N Btags CSVM (TTJets CR)", 5, 0, 5);
  _nBTagsCSVTTTJetsCRHisto = new TH1F("_nBTagsCSVTTTJetsCRHisto", "N Btags CSVT (TTJets CR)", 5, 0, 5);
  _tauPtTTbarCR1Histo = new TH1F("_tauPtTTbarCR1Histo", " Tau p_{T} (TTbarCR1)", 100, 0, 500);
  _tauPtTTbarCR1HistoNoTauId  = new TH1F("_tauPtTTbarCR1HistoNoTauId", " Tau p_{T} (TTbarCR1 No TauId)", 100, 0, 500);
  _tauPtAfterAllCutsHisto = new TH1F("_tauPtAfterAllCutsHisto", "Tau Pt", sizeof(ptAndMetBinning)/sizeof(float)-1, ptAndMetBinning);
  _elecPtAfterAllCutsHisto = new TH1F("_elecPtAfterAllCutsHisto", "Elec P_{T}", sizeof(ptAndMetBinning)/sizeof(float)-1, ptAndMetBinning);
  _metAfterAllCutsHisto = new TH1F("_metAfterAllCutsHisto", "#slash{E}_{T}", sizeof(ptAndMetBinning)/sizeof(float)-1, ptAndMetBinning);
  
  _elecPtNMinus1Histo = new TH1F("_elecPtNMinus1Histo", " Elec P_{T} (N-1)", sizeof(ptAndMetBinning)/sizeof(float)-1, ptAndMetBinning);
  _elecEtaNMinus1Histo = new TH1F("_elecEtaNMinus1Histo", " Elec #eta (N-1)", 100, -5., 5.);
  _elecDPhiInEBNMinus1Histo = new TH1F("_elecDPhiInEBNMinus1Histo", " Elec #Delta#phi_{in} EB (N-1)", 100, -0.5, 0.5);
  _elecDPhiInEENMinus1Histo = new TH1F("_elecDPhiInEENMinus1Histo", " Elec #Delta#phi_{in} EE (N-1)", 100, -0.5, 0.5);
  _elecDEtaInEBNMinus1Histo = new TH1F("_elecDEtaInEBNMinus1Histo", " Elec #Delta#eta_{in} EB (N-1)", 200, -0.1, 0.1);
  _elecDEtaInEENMinus1Histo = new TH1F("_elecDEtaInEENMinus1Histo", " Elec #Delta#eta_{in} EE (N-1)", 200, -0.1, 0.1);
  _elecSigmaIEtaIEtaEBNMinus1Histo = new TH1F("_elecSigmaIEtaIEtaEBNMinus1Histo", " Elec #sigma_{i#eta i#eta} EB (N-1)", 100, 0, 0.1); 
  _elecSigmaIEtaIEtaEENMinus1Histo = new TH1F("_elecSigmaIEtaIEtaEENMinus1Histo", " Elec #sigma_{i#eta i#eta} EE (N-1)", 100, 0, 0.1);
  _elecHOverEEBNMinus1Histo = new TH1F("_elecHOverEEBNMinus1Histo", " Elec H/E EB (N-1)", 200, 0, 2);
  _elecHOverEEENMinus1Histo = new TH1F("_elecHOverEEENMinus1Histo", " Elec H/E EE (N-1)", 200, 0, 2);
  _elecOneOverEMinusOneOverPEBNMinus1Histo = new TH1F("_elecOneOverEMinusOneOverPEBNMinus1Histo", " Elec 1/E - 1/p EB (N-1)", 100, -0.5, 0.5);
  _elecOneOverEMinusOneOverPEENMinus1Histo = new TH1F("_elecOneOverEMinusOneOverPEENMinus1Histo", " Elec 1/E - 1/p EE (N-1)", 100, -0.5, 0.5);
  _elecDxyVtxNMinus1Histo = new TH1F("_elecDxyVtxNMinus1Histo", " Elec d_{xy} (N-1)", 200, -2, 2);
  _elecDzVtxNMinus1Histo = new TH1F("_elecDzVtxNMinus1Histo", " Elec d_{z} (N-1)", 200, -2, 2);
  _elecMissNMinus1Histo = new TH1F("_elecMissNMinus1Histo", " Elec Missing Hits (N-1)", 10, 0, 10);
  _elecPFRelIsoDR03NMinus1Histo = new TH1F("_elecPFRelIsoDR03NMinus1Histo", " Elec PFRelIso #Delta#beta #Delta R = 0.3 (N-1)", 250, 0, 10);
  _elecPFRelIsoDR04NMinus1Histo = new TH1F("_elecPFRelIsoDR04NMinus1Histo", " Elec PFRelIso #Delta#beta #Delta R = 0.4 (N-1)", 250, 0, 10);
  
  _elecTauDRNMinus1Histo = new TH1F("_elecTauDRNMinus1Histo", " #Delta R(e,#tau) (N-1)", 100, 0, 6.28);
  _tauPtNMinus1Histo = new TH1F("_tauPtNMinus1Histo", " Tau p_{T} (N-1)", sizeof(ptAndMetBinning)/sizeof(float)-1, ptAndMetBinning);
  _tauAltLVPtNMinus1Histo = new TH1F("_tauAltLVPtNMinus1Histo", " Tau p_{T} (N-1)", sizeof(ptAndMetBinning)/sizeof(float)-1, ptAndMetBinning);
  _tauEtaNMinus1Histo = new TH1F("_tauEtaNMinus1Histo", " Tau #eta (N-1)",  100, -5., 5.);
  _tauDecayModeNMinus1Histo = new TH1F("_tauDecayModeNMinus1Histo", "Tau Decay Mode Finding (N-1)", 2, 0, 2);
  _tauNProngsNMinus1Histo = new TH1F("_tauNProngsNMinus1Histo", "Tau NProngs (N-1)", 5, 0, 5);
  _tauMuonDiscNMinus1Histo = new TH1F("_tauMuonDiscNMinus1Histo", " Tau Muon Disc (N-1)", 3, 0, 3);
  _tauElecDiscNMinus1Histo = new TH1F("_tauElecDiscNMinus1Histo", " Tau Elec Disc (N-1)", 5, 0, 5);
  _tauElecMVA5DiscNMinus1Histo = new TH1F("_tauElecMVA5DiscNMinus1Histo", " Tau Elec Disc MVA (N-1)", 6, 0, 6);
  _tauIsoDB3HitsNMinus1Histo = new TH1F("_tauIsoDB3HitsNMinus1Histo", " TauIso3Hits #Delta#beta (N-1)", 4, 0, 4);
  _tauIsoMVANMinus1Histo = new TH1F("_tauIsoMVANMinus1Histo", " TauIso MVA (N-1)", 6, 0, 6);
  _elecTauChargeNMinus1Histo = new TH1F("_elecTauChargeNMinus1Histo", " Q_{e} * Q_{#tau} (N-1)", 3, -1, 2);
  _elecTauCosDPhiNMinus1Histo = new TH1F("_elecTauCosDPhiNMinus1Histo", " cos#Delta#phi(e,#tau) (N-1)", 200, -1, 1);
  _metNMinus1Histo = new TH1F("_metNMinus1Histo"," #slash{E}_{T} (N-1)", sizeof(ptAndMetBinning)/sizeof(float)-1, ptAndMetBinning);
  //_metNMinus1Histo = new TH1F("_metNMinus1Histo"," #slash{E}_{T} (N-1)", 100, 0, 500);
  _pZetaVsPZetaVis2DNMinus1Histo = new TH2F("_pZetaVsPZetaVis2DNMinus1Histo", " p_{#zeta} Vs p_{#zeta}^{Vis} (N-1)", 100, 0, 100, 200, -100, 100);
  _pZetaNMinus1Histo = new TH1F("_pZetaNMinus1Histo", " p_{#zeta} (N-1)", 50, -50, 50);
  _pZetaVisNMinus1Histo = new TH1F("_pZetaVisNMinus1Histo", " p_{#zeta}^{Vis} (N-1)", 100, 0, 100);
  stringstream slopeString;
  slopeString << _zeta2DSlope;
  string pzetaHistoTitle = "p_{#zeta} - " + slopeString.str() + "p_{#zeta}^{Vis}";
  _pZetaMinusPZetaVisNMinus1Histo = new TH1F("_pZetaMinusPZetaVisNMinus1Histo", pzetaHistoTitle.c_str(), 200, -100, 100);
  _nBTagsTCHEMNMinus1Histo = new TH1F("_nBTagsTCHENMinus1Histo", " N Btags TCHEM (N-1)", 5, 0, 5);
  _nBTagsTCHPTNMinus1Histo = new TH1F("_nBTagsTCHPTNMinus1Histo", " N Btags TCHPT (N-1)", 5, 0, 5);
  _nBTagsCSVLNMinus1Histo = new TH1F("_nBTagsCSVLNMinus1Histo", " N Btags CSVL (N-1)", 5, 0, 5);
  _nBTagsCSVMNMinus1Histo = new TH1F("_nBTagsCSVMNMinus1Histo", " N Btags CSVM (N-1)", 5, 0, 5);
  _nBTagsCSVTNMinus1Histo = new TH1F("_nBTagsCSVTNMinus1Histo", " N Btags CSVT (N-1)", 5, 0, 5);
  _nJetsNMinus1Histo = new TH1F("_nJetsNMinus1Histo", " NJets (N-1)", 10, 0, 10);
  _jetEtSumNMinus1Histo = new TH1F("_jetEtSumNMinus1Histo", "Jets Sum E_{T}", 100, 0, 1000);
  _nJetsVsNBtagsCSVM2DHisto = new TH2F("_nJetsVsNBtagsCSVM2DHisto", "NJets Vs NBtags CSVM", 5, 0, 5, 10, 0, 10);
  _jetEtSumVsNBtagsCSVM2DHisto = new TH2F("_jetEtSumVsNBtagsCSVM2DHisto", "Jets Sum E_{T} Vs NBtags CSVM", 5, 0, 5, 100, 0, 1000);
  _nVBFJetsNMinus1Histo = new TH1F("_nVBFJetsNMinus1Histo", " N VBF Jets", 10, 0, 10);
  _jet1PtNMinus1Histo = new TH1F("_jet1PtNMinus1Histo", " VBF Jet_{1} p_{T} (N-1)", 100, 0, 1000);
  _jet2PtNMinus1Histo = new TH1F("_jet2PtNMinus1Histo", " VBF Jet_{2} p_{T} (N-1)", 100, 0, 1000);
  _vbfDiJetDEtaNMinus1Histo = new TH1F("_vbfDiJetDEtaNMinus1Histo", " VBF DiJet #Delta#eta (N-1)", 100, 0, 10); 
  _vbfDiJetMassNMinus1Histo = new TH1F("_vbfDiJetMassNMinus1Histo", " VBF DiJet Mass (N-1)", sizeof(vbfMassBinning)/sizeof(float)-1,vbfMassBinning);
  _zeeCRVBFDiJetMassHisto = new TH1F("_zeeCRVBFDiJetMassHisto", " VBF DiJet Mass (N-1)", sizeof(vbfMassBinning)/sizeof(float)-1,vbfMassBinning);
  _ztautauCRVBFDiJetMassHisto = new TH1F("_ztautauCRVBFDiJetMassHisto", " VBF DiJet Mass (N-1)", sizeof(vbfMassBinning)/sizeof(float)-1,vbfMassBinning);
  _wjetsCRVBFDiJetMassHisto = new TH1F("_wjetsCRVBFDiJetMassHisto", " VBF DiJet Mass (N-1)", sizeof(vbfMassBinning)/sizeof(float)-1,vbfMassBinning);
  _ttjetsCRVBFDiJetMassHisto = new TH1F("_ttjetsCRVBFDiJetMassHisto", " VBF DiJet Mass (N-1)", sizeof(vbfMassBinning)/sizeof(float)-1,vbfMassBinning);

}

void ElecTauAnalysis::GetPUWeights(){
  //fisrt get the data distribution
  TFile* dataPUFile = new TFile(_thePileUpRootFileName.c_str());
  dataPUHisto = (TH1F*)dataPUFile->Get("pileup");
  dataPUHisto->Scale(1./dataPUHisto->Integral());
  
  //get MC distribution; copy histo from data
  TH1F* mcPUHisto = (TH1F*)dataPUHisto->Clone("mcPUHisto");
  fChain->Project("mcPUHisto", "TrueNumPU");
  mcPUHisto->Scale(1./mcPUHisto->Integral());
  TH1F* theWeights = (TH1F*)dataPUHisto->Clone("theWeights");
  theWeights->SetName("PUWeights");
  theWeights->Divide(mcPUHisto);
  for(int iBin = 1; iBin < theWeights->GetNbinsX() + 1; iBin++){
    _thePUWeights.push_back(theWeights->GetBinContent(iBin));
  }
}

pair<double, double> ElecTauAnalysis::GetFactorizedEvents(pair<double, double> nLoose, pair<double,double> effLoose, 
                                                      pair<double, double> effTight){
  float normFactor = effTight.first/effLoose.first;
  double normFactorSqrd = normFactor*normFactor;
  double sigmaEffLooseSqrd = effLoose.second*effLoose.second;
  double sigmaEffTightSqrd = effTight.second*effTight.second;
  
  double nLooseSqrd = nLoose.first*nLoose.first;
  double sigmaNLooseSqrd = nLoose.second*nLoose.second;
  
  double sigmaNormFactor = sqrt(sigmaEffTightSqrd + (normFactorSqrd*sigmaEffLooseSqrd)/effLoose.first);
  double normNEvents = nLoose.first * normFactor;
  
  double sigmaNormEvents = sqrt((normFactorSqrd*sigmaNLooseSqrd) +(nLooseSqrd*sigmaNormFactor*sigmaNormFactor));
  pair<double, double> theFactorizedNEvents = make_pair<double, double>(normNEvents, sigmaNormEvents);
  
  //std::cout << "Norm NEvents " << normNEvents << "\t" << sigmaNormEvents << "\n";

  return theFactorizedNEvents;

}