cutClass.cpp~

#include "cutClass.h"
#include <cmath>
#include <math.h>
#include "TLorentzVector.h"
#include <iostream>
#include <libconfig.h++>
#include <iomanip>

Cuts::Cuts(bool doPlots, bool fillCutFlows,bool invertIsoCut, bool lepCutFlow, bool moreMuonIdCuts, bool dumpEventNumber):
  // Set all default parameters. These will be editable later on, probably.
  numTightEle_(3),
  tightElePt_(20.),
  tightEleEta_(2.5),
  tightEled0_(0.04),
  tightEleMissLayers_(0),
  tightEleCheckPhotonVeto_(true),
  tightEleMVA_(0.5),
  tightEleRelIso_(0.15),
  //Loose electron initialisation
  numLooseEle_(3),
  looseElePt_(10),
  looseEleEta_(2.5),
  looseEleMVA_(0.5),
  looseEleRelIso_(0.15),
  //Tight muon initialisation
  numTightMu_(0),
  tightMuonPt_(20),
  tightMuonEta_(2.4),
  tightMuonRelIso_(0.2),
  //Loose muons
  numLooseMu_(0),
  looseMuonPt_(10),
  looseMuonEta_(2.4),
  looseMuonRelIso_(0.2),
  //Jet initialisation
  numJets_(2),
  jetPt_(30.),
  jetEta_(2.5),
  jetNConsts_(2),
  jetIDDo_(true),
  //B-discriminator cut
  numbJets_(1),
  //bDiscCut_(0.679), // Tight level
  bDiscCut_(0.244), // Loose cut
  //Do plots?
  doPlots_(doPlots),
  fillCutFlow_(fillCutFlows),
  //background estimation. May not be possible
  invertIsoCut_(invertIsoCut),
  //Synchronisation cut flow.
  synchCutFlow_(lepCutFlow),
  //Dump event numbers?
  makeEventDump_(dumpEventNumber),
  singleEventInfoDump_(false),
  //Use the tighter muon ID cuts?
  tighterMuonID_(moreMuonIdCuts),
  //Set isMC. Default is true, but it's called everytime a new dataset is processed anyway.
  isMC_(true),
  //Same for trigger flag.
  triggerFlag_(""),
  //Make cloned tree false for now
  postLepSelTree_(0)
{
  //Space here in case other stuff needs to be done.
  //If doing synchronisation., initialise that here.
  if (synchCutFlow_){
    synchCutFlowHist_ = new TH1F("synchCutFlow","synchCutFlow",7,0,7);
    synchNumEles_ = new TH1I("synchNumEles","synchNumEles",10,0,10);
    synchNumMus_ = new TH1I("synchNumMuos","synchNumMuos",10,0,10);
    synchMuonCutFlow_ = new TH1I("synchMuonCutFlow","synchMuonCutFlow",11,0,11);
  }
}

Cuts::~Cuts(){
  if (synchCutFlow_) {
    delete synchCutFlowHist_;
    delete synchNumEles_;
    delete synchNumMus_;
    delete synchMuonCutFlow_;
    if (makeEventDump_){ step0EventDump_.close();
      step2EventDump_.close();
      step4EventDump_.close();
      step6EventDump_.close();
    }
  }
}

bool Cuts::parse_config(std::string confName){
  libconfig::Config cutConf;
  
  //Get the configuration file
  try {
    cutConf.readFile(confName.c_str());
  }
  catch (const libconfig::FileIOException &exep){
    std::cerr << "Error opening cut configuration file" << std::endl;
    return 0;
  }
  catch(const libconfig::ParseException &e){
    std::cerr << "Parse error in cut configuration at " << e.getFile() << ":" << e.getLine() << " - " << e.getError() << std::endl;
    return 0;
  }
  const libconfig::Setting& root = cutConf.getRoot();

  if (! root.exists("cuts")){
    std::cerr << "I don't know what you think you're doing, but there aren't any cuts here. Get back to the drawing board." << std::endl;
    return false;
  }
  if (root.exists("trigLabel")) root.lookupValue("trigLabel",cutConfTrigLabel_);
  if (root.exists("plotPostfix")) root.lookupValue("plotPostfix",postfixName_);
  libconfig::Setting& cuts = root["cuts"];
  if (cuts.exists("tightElectrons")){
    libconfig::Setting& eles = cuts["tightElectrons"];
    eles.lookupValue("pt",tightElePt_);
    eles.lookupValue("eta",tightEleEta_);
    eles.lookupValue("relIso",tightEleRelIso_);
    eles.lookupValue("MVA",tightEleMVA_);
    eles.lookupValue("number",numTightEle_);
  }
  if (cuts.exists("looseElectrons")){
    libconfig::Setting& eles = cuts["looseElectrons"];
    eles.lookupValue("pt",looseElePt_);
    eles.lookupValue("eta",looseEleEta_);
    eles.lookupValue("relIso",looseEleRelIso_);
    eles.lookupValue("MVA",looseEleMVA_);
    eles.lookupValue("number",numLooseEle_);
  }
  if (cuts.exists("tightMuons")){
    libconfig::Setting& mus = cuts["tightMuons"];
    mus.lookupValue("pt",tightMuonPt_);
    mus.lookupValue("eta",tightMuonEta_);
    mus.lookupValue("relIso",tightMuonRelIso_);
    mus.lookupValue("number",numTightMu_);
  }
  if (cuts.exists("looseMuons")){
    libconfig::Setting& mus = cuts["looseMuons"];
    mus.lookupValue("pt",looseMuonPt_);
    mus.lookupValue("eta",looseMuonEta_);
    mus.lookupValue("relIso",looseMuonRelIso_);
    mus.lookupValue("number",numLooseMu_);
  }
  std::cerr << "And so it's looking for " << numTightMu_ << " muons and " << numTightEle_ << " electrons" << std::endl;
  
  if (makeEventDump_){ 
    step0EventDump_.open("step0EventDump"+postfixName_+".txt");
    step2EventDump_.open("step2EventDump"+postfixName_+".txt");
    step4EventDump_.open("step4EventDump"+postfixName_+".txt");
    step6EventDump_.open("step6EventDump"+postfixName_+".txt");
  }

  
  return true;
}

bool Cuts::makeCuts(AnalysisEvent *event, float eventWeight, std::map<std::string,Plots*> plotMap, TH1F* cutFlow){

  //If we're doing synchronisation, do this function.
  if (synchCutFlow_){
    return synchCuts(event);
  }
  
  if (!isMC_) if (!triggerCuts(event)) return false;

  //Make lepton cuts. Does the inverted iso cuts if necessary.
  //  if (!(invertIsoCut_?invertIsoCut(event,eventWeight, plotMap,cutFlow):invertIsoCut(event,eventWeight, plotMap,cutFlow))) return false;
  if (!makeLeptonCuts(event,eventWeight,plotMap,cutFlow)) return false;

  //This is to make some skims for faster running. Do lepSel and save some files.
  if(postLepSelTree_) {
    postLepSelTree_->Fill();
  }

  event->jetIndex = makeJetCuts(event);
  if (doPlots_) plotMap["zMass"]->fillAllPlots(event,eventWeight);
  if (event->jetIndex.size() < numJets_) return false;
  

  if (doPlots_||fillCutFlow_) cutFlow->Fill(2.5,eventWeight);

  event->bTagIndex = makeBCuts(event,event->jetIndex);
  if (doPlots_) plotMap["jetSel"]->fillAllPlots(event,eventWeight);
  if (event->bTagIndex.size() < numbJets_) return false;

  if (doPlots_) plotMap["bTag"]->fillAllPlots(event,eventWeight);
  if (doPlots_||fillCutFlow_) cutFlow->Fill(3.5,eventWeight);


  return true;
}

