NonMCBackground/src/ElectronFakeRateMod.cc

#include "ElectronFakeRateMod.h"
#include "SelectionFuncs.h"


//needed because of externs
vector<bool>   PFnoPUflag;
Float_t computePFMuonIso(const mithep::Muon *mu, 
                        const mithep::Vertex * vtx,
                        const mithep::Array<mithep::PFCandidate> * pfCandidates,
                        const Double_t dRMax);
                         
Bool_t passElectronDenominatorCuts(const mithep::Electron *ele, const mithep::Vertex * vtx, Int_t DenominatorType);

string IntToString(int i) {
  char temp[100];
  sprintf(temp, "%d", i);
  string str = temp;
  return str;
}

class TTree;
class TFile;
class TString;

const Float_t g_electron_mass = 0.51099892e-3;
const Float_t g_muon_mass = 105.658369e-3;


void mithep::ElectronFakeRateMod::SlaveBegin()
{
  

  denominatorType.push_back(1); 
  denominatorType.push_back(2); 
  denominatorType.push_back(3); 
  denominatorType.push_back(4); 
  
  initElectronIDMVAV1();

  ReqBranch("PFMet",           fPFMet);

  ReqBranch(mithep::Names::gkMvfConversionBrn,  fConversions);
  
  fMuonName     = mithep::Names::gkMuonBrn;
  ReqBranch(fMuonName,      fMuons);
  
  fElectronName     = mithep::Names::gkElectronBrn;
  ReqBranch(fElectronName,      fElectrons);
  
  fPFCandidateName = mithep::Names::gkPFCandidatesBrn;  
  ReqBranch(fPFCandidateName,     fPFCandidates);
  
  fPrimVtxName = mithep::Names::gkPVBrn;          
  ReqBranch(fPrimVtxName,         fPrimVerts);
  
  fPUEnergyDensityName = mithep::Names::gkPileupEnergyDensityBrn;         
  ReqBranch(fPUEnergyDensityName, fPileupEnergyDensity);
  
  fTrigMaskName = mithep::Names::gkHltBitBrn;     
  ReqBranch(fTrigMaskName,  fTrigMask);
  
  fPileupName = mithep::Names::gkPileupInfoBrn;   
  ReqBranch(fPileupName,  fPileup);
 
 
  fOutputFile = new TFile(fOutputName, "RECREATE");
  fOutputTrees = vector<TTree*>(denominatorType.size(),0);
  for (UInt_t denominatorTypeIndex = 0; denominatorTypeIndex < denominatorType.size(); ++denominatorTypeIndex) {
    fOutputTrees[denominatorTypeIndex] = new TTree(TString("denominator_v"+IntToString(denominatorType[denominatorTypeIndex])),TString("denominator_v"+IntToString(denominatorType[denominatorTypeIndex])));
  }

  rlrm.AddJSONFile("/afs/cern.ch/cms/CAF/CMSCOMM/COMM_DQM/certification/Collisions12/8TeV/Prompt/Cert_190456-193557_8TeV_PromptReco_Collisions12_JSON.txt"); 

  fake_electron_pt = vector<Float_t>(denominatorType.size(),0);
  fake_electron_eta = vector<Float_t>(denominatorType.size(),0);
  fake_electron_phi = vector<Float_t>(denominatorType.size(),0);
  fake_electron_pass = vector<UInt_t>(denominatorType.size(),0);


  for (UInt_t denominatorTypeIndex = 0; denominatorTypeIndex < denominatorType.size(); ++denominatorTypeIndex) {
    fOutputTrees[denominatorTypeIndex]->Branch("nvertices",&nvertices,"nvertices/i");
    fOutputTrees[denominatorTypeIndex]->Branch("evt_num",&evt_num,"evt_num/i");
    fOutputTrees[denominatorTypeIndex]->Branch("lumi_sec",&lumi_sec,"lumi_sec/i");
    fOutputTrees[denominatorTypeIndex]->Branch("run_num",&run_num,"run_num/i");

    fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_electron_pt"),&fake_electron_pt[denominatorTypeIndex],TString("fake_electron_pt/F"));
    fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_electron_eta"),&fake_electron_eta[denominatorTypeIndex],TString("fake_electron_eta/F"));
    fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_electron_phi"),&fake_electron_phi[denominatorTypeIndex],TString("fake_electron_phi/F"));
    fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_electron_pass"),&fake_electron_pass[denominatorTypeIndex],TString("fake_electron_pass/i"));
  }

}

void mithep::ElectronFakeRateMod::SlaveTerminate()
{       
  fOutputFile->cd();  
  for (UInt_t denominatorTypeIndex = 0; denominatorTypeIndex < denominatorType.size(); ++denominatorTypeIndex) {
      fOutputTrees[denominatorTypeIndex]->Print();

  }
  fOutputFile->Write();
  fOutputFile->Close();
  
}


