LeptonSelection/src/MuonSelection.cc

#include <math.h>
#include <iostream>

#include "MuonSelection.h"
#include "ParseArgs.h"
#include "SelectionStatus.h"

#include "Track.h"
#include "MuonTools.h"
#include "MuonQuality.h"
#include "MuonIDMVA.h"

mithep::MuonIDMVA * muIDMVA;
mithep::MuonTools muTools;

extern vector<bool> PFnoPUflag;
extern Float_t computePFMuonIso(const mithep::Muon *mu, 
                        const mithep::Vertex * vtx,
                        const mithep::Array<mithep::PFCandidate> * pfCandidates,
                        const Double_t dRMax);

//--------------------------------------------------------------------------------------------------
SelectionStatus muonDummyVeto(ControlFlags &ctrl, 
                              const mithep::Muon *muon, 
                              const mithep::Vertex * vtx) 
//--------------------------------------------------------------------------------------------------
{ 
  SelectionStatus status;
  status.setStatus(SelectionStatus::PRESELECTION);
  return status;  
}

//--------------------------------------------------------------------------------------------------
SelectionStatus muonCutBasedVeto(ControlFlags &ctrl, 
                                 const mithep::Muon *muon, 
                                 const mithep::Vertex *vtx) 
//--------------------------------------------------------------------------------------------------
{ 
  //
  // Loose cut-based ID for isolation veto
  //
  bool ret = true;
  
  if(!(muon->IsGlobalMuon() || muon->IsTrackerMuon())) ret = false;
  if( muon->NValidHits() < 11 )                        ret = false;
  if( fabs(muon->Ip3dPVSignificance()) >= 4 )          ret = false;

  SelectionStatus status;
  if( ret ) status.setStatus(SelectionStatus::PRESELECTION);
  return status;
}


//--------------------------------------------------------------------------------------------------
SelectionStatus noPreSelection( ControlFlags &ctrl,  const mithep::Muon * mu ) 

//--------------------------------------------------------------------------------------------------
{
        SelectionStatus status;
        status.setStatus(SelectionStatus::PRESELECTION); 
        if(ctrl.debug) cout << "muon presel returning status : " << status.getStatus() << endl;
        return status;
}

//--------------------------------------------------------------------------------------------------
SelectionStatus muonPreSelection( ControlFlags &ctrl,  
                                  const mithep::Muon * mu, 
                                  const mithep::Vertex * vtx, 
                                  const mithep::Array<mithep::PFCandidate> * pfCandidates ) 
//--------------------------------------------------------------------------------------------------
{ 
  bool ret = true; 

  ret &= ( fabs(mu->Ip3dPVSignificance()) < 100 ); 
  ret &= ( mu->Pt() >= 5 ); 
  ret &= ( fabs(mu->Eta()) <= 2.4 ); 
  ret &= ( mu->IsTrackerMuon() || mu->IsGlobalMuon() ); 
//   if( mu->IsTrackerMuon() && mu->HasTrackerTrk() && !mu->IsGlobalMuon() )  
//     ret &=  ( mu->Quality().QualityMask().Mask() & mithep::MuonQuality::AllArbitrated);
  ret &= (  mu->IsoR03SumPt()/mu->Pt() < 0.7 );

  SelectionStatus status;
  if( ret ) status.setStatus(SelectionStatus::PRESELECTION);
  if(ctrl.debug) cout << "muon presel returning status : " << status.getStatus() << endl;  
  return status;
} 


//--------------------------------------------------------------------------------------------------
SelectionStatus muonReferencePreSelection( ControlFlags &ctrl,  
                                           const mithep::Muon * mu, 
                                           const mithep::Vertex * vtx, 
                                           const mithep::Array<mithep::PFCandidate> * pfCandidates ) 