//Make lepton cuts. Will become customisable in a config later on.
bool Cuts::makeLeptonCuts(AnalysisEvent* event,float eventWeight,std::map<std::string,Plots*> plotMap, TH1F* cutFlow){

  //Do lepton selection. 
  event->electronIndexTight = getTightEles(event);
  if (event->electronIndexTight.size() != numTightEle_) return false;
  event->electronIndexLoose = getLooseEles(event);
  if (event->electronIndexLoose.size() != numLooseEle_) return false;
  event->muonIndexTight = getTightMuons(event);
  if (event->muonIndexTight.size() != numTightMu_) return false;
  event->muonIndexLoose = getLooseMuons(event);
  if (event->muonIndexLoose.size() != numLooseMu_) return false;


  //Should I make it return which leptons are the zMass candidate? Probably.
  float invMass = getZCand(event, event->electronIndexTight, event->muonIndexTight);
  if(doPlots_) plotMap["lepSel"]->fillAllPlots(event,eventWeight);
  if(doPlots_||fillCutFlow_){
    cutFlow->Fill(0.5,eventWeight);
  }

  if (fabs(invMass) > 10.) return false;

  //  plotMap["zMass"]->fillAllPlots(event,eventWeight);
  if(doPlots_||fillCutFlow_) cutFlow->Fill(1.5,eventWeight);

  
  return true;
}

std::vector<int> Cuts::getTightEles(AnalysisEvent* event) {
  std::vector<int> electrons;
  for (int i = 0; i < event->numElePF2PAT; i++){
    if (!event->elePF2PATIsGsf[i]) continue;
    TLorentzVector tempVec(event->elePF2PATGsfPx[i],event->elePF2PATGsfPy[i],event->elePF2PATGsfPz[i],event->elePF2PATGsfE[i]);
    /*    std::cout << tempVec.Pt();
    std::cout << " | " << event->elePF2PATComRelIsoRho[i]/tempVec.Pt();
    std::cout << " | " << event->elePF2PATRhoIso[i];
    std::cout << " | " << event->elePF2PATAEff03[i];
    std::cout << " | " << event->elePF2PATChHadIso[i];
    std::cout << " | " << event->elePF2PATNtHadIso[i];
    std::cout << " | " << event->elePF2PATGammaIso[i] << std::endl;*/

    if (tempVec.Pt() < tightElePt_) continue;
    if (std::abs(tempVec.Eta()) > tightEleEta_)continue;
    if (std::abs(event->elePF2PATD0PV[i]) > tightEled0_)continue;
    if (event->elePF2PATMissingInnerLayers[i] > tightEleMissLayers_) continue;
    if (!event->elePF2PATPhotonConversionVeto[i] && tightEleCheckPhotonVeto_)continue;
    if (event->elePF2PATMVA[i] < tightEleMVA_)continue;
    if (event->elePF2PATComRelIsoRho[i]/tempVec.Pt() > tightEleRelIso_)continue;
    
    electrons.push_back(i);
  }
  return electrons;
}

std::vector<int> Cuts::getLooseEles(AnalysisEvent* event){
  std::vector<int> electrons;
  for (int i = 0; i < event->numElePF2PAT; i++){
    TLorentzVector tempVec(event->elePF2PATGsfPx[i],event->elePF2PATGsfPy[i],event->elePF2PATGsfPz[i],event->elePF2PATGsfE[i]);
    if (tempVec.Pt() < looseElePt_) continue;
    if (std::abs(tempVec.Eta()) > looseEleEta_)continue;
    if (event->elePF2PATMVA[i] < looseEleMVA_)continue;
    if (event->elePF2PATComRelIsoRho[i]/tempVec.Pt() > looseEleRelIso_)continue;    

    electrons.push_back(i);
  }
  return electrons;
}

std::vector<int> Cuts::getTightMuons(AnalysisEvent* event){
  std::vector<int> muons;
  for (int i = 0; i < event->numMuonPF2PAT; i++){
    //if (i == 0) std::cout << "Starts doing first..." << event->muonPF2PATIsPFMuon[i];
    //if (i > 0) std::cout << "   " << event->muonPF2PATIsPFMuon[i];
    //    std::cout << i << "\t" << event->muonPF2PATPt[i] << "\t" << event->muonPF2PATEta[i] << "\t" << event->muonPF2PATComRelIsodBeta[i] << "\t" << event->muonPF2PATChi2[i]/event->muonPF2PATNDOF[i] << "\t" << event->muonPF2PATGlobalID[i] << "\t" << event->muonPF2PATIsPFMuon[i] << "\t" << event->muonPF2PATTrackID[i] << "\t" <<  event->muonPF2PATTkLysWithMeasurements[i] << "\t" << event->muonPF2PATDBPV[i] << "\t" << event->muonPF2PATTrackNHits[i] << "\t" << event->muonPF2PATMuonNHits[i] << "\t" << event->muonPF2PATVldPixHits[i] << "\t" << event->muonPF2PATMatchedStations[i] << "\t" << std::abs(event->pvZ - event->muonPF2PATVertZ[i]) << std::endl;
    if(!event->muonPF2PATIsPFMuon[i]) continue;
    if (!event->muonPF2PATTrackID[i]) continue;
    if (!event->muonPF2PATGlobalID[i]) continue;
    //    if (!event->muonPF2PATIsPFMuon[i]) continue;
    if (event->muonPF2PATPt[i] < tightMuonPt_) continue;
    if (std::abs(event->muonPF2PATEta[i]) > tightMuonEta_) continue;
    if (event->muonPF2PATComRelIsodBeta[i] > tightMuonRelIso_) continue;

    //Do a little test of muon id stuff here.
    if (event->muonPF2PATChi2[i] < 0.) continue;
    if (event->muonPF2PATChi2[i]/event->muonPF2PATNDOF[i] >= 10.) continue;
    //    if (std::abs(event->muonPF2PATDBInnerTrackD0[i]) > 0.2) continue;
    //if (event->muonPF2PATNChambers[i] < 2) continue;
    //if (i == 0) std::cout << "gets to tighter ";
    //if (i == 0) std::cout << "First muon ";
    //if (i > 0) std::cout << "Checking second muon";
    
    if (event->muonPF2PATTkLysWithMeasurements[i] < 6) continue;
    if (fabs(event->muonPF2PATDBPV[i]) >= 0.02) continue;
    //      if (event->muonPF2PATTrackNHits[i] < 11) continue;
    if (event->muonPF2PATMuonNHits[i] < 1) continue;
    if (event->muonPF2PATVldPixHits[i] < 1) continue;
    if (event->muonPF2PATMatchedStations[i] < 2) continue;
    if (std::abs(event->pvZ - event->muonPF2PATVertZ[i]) > 0.5) continue;
    //if(i == 0) std::cout << "does first ";
    //if (i > 0) std::cout << "allows second muon";
    muons.push_back(i);
  }
  //  std::cout << muons.size() << std::endl;
  return muons;
}