void mithep::ElectronFakeRateMod::Process()
{

  vector<UInt_t> nElectrons(denominatorType.size(),0);

  gDebugMask  = mithep::Debug::kAnalysis;  // debug message category
  gDebugLevel = 1;                 // higher level allows more messages to print
  
  LoadBranch(fMuonName);
  LoadBranch(fElectronName);
  LoadBranch(fPFCandidateName);
  LoadBranch(fPrimVtxName);      
  LoadBranch(fPUEnergyDensityName);      
  LoadBranch(fTrigMaskName);
  LoadBranch("PFMet"); 
  LoadBranch(mithep::Names::gkMvfConversionBrn);

  if(fPFMet->At(0)->Et() >25) return;
  
  if(store_npu)
    LoadBranch(fPileupName);
    
  mithep::RunLumiRangeMap::RunLumiPairType rl(GetEventHeader()->RunNum(), GetEventHeader()->LumiSec());      
  if(useJSON && !rlrm.HasRunLumi(rl)) return;    
  
  fillTriggerBits( GetHLTTable(), fTrigMask, fTriggerBits );
  
  PFnoPUflag.clear();
  makePFnoPUArray(fPFCandidates, PFnoPUflag, fPrimVerts );
    
  vector<const Muon*> goodMuons;
  vector<const Electron*> goodElectrons;

  for(Int_t i=0; i<fMuons->GetEntries(); i++) {
    const Muon *mu = fMuons->At(i);     
    if( muonPOG2012CutBasedIDTight(mu,fPrimVerts->At(0),fPFCandidates,fPileupEnergyDensity) )           
      goodMuons.push_back(mu);
  }
  for(Int_t i=0; i<fElectrons->GetEntries(); i++) {
    const Electron *ele = fElectrons->At(i); 
    Bool_t fails_some_denominator = kFALSE;
    for ( UInt_t denominatorTypeIndex = 0 ; denominatorTypeIndex < denominatorType.size() ; ++denominatorTypeIndex ) {
      if (passElectronDenominatorCuts(ele, fPrimVerts->At(0), denominatorType[denominatorTypeIndex]))  
        nElectrons[denominatorTypeIndex]++; 
      else
        fails_some_denominator = kTRUE;
    }

    // if(!fails_some_denominator && electronPOG2012CutBasedIDMedium(ele, fPrimVerts->At(0), fPFCandidates, fConversions, fPileupEnergyDensity->At(0)->RhoLowEta(), mithep::ElectronTools::kEleEAData2011))
    assert(0); // doesn't compile
      goodElectrons.push_back(ele);
  }     

  Int_t NZCandidates = 0;
  const mithep::Muon* ZMuon1 = 0;
  const mithep::Muon* ZMuon2 = 0;
  const mithep::Electron *ZEle1 = 0;
  const mithep::Electron *ZEle2 = 0;
  Double_t ZPt = 0;
  Double_t ZMass = 0;
  for(UInt_t i=0; i<goodMuons.size(); i++) {
    FourVectorM mu1;
    mu1.SetCoordinates(goodMuons[i]->Pt(), goodMuons[i]->Eta(), goodMuons[i]->Phi(), g_muon_mass );
    for(UInt_t j=i+1; j<goodMuons.size(); j++) {
      FourVectorM mu2;
      mu2.SetCoordinates(goodMuons[j]->Pt(), goodMuons[j]->Eta(), goodMuons[j]->Phi(), g_muon_mass );
      FourVectorM dilepton = mu1+mu2;
      if (dilepton.M() > massLo && dilepton.M() < massHi) {
        ZMuon1 = goodMuons[i];
        ZMuon2 = goodMuons[j];
        ZPt = dilepton.Pt();
        ZMass = dilepton.M();
        NZCandidates++;
      }
    }
  }
  
  for(UInt_t i=0; i<goodElectrons.size(); i++) {
    FourVectorM ele1;
    ele1.SetCoordinates(goodElectrons[i]->Pt(), goodElectrons[i]->Eta(), goodElectrons[i]->Phi(),g_electron_mass );
    for(UInt_t j=i+1; j<goodElectrons.size(); j++) {
      FourVectorM ele2;
      ele2.SetCoordinates(goodElectrons[j]->Pt(), goodElectrons[j]->Eta(), goodElectrons[j]->Phi(), g_electron_mass );
      FourVectorM dilepton = ele1+ele2;
      if (dilepton.M() > massLo && dilepton.M() < massHi) {
        ZEle1 = goodElectrons[i];
        ZEle2 = goodElectrons[j];
        ZPt = dilepton.Pt();
        ZMass = dilepton.M();
        NZCandidates++;
      }
    }
  }

  if (NZCandidates != 1) return;
  Int_t ZDecayType = 0;
  if (ZMuon1 && ZMuon2) ZDecayType = 13; 
  else if (ZEle1 && ZEle2) ZDecayType = 11; 
  else cout << "Error. Z Leptons not properly found.\n";              
  
  if(ZDecayType == 13){
    
    if(!(ZMuon1->Pt() > 20 && ZMuon2->Pt() > 10)  && !(ZMuon2->Pt() > 20 && ZMuon1->Pt() > 10)) return;

    std::bitset<1024> zmuon1_matchbits = mithep::ElectronFakeRateMod::GetHLTMatchBits(ZMuon1->Pt(), ZMuon1->Eta(), ZMuon1->Phi());
    std::bitset<1024> zmuon2_matchbits = mithep::ElectronFakeRateMod::GetHLTMatchBits(ZMuon2->Pt(), ZMuon2->Eta(), ZMuon2->Phi());
    
    Bool_t pass_dimuon_trigger1 = (ZDecayType == 13) 
      && (fTriggerBits.test(kHLT_Mu17_Mu8)) 
      &&
      (
       ((zmuon1_matchbits.test(kHLT_Mu17_Mu8_Mu1Obj)) 
        && (zmuon2_matchbits.test(kHLT_Mu17_Mu8_Mu2Obj)))
       ||
       ((zmuon1_matchbits.test(kHLT_Mu17_Mu8_Mu2Obj)) 
        && (zmuon2_matchbits.test(kHLT_Mu17_Mu8_Mu1Obj)))
       );
    
    Bool_t pass_dimuon_trigger2 = (ZDecayType == 13)
      && (fTriggerBits.test(kHLT_Mu17_TkMu8)) 
      &&
      (
       ((zmuon1_matchbits.test(kHLT_Mu17_TkMu8_Mu1Obj))
        && (zmuon2_matchbits.test(kHLT_Mu17_TkMu8_Mu2Obj)))
       ||
       ((zmuon1_matchbits.test(kHLT_Mu17_TkMu8_Mu2Obj)) 
        && (zmuon2_matchbits.test(kHLT_Mu17_TkMu8_Mu1Obj)))
       );
    
    
    if(!pass_dimuon_trigger1 && !pass_dimuon_trigger2) return;
    if(!(die_or_dimu == e_dimuon_data)) return;                         
  }     
  else if (ZDecayType == 11){
    
    if(!(ZEle1->Pt() > 20 && ZEle2->Pt() > 10)  && !(ZEle2->Pt() > 20 && ZEle1->Pt() > 10)) return;

    std::bitset<1024> zele1_matchbits = mithep::ElectronFakeRateMod::GetHLTMatchBits(ZEle1->Pt(), ZEle1->Eta(), ZEle1->Phi());
    std::bitset<1024> zele2_matchbits = mithep::ElectronFakeRateMod::GetHLTMatchBits(ZEle2->Pt(), ZEle2->Eta(), ZEle2->Phi());
    
    Bool_t pass_dielectron_trigger =
      (fTriggerBits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL)) 
      &&
      (
       ((zele1_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele1Obj)) 
        && (zele2_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele2Obj)))
       ||
       ((zele1_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele2Obj)) 
        && (zele2_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele1Obj)))
       );               
    
    
    if(!pass_dielectron_trigger) return;        
    if(!(die_or_dimu == e_dielectron_data))     return;
    
  }
  else assert(0);

  nvertices = fPrimVerts->GetEntries();
  run_num = GetEventHeader()->RunNum();
  lumi_sec =  GetEventHeader()->LumiSec();
  evt_num = GetEventHeader()->EvtNum();

  
  for(Int_t i=0; i<fElectrons->GetEntries(); i++) {
    const Electron *ele = fElectrons->At(i);                    
    
    if (ZDecayType == 11 && (ele == ZEle1 || ele == ZEle2)) continue;
    if (ZEle1 && mithep::MathUtils::DeltaR(ele->Eta(), ele->Phi(), ZEle1->Eta(), ZEle1->Phi()) < 0.3) continue;
    if (ZEle2 && mithep::MathUtils::DeltaR(ele->Eta(), ele->Phi(), ZEle2->Eta(), ZEle2->Phi()) < 0.3) continue;
    if (ZMuon1 && mithep::MathUtils::DeltaR(ele->Eta(), ele->Phi(), ZMuon1->Eta(), ZMuon1->Phi()) < 0.3) continue;
    if (ZMuon2 && mithep::MathUtils::DeltaR(ele->Eta(), ele->Phi(), ZMuon2->Eta(), ZMuon2->Phi()) < 0.3) continue;  
    
    for ( UInt_t denominatorTypeIndex = 0 ; denominatorTypeIndex < denominatorType.size() ; ++denominatorTypeIndex ) {
      
      if ((ZDecayType == 11 &&  nElectrons[denominatorTypeIndex] > 3) 
          ||
          (ZDecayType == 13 &&  nElectrons[denominatorTypeIndex] > 1) 
          ) 
        continue;

      if (passElectronDenominatorCuts(ele, fPrimVerts->At(0), denominatorType[denominatorTypeIndex])) {
        
        fake_electron_pt[denominatorTypeIndex] = ele->Pt();
        fake_electron_eta[denominatorTypeIndex] = ele->Eta();
        fake_electron_phi[denominatorTypeIndex] = ele->Phi();

        Bool_t passNumerator = kTRUE;
        
        ControlFlags ctrl;
        mithep::ElectronTools::EElectronEffectiveAreaTarget eraEle = mithep::ElectronTools::kEleEAData2012;
        vector<const mithep::PFCandidate*> photonsToVeto;
        
        if(!electronReferencePreSelection(ctrl,  ele, fPrimVerts->At(0)).passPre()) passNumerator = kFALSE;
        if(fabs(ele->Ip3dPVSignificance()) > 4) passNumerator = kFALSE;
        if(!electronReferenceIDMVASelectionV1(ctrl,  ele, fPrimVerts->At(0)).looseID()) passNumerator = kFALSE;
        if(!electronReferenceIsoSelection(ctrl,ele,fPrimVerts->At(0),fPFCandidates,fPileupEnergyDensity,eraEle,photonsToVeto).passLooseIso()) passNumerator = kFALSE;   
        
        fake_electron_pass[denominatorTypeIndex] = passNumerator ? 1 : 0;

        fOutputTrees[denominatorTypeIndex]->Fill();
      }//if pass denominator cut
    } //loop over denominator types
  } //loop over muons
}