//--------------------------------------------------------------------------------------------------
{ 
  bool ret = true; 

  if(ctrl.debug) cout << "inside muonpresel ... " << endl;
  ret &= ( fabs(mu->Ip3dPVSignificance()) < 100 ); 
  if(ctrl.debug) cout << "and pass IP (" << mu->Ip3dPVSignificance() << ") ? ... " << ret << endl;  
  ret &= ( mu->Pt() >= 5 ); 
  if(ctrl.debug) cout << "and >5 GeV (" << mu->Pt() << ") ? ... " << ret << endl;  
  ret &= ( fabs(mu->Eta()) <= 2.4 ); 
  if(ctrl.debug) cout << "and < 2.4 eta (" << mu->Eta() << ")? ... " << ret << endl;  
  ret &= ( mu->IsTrackerMuon() || mu->IsGlobalMuon() ); 
  if(ctrl.debug) cout << "is Tracker or Global? ... " << ret << endl;
        ret &= (mu->HasTrackerTrk() && mu->TrackerTrk()->D0Corrected(*vtx) < 0.5 && mu->TrackerTrk()->DzCorrected(*vtx) < 1.0); 
  if(ctrl.debug) cout << "d0 and dz? ... " << ret << endl;  
          
  //  if( mu->IsTrackerMuon() && mu->HasTrackerTrk() && !mu->IsGlobalMuon() )  
  //    ret &=  ( mu->Quality().QualityMask().Mask() & mithep::MuonQuality::AllArbitrated);
  //  if(ctrl.debug) cout << "if isTrackerMuon, arbitrated ? ... " << ret << endl;  
  //  ret &= ( mu->BestTrk()->DzCorrected(vtx) < 0.1);   
  //  if( ctrl.debug ) cout << "elepresel : ret after dZcorr : " << ret <<  endl;


  SelectionStatus status;
  if( ret ) status.setStatus(SelectionStatus::PRESELECTION);
  if(ctrl.debug) cout << "muon presel returning status : " << status.getStatus() << endl;  
  return status;
} 


//--------------------------------------------------------------------------------------------------
SelectionStatus muonIDMVASelection(ControlFlags &ctrl, 
                                 const mithep::Muon *mu,
                                 const mithep::Vertex * vtx   )
//--------------------------------------------------------------------------------------------------
{
  double global_rchi2 = (mu->IsGlobalMuon()) ? mu->GlobalTrk()->RChi2() : 0.;

  // 
  // WARNING !!! KH hacked out to sync
  //
  /*
  double mvaval = muIDMVA->MVAValue_ID(mu->Pt(),
                                    mu->Eta(),
                                    mu->IsGlobalMuon(),
                                    mu->IsTrackerMuon(),
                                    mu->TrackerTrk()->RChi2(),
                                    global_rchi2,
                                    (Double_t)(mu->NValidHits()),
                                    (Double_t)(mu->TrackerTrk()->NHits()),
                                    (Double_t)(mu->TrackerTrk()->NPixelHits()),
                                    (Double_t)(mu->NMatches()),
                                    mu->TrkKink(),
                                    muTools.GetSegmentCompatability(mu),
                                    muTools.GetCaloCompatability(mu,true,true),
                                    mu->HadEnergy(),
                                    mu->EmEnergy(),
                                    mu->HadS9Energy(),
                                    mu->EmS9Energy(),
                                    (Bool_t)(ctrl.debug) );
  */
  double mvaval = 0.0;

  SelectionStatus status;
  bool pass = false;
  
  if( mu->IsGlobalMuon() && mu->IsTrackerMuon() 
      && fabs(mu->Eta()) <= 1.5 && mu->Pt() <= 10 && mvaval >= MUON_IDMVA_CUT_BIN0)  pass = true;
  else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
           && fabs(mu->Eta()) <= 1.5 && mu->Pt() > 10 && mvaval >= MUON_IDMVA_CUT_BIN1)  pass = true;
  else if( mu->IsGlobalMuon() && mu->IsTrackerMuon() 
           && fabs(mu->Eta()) > 1.5 && mu->Pt() <= 10 && mvaval >= MUON_IDMVA_CUT_BIN2)  pass = true;
  else if( mu->IsGlobalMuon() && mu->IsTrackerMuon() 
           && fabs(mu->Eta()) > 1.5 && mu->Pt() > 10 && mvaval >= MUON_IDMVA_CUT_BIN3)  pass = true;
  else if( !(mu->IsGlobalMuon()) && mu->IsTrackerMuon() && mvaval >= MUON_IDMVA_CUT_BIN4 )  pass = true;
  else if( mu->IsGlobalMuon() && !(mu->IsTrackerMuon()) && mvaval >= MUON_IDMVA_CUT_BIN5 )  pass = true;
   

  if( pass ) {
    status.orStatus(SelectionStatus::LOOSEID);
    status.orStatus(SelectionStatus::TIGHTID);
  }

  return status;
}