std::vector<int> Cuts::getLooseMuons(AnalysisEvent* event){
  std::vector<int> muons;
  for (int i = 0; i < event->numMuonPF2PAT; i++){
    if (!event->muonPF2PATGlobalID[i] && !event->muonPF2PATTrackID[i]) continue;
    if (event->muonPF2PATPt[i] < looseMuonPt_) continue;
    if (std::abs(event->muonPF2PATEta[i]) > looseMuonEta_) continue;
    if (event->muonPF2PATComRelIsodBeta[i] > looseMuonRelIso_) continue;

    muons.push_back(i);
  }
  return muons;
}

float Cuts::getZCand(AnalysisEvent *event, std::vector<int> electrons, std::vector<int> muons){
  float closestMass = 9999.;
  //Use electrons if there are at least 2. Otherwise use muons.
  if (electrons.size() > 1){
    for (unsigned int i = 0; i < electrons.size(); i++){
      for (unsigned int j = i + 1; j < electrons.size(); j++) {
	if (event->elePF2PATCharge[electrons[i]] * event->elePF2PATCharge[electrons[j]] > 0) continue;
	TLorentzVector lepton1 = TLorentzVector(event->elePF2PATGsfPx[electrons[i]],event->elePF2PATGsfPy[electrons[i]],event->elePF2PATGsfPz[electrons[i]],event->elePF2PATGsfE[electrons[i]]);
	TLorentzVector lepton2 = TLorentzVector(event->elePF2PATGsfPx[electrons[j]],event->elePF2PATGsfPy[electrons[j]],event->elePF2PATGsfPz[electrons[j]],event->elePF2PATGsfE[electrons[j]]);
	float invMass = (lepton1 + lepton2).M() -91;
	if (fabs(invMass) < fabs(closestMass)){
	  // set up the tlorentz vectors in the event. For plotting and jazz.
	  event->zPairLeptons.first = lepton1.Pt() > lepton2.Pt()?lepton1:lepton2;
	  event->zPairRelIso.first = lepton1.Pt() > lepton2.Pt()?event->elePF2PATComRelIsoRho[electrons[i]]/lepton1.Pt():event->elePF2PATComRelIsoRho[electrons[j]]/lepton2.Pt();
	  event->zPairRelIso.second = lepton1.Pt() > lepton2.Pt()?event->elePF2PATComRelIsoRho[electrons[j]]/lepton2.Pt():event->elePF2PATComRelIsoRho[electrons[i]]/lepton1.Pt();
	  event->zPairLeptons.second = lepton1.Pt() > lepton2.Pt()?lepton2:lepton1;
	  closestMass = invMass;
	  //Now set up W lepton.
	  if (electrons.size() ==2) {
	    event->wLepton = TLorentzVector(event->muonPF2PATPX[muons[0]],event->muonPF2PATPY[muons[0]],event->muonPF2PATPZ[muons[0]],event->muonPF2PATE[muons[0]]);
	    event->wLeptonRelIso = event->muonPF2PATComRelIsodBeta[muons[0]];
	  }
	  else {
	    for (unsigned int k = 0; k < electrons.size(); k++){
	      if (k == i || k == j) continue;
	      event->wLepton = TLorentzVector(event->elePF2PATGsfPx[electrons[k]],event->elePF2PATGsfPy[electrons[k]],event->elePF2PATGsfPz[electrons[k]],event->elePF2PATGsfE[electrons[k]]);
	      event->wLeptonRelIso = event->elePF2PATComRelIsoRho[electrons[k]]/event->wLepton.Pt();
	    }
	  }
	}
      }
    }
  } else {
    for (unsigned int i = 0; i < muons.size(); i++){
      for (unsigned int j = i + 1; j < muons.size(); j++) {
	if (event->muonPF2PATCharge[muons[i]] * event->muonPF2PATCharge[muons[j]] > 0) continue;
	TLorentzVector lepton1 = TLorentzVector(event->muonPF2PATPX[muons[i]],event->muonPF2PATPY[muons[i]],event->muonPF2PATPZ[muons[i]],event->muonPF2PATE[muons[i]]);
	TLorentzVector lepton2 = TLorentzVector(event->muonPF2PATPX[muons[j]],event->muonPF2PATPY[muons[j]],event->muonPF2PATPZ[muons[j]],event->muonPF2PATE[muons[j]]);
	float invMass = (lepton1 + lepton2).M() -91;
	if (fabs(invMass) < fabs(closestMass)){
	  // set up the tlorentz vectors in the event. For plotting and jazz.
	  event->zPairLeptons.first = lepton1;
	  event->zPairLeptons.second = lepton2;
	  event->zPairRelIso.first = event->muonPF2PATComRelIsodBeta[muons[i]];
	  event->zPairRelIso.second = event->muonPF2PATComRelIsodBeta[muons[j]];
	  closestMass = invMass;
	  //Now set up W lepton.
	  if (muons.size() ==2){
	    event->wLepton = TLorentzVector(event->elePF2PATGsfPx[electrons[0]],event->elePF2PATGsfPy[electrons[0]],event->elePF2PATGsfPz[electrons[0]],event->elePF2PATGsfE[electrons[0]]);
	    event->wLeptonRelIso = event->elePF2PATComRelIsoRho[electrons[0]]/event->wLepton.Pt();
	  }
	  else {
	    for (unsigned int k = 0; k < muons.size(); k++){
	      if (k == i || k == j) continue;
	      event->wLepton = TLorentzVector(event->muonPF2PATPX[muons[k]],event->muonPF2PATPY[muons[k]],event->muonPF2PATPZ[muons[k]],event->muonPF2PATE[muons[k]]);
	      event->wLeptonRelIso = event->muonPF2PATComRelIsodBeta[muons[k]];
	    }
	  }
	}
      }
    }
  }
  return closestMass;
}