std::bitset<1024> mithep::ElectronFakeRateMod::GetHLTMatchBits(const Double_t pt, const Double_t eta, const Double_t phi)
{
  std::bitset<1024> object_bitset;
  const Double_t hltMatchR = 0.2;
          
  if(HasHLTInfo()) {
    const TriggerTable *hltTable = GetHLTTable();
    assert(hltTable);
    for(UInt_t itrig=0; itrig<triggerNames.size(); itrig++) {
      const TriggerName *trigname = hltTable->Get(triggerNames[itrig].c_str());
      if(!trigname) continue;
                              
      const TList *list = GetHLTObjectsTable()->GetList(triggerNames[itrig].c_str());
      if(!list) continue;
      TIter iter(list->MakeIterator());
      const TriggerObject *to = dynamic_cast<const TriggerObject*>(iter.Next()); 
                                
      while(to) {         
        if(to->IsHLT()) {
          
          //assert(triggerObjNames1[itrig].length()>0);
          
          if(triggerObjNames1[itrig].length()>0 && (triggerObjNames1[itrig] == string(to->ModuleName()))) {
            if(MathUtils::DeltaR(phi,eta,to->Phi(),to->Eta()) < hltMatchR){
              object_bitset.set(triggerObjIds1[itrig]);
            }   
          }
          
          if(triggerObjNames2[itrig].length()>0 && (triggerObjNames2[itrig] == string(to->ModuleName()))) {
            if(MathUtils::DeltaR(phi,eta,to->Phi(),to->Eta()) < hltMatchR){
              object_bitset.set(triggerObjIds2[itrig]);
            }   
          }
          
        }
        to = dynamic_cast<const TriggerObject*>(iter.Next());
      }    
    }
  }
  
  return object_bitset;
}

