NonMCBackground/src/MuonFakeRateMod.cc

#include "MuonFakeRateMod.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);
                        
unsigned makePFnoPUArray(const mithep::Array<mithep::PFCandidate> * fPFCandidates,
                         vector<bool> &pfNoPileUpFlag,
                         const mithep::Array<mithep::Vertex>     * vtxArr );
                         
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 );
  
  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(electron2012CutBasedIDMedium(ele, fPrimVerts->At(0), fPFCandidates, fConversions, fPileupEnergyDensity)) 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(), 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
}

unsigned makePFnoPUArray(const mithep::Array<mithep::PFCandidate> * fPFCandidates,
                         vector<bool> &pfNoPileUpflag,
                         const mithep::Array<mithep::Vertex>      * vtxArr )
{
  for(UInt_t i = 0; i < fPFCandidates->GetEntries(); i++) {
    const mithep::PFCandidate *pf = fPFCandidates->At(i);
    assert(pf);

    if(pf->PFType() == mithep::PFCandidate::eHadron) {
      if(pf->HasTrackerTrk() && 
         vtxArr->At(0)->HasTrack(pf->TrackerTrk()) &&
         vtxArr->At(0)->TrackWeight(pf->TrackerTrk()) > 0) {

        pfNoPileUpflag.push_back(1);

      } else { 

        Bool_t vertexFound = kFALSE;
        const mithep::Vertex *closestVtx = 0;
        Double_t dzmin = 10000;
        
        // loop over vertices
        for(UInt_t j = 0; j < vtxArr->GetEntries(); j++) {
          const mithep::Vertex *vtx = vtxArr->At(j);
          assert(vtx);
          
          if(pf->HasTrackerTrk() && 
             vtx->HasTrack(pf->TrackerTrk()) &&
             vtx->TrackWeight(pf->TrackerTrk()) > 0) {
            vertexFound = kTRUE;
            closestVtx = vtx;
            break;
          }
          Double_t dz = fabs(pf->SourceVertex().Z() - vtx->Z());
          if(dz < dzmin) {
            closestVtx = vtx;
            dzmin = dz;
          }
        }

        //      if(fCheckClosestZVertex) {
        if(1) {
          // Fallback: if track is not associated with any vertex,
          // associate it with the vertex closest in z
          if(vertexFound || closestVtx != vtxArr->At(0)) {
            //      pfPileUp->Add(pf);
            pfNoPileUpflag.push_back(0);
          } else {
            pfNoPileUpflag.push_back(1);
          }
        } else {
          if(vertexFound && closestVtx != vtxArr->At(0)) {
            //      pfPileUp->Add(pf);
            pfNoPileUpflag.push_back(0);
          } else {
            //      PFCandidate * pfNoPileUp->AddNew(); // Ridiculous but that's how it is
            pfNoPileUpflag.push_back(1);
          }
        }
      } //hadron & trk stuff
    } else { // hadron
      //      PFCandidate * ptr = pfNoPileUp->AddNew();
      pfNoPileUpflag.push_back(1);
    }
  } // Loop over PF candidates

  return pfNoPileUpflag.size();
}

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.5
Committed:	Sat Jul 7 13:33:32 2012 UTC (12 years, 10 months ago) by anlevin
Content type:	text/plain
Branch:	MAIN
Changes since 1.4:	+4 -3 lines
Log Message:	updated muon fake rate mod
#	User	Rev	Content
1	anlevin	1.1	#include "MuonFakeRateMod.h"
2
3
4			//needed because of externs
5			vector<bool> PFnoPUflag;
6			Float_t computePFMuonIso(const mithep::Muon *mu,
7			const mithep::Vertex * vtx,
8			const mithep::Array<mithep::PFCandidate> * pfCandidates,
9			const Double_t dRMax);
10
11			unsigned makePFnoPUArray(const mithep::Array<mithep::PFCandidate> * fPFCandidates,
12			vector<bool> &pfNoPileUpFlag,
13			const mithep::Array<mithep::Vertex> * vtxArr );
14
15			Bool_t passMuonDenominatorCuts(const mithep::Muon mu, const mithep::Vertex vtx, Int_t DenominatorType);
16
17			string IntToString(int i) {
18			char temp[100];
19			sprintf(temp, "%d", i);
20			string str = temp;
21			return str;
22			}
23
24			class TTree;
25			class TFile;
26			class TString;
27
28			const Float_t g_muon_mass = 105.658369e-3;
29
30
31			void mithep::MuonFakeRateMod::SlaveBegin()
32			{
33	anlevin	1.5
34	anlevin	1.1	denominatorType.push_back(1);
35			denominatorType.push_back(2);
36
37			nProbes = 0;
38
39			ReqBranch("PFMet", fPFMet);
40
41			ReqBranch(mithep::Names::gkMvfConversionBrn, fConversions);
42
43			fMuonName = mithep::Names::gkMuonBrn;
44			ReqBranch(fMuonName, fMuons);
45
46			fElectronName = mithep::Names::gkElectronBrn;
47			ReqBranch(fElectronName, fElectrons);
48
49			fPFCandidateName = mithep::Names::gkPFCandidatesBrn;
50			ReqBranch(fPFCandidateName, fPFCandidates);
51
52			fPrimVtxName = mithep::Names::gkPVBrn;
53			ReqBranch(fPrimVtxName, fPrimVerts);
54
55			fPUEnergyDensityName = mithep::Names::gkPileupEnergyDensityBrn;
56			ReqBranch(fPUEnergyDensityName, fPileupEnergyDensity);
57
58			fTrigMaskName = mithep::Names::gkHltBitBrn;
59			ReqBranch(fTrigMaskName, fTrigMask);
60
61			fPileupName = mithep::Names::gkPileupInfoBrn;
62			ReqBranch(fPileupName, fPileup);
63
64			fOutputFile = new TFile(fOutputName, "RECREATE");
65			fOutputTrees = vector<TTree*>(denominatorType.size(),0);
66			for (UInt_t denominatorTypeIndex = 0; denominatorTypeIndex < denominatorType.size(); ++denominatorTypeIndex) {
67			fOutputTrees[denominatorTypeIndex] = new TTree(TString("denominator_v"+IntToString(denominatorType[denominatorTypeIndex])),TString("denominator_v"+IntToString(denominatorType[denominatorTypeIndex])));
68			}
69
70			fake_muon_pt = vector<Float_t>(denominatorType.size(),0);
71			fake_muon_eta = vector<Float_t>(denominatorType.size(),0);
72			fake_muon_phi = vector<Float_t>(denominatorType.size(),0);
73			fake_muon_pass = vector<UInt_t>(denominatorType.size(),0);
74
75
76			for (UInt_t denominatorTypeIndex = 0; denominatorTypeIndex < denominatorType.size(); ++denominatorTypeIndex) {
77			fOutputTrees[denominatorTypeIndex]->Branch("nvertices",&nvertices,"nvertices/i");
78			fOutputTrees[denominatorTypeIndex]->Branch("evt_num",&evt_num,"evt_num/i");
79			fOutputTrees[denominatorTypeIndex]->Branch("lumi_sec",&lumi_sec,"lumi_sec/i");
80			fOutputTrees[denominatorTypeIndex]->Branch("run_num",&run_num,"run_num/i");
81
82			fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_muon_pt"),&fake_muon_pt[denominatorTypeIndex],TString("fake_muon_pt/F"));
83			fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_muon_eta"),&fake_muon_eta[denominatorTypeIndex],TString("fake_muon_eta/F"));
84			fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_muon_phi"),&fake_muon_phi[denominatorTypeIndex],TString("fake_muon_phi/F"));
85			fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_muon_pass"),&fake_muon_pass[denominatorTypeIndex],TString("fake_muon_pass/i"));
86			}
87
88			rlrm.AddJSONFile("/afs/cern.ch/cms/CAF/CMSCOMM/COMM_DQM/certification/Collisions12/8TeV/Prompt/Cert_190456-194479_8TeV_PromptReco_Collisions12_JSON.txt");
89	anlevin	1.5
90	anlevin	1.1	}
91
92			void mithep::MuonFakeRateMod::SlaveTerminate()
93			{
94			fOutputFile->cd();
95			for (UInt_t denominatorTypeIndex = 0; denominatorTypeIndex < denominatorType.size(); ++denominatorTypeIndex) {
96			fOutputTrees[denominatorTypeIndex]->Print();
97
98			}
99			fOutputFile->Write();
100			fOutputFile->Close();
101
102			}
103
104
105			void mithep::MuonFakeRateMod::Process()
106			{
107			vector<UInt_t> nLeptons(denominatorType.size(),0);
108			vector<UInt_t> nMuons(denominatorType.size(),0);
109
110
111			gDebugMask = mithep::Debug::kAnalysis; // debug message category
112			gDebugLevel = 1; // higher level allows more messages to print
113
114			LoadBranch(fMuonName);
115			LoadBranch(fElectronName);
116			LoadBranch(fPFCandidateName);
117			LoadBranch(fPrimVtxName);
118			LoadBranch(fPUEnergyDensityName);
119			LoadBranch(fTrigMaskName);
120			LoadBranch("PFMet");
121
122			LoadBranch(mithep::Names::gkMvfConversionBrn);
123			if(store_npu)
124			LoadBranch(fPileupName);
125
126			mithep::RunLumiRangeMap::RunLumiPairType rl(GetEventHeader()->RunNum(), GetEventHeader()->LumiSec());
127			if(useJSON && !rlrm.HasRunLumi(rl)) return;
128
129			fillTriggerBits( GetHLTTable(), fTrigMask, fTriggerBits );
130
131			makePFnoPUArray(fPFCandidates, PFnoPUflag, fPrimVerts );
132
133			if(fPFMet->At(0)->Et() >25) return;
134
135			vector<const Muon*> goodMuons;
136			vector<const Electron*> goodElectrons;
137
138			for(Int_t i=0; i<fMuons->GetEntries(); i++) {
139			const Muon *mu = fMuons->At(i);
140			for ( UInt_t denominatorTypeIndex = 0 ; denominatorTypeIndex < denominatorType.size() ; ++denominatorTypeIndex ) {
141			if (passMuonDenominatorCuts(mu, fPrimVerts->At(0), denominatorType[denominatorTypeIndex]))
142			nMuons[denominatorTypeIndex]++;
143			}
144	anlevin	1.3
145	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...
146			passMuonDenominatorCuts(mu, fPrimVerts->At(0), denominatorType[1])
147			&&
148	anlevin	1.4	muonPOG2012CutBasedIDTight(mu,fPrimVerts->At(0),fPFCandidates,fPileupEnergyDensity)
149	anlevin	1.1	)
150			goodMuons.push_back(mu);
151
152			}
153			for(Int_t i=0; i<fElectrons->GetEntries(); i++) {
154			const Electron *ele = fElectrons->At(i);
155	anlevin	1.3	if(electron2012CutBasedIDMedium(ele, fPrimVerts->At(0), fPFCandidates, fConversions, fPileupEnergyDensity)) goodElectrons.push_back(ele);
156	anlevin	1.1
157			}
158
159			Int_t NZCandidates = 0;
160			const mithep::Muon* ZMuon1 = 0;
161			const mithep::Muon* ZMuon2 = 0;
162			const mithep::Electron *ZEle1 = 0;
163			const mithep::Electron *ZEle2 = 0;
164			Double_t ZPt = 0;
165			Double_t ZMass = 0;
166			for(UInt_t i=0; i<goodMuons.size(); i++) {
167			FourVectorM mu1;
168			mu1.SetCoordinates(goodMuons[i]->Pt(), goodMuons[i]->Eta(), goodMuons[i]->Phi(), 105.658369e-3 );
169			for(UInt_t j=i+1; j<goodMuons.size(); j++) {
170			FourVectorM mu2;
171			mu2.SetCoordinates(goodMuons[j]->Pt(), goodMuons[j]->Eta(), goodMuons[j]->Phi(), 105.658369e-3 );
172			FourVectorM dilepton = mu1+mu2;
173			if (dilepton.M() > massLo && dilepton.M() < massHi) {
174			ZMuon1 = goodMuons[i];
175			ZMuon2 = goodMuons[j];
176			ZPt = dilepton.Pt();
177			ZMass = dilepton.M();
178			NZCandidates++;
179			}
180			}
181			}
182
183			for(UInt_t i=0; i<goodElectrons.size(); i++) {
184			FourVectorM ele1;
185			ele1.SetCoordinates(goodElectrons[i]->Pt(), goodElectrons[i]->Eta(), goodElectrons[i]->Phi(), 0.51099892e-3 );
186			for(UInt_t j=i+1; j<goodElectrons.size(); j++) {
187			FourVectorM ele2;
188			ele2.SetCoordinates(goodElectrons[j]->Pt(), goodElectrons[j]->Eta(), goodElectrons[j]->Phi(), 0.51099892e-3 );
189			FourVectorM dilepton = ele1+ele2;
190			if (dilepton.M() > massLo && dilepton.M() < massHi) {
191			ZEle1 = goodElectrons[i];
192			ZEle2 = goodElectrons[j];
193			ZPt = dilepton.Pt();
194			ZMass = dilepton.M();
195			NZCandidates++;
196			}
197			}
198			}
199
200			if (NZCandidates != 1) return;
201			Int_t ZDecayType = 0;
202			if (ZMuon1 && ZMuon2) ZDecayType = 13;
203			else if (ZEle1 && ZEle2) ZDecayType = 11;
204			else cout << "Error. Z Leptons not properly found.\n";
205
206			if(ZDecayType == 13){
207
208			std::bitset<1024> zmuon1_matchbits = mithep::MuonFakeRateMod::GetHLTMatchBits(ZMuon1->Pt(), ZMuon1->Eta(), ZMuon1->Phi());
209			std::bitset<1024> zmuon2_matchbits = mithep::MuonFakeRateMod::GetHLTMatchBits(ZMuon2->Pt(), ZMuon2->Eta(), ZMuon2->Phi());
210
211			Bool_t pass_dimuon_trigger1 = (ZDecayType == 13)
212			&& (fTriggerBits.test(kHLT_Mu17_Mu8))
213			&&
214			(
215			((zmuon1_matchbits.test(kHLT_Mu17_Mu8_Mu1Obj))
216			&& (zmuon2_matchbits.test(kHLT_Mu17_Mu8_Mu2Obj)))
217			\|\|
218			((zmuon1_matchbits.test(kHLT_Mu17_Mu8_Mu2Obj))
219			&& (zmuon2_matchbits.test(kHLT_Mu17_Mu8_Mu1Obj)))
220			);
221
222			Bool_t pass_dimuon_trigger2 = (ZDecayType == 13)
223			&& (fTriggerBits.test(kHLT_Mu17_TkMu8))
224			&&
225			(
226			((zmuon1_matchbits.test(kHLT_Mu17_TkMu8_Mu1Obj))
227			&& (zmuon2_matchbits.test(kHLT_Mu17_TkMu8_Mu2Obj)))
228			\|\|
229			((zmuon1_matchbits.test(kHLT_Mu17_TkMu8_Mu2Obj))
230			&& (zmuon2_matchbits.test(kHLT_Mu17_TkMu8_Mu1Obj)))
231			);
232
233
234			if(!pass_dimuon_trigger1 && !pass_dimuon_trigger2) return;
235			if(!(die_or_dimu == e_dimuon_data)) return;
236			}
237			else if (ZDecayType == 11){
238
239			std::bitset<1024> zele1_matchbits = mithep::MuonFakeRateMod::GetHLTMatchBits(ZEle1->Pt(), ZEle1->Eta(), ZEle1->Phi());
240			std::bitset<1024> zele2_matchbits = mithep::MuonFakeRateMod::GetHLTMatchBits(ZEle2->Pt(), ZEle2->Eta(), ZEle2->Phi());
241
242			Bool_t pass_dielectron_trigger =
243			(fTriggerBits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL))
244			&&
245			(
246			((zele1_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele1Obj))
247			&& (zele2_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele2Obj)))
248			\|\|
249			((zele1_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele2Obj))
250			&& (zele2_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele1Obj)))
251			);
252
253			if(!pass_dielectron_trigger) return;
254			if(!(die_or_dimu == e_dielectron_data)) return;
255
256			}
257			else assert(0);
258
259			nvertices = fPrimVerts->GetEntries();
260			run_num = GetEventHeader()->RunNum();
261			lumi_sec = GetEventHeader()->LumiSec();
262			evt_num = GetEventHeader()->EvtNum();
263
264
265			for(Int_t i=0; i<fMuons->GetEntries(); i++) {
266			const Muon *mu = fMuons->At(i);
267
268			if (ZDecayType == 13 && (mu == ZMuon1 \|\| mu == ZMuon2)) continue;
269			if (ZEle1 && mithep::MathUtils::DeltaR(mu->Eta(), mu->Phi(), ZEle1->Eta(), ZEle1->Phi()) < 0.3) continue;
270			if (ZEle2 && mithep::MathUtils::DeltaR(mu->Eta(), mu->Phi(), ZEle2->Eta(), ZEle2->Phi()) < 0.3) continue;
271			if (ZMuon1 && mithep::MathUtils::DeltaR(mu->Eta(), mu->Phi(), ZMuon1->Eta(), ZMuon1->Phi()) < 0.3) continue;
272			if (ZMuon2 && mithep::MathUtils::DeltaR(mu->Eta(), mu->Phi(), ZMuon2->Eta(), ZMuon2->Phi()) < 0.3) continue;
273
274			for ( UInt_t denominatorTypeIndex = 0 ; denominatorTypeIndex < denominatorType.size() ; ++denominatorTypeIndex ) {
275			if ((ZDecayType == 11 && nMuons[denominatorTypeIndex] > 1)
276			\|\|
277			(ZDecayType == 13 && nMuons[denominatorTypeIndex] > 3)
278			) continue;
279
280			if (passMuonDenominatorCuts(mu, fPrimVerts->At(0), denominatorType[denominatorTypeIndex])) {
281
282			fake_muon_pt[denominatorTypeIndex] = mu->Pt();
283			fake_muon_eta[denominatorTypeIndex] = mu->Eta();
284			fake_muon_phi[denominatorTypeIndex] = mu->Phi();
285
286			Bool_t passNumerator = kTRUE;
287			vector<const mithep::PFCandidate*> photonsToVeto;
288
289			ControlFlags ctrl;
290			mithep::MuonTools::EMuonEffectiveAreaTarget eraMu = mithep::MuonTools::kMuEAData2012;
291			if(!muonReferencePreSelection(ctrl, mu, fPrimVerts->At(0), fPFCandidates ).passPre()) passNumerator = kFALSE;
292			if(!muonIDPFSelection(ctrl, mu, fPrimVerts->At(0), fPFCandidates ).looseID()) passNumerator = kFALSE;
293			if(fabs(mu->Ip3dPVSignificance()) > 4) passNumerator = kFALSE;
294			if(!muonReferenceIsoSelection(ctrl,mu,fPrimVerts->At(0),fPFCandidates,fPileupEnergyDensity,eraMu,photonsToVeto).passLooseIso()) passNumerator = kFALSE;
295	anlevin	1.5
296			fake_muon_pass[denominatorTypeIndex] = passNumerator ? 1 : 0;
297	anlevin	1.1
298			fOutputTrees[denominatorTypeIndex]->Fill();
299
300			}//if pass denominator cut
301			} //loop over denominator types
302			} //loop over muons
303			}
304
305			unsigned makePFnoPUArray(const mithep::Array<mithep::PFCandidate> * fPFCandidates,
306			vector<bool> &pfNoPileUpflag,
307			const mithep::Array<mithep::Vertex> * vtxArr )
308			{
309			for(UInt_t i = 0; i < fPFCandidates->GetEntries(); i++) {
310			const mithep::PFCandidate *pf = fPFCandidates->At(i);
311			assert(pf);
312
313			if(pf->PFType() == mithep::PFCandidate::eHadron) {
314			if(pf->HasTrackerTrk() &&
315			vtxArr->At(0)->HasTrack(pf->TrackerTrk()) &&
316			vtxArr->At(0)->TrackWeight(pf->TrackerTrk()) > 0) {
317
318			pfNoPileUpflag.push_back(1);
319
320			} else {
321
322			Bool_t vertexFound = kFALSE;
323			const mithep::Vertex *closestVtx = 0;
324			Double_t dzmin = 10000;
325
326			// loop over vertices
327			for(UInt_t j = 0; j < vtxArr->GetEntries(); j++) {
328			const mithep::Vertex *vtx = vtxArr->At(j);
329			assert(vtx);
330
331			if(pf->HasTrackerTrk() &&
332			vtx->HasTrack(pf->TrackerTrk()) &&
333			vtx->TrackWeight(pf->TrackerTrk()) > 0) {
334			vertexFound = kTRUE;
335			closestVtx = vtx;
336			break;
337			}
338			Double_t dz = fabs(pf->SourceVertex().Z() - vtx->Z());
339			if(dz < dzmin) {
340			closestVtx = vtx;
341			dzmin = dz;
342			}
343			}
344
345			// if(fCheckClosestZVertex) {
346			if(1) {
347			// Fallback: if track is not associated with any vertex,
348			// associate it with the vertex closest in z
349			if(vertexFound \|\| closestVtx != vtxArr->At(0)) {
350			// pfPileUp->Add(pf);
351			pfNoPileUpflag.push_back(0);
352			} else {
353			pfNoPileUpflag.push_back(1);
354			}
355			} else {
356			if(vertexFound && closestVtx != vtxArr->At(0)) {
357			// pfPileUp->Add(pf);
358			pfNoPileUpflag.push_back(0);
359			} else {
360			// PFCandidate * pfNoPileUp->AddNew(); // Ridiculous but that's how it is
361			pfNoPileUpflag.push_back(1);
362			}
363			}
364			} //hadron & trk stuff
365			} else { // hadron
366			// PFCandidate * ptr = pfNoPileUp->AddNew();
367			pfNoPileUpflag.push_back(1);
368			}
369			} // Loop over PF candidates
370
371			return pfNoPileUpflag.size();
372			}
373
374			std::bitset<1024> mithep::MuonFakeRateMod::GetHLTMatchBits(const Double_t pt, const Double_t eta, const Double_t phi)
375			{
376			std::bitset<1024> object_bitset;
377			const Double_t hltMatchR = 0.2;
378
379			if(HasHLTInfo()) {
380			const TriggerTable *hltTable = GetHLTTable();
381			assert(hltTable);
382			for(UInt_t itrig=0; itrig<triggerNames.size(); itrig++) {
383			const TriggerName *trigname = hltTable->Get(triggerNames[itrig].c_str());
384			if(!trigname) continue;
385
386			const TList *list = GetHLTObjectsTable()->GetList(triggerNames[itrig].c_str());
387			if(!list) continue;
388			TIter iter(list->MakeIterator());
389			const TriggerObject to = dynamic_cast<const TriggerObject>(iter.Next());
390
391			while(to) {
392			if(to->IsHLT()) {
393
394			//assert(triggerObjNames1[itrig].length()>0);
395
396			if(triggerObjNames1[itrig].length()>0 && (triggerObjNames1[itrig] == string(to->ModuleName()))) {
397			if(MathUtils::DeltaR(phi,eta,to->Phi(),to->Eta()) < hltMatchR){
398			object_bitset.set(triggerObjIds1[itrig]);
399			}
400			}
401
402			if(triggerObjNames2[itrig].length()>0 && (triggerObjNames2[itrig] == string(to->ModuleName()))) {
403			if(MathUtils::DeltaR(phi,eta,to->Phi(),to->Eta()) < hltMatchR){
404			object_bitset.set(triggerObjIds2[itrig]);
405			}
406			}
407
408			}
409			to = dynamic_cast<const TriggerObject*>(iter.Next());
410			}
411			}
412			}
413
414			return object_bitset;
415			}
416
417			Bool_t passMuonDenominatorCuts(const mithep::Muon mu, const mithep::Vertex vtx, Int_t DenominatorType) {
418
419			Bool_t pass = kTRUE;
420
421			if (DenominatorType == 1) {
422
423			if (fabs(mu->Eta()) > 2.4) pass = kFALSE;
424
425			if (mu->Pt() < 5) pass = kFALSE;
426
427			if(!mu->IsTrackerMuon() && !mu->IsGlobalMuon()) pass = kFALSE;
428
429			if(fabs(mu->Ip3dPVSignificance()) > 100) pass = kFALSE;
430
431			if(!(mu->HasTrackerTrk() && fabs(mu->TrackerTrk()->D0Corrected(vtx)) < 0.5 && fabs(mu->TrackerTrk()->DzCorrected(vtx)) < 1.0)) pass = kFALSE;
432
433			}
434			else if (DenominatorType == 2) {
435
436			if (fabs(mu->Eta()) > 2.4) pass = kFALSE;
437
438			if (mu->Pt() < 5) pass = kFALSE;
439
440			if(!mu->IsTrackerMuon() && !mu->IsGlobalMuon())
441			pass = kFALSE;
442
443			if(fabs(mu->Ip3dPVSignificance()) > 4) pass = kFALSE;
444
445			if(!(mu->HasTrackerTrk() && fabs(mu->TrackerTrk()->D0Corrected(vtx)) < 0.5 && fabs(mu->TrackerTrk()->DzCorrected(vtx)) < 1.0)) pass = kFALSE;
446
447			}
448			else
449			assert(0);
450
451			return pass;
452			}
453
454