std::vector<int> Cuts::makeJetCuts(AnalysisEvent *event){
  
  std::vector<int> jets;
  //  std::cout << event->eventNum << std::endl << "Jets: " << std::endl;
  for (int i = 0; i < event->numJetPF2PAT; i++){
    //if (std::sqrt(event->jetPF2PATPx[i] * event->jetPF2PATPx[i] + event->jetPF2PATPy[i] * event->jetPF2PATPy[i]) < jetPt_) continue;
    if (event->jetPF2PATPtRaw[i] < jetPt_) continue;
    if (std::abs(event->jetPF2PATEta[i]) > jetEta_) continue;
    if (event->jetPF2PATNConstituents[i] < jetNConsts_) continue;
    //std::cerr << ((event->jetPF2PATNeutralHadronEnergyFractionCorr[i] < 0.99 && event->jetPF2PATNeutralEmEnergyFractionCorr[i] < 0.99)) << "\t" << std::abs(event->jetPF2PATEta[i]) << "\t" <<  (event->jetPF2PATChargedEmEnergyFraction[i] < 0.99) <<  "\t" << (event->jetPF2PATChargedHadronEnergyFraction[i] > 0.) << "\t" << (event->jetPF2PATChargedMultiplicity[i] > 0.)  << "\t" << ((event->jetPF2PATNeutralHadronEnergyFractionCorr[i] < 0.99 && event->jetPF2PATNeutralEmEnergyFractionCorr[i] < 0.99) && ((std::abs(event->jetPF2PATEta[i]) > 2.4) || (event->jetPF2PATChargedEmEnergyFraction[i] < 0.99 && event->jetPF2PATChargedHadronEnergyFraction[i] > 0. && event->jetPF2PATChargedMultiplicity[i] > 0.))) << std::endl;
    if (jetIDDo_ && !((event->jetPF2PATNeutralHadronEnergyFractionCorr[i] < 0.99 && event->jetPF2PATNeutralEmEnergyFractionCorr[i] < 0.99) && ((std::abs(event->jetPF2PATEta[i]) > 2.4) || (event->jetPF2PATChargedEmEnergyFraction[i] < 0.99 && event->jetPF2PATChargedHadronEnergyFraction[i] > 0. && event->jetPF2PATChargedMultiplicity[i] > 0.)))) continue;
    double deltaLep = 10000.;
    //Testing out if electron only cleaning works.
    /*    if (event->electronIndexTight.size() > 1){
      if (deltaLep > deltaR(event->zPairLeptons.first.Eta(),event->zPairLeptons.first.Phi(),event->jetPF2PATEta[i],event->jetPF2PATPhi[i]))
      deltaLep = deltaR(event->zPairLeptons.first.Eta(),event->zPairLeptons.first.Phi(),event->jetPF2PATEta[i],event->jetPF2PATPhi[i]);
    if (deltaLep > deltaR(event->zPairLeptons.second.Eta(),event->zPairLeptons.second.Phi(),event->jetPF2PATEta[i],event->jetPF2PATPhi[i]))
      deltaLep = deltaR(event->zPairLeptons.second.Eta(),event->zPairLeptons.second.Phi(),event->jetPF2PATEta[i],event->jetPF2PATPhi[i]);
    }
    else if (event->electronIndexTight.size() > 0){
      if (deltaLep > deltaR(event->wLepton.Eta(),event->wLepton.Phi(),event->jetPF2PATEta[i],event->jetPF2PATPhi[i]))
	deltaLep = deltaR(event->wLepton.Eta(),event->wLepton.Phi(),event->jetPF2PATEta[i],event->jetPF2PATPhi[i]);
	}*/
    //    else { deltaLep = 10000.;}
    if (deltaLep > deltaR(event->zPairLeptons.first.Eta(),event->zPairLeptons.first.Phi(),event->jetPF2PATEta[i],event->jetPF2PATPhi[i]))
      deltaLep = deltaR(event->zPairLeptons.first.Eta(),event->zPairLeptons.first.Phi(),event->jetPF2PATEta[i],event->jetPF2PATPhi[i]);
    if (deltaLep > deltaR(event->zPairLeptons.second.Eta(),event->zPairLeptons.second.Phi(),event->jetPF2PATEta[i],event->jetPF2PATPhi[i]))
      deltaLep = deltaR(event->zPairLeptons.second.Eta(),event->zPairLeptons.second.Phi(),event->jetPF2PATEta[i],event->jetPF2PATPhi[i]);
    if (deltaLep > deltaR(event->wLepton.Eta(),event->wLepton.Phi(),event->jetPF2PATEta[i],event->jetPF2PATPhi[i]))
      deltaLep = deltaR(event->wLepton.Eta(),event->wLepton.Phi(),event->jetPF2PATEta[i],event->jetPF2PATPhi[i]);
    //std::cout << event->jetPF2PATPtRaw[i] << " " << deltaLep << std::endl;
    if (deltaLep < 0.5) continue;
    //    if (event->jetPF2PATdRClosestLepton[i] < 0.5) continue;
    jets.push_back(i);
  }
  return jets;
}

std::vector<int> Cuts::makeBCuts(AnalysisEvent* event, std::vector<int> jets){
  
  std::vector<int> bJets;

  for (unsigned int i = 0; i < jets.size(); i++){
    if (singleEventInfoDump_) std::cout << event->jetPF2PATPtRaw[jets[i]] << " " << event->jetPF2PATBDiscriminator[jets[i]] << std::endl;
    if (event->jetPF2PATBDiscriminator[jets[i]] < bDiscCut_) continue;

    bJets.push_back(i);
  }
  return bJets;
}

void Cuts::setTightEle(float,float,float)
{
}

//This is only called if the thing is data. There's also a little clause to run over certain triggers if they exist. Because I failed miserably to get them all first time through processing...
bool Cuts::triggerCuts(AnalysisEvent* event){
  //MuEG triggers 
  bool muEGTrig = false;
  if (event->HLT_Mu17_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_v4 > 0 || event->HLT_Mu17_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_v5 > 0 || event->HLT_Mu17_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_v6 > 0 || event->HLT_Mu17_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_v7 > 0 || event->HLT_Mu17_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_v8 > 0 || event->HLT_Mu17_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_v9 > 0 ) muEGTrig = true;
  //else if (cutConfTrigLabel_.find("d") != std::string::npos) muEGTrig = true;
  //double electron triggers
  bool eeTrig = false;
  if (event->HLT_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_v15 > 0 || event->HLT_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_v16 > 0 || event->HLT_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_v17 > 0 || event->HLT_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_v18 ==1 || event->HLT_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_v19 > 0) eeTrig = true;
  //  else if (event->HLT_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_v11 > 0) eeTrig = true;

  //double muon triggers
  bool mumuTrig = false;
  if (!isMC_) {if (event->HLT_Mu17_Mu8_v14 > 0 || event->HLT_Mu17_Mu8_v15 > 0 || event->HLT_Mu17_Mu8_v16 > 0 || event->HLT_Mu17_Mu8_v17 > 0 || event->HLT_Mu17_Mu8_v18 > 0 || event->HLT_Mu17_Mu8_v19 > 0 || event->HLT_Mu17_Mu8_v20 > 0 || event->HLT_Mu17_Mu8_v21 > 0 || event->HLT_Mu17_Mu8_v22 > 0) mumuTrig = true;}
  else if (event->HLT_Mu17_Mu8_v12 > 0 || event->HLT_Mu17_TkMu8_v5 > 0) mumuTrig = true;
  
  if (cutConfTrigLabel_.find("d") != std::string::npos){if (muEGTrig) return true;}
  if (cutConfTrigLabel_.find("e") != std::string::npos){if (eeTrig && !(muEGTrig || mumuTrig)) return true;}
  if (cutConfTrigLabel_.find("m") != std::string::npos){if (mumuTrig && !(eeTrig || muEGTrig)) return true;}
  return false;
}