Bool_t passElectronDenominatorCuts(const mithep::Electron *ele, const mithep::Vertex * vtx, Int_t DenominatorType) {
  
  Bool_t pass = kTRUE;  

  if (DenominatorType == 1) {           
    if (fabs(ele->Ip3dPVSignificance()) > 100) pass = kFALSE;
    if (fabs(ele->Eta()) > 2.5) pass = kFALSE;
    if (ele->Pt()  < 7) pass = kFALSE;
    if (ele->CorrectedNExpectedHitsInner() > 1) pass = kFALSE;  
    if (fabs(ele->BestTrk()->D0Corrected(*vtx)) > 0.5 ||
        fabs(ele->BestTrk()->DzCorrected(*vtx)) > 1.0 ) pass = kFALSE;
  }
  else if (DenominatorType == 2) {              
    if (fabs(ele->Ip3dPVSignificance()) > 4) pass = kFALSE;
    if (fabs(ele->Eta()) > 2.5) pass = kFALSE;
    if (ele->Pt()  < 7) pass = kFALSE;
    if (ele->CorrectedNExpectedHitsInner() > 1) pass = kFALSE;  
    if (fabs(ele->BestTrk()->D0Corrected(*vtx)) > 0.5 ||
        fabs(ele->BestTrk()->DzCorrected(*vtx)) > 1.0 ) pass = kFALSE;
  }
  else if (DenominatorType == 3) {
    if (fabs(ele->Ip3dPVSignificance()) > 100) pass = kFALSE;
    if (fabs(ele->Eta()) > 2.5) pass = kFALSE;
    if (ele->Pt()  < 7) pass = kFALSE;
    if (ele->CorrectedNExpectedHitsInner() > 1) pass = kFALSE;  
  }
  else if (DenominatorType == 4) {      
    if (fabs(ele->Ip3dPVSignificance()) > 4) pass = kFALSE;
    if (fabs(ele->Eta()) > 2.5) pass = kFALSE;
    if (ele->Pt()  < 7) pass = kFALSE;
    if (ele->CorrectedNExpectedHitsInner() > 1) pass = kFALSE;  
  }
  else
    assert(0);
  
  return pass;
}


