NonMCBackground/src/MuonFakeRateMod.cc

#include "MuonFakeRateMod.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 passMuonDenominatorCuts(const mithep::Muon *mu, 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_muon_mass = 105.658369e-3;


void mithep::MuonFakeRateMod::SlaveBegin()
{

  denominatorType.push_back(1); 
  denominatorType.push_back(2); 

  nProbes = 0;

  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])));
  }

  fake_muon_pt = vector<Float_t>(denominatorType.size(),0);
  fake_muon_eta = vector<Float_t>(denominatorType.size(),0);
  fake_muon_phi = vector<Float_t>(denominatorType.size(),0);
  fake_muon_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_muon_pt"),&fake_muon_pt[denominatorTypeIndex],TString("fake_muon_pt/F"));
    fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_muon_eta"),&fake_muon_eta[denominatorTypeIndex],TString("fake_muon_eta/F"));
    fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_muon_phi"),&fake_muon_phi[denominatorTypeIndex],TString("fake_muon_phi/F"));
    fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_muon_pass"),&fake_muon_pass[denominatorTypeIndex],TString("fake_muon_pass/i"));
  }
  
  rlrm.AddJSONFile("/afs/cern.ch/cms/CAF/CMSCOMM/COMM_DQM/certification/Collisions12/8TeV/Prompt/Cert_190456-194479_8TeV_PromptReco_Collisions12_JSON.txt"); 

}

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

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


void mithep::MuonFakeRateMod::Process()
{
  vector<UInt_t> nLeptons(denominatorType.size(),0);
  vector<UInt_t> nMuons(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(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 );
    
  if(fPFMet->At(0)->Et() >25) return;
  
  vector<const Muon*> goodMuons;
  vector<const Electron*> goodElectrons;
  
  for(Int_t i=0; i<fMuons->GetEntries(); i++) {
    const Muon *mu = fMuons->At(i); 
    for ( UInt_t denominatorTypeIndex = 0 ; denominatorTypeIndex < denominatorType.size() ; ++denominatorTypeIndex ) {
      if (passMuonDenominatorCuts(mu, fPrimVerts->At(0), denominatorType[denominatorTypeIndex]))  
        nMuons[denominatorTypeIndex]++; 
    }
        
    if ( //maybe you should do a loop to make sure the muon passes all of the denominators? right now this is ok because 2 is tighter...
        passMuonDenominatorCuts(mu, fPrimVerts->At(0), denominatorType[1])
        &&
        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);
    // if(electronPOG2012CutBasedIDMedium(ele, fPrimVerts->At(0), fPFCandidates, fConversions, fPileupEnergyDensity->At(0)->RhoLowEta(), mithep::ElectronTools::kEleEAData2011)) goodElectrons.push_back(ele);
    assert(0); // doesn't compile

  }     
  
  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(), 105.658369e-3 );
    for(UInt_t j=i+1; j<goodMuons.size(); j++) {
      FourVectorM mu2;
      mu2.SetCoordinates(goodMuons[j]->Pt(), goodMuons[j]->Eta(), goodMuons[j]->Phi(), 105.658369e-3 );
      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(), 0.51099892e-3 );
    for(UInt_t j=i+1; j<goodElectrons.size(); j++) {
      FourVectorM ele2;
      ele2.SetCoordinates(goodElectrons[j]->Pt(), goodElectrons[j]->Eta(), goodElectrons[j]->Phi(), 0.51099892e-3 );
      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){

    std::bitset<1024> zmuon1_matchbits = mithep::MuonFakeRateMod::GetHLTMatchBits(ZMuon1->Pt(), ZMuon1->Eta(), ZMuon1->Phi());
    std::bitset<1024> zmuon2_matchbits = mithep::MuonFakeRateMod::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){
    
    std::bitset<1024> zele1_matchbits = mithep::MuonFakeRateMod::GetHLTMatchBits(ZEle1->Pt(), ZEle1->Eta(), ZEle1->Phi());
    std::bitset<1024> zele2_matchbits = mithep::MuonFakeRateMod::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<fMuons->GetEntries(); i++) {
    const Muon *mu =    fMuons->At(i);                  
    
    if (ZDecayType == 13 && (mu == ZMuon1 || mu == ZMuon2)) continue;
    if (ZEle1 && mithep::MathUtils::DeltaR(mu->Eta(), mu->Phi(), ZEle1->Eta(), ZEle1->Phi()) < 0.3) continue;
    if (ZEle2 && mithep::MathUtils::DeltaR(mu->Eta(), mu->Phi(), ZEle2->Eta(), ZEle2->Phi()) < 0.3) continue;
    if (ZMuon1 && mithep::MathUtils::DeltaR(mu->Eta(), mu->Phi(), ZMuon1->Eta(), ZMuon1->Phi()) < 0.3) continue;
    if (ZMuon2 && mithep::MathUtils::DeltaR(mu->Eta(), mu->Phi(), ZMuon2->Eta(), ZMuon2->Phi()) < 0.3) continue;  
    
    for ( UInt_t denominatorTypeIndex = 0 ; denominatorTypeIndex < denominatorType.size() ; ++denominatorTypeIndex ) {
      if ((ZDecayType == 11 &&  nMuons[denominatorTypeIndex] > 1) 
          ||
          (ZDecayType == 13 &&  nMuons[denominatorTypeIndex] > 3) 
          ) continue;
      
      if (passMuonDenominatorCuts(mu, fPrimVerts->At(0), denominatorType[denominatorTypeIndex])) {
        
        fake_muon_pt[denominatorTypeIndex] = mu->Pt();
        fake_muon_eta[denominatorTypeIndex] = mu->Eta();
        fake_muon_phi[denominatorTypeIndex] = mu->Phi();
        
        Bool_t passNumerator = kTRUE;
        vector<const mithep::PFCandidate*> photonsToVeto;
        
        ControlFlags ctrl;
        mithep::MuonTools::EMuonEffectiveAreaTarget eraMu = mithep::MuonTools::kMuEAData2012;
        if(!muonReferencePreSelection(ctrl,  mu, fPrimVerts->At(0), fPFCandidates ).passPre()) passNumerator = kFALSE;
        if(!muonIDPFSelection(ctrl,  mu, fPrimVerts->At(0), fPFCandidates ).looseID()) passNumerator = kFALSE;
        if(fabs(mu->Ip3dPVSignificance()) > 4) passNumerator = kFALSE;  
        if(!muonReferenceIsoSelection(ctrl,mu,fPrimVerts->At(0),fPFCandidates,fPileupEnergyDensity,eraMu,photonsToVeto).passLooseIso()) passNumerator = kFALSE;                 
       
        fake_muon_pass[denominatorTypeIndex] = passNumerator ? 1 : 0;
        
        fOutputTrees[denominatorTypeIndex]->Fill();
        
      }//if pass denominator cut
    } //loop over denominator types
  } //loop over muons
}