//--------------------------------------------------------------------------------------------------
void initMuonIDMVA() 
//--------------------------------------------------------------------------------------------------
{ 
  muIDMVA = new mithep::MuonIDMVA();
  vector<string> weightFiles;
  weightFiles.push_back("./data/MuonIDMVAWeights/MuonIDMVA_BDTG_barrel_lowpt_V2.weights.xml");
  weightFiles.push_back("./data/MuonIDMVAWeights/MuonIDMVA_BDTG_barrel_highpt_V2.weights.xml");
  weightFiles.push_back("./data/MuonIDMVAWeights/MuonIDMVA_BDTG_endcap_lowpt_V2.weights.xml");
  weightFiles.push_back("./data/MuonIDMVAWeights/MuonIDMVA_BDTG_endcap_highpt_V2.weights.xml");
  weightFiles.push_back("./data/MuonIDMVAWeights/MuonIDMVA_BDTG_tracker_V2.weights.xml");
  weightFiles.push_back("./data/MuonIDMVAWeights/MuonIDMVA_BDTG_global_V2.weights.xml");
  muIDMVA->Initialize( "MuonIDMVA",
                       mithep::MuonIDMVA::kIDV0,
                       kTRUE, weightFiles);
}

//--------------------------------------------------------------------------------------------------
SelectionStatus muonIDPFSelection(ControlFlags &ctrl, 
                                  const mithep::Muon *mu,
                                  const mithep::Vertex * vtx,
                                  const mithep::Array<mithep::PFCandidate> * pfCandidates )
//--------------------------------------------------------------------------------------------------
{
  bool pass = false;
  SelectionStatus status; // init'ed to FAIL

  // check that it matches to a PF muon
  for( int i=0; i<pfCandidates->GetEntries(); i++ ) {
    // tmp 
    if( !(PFnoPUflag[i]) ) continue; // my PF no PU hack

    const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*pfCandidates)[i]);
    if( (pf->TrackerTrk() == mu->TrackerTrk()) && abs(pf->PFType()) == mithep::PFCandidate::eMuon ) {
      pass = true;
      break;
    }
  }

  //
  // NB : ipsig here for reference synch ...
  //
#ifdef SYNC
  if( pass && fabs(mu->Ip3dPVSignificance()) < 4 ) {
    status.orStatus(SelectionStatus::LOOSEID);
    status.orStatus(SelectionStatus::TIGHTID);
  }
#else
  if( pass ) {
    status.orStatus(SelectionStatus::LOOSEID);
    status.orStatus(SelectionStatus::TIGHTID);
  }