Revision:	1.9
Committed:	Tue Oct 23 12:43:04 2012 UTC (12 years, 6 months ago) by dkralph
Content type:	text/plain
Branch:	MAIN
CVS Tags:	compiled, HEAD
Changes since 1.8:	+2 -1 lines
Log Message:	* empty log message *
#	User	Rev	Content
1	anlevin	1.1	#include "ElectronFakeRateMod.h"
2	anlevin	1.7	#include "SelectionFuncs.h"
3	anlevin	1.1
4
5			//needed because of externs
6			vector<bool> PFnoPUflag;
7			Float_t computePFMuonIso(const mithep::Muon *mu,
8			const mithep::Vertex * vtx,
9			const mithep::Array<mithep::PFCandidate> * pfCandidates,
10			const Double_t dRMax);
11
12			Bool_t passElectronDenominatorCuts(const mithep::Electron ele, const mithep::Vertex vtx, Int_t DenominatorType);
13
14			string IntToString(int i) {
15			char temp[100];
16			sprintf(temp, "%d", i);
17			string str = temp;
18			return str;
19			}
20
21			class TTree;
22			class TFile;
23			class TString;
24
25			const Float_t g_electron_mass = 0.51099892e-3;
26			const Float_t g_muon_mass = 105.658369e-3;
27
28
29			void mithep::ElectronFakeRateMod::SlaveBegin()
30			{
31
32
33			denominatorType.push_back(1);
34			denominatorType.push_back(2);
35			denominatorType.push_back(3);
36			denominatorType.push_back(4);
37
38			initElectronIDMVAV1();
39
40			ReqBranch("PFMet", fPFMet);
41
42			ReqBranch(mithep::Names::gkMvfConversionBrn, fConversions);
43
44			fMuonName = mithep::Names::gkMuonBrn;
45			ReqBranch(fMuonName, fMuons);
46
47			fElectronName = mithep::Names::gkElectronBrn;
48			ReqBranch(fElectronName, fElectrons);
49
50			fPFCandidateName = mithep::Names::gkPFCandidatesBrn;
51			ReqBranch(fPFCandidateName, fPFCandidates);
52
53			fPrimVtxName = mithep::Names::gkPVBrn;
54			ReqBranch(fPrimVtxName, fPrimVerts);
55
56			fPUEnergyDensityName = mithep::Names::gkPileupEnergyDensityBrn;
57			ReqBranch(fPUEnergyDensityName, fPileupEnergyDensity);
58
59			fTrigMaskName = mithep::Names::gkHltBitBrn;
60			ReqBranch(fTrigMaskName, fTrigMask);
61
62			fPileupName = mithep::Names::gkPileupInfoBrn;
63			ReqBranch(fPileupName, fPileup);
64
65
66			fOutputFile = new TFile(fOutputName, "RECREATE");
67			fOutputTrees = vector<TTree*>(denominatorType.size(),0);
68			for (UInt_t denominatorTypeIndex = 0; denominatorTypeIndex < denominatorType.size(); ++denominatorTypeIndex) {
69			fOutputTrees[denominatorTypeIndex] = new TTree(TString("denominator_v"+IntToString(denominatorType[denominatorTypeIndex])),TString("denominator_v"+IntToString(denominatorType[denominatorTypeIndex])));
70			}
71
72			rlrm.AddJSONFile("/afs/cern.ch/cms/CAF/CMSCOMM/COMM_DQM/certification/Collisions12/8TeV/Prompt/Cert_190456-193557_8TeV_PromptReco_Collisions12_JSON.txt");
73
74			fake_electron_pt = vector<Float_t>(denominatorType.size(),0);
75			fake_electron_eta = vector<Float_t>(denominatorType.size(),0);
76			fake_electron_phi = vector<Float_t>(denominatorType.size(),0);
77			fake_electron_pass = vector<UInt_t>(denominatorType.size(),0);
78
79
80			for (UInt_t denominatorTypeIndex = 0; denominatorTypeIndex < denominatorType.size(); ++denominatorTypeIndex) {
81			fOutputTrees[denominatorTypeIndex]->Branch("nvertices",&nvertices,"nvertices/i");
82			fOutputTrees[denominatorTypeIndex]->Branch("evt_num",&evt_num,"evt_num/i");
83			fOutputTrees[denominatorTypeIndex]->Branch("lumi_sec",&lumi_sec,"lumi_sec/i");
84			fOutputTrees[denominatorTypeIndex]->Branch("run_num",&run_num,"run_num/i");
85
86			fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_electron_pt"),&fake_electron_pt[denominatorTypeIndex],TString("fake_electron_pt/F"));
87			fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_electron_eta"),&fake_electron_eta[denominatorTypeIndex],TString("fake_electron_eta/F"));
88			fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_electron_phi"),&fake_electron_phi[denominatorTypeIndex],TString("fake_electron_phi/F"));
89			fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_electron_pass"),&fake_electron_pass[denominatorTypeIndex],TString("fake_electron_pass/i"));
90			}
91
92			}
93
94			void mithep::ElectronFakeRateMod::SlaveTerminate()
95			{
96			fOutputFile->cd();
97			for (UInt_t denominatorTypeIndex = 0; denominatorTypeIndex < denominatorType.size(); ++denominatorTypeIndex) {
98			fOutputTrees[denominatorTypeIndex]->Print();
99
100			}
101			fOutputFile->Write();
102			fOutputFile->Close();
103
104			}
105
106
107			void mithep::ElectronFakeRateMod::Process()
108			{
109
110			vector<UInt_t> nElectrons(denominatorType.size(),0);
111
112			gDebugMask = mithep::Debug::kAnalysis; // debug message category
113			gDebugLevel = 1; // higher level allows more messages to print
114
115			LoadBranch(fMuonName);
116			LoadBranch(fElectronName);
117			LoadBranch(fPFCandidateName);
118			LoadBranch(fPrimVtxName);
119			LoadBranch(fPUEnergyDensityName);
120			LoadBranch(fTrigMaskName);
121			LoadBranch("PFMet");
122			LoadBranch(mithep::Names::gkMvfConversionBrn);
123
124			if(fPFMet->At(0)->Et() >25) return;
125
126			if(store_npu)
127			LoadBranch(fPileupName);
128
129			mithep::RunLumiRangeMap::RunLumiPairType rl(GetEventHeader()->RunNum(), GetEventHeader()->LumiSec());
130			if(useJSON && !rlrm.HasRunLumi(rl)) return;
131
132			fillTriggerBits( GetHLTTable(), fTrigMask, fTriggerBits );
133
134	anlevin	1.6	PFnoPUflag.clear();
135	anlevin	1.1	makePFnoPUArray(fPFCandidates, PFnoPUflag, fPrimVerts );
136
137			vector<const Muon*> goodMuons;
138			vector<const Electron*> goodElectrons;
139
140			for(Int_t i=0; i<fMuons->GetEntries(); i++) {
141	anlevin	1.3	const Muon *mu = fMuons->At(i);
142	anlevin	1.4	if( muonPOG2012CutBasedIDTight(mu,fPrimVerts->At(0),fPFCandidates,fPileupEnergyDensity) )
143	anlevin	1.1	goodMuons.push_back(mu);
144			}
145			for(Int_t i=0; i<fElectrons->GetEntries(); i++) {
146			const Electron *ele = fElectrons->At(i);
147			Bool_t fails_some_denominator = kFALSE;
148			for ( UInt_t denominatorTypeIndex = 0 ; denominatorTypeIndex < denominatorType.size() ; ++denominatorTypeIndex ) {
149			if (passElectronDenominatorCuts(ele, fPrimVerts->At(0), denominatorType[denominatorTypeIndex]))
150			nElectrons[denominatorTypeIndex]++;
151			else
152			fails_some_denominator = kTRUE;
153			}
154
155	dkralph	1.9	// if(!fails_some_denominator && electronPOG2012CutBasedIDMedium(ele, fPrimVerts->At(0), fPFCandidates, fConversions, fPileupEnergyDensity->At(0)->RhoLowEta(), mithep::ElectronTools::kEleEAData2011))
156			assert(0); // doesn't compile
157	anlevin	1.1	goodElectrons.push_back(ele);
158			}
159
160			Int_t NZCandidates = 0;
161			const mithep::Muon* ZMuon1 = 0;
162			const mithep::Muon* ZMuon2 = 0;
163			const mithep::Electron *ZEle1 = 0;
164			const mithep::Electron *ZEle2 = 0;
165			Double_t ZPt = 0;
166			Double_t ZMass = 0;
167			for(UInt_t i=0; i<goodMuons.size(); i++) {
168			FourVectorM mu1;
169			mu1.SetCoordinates(goodMuons[i]->Pt(), goodMuons[i]->Eta(), goodMuons[i]->Phi(), g_muon_mass );
170			for(UInt_t j=i+1; j<goodMuons.size(); j++) {
171			FourVectorM mu2;
172			mu2.SetCoordinates(goodMuons[j]->Pt(), goodMuons[j]->Eta(), goodMuons[j]->Phi(), g_muon_mass );
173			FourVectorM dilepton = mu1+mu2;
174			if (dilepton.M() > massLo && dilepton.M() < massHi) {
175			ZMuon1 = goodMuons[i];
176			ZMuon2 = goodMuons[j];
177			ZPt = dilepton.Pt();
178			ZMass = dilepton.M();
179			NZCandidates++;
180			}
181			}
182			}
183
184			for(UInt_t i=0; i<goodElectrons.size(); i++) {
185			FourVectorM ele1;
186			ele1.SetCoordinates(goodElectrons[i]->Pt(), goodElectrons[i]->Eta(), goodElectrons[i]->Phi(),g_electron_mass );
187			for(UInt_t j=i+1; j<goodElectrons.size(); j++) {
188			FourVectorM ele2;
189			ele2.SetCoordinates(goodElectrons[j]->Pt(), goodElectrons[j]->Eta(), goodElectrons[j]->Phi(), g_electron_mass );
190			FourVectorM dilepton = ele1+ele2;
191			if (dilepton.M() > massLo && dilepton.M() < massHi) {
192			ZEle1 = goodElectrons[i];
193			ZEle2 = goodElectrons[j];
194			ZPt = dilepton.Pt();
195			ZMass = dilepton.M();
196			NZCandidates++;
197			}
198			}
199			}
200
201			if (NZCandidates != 1) return;
202			Int_t ZDecayType = 0;
203			if (ZMuon1 && ZMuon2) ZDecayType = 13;
204			else if (ZEle1 && ZEle2) ZDecayType = 11;
205			else cout << "Error. Z Leptons not properly found.\n";
206
207			if(ZDecayType == 13){
208
209			if(!(ZMuon1->Pt() > 20 && ZMuon2->Pt() > 10) && !(ZMuon2->Pt() > 20 && ZMuon1->Pt() > 10)) return;
210
211			std::bitset<1024> zmuon1_matchbits = mithep::ElectronFakeRateMod::GetHLTMatchBits(ZMuon1->Pt(), ZMuon1->Eta(), ZMuon1->Phi());
212			std::bitset<1024> zmuon2_matchbits = mithep::ElectronFakeRateMod::GetHLTMatchBits(ZMuon2->Pt(), ZMuon2->Eta(), ZMuon2->Phi());
213
214			Bool_t pass_dimuon_trigger1 = (ZDecayType == 13)
215			&& (fTriggerBits.test(kHLT_Mu17_Mu8))
216			&&
217			(
218			((zmuon1_matchbits.test(kHLT_Mu17_Mu8_Mu1Obj))
219			&& (zmuon2_matchbits.test(kHLT_Mu17_Mu8_Mu2Obj)))
220			\|\|
221			((zmuon1_matchbits.test(kHLT_Mu17_Mu8_Mu2Obj))
222			&& (zmuon2_matchbits.test(kHLT_Mu17_Mu8_Mu1Obj)))
223			);
224
225			Bool_t pass_dimuon_trigger2 = (ZDecayType == 13)
226			&& (fTriggerBits.test(kHLT_Mu17_TkMu8))
227			&&
228			(
229			((zmuon1_matchbits.test(kHLT_Mu17_TkMu8_Mu1Obj))
230			&& (zmuon2_matchbits.test(kHLT_Mu17_TkMu8_Mu2Obj)))
231			\|\|
232			((zmuon1_matchbits.test(kHLT_Mu17_TkMu8_Mu2Obj))
233			&& (zmuon2_matchbits.test(kHLT_Mu17_TkMu8_Mu1Obj)))
234			);
235
236
237			if(!pass_dimuon_trigger1 && !pass_dimuon_trigger2) return;
238			if(!(die_or_dimu == e_dimuon_data)) return;
239			}
240			else if (ZDecayType == 11){
241
242			if(!(ZEle1->Pt() > 20 && ZEle2->Pt() > 10) && !(ZEle2->Pt() > 20 && ZEle1->Pt() > 10)) return;
243
244			std::bitset<1024> zele1_matchbits = mithep::ElectronFakeRateMod::GetHLTMatchBits(ZEle1->Pt(), ZEle1->Eta(), ZEle1->Phi());
245			std::bitset<1024> zele2_matchbits = mithep::ElectronFakeRateMod::GetHLTMatchBits(ZEle2->Pt(), ZEle2->Eta(), ZEle2->Phi());
246
247			Bool_t pass_dielectron_trigger =
248			(fTriggerBits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL))
249			&&
250			(
251			((zele1_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele1Obj))
252			&& (zele2_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele2Obj)))
253			\|\|
254			((zele1_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele2Obj))
255			&& (zele2_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele1Obj)))
256			);
257
258
259			if(!pass_dielectron_trigger) return;
260			if(!(die_or_dimu == e_dielectron_data)) return;
261
262			}
263			else assert(0);
264
265			nvertices = fPrimVerts->GetEntries();
266			run_num = GetEventHeader()->RunNum();
267			lumi_sec = GetEventHeader()->LumiSec();
268			evt_num = GetEventHeader()->EvtNum();
269
270
271			for(Int_t i=0; i<fElectrons->GetEntries(); i++) {
272			const Electron *ele = fElectrons->At(i);
273
274			if (ZDecayType == 11 && (ele == ZEle1 \|\| ele == ZEle2)) continue;
275			if (ZEle1 && mithep::MathUtils::DeltaR(ele->Eta(), ele->Phi(), ZEle1->Eta(), ZEle1->Phi()) < 0.3) continue;
276			if (ZEle2 && mithep::MathUtils::DeltaR(ele->Eta(), ele->Phi(), ZEle2->Eta(), ZEle2->Phi()) < 0.3) continue;
277			if (ZMuon1 && mithep::MathUtils::DeltaR(ele->Eta(), ele->Phi(), ZMuon1->Eta(), ZMuon1->Phi()) < 0.3) continue;
278			if (ZMuon2 && mithep::MathUtils::DeltaR(ele->Eta(), ele->Phi(), ZMuon2->Eta(), ZMuon2->Phi()) < 0.3) continue;
279
280			for ( UInt_t denominatorTypeIndex = 0 ; denominatorTypeIndex < denominatorType.size() ; ++denominatorTypeIndex ) {
281
282			if ((ZDecayType == 11 && nElectrons[denominatorTypeIndex] > 3)
283			\|\|
284			(ZDecayType == 13 && nElectrons[denominatorTypeIndex] > 1)
285			)
286			continue;
287
288			if (passElectronDenominatorCuts(ele, fPrimVerts->At(0), denominatorType[denominatorTypeIndex])) {
289
290			fake_electron_pt[denominatorTypeIndex] = ele->Pt();
291			fake_electron_eta[denominatorTypeIndex] = ele->Eta();
292			fake_electron_phi[denominatorTypeIndex] = ele->Phi();
293
294			Bool_t passNumerator = kTRUE;
295
296			ControlFlags ctrl;
297			mithep::ElectronTools::EElectronEffectiveAreaTarget eraEle = mithep::ElectronTools::kEleEAData2012;
298			vector<const mithep::PFCandidate*> photonsToVeto;
299
300			if(!electronReferencePreSelection(ctrl, ele, fPrimVerts->At(0)).passPre()) passNumerator = kFALSE;
301			if(fabs(ele->Ip3dPVSignificance()) > 4) passNumerator = kFALSE;
302			if(!electronReferenceIDMVASelectionV1(ctrl, ele, fPrimVerts->At(0)).looseID()) passNumerator = kFALSE;
303			if(!electronReferenceIsoSelection(ctrl,ele,fPrimVerts->At(0),fPFCandidates,fPileupEnergyDensity,eraEle,photonsToVeto).passLooseIso()) passNumerator = kFALSE;
304
305			fake_electron_pass[denominatorTypeIndex] = passNumerator ? 1 : 0;
306
307			fOutputTrees[denominatorTypeIndex]->Fill();
308			}//if pass denominator cut
309			} //loop over denominator types
310			} //loop over muons
311			}
312
313
314			std::bitset<1024> mithep::ElectronFakeRateMod::GetHLTMatchBits(const Double_t pt, const Double_t eta, const Double_t phi)
315			{
316			std::bitset<1024> object_bitset;
317			const Double_t hltMatchR = 0.2;
318
319			if(HasHLTInfo()) {
320			const TriggerTable *hltTable = GetHLTTable();
321			assert(hltTable);
322			for(UInt_t itrig=0; itrig<triggerNames.size(); itrig++) {
323			const TriggerName *trigname = hltTable->Get(triggerNames[itrig].c_str());
324			if(!trigname) continue;
325
326			const TList *list = GetHLTObjectsTable()->GetList(triggerNames[itrig].c_str());
327			if(!list) continue;
328			TIter iter(list->MakeIterator());
329			const TriggerObject to = dynamic_cast<const TriggerObject>(iter.Next());
330
331			while(to) {
332			if(to->IsHLT()) {
333
334			//assert(triggerObjNames1[itrig].length()>0);
335
336			if(triggerObjNames1[itrig].length()>0 && (triggerObjNames1[itrig] == string(to->ModuleName()))) {
337			if(MathUtils::DeltaR(phi,eta,to->Phi(),to->Eta()) < hltMatchR){
338			object_bitset.set(triggerObjIds1[itrig]);
339			}
340			}
341
342			if(triggerObjNames2[itrig].length()>0 && (triggerObjNames2[itrig] == string(to->ModuleName()))) {
343			if(MathUtils::DeltaR(phi,eta,to->Phi(),to->Eta()) < hltMatchR){
344			object_bitset.set(triggerObjIds2[itrig]);
345			}
346			}
347
348			}
349			to = dynamic_cast<const TriggerObject*>(iter.Next());
350			}
351			}
352			}
353
354			return object_bitset;
355			}
356
357			Bool_t passElectronDenominatorCuts(const mithep::Electron ele, const mithep::Vertex vtx, Int_t DenominatorType) {
358
359			Bool_t pass = kTRUE;
360
361			if (DenominatorType == 1) {
362			if (fabs(ele->Ip3dPVSignificance()) > 100) pass = kFALSE;
363			if (fabs(ele->Eta()) > 2.5) pass = kFALSE;
364			if (ele->Pt() < 7) pass = kFALSE;
365			if (ele->CorrectedNExpectedHitsInner() > 1) pass = kFALSE;
366			if (fabs(ele->BestTrk()->D0Corrected(*vtx)) > 0.5 \|\|
367			fabs(ele->BestTrk()->DzCorrected(*vtx)) > 1.0 ) pass = kFALSE;
368			}
369			else if (DenominatorType == 2) {
370			if (fabs(ele->Ip3dPVSignificance()) > 4) pass = kFALSE;
371			if (fabs(ele->Eta()) > 2.5) pass = kFALSE;
372			if (ele->Pt() < 7) pass = kFALSE;
373			if (ele->CorrectedNExpectedHitsInner() > 1) pass = kFALSE;
374			if (fabs(ele->BestTrk()->D0Corrected(*vtx)) > 0.5 \|\|
375			fabs(ele->BestTrk()->DzCorrected(*vtx)) > 1.0 ) pass = kFALSE;
376			}
377			else if (DenominatorType == 3) {
378			if (fabs(ele->Ip3dPVSignificance()) > 100) pass = kFALSE;
379			if (fabs(ele->Eta()) > 2.5) pass = kFALSE;
380			if (ele->Pt() < 7) pass = kFALSE;
381			if (ele->CorrectedNExpectedHitsInner() > 1) pass = kFALSE;
382			}
383			else if (DenominatorType == 4) {
384			if (fabs(ele->Ip3dPVSignificance()) > 4) pass = kFALSE;
385			if (fabs(ele->Eta()) > 2.5) pass = kFALSE;
386			if (ele->Pt() < 7) pass = kFALSE;
387			if (ele->CorrectedNExpectedHitsInner() > 1) pass = kFALSE;
388			}
389			else
390			assert(0);
391
392			return pass;
393			}
394
395
396