std::bitset<1024> mithep::MuonFakeRateMod::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 passMuonDenominatorCuts(const mithep::Muon *mu, const mithep::Vertex * vtx, Int_t DenominatorType) {
  
  Bool_t pass = kTRUE;  

  if (DenominatorType == 1) {
                
    if (fabs(mu->Eta()) > 2.4) pass = kFALSE;

    if (mu->Pt() < 5) pass = kFALSE;
        
    if(!mu->IsTrackerMuon() && !mu->IsGlobalMuon()) pass = kFALSE;      
                        
    if(fabs(mu->Ip3dPVSignificance()) > 100) pass = kFALSE;
                
    if(!(mu->HasTrackerTrk() && fabs(mu->TrackerTrk()->D0Corrected(*vtx)) < 0.5 && fabs(mu->TrackerTrk()->DzCorrected(*vtx)) < 1.0)) pass = kFALSE; 
    
  }
  else if (DenominatorType == 2) {
                
    if (fabs(mu->Eta()) > 2.4) pass = kFALSE;

    if (mu->Pt() < 5) pass = kFALSE;
        
    if(!mu->IsTrackerMuon() && !mu->IsGlobalMuon())
      pass = kFALSE;    
                        
    if(fabs(mu->Ip3dPVSignificance()) > 4) pass = kFALSE;
    
    if(!(mu->HasTrackerTrk() && fabs(mu->TrackerTrk()->D0Corrected(*vtx)) < 0.5 && fabs(mu->TrackerTrk()->DzCorrected(*vtx)) < 1.0)) pass = kFALSE;             
                
  }
  else
    assert(0);
  
  return pass;
}


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