//Does the cuts required for the synchronisation.
bool Cuts::synchCuts(AnalysisEvent* event){
  //Trigger stuff would go first, but in MC at least (what I'm starting with) I don't have that saved. Idiot.
  
  if (singleEventInfoDump_){
    std::cout << std::setprecision(6) << std::fixed;
  }
  
  int looseLeps = getLooseLepsNum(event);
  if (isMC_ && looseLeps < 2) return false;
  if (!isMC_ && looseLeps < 3) return false;
  if (!isMC_ && !triggerCuts(event)) return false;
  synchCutFlowHist_->Fill(0.5);
  if (makeEventDump_) {step0EventDump_ << event->eventNum << std::endl;}
  
  event->electronIndexTight = getTightEles(event);
  event->muonIndexTight = getTightMuons(event);
  synchNumEles_->Fill(event->electronIndexTight.size());
  synchNumMus_->Fill(event->muonIndexTight.size());
  //If electrons are expected to be the Z, check there are an oppositely charged pair.
  bool zCand = false;
  if (numTightEle_ > 1){
    for (unsigned int i = 0; i < event->electronIndexTight.size();++i){
      for (unsigned int j = i +1; j < event->electronIndexTight.size();++j){
	if (event->elePF2PATCharge[event->electronIndexTight[i]] * event->elePF2PATCharge[event->electronIndexTight[j]] > 0) continue;
	zCand = true;
      }
    }

  }
  else if (numTightMu_ > 1){
    for (unsigned int i = 0; i < event->muonIndexTight.size();++i){
      for (unsigned int j = i +1; j < event->muonIndexTight.size();++j){
	if (event->muonPF2PATCharge[event->muonIndexTight[i]] * event->muonPF2PATCharge[event->muonIndexTight[j]] > 0) continue;
	zCand = true;
      }
    }
  }
  if (!zCand) return false;
  synchCutFlowHist_->Fill(1.5); 
  if (singleEventInfoDump_) std::cout << "Gets past z cand stuff with " << event->electronIndexTight.size() << " tight electrons and " << event->muonIndexTight.size() << " tight muons" << std::endl;
  if (event->electronIndexTight.size() != numTightEle_) return false;
  if (event->muonIndexTight.size() != numTightMu_) return false;
  //loose lepton veto
  if (singleEventInfoDump_) std::cout << "Correct number of leptons and loose: " << getLooseEles(event).size() << " " << getLooseMuons(event).size() << std::endl;
  if (event->electronIndexTight.size() != getLooseEles(event).size()) return false;
  if (event->muonIndexTight.size() != getLooseMuons(event).size()) return false;
  if (singleEventInfoDump_) std::cout << " and passes veto too." << std::endl;
  synchCutFlowHist_->Fill(2.5);
  if (makeEventDump_){dumpToFile(event,2);}
  if (singleEventInfoDump_) std::cout << "Z mass: " << getZCand(event,event->electronIndexTight,event->muonIndexTight) << std::endl;
  if (fabs(getZCand(event,event->electronIndexTight,event->muonIndexTight)) > 15.) return false;
  synchCutFlowHist_->Fill(3.5);

  //Add in extra steps here.
  event->jetIndex = makeJetCuts(event);
  if (singleEventInfoDump_) std::cout << "Number of jets: " << event->jetIndex.size() << std::endl;
  if (event->jetIndex.size() < 1) return false;
  if (makeEventDump_){
    dumpToFile(event,4);
  }
  //  std::cout << event->jetIndex.size() << std::endl;
  synchCutFlowHist_->Fill(4.5);
  event->bTagIndex = makeBCuts(event,event->jetIndex);
  if (singleEventInfoDump_) std::cout << "Number of bJets: " << event->bTagIndex.size() << std::endl;
  if (event->bTagIndex.size() < 1) return false;
  synchCutFlowHist_->Fill(5.5);
  if (singleEventInfoDump_) std::cout << "MET: " << event->metPF2PATPt << std::endl;
  if (event->metPF2PATPt < 20.) return false;
  synchCutFlowHist_->Fill(6.5);
  if (singleEventInfoDump_) std::cout << "Passes all cuts! Yay!" << std::endl;
  if (makeEventDump_){
    dumpToFile(event,6);
  }
  if (singleEventInfoDump_) std::cout << std::setprecision(1) << std::fixed;
  return true;
}

TH1F* Cuts::getSynchCutFlow(){
  std::cout << "Eles: " << numTightEle_ << " Muons: " << numTightMu_ << std::endl;
  char const *names[] = {"Trigger","Lep Pair", "3 Leptons", "zMass","1 jet","b-tag","MET"};
  for (unsigned int i = 1; i < 8; ++i){
    std::cout << names[i-1] << ": " << synchCutFlowHist_->GetBinContent(i) << std::endl;
  }
  std::cout << "The number of leptons in the passing trigger events:" << std::endl;
  std::cout << "Leps\tEle\tMuon" << std::endl;
  for (unsigned int i = 0; i < 10; i++){
    std::cout << i << "\t" << synchNumEles_->GetBinContent(i) << "\t" << synchNumMus_->GetBinContent(i) << std::endl;
  }
  char const *labels[] = {"In", "ID", "PtEtaIso", "chi2","tklay","DBPV","TrackHit","MuHits","PixHits","MtStats","DZPV"};
  for (unsigned int i = 1; i < 12; ++i){
    std::cout << labels[i-1] << ": \t" << synchMuonCutFlow_->GetBinContent(i) << std::endl;
  }
  return synchCutFlowHist_;
}

//Method used for the synchronisation. Mimics the IPHC preselection skims.
int Cuts::getLooseLepsNum(AnalysisEvent* event){
  return getLooseElecs(event) + getLooseMus(event);
}

int Cuts::getLooseElecs(AnalysisEvent* event){
  int looseLeps = 0.;
  for (int i = 0; i < event->numElePF2PAT; i++){
    if (!event->elePF2PATIsGsf[i]) continue;
    TLorentzVector tempVec(event->elePF2PATGsfPx[i],event->elePF2PATGsfPy[i],event->elePF2PATGsfPz[i],event->elePF2PATGsfE[i]);
    if (tempVec.Pt() < 15) continue;
    if (std::abs(tempVec.Eta()) > 2.5)continue;
    looseLeps++;
  }
  return looseLeps;
}

int Cuts::getLooseMus(AnalysisEvent* event){
  int looseLeps = 0.;
  for (int i = 0; i < event->numMuonPF2PAT; i++){
    if (!event->muonPF2PATGlobalID[i]||!event->muonPF2PATTrackID[i]) continue;
    if (event->muonPF2PATPt[i] < 15) continue;
    if (std::abs(event->muonPF2PATEta[i]) > 2.5) continue;
    looseLeps++;
  }
  return looseLeps;
}