#endif

  return status;
}
Revision:	1.21
Committed:	Tue May 22 22:31:17 2012 UTC (12 years, 11 months ago) by anlevin
Content type:	text/plain
Branch:	MAIN
Changes since 1.20:	+11 -8 lines
Log Message:	made changes for new round of synchronization
#	User	Rev	Content
1	khahn	1.1	#include <math.h>
2			#include <iostream>
3
4			#include "MuonSelection.h"
5			#include "ParseArgs.h"
6			#include "SelectionStatus.h"
7
8	khahn	1.4	#include "Track.h"
9			#include "MuonTools.h"
10			#include "MuonQuality.h"
11			#include "MuonIDMVA.h"
12
13			mithep::MuonIDMVA * muIDMVA;
14			mithep::MuonTools muTools;
15
16	khahn	1.19	extern vector<bool> PFnoPUflag;
17	khahn	1.4	extern Float_t computePFMuonIso(const mithep::Muon *mu,
18	anlevin	1.21	const mithep::Vertex * vtx,
19	khahn	1.4	const mithep::Array<mithep::PFCandidate> * pfCandidates,
20			const Double_t dRMax);
21
22			//--------------------------------------------------------------------------------------------------
23	khahn	1.8	SelectionStatus muonDummyVeto(ControlFlags &ctrl,
24			const mithep::Muon *muon,
25	anlevin	1.21	const mithep::Vertex * vtx)
26	khahn	1.8	//--------------------------------------------------------------------------------------------------
27			{
28			SelectionStatus status;
29			status.setStatus(SelectionStatus::PRESELECTION);
30			return status;
31			}
32
33			//--------------------------------------------------------------------------------------------------
34			SelectionStatus muonCutBasedVeto(ControlFlags &ctrl,
35			const mithep::Muon *muon,
36	anlevin	1.21	const mithep::Vertex *vtx)
37	khahn	1.8	//--------------------------------------------------------------------------------------------------
38			{
39			//
40			// Loose cut-based ID for isolation veto
41			//
42			bool ret = true;
43
44			if(!(muon->IsGlobalMuon() \|\| muon->IsTrackerMuon())) ret = false;
45			if( muon->NValidHits() < 11 ) ret = false;
46			if( fabs(muon->Ip3dPVSignificance()) >= 4 ) ret = false;
47
48			SelectionStatus status;
49			if( ret ) status.setStatus(SelectionStatus::PRESELECTION);
50			return status;
51			}
52
53
54			//--------------------------------------------------------------------------------------------------
55	khahn	1.4	SelectionStatus noPreSelection( ControlFlags &ctrl, const mithep::Muon * mu )
56
57			//--------------------------------------------------------------------------------------------------
58			{
59	anlevin	1.3	SelectionStatus status;
60			status.setStatus(SelectionStatus::PRESELECTION);
61			if(ctrl.debug) cout << "muon presel returning status : " << status.getStatus() << endl;
62			return status;
63			}
64	khahn	1.1
65	khahn	1.4	//--------------------------------------------------------------------------------------------------
66			SelectionStatus muonPreSelection( ControlFlags &ctrl,
67			const mithep::Muon * mu,
68	anlevin	1.21	const mithep::Vertex * vtx,
69	khahn	1.4	const mithep::Array<mithep::PFCandidate> * pfCandidates )
70			//--------------------------------------------------------------------------------------------------
71			{
72			bool ret = true;
73	khahn	1.20
74			ret &= ( fabs(mu->Ip3dPVSignificance()) < 100 );
75	khahn	1.14	ret &= ( mu->Pt() >= 5 );
76	khahn	1.7	ret &= ( fabs(mu->Eta()) <= 2.4 );
77	khahn	1.11	ret &= ( mu->IsTrackerMuon() \|\| mu->IsGlobalMuon() );
78	khahn	1.20	// if( mu->IsTrackerMuon() && mu->HasTrackerTrk() && !mu->IsGlobalMuon() )
79			// ret &= ( mu->Quality().QualityMask().Mask() & mithep::MuonQuality::AllArbitrated);
80			ret &= ( mu->IsoR03SumPt()/mu->Pt() < 0.7 );
81	khahn	1.1
82			SelectionStatus status;
83			if( ret ) status.setStatus(SelectionStatus::PRESELECTION);
84			if(ctrl.debug) cout << "muon presel returning status : " << status.getStatus() << endl;
85			return status;
86			}
87
88	khahn	1.4
89	khahn	1.1
90	khahn	1.13	//--------------------------------------------------------------------------------------------------
91			SelectionStatus muonReferencePreSelection( ControlFlags &ctrl,
92			const mithep::Muon * mu,
93	anlevin	1.21	const mithep::Vertex * vtx,
94	khahn	1.13	const mithep::Array<mithep::PFCandidate> * pfCandidates )
95			//--------------------------------------------------------------------------------------------------
96			{
97			bool ret = true;
98	khahn	1.20
99	khahn	1.13	if(ctrl.debug) cout << "inside muonpresel ... " << endl;
100	khahn	1.16	ret &= ( fabs(mu->Ip3dPVSignificance()) < 100 );
101	khahn	1.13	if(ctrl.debug) cout << "and pass IP (" << mu->Ip3dPVSignificance() << ") ? ... " << ret << endl;
102	khahn	1.14	ret &= ( mu->Pt() >= 5 );
103	khahn	1.13	if(ctrl.debug) cout << "and >5 GeV (" << mu->Pt() << ") ? ... " << ret << endl;
104			ret &= ( fabs(mu->Eta()) <= 2.4 );
105			if(ctrl.debug) cout << "and < 2.4 eta (" << mu->Eta() << ")? ... " << ret << endl;
106			ret &= ( mu->IsTrackerMuon() \|\| mu->IsGlobalMuon() );
107	anlevin	1.21	if(ctrl.debug) cout << "is Tracker or Global? ... " << ret << endl;
108			ret &= (mu->HasTrackerTrk() && mu->TrackerTrk()->D0Corrected(vtx) < 0.5 && mu->TrackerTrk()->DzCorrected(vtx) < 1.0);
109			if(ctrl.debug) cout << "d0 and dz? ... " << ret << endl;
110
111	khahn	1.13	// if( mu->IsTrackerMuon() && mu->HasTrackerTrk() && !mu->IsGlobalMuon() )
112			// ret &= ( mu->Quality().QualityMask().Mask() & mithep::MuonQuality::AllArbitrated);
113			// if(ctrl.debug) cout << "if isTrackerMuon, arbitrated ? ... " << ret << endl;
114			// ret &= ( mu->BestTrk()->DzCorrected(vtx) < 0.1);
115			// if( ctrl.debug ) cout << "elepresel : ret after dZcorr : " << ret << endl;
116
117
118			SelectionStatus status;
119			if( ret ) status.setStatus(SelectionStatus::PRESELECTION);
120			if(ctrl.debug) cout << "muon presel returning status : " << status.getStatus() << endl;
121			return status;
122			}
123
124
125
126	anlevin	1.3
127	khahn	1.4	//--------------------------------------------------------------------------------------------------
128			SelectionStatus muonIDMVASelection(ControlFlags &ctrl,
129			const mithep::Muon *mu,
130	anlevin	1.21	const mithep::Vertex * vtx )
131	khahn	1.4	//--------------------------------------------------------------------------------------------------
132			{
133			double global_rchi2 = (mu->IsGlobalMuon()) ? mu->GlobalTrk()->RChi2() : 0.;
134	anlevin	1.3
135	khahn	1.13	//
136			// WARNING !!! KH hacked out to sync
137			//
138			/*
139	khahn	1.4	double mvaval = muIDMVA->MVAValue_ID(mu->Pt(),
140			mu->Eta(),
141			mu->IsGlobalMuon(),
142			mu->IsTrackerMuon(),
143			mu->TrackerTrk()->RChi2(),
144			global_rchi2,
145			(Double_t)(mu->NValidHits()),
146	khahn	1.11	(Double_t)(mu->TrackerTrk()->NHits()),
147	khahn	1.4	(Double_t)(mu->TrackerTrk()->NPixelHits()),
148			(Double_t)(mu->NMatches()),
149			mu->TrkKink(),
150			muTools.GetSegmentCompatability(mu),
151			muTools.GetCaloCompatability(mu,true,true),
152			mu->HadEnergy(),
153			mu->EmEnergy(),
154			mu->HadS9Energy(),
155			mu->EmS9Energy(),
156			(Bool_t)(ctrl.debug) );
157	khahn	1.13	*/
158			double mvaval = 0.0;
159	anlevin	1.3
160	khahn	1.4	SelectionStatus status;
161			bool pass = false;
162
163			if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
164	khahn	1.7	&& fabs(mu->Eta()) <= 1.5 && mu->Pt() <= 10 && mvaval >= MUON_IDMVA_CUT_BIN0) pass = true;
165	khahn	1.4	else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
166	khahn	1.7	&& fabs(mu->Eta()) <= 1.5 && mu->Pt() > 10 && mvaval >= MUON_IDMVA_CUT_BIN1) pass = true;
167	khahn	1.4	else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
168	khahn	1.7	&& fabs(mu->Eta()) > 1.5 && mu->Pt() <= 10 && mvaval >= MUON_IDMVA_CUT_BIN2) pass = true;
169	khahn	1.4	else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
170			&& fabs(mu->Eta()) > 1.5 && mu->Pt() > 10 && mvaval >= MUON_IDMVA_CUT_BIN3) pass = true;
171	khahn	1.11	else if( !(mu->IsGlobalMuon()) && mu->IsTrackerMuon() && mvaval >= MUON_IDMVA_CUT_BIN4 ) pass = true;
172			else if( mu->IsGlobalMuon() && !(mu->IsTrackerMuon()) && mvaval >= MUON_IDMVA_CUT_BIN5 ) pass = true;
173
174	khahn	1.4
175
176			if( pass ) {
177			status.orStatus(SelectionStatus::LOOSEID);
178			status.orStatus(SelectionStatus::TIGHTID);
179	anlevin	1.3	}
180
181	khahn	1.4	return status;
182			}
183	anlevin	1.3
184	khahn	1.4	//--------------------------------------------------------------------------------------------------
185			void initMuonIDMVA()
186			//--------------------------------------------------------------------------------------------------
187			{
188			muIDMVA = new mithep::MuonIDMVA();
189			vector<string> weightFiles;
190	khahn	1.6	weightFiles.push_back("./data/MuonIDMVAWeights/MuonIDMVA_BDTG_barrel_lowpt_V2.weights.xml");
191			weightFiles.push_back("./data/MuonIDMVAWeights/MuonIDMVA_BDTG_barrel_highpt_V2.weights.xml");
192			weightFiles.push_back("./data/MuonIDMVAWeights/MuonIDMVA_BDTG_endcap_lowpt_V2.weights.xml");
193			weightFiles.push_back("./data/MuonIDMVAWeights/MuonIDMVA_BDTG_endcap_highpt_V2.weights.xml");
194			weightFiles.push_back("./data/MuonIDMVAWeights/MuonIDMVA_BDTG_tracker_V2.weights.xml");
195			weightFiles.push_back("./data/MuonIDMVAWeights/MuonIDMVA_BDTG_global_V2.weights.xml");
196	khahn	1.4	muIDMVA->Initialize( "MuonIDMVA",
197			mithep::MuonIDMVA::kIDV0,
198			kTRUE, weightFiles);
199			}
200	khahn	1.9
201			//--------------------------------------------------------------------------------------------------
202			SelectionStatus muonIDPFSelection(ControlFlags &ctrl,
203			const mithep::Muon *mu,
204	anlevin	1.21	const mithep::Vertex * vtx,
205	khahn	1.9	const mithep::Array<mithep::PFCandidate> * pfCandidates )
206			//--------------------------------------------------------------------------------------------------
207			{
208			bool pass = false;
209			SelectionStatus status; // init'ed to FAIL
210
211			// check that it matches to a PF muon
212			for( int i=0; i<pfCandidates->GetEntries(); i++ ) {
213	khahn	1.19	// tmp
214			if( !(PFnoPUflag[i]) ) continue; // my PF no PU hack
215
216	khahn	1.9	const mithep::PFCandidate pf = (mithep::PFCandidate)((*pfCandidates)[i]);
217	khahn	1.19	if( (pf->TrackerTrk() == mu->TrackerTrk()) && abs(pf->PFType()) == mithep::PFCandidate::eMuon ) {
218	khahn	1.17	pass = true;
219			break;
220	khahn	1.9	}
221			}
222
223	khahn	1.13	//
224			// NB : ipsig here for reference synch ...
225			//
226	khahn	1.19	#ifdef SYNC
227	khahn	1.16	if( pass && fabs(mu->Ip3dPVSignificance()) < 4 ) {
228	khahn	1.9	status.orStatus(SelectionStatus::LOOSEID);
229			status.orStatus(SelectionStatus::TIGHTID);
230			}
231	khahn	1.19	#else
232			if( pass ) {
233			status.orStatus(SelectionStatus::LOOSEID);
234			status.orStatus(SelectionStatus::TIGHTID);
235			}
236			#endif
237	khahn	1.9
238			return status;
239			}