//First tentative attempt at doing the background isolation. 
bool Cuts::invertIsoCut(AnalysisEvent* event,float eventWeight,std::map<std::string,Plots*> plotMap, TH1F* cutFlow){
  //Check there are exactly 2 tight leptons with the correct isolation cut.
  event->electronIndexTight = getTightEles(event);
  event->muonIndexTight = getTightMuons(event);

  if ((event->electronIndexTight.size() + event->muonIndexTight.size()) != 2) return false;
  if (event->electronIndexTight.size() == 1) return false;
  
  //Check they are a valid zCandidate. (Just call the zCand method here? - No, that won't actually work.)
  float invMass = 100.;
  if (numTightEle_ > 1){
    if (event->electronIndexTight.size() < 2) return false;
    if (event->elePF2PATCharge[event->electronIndexTight[0]] * event->elePF2PATCharge[event->electronIndexTight[1]] > 0) return false;
    event->zPairLeptons.first = TLorentzVector(event->elePF2PATPX[event->electronIndexTight[0]],event->elePF2PATPY[event->electronIndexTight[0]],event->elePF2PATPZ[event->electronIndexTight[0]],event->elePF2PATE[event->electronIndexTight[0]]);
    event->zPairLeptons.second = TLorentzVector(event->elePF2PATPX[event->electronIndexTight[1]],event->elePF2PATPY[event->electronIndexTight[1]],event->elePF2PATPZ[event->electronIndexTight[1]],event->elePF2PATE[event->electronIndexTight[1]]); 
    invMass = (event->zPairLeptons.first + event->zPairLeptons.second).M() -91.2;
  }
  else{
    if (event->muonIndexTight.size() < 2) return false;
    if (event->muonPF2PATCharge[event->muonIndexTight[0]] * event->muonPF2PATCharge[event->muonIndexTight[1]] > 0) return false;
    event->zPairLeptons.first = TLorentzVector(event->muonPF2PATPX[event->muonIndexTight[0]],event->muonPF2PATPY[event->muonIndexTight[0]],event->muonPF2PATPZ[event->muonIndexTight[0]],event->muonPF2PATE[event->muonIndexTight[0]]);
    event->zPairLeptons.second = TLorentzVector(event->muonPF2PATPX[event->muonIndexTight[1]],event->muonPF2PATPY[event->muonIndexTight[1]],event->muonPF2PATPZ[event->muonIndexTight[1]],event->muonPF2PATE[event->muonIndexTight[1]]);
    invMass = (event->zPairLeptons.first + event->zPairLeptons.second).M() -91.2;
  }
  
  //Get rev iso candidates.
  std::vector<int> invIsoEle = getInvIsoEles(event);
  std::vector<int> invIsoMus = getInvIsoMuons(event);
  
  //Check we have the right number of leptons.
  if (event->electronIndexTight.size() + invIsoEle.size() != numTightEle_) return false;
  if (event->muonIndexTight.size() + invIsoMus.size() != numTightMu_) return false;

  //Put extra lepton into W boson thing.
  if (invIsoEle.size() == 1){
    event->wLepton = TLorentzVector(event->elePF2PATPX[invIsoEle[0]],event->elePF2PATPY[invIsoEle[0]],event->elePF2PATPZ[invIsoEle[0]],event->elePF2PATE[invIsoEle[0]]);
    event->wLeptonRelIso = event->elePF2PATComRelIsoRho[invIsoEle[0]]/event->elePF2PATPT[invIsoEle[0]];
  }
  else{
    event->wLepton = TLorentzVector(event->muonPF2PATPX[invIsoMus[0]],event->muonPF2PATPY[invIsoMus[0]],event->muonPF2PATPZ[invIsoMus[0]],event->muonPF2PATE[invIsoMus[0]]);
    event->wLeptonRelIso = event->muonPF2PATComRelIsodBeta[invIsoMus[0]];
  }

  if(doPlots_) plotMap["lepSel"]->fillAllPlots(event,eventWeight);
  if(doPlots_||fillCutFlow_){
    cutFlow->Fill(0.5,eventWeight);
  }

  if (std::abs(invMass) > 10.) return false;
  if(doPlots_||fillCutFlow_) cutFlow->Fill(1.5,eventWeight);
  return true;

}

std::vector<int> Cuts::getInvIsoEles(AnalysisEvent* event) {                                                                                                                                                                                  
  std::vector<int> electrons;
  for (int i = 0; i < event->numElePF2PAT; i++){
    if (event->elePF2PATPT[i] < tightElePt_) continue;
    if (std::abs(event->elePF2PATEta[i]) > tightEleEta_)continue;
    if (std::abs(event->elePF2PATBeamSpotCorrectedTrackD0[i]) > tightEled0_)continue;
    if (event->elePF2PATMissingInnerLayers[i] > tightEleMissLayers_) continue;
    if (!event->elePF2PATPhotonConversionVeto[i] && tightEleCheckPhotonVeto_)continue;
    if (event->elePF2PATMVA[i] < tightEleMVA_)continue;
    if (event->elePF2PATComRelIsoRho[i]/event->elePF2PATPT[i] < tightEleRelIso_)continue;

    electrons.push_back(i);
  }
  return electrons;
}

std::vector<int> Cuts::getInvIsoMuons(AnalysisEvent* event){
  std::vector<int> muons;
  for (int i = 0; i < event->numMuonPF2PAT; i++){
    if (!event->muonPF2PATGlobalID[i] && !event->muonPF2PATTrackID[i]) continue;
    if (event->muonPF2PATPt[i] < tightMuonPt_) continue;
    if (std::abs(event->muonPF2PATEta[i]) > tightMuonEta_) continue;
    if (event->muonPF2PATComRelIsodBeta[i] < tightMuonRelIso_) continue;

    //Do a little test of muon id stuff here.
    if (event->muonPF2PATChi2[i]/event->muonPF2PATNDOF[i] > 10.) continue;
    if (std::abs(event->muonPF2PATDBInnerTrackD0[i]) > 0.2) continue;
    if (event->muonPF2PATNChambers[i] < 2) continue;

    muons.push_back(i);
  }
  return muons;
}

std::vector<int> Cuts::synchTightMuons(AnalysisEvent* event){
  std::vector<int> muons;
  for (int i = 0; i < event->numMuonPF2PAT; i++){
    synchMuonCutFlow_->Fill(0.5);
    //if (i == 0) std::cout << "Starts doing first..." << event->muonPF2PATIsPFMuon[i];
    //if (i > 0) std::cout << "   " << event->muonPF2PATIsPFMuon[i];
    if (!event->muonPF2PATGlobalID[i] && !event->muonPF2PATTrackID[i]) continue;
    synchMuonCutFlow_->Fill(1.5);
    //    if (!event->muonPF2PATIsPFMuon[i]) continue;
    if (event->muonPF2PATPt[i] < tightMuonPt_) continue;
    if (std::abs(event->muonPF2PATEta[i]) > tightMuonEta_) continue;
    if (event->muonPF2PATComRelIsodBeta[i] > tightMuonRelIso_) continue;
    synchMuonCutFlow_->Fill(2.5);
    //Do a little test of muon id stuff here.
    if (event->muonPF2PATChi2[i] < 0.) continue;
    if (event->muonPF2PATChi2[i]/event->muonPF2PATNDOF[i] >= 10.) continue;
    //    if (std::abs(event->muonPF2PATDBInnerTrackD0[i]) > 0.2) continue;
    //if (event->muonPF2PATNChambers[i] < 2) continue;
    //if (i == 0) std::cout << "gets to tighter ";
    synchMuonCutFlow_->Fill(3.5);
    //if (i == 0) std::cout << "First muon ";
    //if (i > 0) std::cout << "Checking second muon";
    //      std::cout << event->muonPF2PATTkLysWithMeasurements[i] << "\t" << event->muonPF2PATDBPV[i] << "\t" << event->muonPF2PATTrackNHits[i] << "\t" << event->muonPF2PATMuonNHits[i] << "\t" << event->muonPF2PATVldPixHits[i] << "\t" << event->muonPF2PATMatchedStations[i] << "\t" << std::abs(event->pvZ - event->muonPF2PATVertZ[i]) << std::endl;
    if (event->muonPF2PATTkLysWithMeasurements[i] < 6) continue;
    synchMuonCutFlow_->Fill(4.5);
    if (event->muonPF2PATDBPV[i] > 0.02) continue;
    synchMuonCutFlow_->Fill(5.5);
    if (event->muonPF2PATTrackNHits[i] < 11) continue;
    synchMuonCutFlow_->Fill(6.5);
    if (event->muonPF2PATMuonNHits[i] < 1) continue;
    synchMuonCutFlow_->Fill(7.5);
    if (event->muonPF2PATVldPixHits[i] < 1) continue;
    synchMuonCutFlow_->Fill(8.5);
    if (event->muonPF2PATMatchedStations[i] < 2) continue;
    synchMuonCutFlow_->Fill(9.5);
    if (std::abs(event->pvZ - event->muonPF2PATVertZ[i]) > 0.5) continue;
    synchMuonCutFlow_->Fill(10.5);
    //if(i == 0) std::cout << "does first ";
    //if (i > 0) std::cout << "allows second muon";
    muons.push_back(i);
  }
  //  std::cout << muons.size() << std::endl;
  return muons;
}

//For synchronisation I am dumping the lepton information here.
void Cuts::dumpLeptonInfo(AnalysisEvent* event){
  std::cout << std::setprecision(6) << std::fixed; 
  std::cerr << std::setprecision(6) << std::fixed; 
  //Dump electron info first
  event->electronIndexTight = getTightEles(event);
  event->muonIndexTight = getTightMuons(event); 
  std::cout << "Electrons: found " << event->electronIndexTight.size() << std::endl;
  for (unsigned int i = 0; i < event->electronIndexTight.size(); i++){
    TLorentzVector tempVec(event->elePF2PATGsfPx[event->electronIndexTight[i]],event->elePF2PATGsfPy[event->electronIndexTight[i]],event->elePF2PATGsfPz[event->electronIndexTight[i]],event->elePF2PATGsfE[event->electronIndexTight[i]]);
    //    std::cout << i << " | " << event->elePF2PATPT[event->electronIndexTight[i]];
    std::cout << " | " << tempVec.Pt();
    //std::cout << " | " << event->elePF2PATEta[event->electronIndexTight[i]];
    std::cout << " | " << tempVec.Eta();
    std::cout << " | " << tempVec.Phi();
    std::cout << " | " << event->elePF2PATPhi[event->electronIndexTight[i]];
    std::cout << " | " << event->elePF2PATD0PV[event->electronIndexTight[i]];
    std::cout << " | " << event->elePF2PATMVA[event->electronIndexTight[i]];
    std::cout << " | " << event->elePF2PATMissingInnerLayers[event->electronIndexTight[i]];
    std::cout << " | " << event->elePF2PATComRelIsoRho[event->electronIndexTight[i]]/tempVec.Pt();
    std::cout << " | " << event->elePF2PATPhotonConversionVeto[event->electronIndexTight[i]];
    std::cout << " | " << event->elePF2PATCharge[event->electronIndexTight[i]];
    /*std::cout << " | " << event->elePF2PATRhoIso[event->electronIndexTight[i]];
    std::cout << " | " << event->elePF2PATAEff03[event->electronIndexTight[i]];
    std::cout << " | " << event->elePF2PATChHadIso[event->electronIndexTight[i]];
    std::cout << " | " << event->elePF2PATNtHadIso[event->electronIndexTight[i]];
    std::cout << " | " << event->elePF2PATGammaIso[event->electronIndexTight[i]];
    std::cout << " | " << event->elePF2PATComRelIsodBeta[event->electronIndexTight[i]];
    std::cout << " | " << event->elePF2PATComRelIso[event->electronIndexTight[i]]; */
    std::cout << std::endl;
  }

  std::cout << "Muons: found " << event->muonIndexTight.size() << std::endl;
  for (unsigned int i = 0; i <  event->muonIndexTight.size(); i++){
    std::cout << i;
    std::cout << " | " << event->muonPF2PATPt[event->muonIndexTight[i]];
    std::cout << " | " << event->muonPF2PATEta[event->muonIndexTight[i]];
    std::cout << " | " << event->muonPF2PATPhi[event->muonIndexTight[i]];
    std::cout << " | " << event->muonPF2PATChi2[event->muonIndexTight[i]]/event->muonPF2PATNDOF[event->muonIndexTight[i]];
    std::cout << " | " << event->muonPF2PATTkLysWithMeasurements[event->muonIndexTight[i]];
    std::cout << " | " << event->muonPF2PATMuonNHits[event->muonIndexTight[i]];
    std::cout << " | " << event->muonPF2PATDBPV[event->muonIndexTight[i]];
    std::cout << " | " << event->pvZ - event->muonPF2PATVertZ[event->muonIndexTight[i]];
    std::cout << " | " << event->muonPF2PATVldPixHits[event->muonIndexTight[i]];
    std::cout << " | " << event->muonPF2PATMatchedStations[event->muonIndexTight[i]];
    std::cout << " | " << event->muonPF2PATComRelIsodBeta[event->muonIndexTight[i]];
    std::cout << " | " << event->muonPF2PATCharge[event->muonIndexTight[i]];
    std::cout << std::endl;
  }
  int numbJets = event->numJetPF2PAT > 4?4:event->numJetPF2PAT;
  std::cout << "Jets: " << event->numJetPF2PAT << std::endl;
  for (int i = 0; i < numbJets; i++){
    std::cout << i;
    //    std::cout << " | " << event->jetPF2PATPt[i];
    //std::cout << " | " << std::sqrt(event->jetPF2PATPx[i] * event->jetPF2PATPx[i] + event->jetPF2PATPy[i] * event->jetPF2PATPy[i]);
    // std::cout << " | " << event->jetPF2PATEt[i];
    std::cout << " | " << event->jetPF2PATPtRaw[i];
    //std::cout << " | " << event->jetPF2PATUnCorPt[i];
    std::cout << " | " << event->jetPF2PATEta[i];
    std::cout << " | " << event->jetPF2PATPhi[i];
    std::cout << " | " << event->jetPF2PATNConstituents[i];
    std::cout << " | " << (event->jetPF2PATNeutralHadronEnergyFractionCorr[i] < 0.99 && event->jetPF2PATNeutralEmEnergyFractionCorr[i] < 0.99) && ((std::abs(event->jetPF2PATEta[i]) > 2.4) || (event->jetPF2PATChargedEmEnergyFraction[i] < 0.99 && event->jetPF2PATChargedHadronEnergyFraction[i] > 0. && event->jetPF2PATChargedMultiplicity[i] > 0.));
    std::cout << " | " << event->jetPF2PATdRClosestLepton[i];
    std::cout << std::endl;
  }
  std::cout << "MET: " << event->metPF2PATEt << " | " << event->metPF2PATPt << std::endl;

  std::cout << std::setprecision(1) << std::fixed; 
  std::cerr << std::setprecision(1) << std::fixed; 
}

void Cuts::dumpLooseLepInfo(AnalysisEvent* event){
  std::cout << std::setprecision(6) << std::fixed; 
  std::cerr << std::setprecision(6) << std::fixed; 

  std::cout << "Electrons: " << event->numElePF2PAT << std::endl;
  for (int i = 0; i < event->numElePF2PAT; i++){

    TLorentzVector tempVec(event->elePF2PATGsfPx[i],event->elePF2PATGsfPy[i],event->elePF2PATGsfPz[i],event->elePF2PATGsfE[i]);
    std::cout << i;
    std::cout << " | " << tempVec.Pt();
    std::cout << " | " << tempVec.Eta();
    std::cout << " | " << event->elePF2PATSCEta[i];
    std::cout << " | " << event->elePF2PATComRelIsoRho[i]/tempVec.Pt();
    std::cout << std::endl;
  }
  std::cout << "Muons: " << event->numMuonPF2PAT << std::endl;
  for (int i = 0; i < event->numMuonPF2PAT; i++){
    std::cout << i;
    std::cout << " | " << event->muonPF2PATPt[i];
    std::cout << " | " << event->muonPF2PATEta[i];
    std::cout << " | " << event->muonPF2PATComRelIsodBeta[i];
    std::cout << " | " << event->muonPF2PATGlobalID[i];
    std::cout << " | " << event->muonPF2PATTrackID[i];
    std::cout << std::endl;
  }

  std::cout << std::setprecision(1) << std::fixed; 
  std::cerr << std::setprecision(1) << std::fixed; 
}

double Cuts::deltaR(float eta1, float phi1, float eta2, float phi2){
  double dEta = eta1-eta2;
  double dPhi = phi1-phi2;
  while (fabs(dPhi) > 3.14159265359){
    dPhi += (dPhi > 0.? -2*3.14159265359:2*3.14159265359);
  }
  //  if(singleEventInfoDump_)  std::cout << eta1 << " " << eta2 << " phi " << phi1 << " " << phi2 << " ds: " << eta1-eta2 << " " << phi1-phi2 << " dR: " << std::sqrt((dEta*dEta)+(dPhi*dPhi)) << std::endl;
  return std::sqrt((dEta*dEta)+(dPhi*dPhi));
}

//For dumping contents of step 4. Bit more complicated than old, so doing it elsewhere.
void Cuts::dumpToFile(AnalysisEvent* event, int step){
  std::vector<TLorentzVector> tempLepVec;
  if (step > 2){
    if (event->zPairLeptons.first.Pt() < event->wLepton.Pt()){
      tempLepVec.push_back(event->wLepton);
      tempLepVec.push_back(event->zPairLeptons.first);
      tempLepVec.push_back(event->zPairLeptons.second);
    }
    else if (event->wLepton.Pt() > event->zPairLeptons.second.Pt()){
      tempLepVec.push_back(event->zPairLeptons.first);
      tempLepVec.push_back(event->wLepton);
      tempLepVec.push_back(event->zPairLeptons.second);
    }
    else {
      tempLepVec.push_back(event->zPairLeptons.first);
      tempLepVec.push_back(event->zPairLeptons.second);
      tempLepVec.push_back(event->wLepton);
    }
  }
  if (step==2){
    int muUsed = 0;
    for (int i = 0; i < event->numElePF2PAT; i++){
      TLorentzVector tempVec(event->elePF2PATGsfPx[i],event->elePF2PATGsfPy[i],event->elePF2PATGsfPz[i],event->elePF2PATGsfE[i]);
      for (int j = muUsed; j < event->numMuonPF2PAT; j++){
	if (tempVec.Pt() > event->muonPF2PATPt[i]){
	  tempLepVec.push_back(tempVec);
	  break;
	}
	else {
	  tempLepVec.push_back(TLorentzVector(event->muonPF2PATPX[j],event->muonPF2PATPY[j],event->muonPF2PATPZ[j],event->muonPF2PATE[j]));
	  muUsed++;
	  if (tempLepVec.size() > 2) break;
	}
	if (tempLepVec.size() > 2) break;
      }
    }
    if (tempLepVec.size() < 3){
      for (int j = 0; j < 3; j++){
	tempLepVec.push_back(TLorentzVector(event->muonPF2PATPX[j],event->muonPF2PATPY[j],event->muonPF2PATPZ[j],event->muonPF2PATE[j]));
      }
    }
  }
  switch (step) {
  case 2:
    step2EventDump_ << event->eventRun << " " << event->eventNum << " ";
    break;
  case 4:
    step4EventDump_ << event->eventRun << " " << event->eventNum << " ";
    break;
  case 6:
    step6EventDump_ << event->eventRun << " " << event->eventNum << " ";
    break;
  }
  for (unsigned int i = 0; i < 3; i++){
    switch (step) {
    case 2:
      step2EventDump_ << tempLepVec[i].Pt() << " " << tempLepVec[i].Eta() << " ";
      break;
    case 4:
      step4EventDump_ << tempLepVec[i].Pt() << " " << tempLepVec[i].Eta() << " ";
      break;
    case 6:
      step6EventDump_ << tempLepVec[i].Pt() << " " << tempLepVec[i].Eta() << " ";
      break;
    }
  }
  event->jetIndex = makeJetCuts(event);
  for (unsigned int i = 0; i < 4; i++){
    switch (step) {
    case 2:
      step2EventDump_ << ((i < event->jetIndex.size())?event->jetPF2PATPtRaw[event->jetIndex[i]]:-666) << " " << ((i < event->jetIndex.size())?event->jetPF2PATEta[event->jetIndex[i]]:-666) << " ";
      break;
    case 4:
      step4EventDump_ << ((i < event->jetIndex.size())?event->jetPF2PATPtRaw[event->jetIndex[i]]:-666) << " " << ((i < event->jetIndex.size())?event->jetPF2PATEta[event->jetIndex[i]]:-666) << " ";
      break;
    case 6:
      step6EventDump_ << ((i < event->jetIndex.size())?event->jetPF2PATPtRaw[event->jetIndex[i]]:-666) << " " << ((i < event->jetIndex.size())?event->jetPF2PATEta[event->jetIndex[i]]:-666) << " ";
      break;
    }
  }
  switch(step){
  case 2:
    step2EventDump_ << event->metPF2PATPt;
    step2EventDump_ << std::endl;
    break;
  case 4:
    step4EventDump_ << event->metPF2PATPt;
    step4EventDump_ << std::endl;
    break;
  case 6:
    step6EventDump_ << event->metPF2PATPt;
    step6EventDump_ << std::endl;
    break;
  }

}