CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MuonMVAEstimator.cc
Go to the documentation of this file.
1 #include <TFile.h>
2 #include <cmath>
3 #include <vector>
4 
5 #include "DQM/PhysicsHWW/interface/MuonMVAEstimator.h"
6 #include "DQM/PhysicsHWW/interface/trackSelections.h"
7 
8 using namespace std;
9 using namespace HWWFunctions;
10 
11 //
12 // dz of muon track : copied from muonSelections.cc
13 //
14 double dzPVmu(const LorentzVector& vtx, const LorentzVector& p4, const LorentzVector& pv){
15  return (vtx.z()-pv.z()) - ((vtx.x()-pv.x())*p4.x()+(vtx.y()-pv.y())*p4.y())/p4.pt() * p4.z()/p4.pt();
16 }
17 
18 //
19 // Definition of particle types (PFType)
20 //
21 enum ParticleType {
22  kX=0, // undefined
23  kh, // charged hadron
24  ke, // electron
25  kmu, // muon
26  kgamma, // photon
27  kh0, // neutral hadron
28  kh_HF, // HF tower identified as a hadron
29  kegamma_HF // HF tower identified as an EM particle
30 };
31 
32 //
33 // This function converts pdgId to PFType
34 //
35 int translatePdgIdToType(int pdgid)
36 {
37  switch ( std::abs(pdgid) )
38  {
39  case 211: return kh;
40  case 11: return ke;
41  case 13: return kmu;
42  case 22: return kgamma;
43  case 130: return kh0;
44  case 1: return kh_HF;
45  case 2: return kegamma_HF;
46  case 0: return kX;
47  default: return kX;
48  }
49 }
50 
51 
52 
53 //--------------------------------------------------------------------------------------------------
54 MuonMVAEstimator::MuonMVAEstimator() :
55  fMethodname("BDTG method"),
56  fisInitialized(kFALSE),
57  fPrintMVADebug(kFALSE),
58  fMVAType(kIsoRings),
59  fUseBinnedVersion(kTRUE),
60  fNMVABins(0)
61 {
62  // Constructor.
63  fTMVAReader = std::vector<TMVA::Reader*>(0);
64 }
65 
66 //--------------------------------------------------------------------------------------------------
67 MuonMVAEstimator::~MuonMVAEstimator(){
68  for (unsigned int i=0;i<fTMVAReader.size(); ++i) {
69  if (fTMVAReader[i]) delete fTMVAReader[i];
70  }
71 }
72 
73 //--------------------------------------------------------------------------------------------------
74 void MuonMVAEstimator::initialize( std::string methodName,
75  std::string weightsfile,
76  MuonMVAEstimator::MVAType type){
77 
78  std::vector<std::string> tempWeightFileVector;
79  tempWeightFileVector.push_back(weightsfile);
80  initialize(methodName,type,kFALSE,tempWeightFileVector);
81 }
82 
83 
84 //--------------------------------------------------------------------------------------------------
85 void MuonMVAEstimator::initialize( std::string methodName,
86  MuonMVAEstimator::MVAType type,
87  Bool_t useBinnedVersion,
88  std::vector<std::string> weightsfiles
89  ) {
90 
91  //clean up first
92  for (unsigned int i=0;i<fTMVAReader.size(); ++i) {
93  if (fTMVAReader[i]) delete fTMVAReader[i];
94  }
95  fTMVAReader.clear();
96 
97  //initialize
98  fisInitialized = kTRUE;
99  fMVAType = type;
100  fMethodname = methodName;
101  fUseBinnedVersion = useBinnedVersion;
102 
103  //Define expected number of bins
104  unsigned int ExpectedNBins = 0;
105  if (!fUseBinnedVersion) {
106  ExpectedNBins = 1;
107  } else if (type == kIDIsoRingsCombined) {
108  ExpectedNBins = 5;
109  } else if (type == kIsoRings) {
110  ExpectedNBins = 6;
111  } else if (type == kIsoDeltaR) {
112  ExpectedNBins = 4;
113  } else if (type == kID) {
114  ExpectedNBins = 6;
115  }
116  fNMVABins = ExpectedNBins;
117 
118  //Check number of weight files given
119  if (fNMVABins != weightsfiles.size() ) {
120  edm::LogError("InvalidInput") << "Error: Expected Number of bins = " << fNMVABins << " does not equal to weightsfiles.size() = "
121  << weightsfiles.size();
122  assert(fNMVABins == weightsfiles.size());
123  }
124 
125 
126  //Loop over all bins
127  for (unsigned int i=0;i<fNMVABins; ++i) {
128 
129  //TMVA::Reader *tmpTMVAReader = new TMVA::Reader( "!Color:!Silent:Error" );
130  TMVA::Reader *tmpTMVAReader = new TMVA::Reader( "!Color:Silent:Error" );
131  //tmpTMVAReader->SetVerbose(kTRUE);
132  tmpTMVAReader->SetVerbose(kFALSE);
133 
134  if (type == kIDIsoRingsCombined) {
135  // Pure tracking variables
136  tmpTMVAReader->AddVariable( "ChargedIso_DR0p0To0p1", &fMVAVar_ChargedIso_DR0p0To0p1 );
137  tmpTMVAReader->AddVariable( "ChargedIso_DR0p1To0p2", &fMVAVar_ChargedIso_DR0p1To0p2 );
138  tmpTMVAReader->AddVariable( "ChargedIso_DR0p2To0p3", &fMVAVar_ChargedIso_DR0p2To0p3 );
139  tmpTMVAReader->AddVariable( "ChargedIso_DR0p3To0p4", &fMVAVar_ChargedIso_DR0p3To0p4 );
140  tmpTMVAReader->AddVariable( "ChargedIso_DR0p4To0p5", &fMVAVar_ChargedIso_DR0p4To0p5 );
141  tmpTMVAReader->AddVariable( "GammaIso_DR0p0To0p1", &fMVAVar_GammaIso_DR0p0To0p1 );
142  tmpTMVAReader->AddVariable( "GammaIso_DR0p1To0p2", &fMVAVar_GammaIso_DR0p1To0p2 );
143  tmpTMVAReader->AddVariable( "GammaIso_DR0p2To0p3", &fMVAVar_GammaIso_DR0p2To0p3 );
144  tmpTMVAReader->AddVariable( "GammaIso_DR0p3To0p4", &fMVAVar_GammaIso_DR0p3To0p4 );
145  tmpTMVAReader->AddVariable( "GammaIso_DR0p4To0p5", &fMVAVar_GammaIso_DR0p4To0p5 );
146  tmpTMVAReader->AddVariable( "NeutralHadronIso_DR0p0To0p1", &fMVAVar_NeutralHadronIso_DR0p0To0p1 );
147  tmpTMVAReader->AddVariable( "NeutralHadronIso_DR0p1To0p2", &fMVAVar_NeutralHadronIso_DR0p1To0p2 );
148  tmpTMVAReader->AddVariable( "NeutralHadronIso_DR0p2To0p3", &fMVAVar_NeutralHadronIso_DR0p2To0p3 );
149  tmpTMVAReader->AddVariable( "NeutralHadronIso_DR0p3To0p4", &fMVAVar_NeutralHadronIso_DR0p3To0p4 );
150  tmpTMVAReader->AddVariable( "NeutralHadronIso_DR0p4To0p5", &fMVAVar_NeutralHadronIso_DR0p4To0p5 );
151  tmpTMVAReader->AddVariable( "TkNchi2", &fMVAVar_MuTkNchi2 );
152  if (i != 4) {
153  tmpTMVAReader->AddVariable( "GlobalNchi2", &fMVAVar_MuGlobalNchi2 );
154  tmpTMVAReader->AddVariable( "NValidHits", &fMVAVar_MuNValidHits );
155  }
156  tmpTMVAReader->AddVariable( "NTrackerHits", &fMVAVar_MuNTrackerHits );
157  tmpTMVAReader->AddVariable( "NPixelHits", &fMVAVar_MuNPixelHits );
158  tmpTMVAReader->AddVariable( "NMatches", &fMVAVar_MuNMatches );
159  tmpTMVAReader->AddVariable( "TrkKink", &fMVAVar_MuTrkKink );
160  tmpTMVAReader->AddVariable( "SegmentCompatibility", &fMVAVar_MuSegmentCompatibility );
161  tmpTMVAReader->AddVariable( "CaloCompatibility", &fMVAVar_MuCaloCompatibility );
162  tmpTMVAReader->AddVariable( "HadEnergy", &fMVAVar_MuHadEnergy );
163  tmpTMVAReader->AddVariable( "EmEnergy", &fMVAVar_MuEmEnergy );
164  tmpTMVAReader->AddVariable( "HadS9Energy", &fMVAVar_MuHadS9Energy );
165  tmpTMVAReader->AddVariable( "EmS9Energy", &fMVAVar_MuEmS9Energy );
166  tmpTMVAReader->AddSpectator("eta", &fMVAVar_MuEta);
167  tmpTMVAReader->AddSpectator("pt", &fMVAVar_MuPt);
168  tmpTMVAReader->AddSpectator("typeBits", &fMVAVar_MuTypeBits);
169  }
170 
171  if (type == kIsoRings) {
172  tmpTMVAReader->AddVariable( "ChargedIso_DR0p0To0p1", &fMVAVar_ChargedIso_DR0p0To0p1 );
173  tmpTMVAReader->AddVariable( "ChargedIso_DR0p1To0p2", &fMVAVar_ChargedIso_DR0p1To0p2 );
174  tmpTMVAReader->AddVariable( "ChargedIso_DR0p2To0p3", &fMVAVar_ChargedIso_DR0p2To0p3 );
175  tmpTMVAReader->AddVariable( "ChargedIso_DR0p3To0p4", &fMVAVar_ChargedIso_DR0p3To0p4 );
176  tmpTMVAReader->AddVariable( "ChargedIso_DR0p4To0p5", &fMVAVar_ChargedIso_DR0p4To0p5 );
177  tmpTMVAReader->AddVariable( "GammaIso_DR0p0To0p1", &fMVAVar_GammaIso_DR0p0To0p1 );
178  tmpTMVAReader->AddVariable( "GammaIso_DR0p1To0p2", &fMVAVar_GammaIso_DR0p1To0p2 );
179  tmpTMVAReader->AddVariable( "GammaIso_DR0p2To0p3", &fMVAVar_GammaIso_DR0p2To0p3 );
180  tmpTMVAReader->AddVariable( "GammaIso_DR0p3To0p4", &fMVAVar_GammaIso_DR0p3To0p4 );
181  tmpTMVAReader->AddVariable( "GammaIso_DR0p4To0p5", &fMVAVar_GammaIso_DR0p4To0p5 );
182  tmpTMVAReader->AddVariable( "NeutralHadronIso_DR0p0To0p1", &fMVAVar_NeutralHadronIso_DR0p0To0p1 );
183  tmpTMVAReader->AddVariable( "NeutralHadronIso_DR0p1To0p2", &fMVAVar_NeutralHadronIso_DR0p1To0p2 );
184  tmpTMVAReader->AddVariable( "NeutralHadronIso_DR0p2To0p3", &fMVAVar_NeutralHadronIso_DR0p2To0p3 );
185  tmpTMVAReader->AddVariable( "NeutralHadronIso_DR0p3To0p4", &fMVAVar_NeutralHadronIso_DR0p3To0p4 );
186  tmpTMVAReader->AddVariable( "NeutralHadronIso_DR0p4To0p5", &fMVAVar_NeutralHadronIso_DR0p4To0p5 );
187 
188  }
189 
190  if (type == kIsoDeltaR) {
191  tmpTMVAReader->AddVariable("PFCharged", &fMVAVar_MuRelIsoPFCharged );
192  tmpTMVAReader->AddVariable("PFNeutral", &fMVAVar_MuRelIsoPFNeutral );
193  tmpTMVAReader->AddVariable("PFPhotons", &fMVAVar_MuRelIsoPFPhotons );
194  tmpTMVAReader->AddVariable("SumDeltaR", &fMVAVar_MuDeltaRSum );
195  tmpTMVAReader->AddVariable("DeltaRMean", &fMVAVar_MuDeltaRMean );
196  tmpTMVAReader->AddVariable("Density", &fMVAVar_MuDensity );
197  }
198  if (type == kID) {
199  tmpTMVAReader->AddVariable( "TkNchi2", &fMVAVar_MuTkNchi2 );
200  if (i != 4) {
201  tmpTMVAReader->AddVariable( "GlobalNchi2", &fMVAVar_MuGlobalNchi2 );
202  tmpTMVAReader->AddVariable( "NValidHits", &fMVAVar_MuNValidHits );
203  }
204  tmpTMVAReader->AddVariable( "NTrackerHits", &fMVAVar_MuNTrackerHits );
205  tmpTMVAReader->AddVariable( "NPixelHits", &fMVAVar_MuNPixelHits );
206  if (i != 5) tmpTMVAReader->AddVariable( "NMatches", &fMVAVar_MuNMatches );
207  tmpTMVAReader->AddVariable( "TrkKink", &fMVAVar_MuTrkKink );
208  tmpTMVAReader->AddVariable( "SegmentCompatibility", &fMVAVar_MuSegmentCompatibility );
209  tmpTMVAReader->AddVariable( "CaloCompatibility", &fMVAVar_MuCaloCompatibility );
210  tmpTMVAReader->AddVariable( "HadEnergy", &fMVAVar_MuHadEnergy );
211  tmpTMVAReader->AddVariable( "EmEnergy", &fMVAVar_MuEmEnergy );
212  tmpTMVAReader->AddVariable( "HadS9Energy", &fMVAVar_MuHadS9Energy );
213  tmpTMVAReader->AddVariable( "EmS9Energy", &fMVAVar_MuEmS9Energy );
214  }
215 
216  tmpTMVAReader->BookMVA(fMethodname , weightsfiles[i]);
217  fTMVAReader.push_back(tmpTMVAReader);
218  }
219 }
220 
221 
222 //--------------------------------------------------------------------------------------------------
223 unsigned int MuonMVAEstimator::GetMVABin( double eta, double pt, Bool_t isGlobal, Bool_t isTrackerMuon) const {
224 
225  //Default is to return the first bin
226  unsigned int bin = 0;
227 
228  if (fMVAType == MuonMVAEstimator::kIsoRings) {
229  if (isGlobal && isTrackerMuon) {
230  if (pt < 10 && fabs(eta) < 1.479) bin = 0;
231  if (pt < 10 && fabs(eta) >= 1.479) bin = 1;
232  if (pt >= 10 && fabs(eta) < 1.479) bin = 2;
233  if (pt >= 10 && fabs(eta) >= 1.479) bin = 3;
234  }
235  else if (isTrackerMuon) bin = 4;
236  else if (isGlobal) bin = 5;
237 
238  }
239 
240  if (fMVAType == MuonMVAEstimator::kID) {
241  if (isGlobal && isTrackerMuon) {
242  if (pt < 10 && fabs(eta) < 1.479) bin = 0;
243  if (pt < 10 && fabs(eta) >= 1.479) bin = 1;
244  if (pt >= 10 && fabs(eta) < 1.479) bin = 2;
245  if (pt >= 10 && fabs(eta) >= 1.479) bin = 3;
246  }
247  else if (isTrackerMuon) bin = 4;
248  else if (isGlobal) bin = 5;
249 
250  }
251 
252  if (fMVAType == MuonMVAEstimator::kIDIsoRingsCombined ) {
253  bin = 0;
254  if (isGlobal && isTrackerMuon) {
255  if (pt < 10 && fabs(eta) < 1.479) bin = 0;
256  if (pt < 10 && fabs(eta) >= 1.479) bin = 1;
257  if (pt >= 10 && fabs(eta) < 1.479) bin = 2;
258  if (pt >= 10 && fabs(eta) >= 1.479) bin = 3;
259  }
260  else if (!isGlobal && isTrackerMuon) {
261  bin = 4;
262  }
263  else {
264  edm::LogWarning("NotTrackerMuon") << "Warning: Muon is not a tracker muon. Such muons are not supported.";
265  bin = 0;
266  }
267  }
268  if (fMVAType == MuonMVAEstimator::kIsoDeltaR){
269  if (pt < 20 && fabs(eta) < 1.479) bin = 0;
270  if (pt < 20 && fabs(eta) >= 1.479) bin = 1;
271  if (pt >= 20 && fabs(eta) < 1.479) bin = 2;
272  if (pt >= 20 && fabs(eta) >= 1.479) bin = 3;
273  }
274  return bin;
275 }
276 
277 void MuonMVAEstimator::bindVariables() {
278 
279  // this binding is needed for variables that sometime diverge.
280 
281  return;
282 }
283 
284 
285 Double_t MuonMVAEstimator::mvaValue_Iso( Double_t Pt,
286  Double_t Eta,
287  Bool_t isGlobalMuon,
288  Bool_t isTrackerMuon,
289  Double_t Rho,
290  MuonEffectiveArea::MuonEffectiveAreaTarget EATarget,
291  Double_t ChargedIso_DR0p0To0p1,
292  Double_t ChargedIso_DR0p1To0p2,
293  Double_t ChargedIso_DR0p2To0p3,
294  Double_t ChargedIso_DR0p3To0p4,
295  Double_t ChargedIso_DR0p4To0p5,
296  Double_t GammaIso_DR0p0To0p1,
297  Double_t GammaIso_DR0p1To0p2,
298  Double_t GammaIso_DR0p2To0p3,
299  Double_t GammaIso_DR0p3To0p4,
300  Double_t GammaIso_DR0p4To0p5,
301  Double_t NeutralHadronIso_DR0p0To0p1,
302  Double_t NeutralHadronIso_DR0p1To0p2,
303  Double_t NeutralHadronIso_DR0p2To0p3,
304  Double_t NeutralHadronIso_DR0p3To0p4,
305  Double_t NeutralHadronIso_DR0p4To0p5,
306  Bool_t printDebug) {
307 
308  if (!fisInitialized) {
309  edm::LogError("NotInitialized") << "Error: MuonMVAEstimator not properly initialized.";
310  return -9999;
311  }
312 
313  fMVAVar_ChargedIso_DR0p0To0p1 = TMath::Min((ChargedIso_DR0p0To0p1)/Pt, 2.5);
314  fMVAVar_ChargedIso_DR0p1To0p2 = TMath::Min((ChargedIso_DR0p1To0p2)/Pt, 2.5);
315  fMVAVar_ChargedIso_DR0p2To0p3 = TMath::Min((ChargedIso_DR0p2To0p3)/Pt, 2.5);
316  fMVAVar_ChargedIso_DR0p3To0p4 = TMath::Min((ChargedIso_DR0p3To0p4)/Pt, 2.5);
317  fMVAVar_ChargedIso_DR0p4To0p5 = TMath::Min((ChargedIso_DR0p4To0p5)/Pt, 2.5);
318  fMVAVar_GammaIso_DR0p0To0p1 = TMath::Max(TMath::Min((GammaIso_DR0p0To0p1 - Rho*MuonEffectiveArea::GetMuonEffectiveArea(MuonEffectiveArea::kMuGammaIsoDR0p0To0p1, Eta, EATarget))/Pt, 2.5), 0.0);
319  fMVAVar_GammaIso_DR0p1To0p2 = TMath::Max(TMath::Min((GammaIso_DR0p1To0p2 - Rho*MuonEffectiveArea::GetMuonEffectiveArea(MuonEffectiveArea::kMuGammaIsoDR0p1To0p2, Eta, EATarget))/Pt, 2.5), 0.0);
320  fMVAVar_GammaIso_DR0p2To0p3 = TMath::Max(TMath::Min((GammaIso_DR0p2To0p3 - Rho*MuonEffectiveArea::GetMuonEffectiveArea(MuonEffectiveArea::kMuGammaIsoDR0p2To0p3, Eta, EATarget))/Pt, 2.5), 0.0);
321  fMVAVar_GammaIso_DR0p3To0p4 = TMath::Max(TMath::Min((GammaIso_DR0p3To0p4 - Rho*MuonEffectiveArea::GetMuonEffectiveArea(MuonEffectiveArea::kMuGammaIsoDR0p3To0p4, Eta, EATarget))/Pt, 2.5), 0.0);
322  fMVAVar_GammaIso_DR0p4To0p5 = TMath::Max(TMath::Min((GammaIso_DR0p4To0p5 - Rho*MuonEffectiveArea::GetMuonEffectiveArea(MuonEffectiveArea::kMuGammaIsoDR0p4To0p5, Eta, EATarget))/Pt, 2.5), 0.0);
323  fMVAVar_NeutralHadronIso_DR0p0To0p1 = TMath::Max(TMath::Min((NeutralHadronIso_DR0p0To0p1 - Rho*MuonEffectiveArea::GetMuonEffectiveArea(MuonEffectiveArea::kMuNeutralHadronIsoDR0p0To0p1, Eta, EATarget))/Pt, 2.5), 0.0);
324  fMVAVar_NeutralHadronIso_DR0p1To0p2 = TMath::Max(TMath::Min((NeutralHadronIso_DR0p1To0p2 - Rho*MuonEffectiveArea::GetMuonEffectiveArea(MuonEffectiveArea::kMuNeutralHadronIsoDR0p1To0p2, Eta, EATarget))/Pt, 2.5), 0.0);
325  fMVAVar_NeutralHadronIso_DR0p2To0p3 = TMath::Max(TMath::Min((NeutralHadronIso_DR0p2To0p3 - Rho*MuonEffectiveArea::GetMuonEffectiveArea(MuonEffectiveArea::kMuNeutralHadronIsoDR0p2To0p3, Eta, EATarget))/Pt, 2.5), 0.0);
326  fMVAVar_NeutralHadronIso_DR0p3To0p4 = TMath::Max(TMath::Min((NeutralHadronIso_DR0p3To0p4 - Rho*MuonEffectiveArea::GetMuonEffectiveArea(MuonEffectiveArea::kMuNeutralHadronIsoDR0p3To0p4, Eta, EATarget))/Pt, 2.5), 0.0);
327  fMVAVar_NeutralHadronIso_DR0p4To0p5 = TMath::Max(TMath::Min((NeutralHadronIso_DR0p4To0p5 - Rho*MuonEffectiveArea::GetMuonEffectiveArea(MuonEffectiveArea::kMuNeutralHadronIsoDR0p4To0p5, Eta, EATarget))/Pt, 2.5), 0.0);
328 
329  // evaluate
330  Double_t mva = fTMVAReader[GetMVABin(Eta,Pt,isGlobalMuon,isTrackerMuon)]->EvaluateMVA(fMethodname);
331 
332  if(printDebug) {
333  LogDebug("MuonMVAEstimator") << " *** Inside the class fMethodname " << fMethodname << " fMVAType " << fMVAType << "\n"
334  << "ChargedIso ( 0.0 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 ): "
335  << fMVAVar_ChargedIso_DR0p0To0p1 << " "
336  << fMVAVar_ChargedIso_DR0p1To0p2 << " "
337  << fMVAVar_ChargedIso_DR0p2To0p3 << " "
338  << fMVAVar_ChargedIso_DR0p3To0p4 << " "
339  << fMVAVar_ChargedIso_DR0p4To0p5 << "\n"
340  << "PF Gamma Iso ( 0.0 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 ): "
341  << fMVAVar_GammaIso_DR0p0To0p1 << " "
342  << fMVAVar_GammaIso_DR0p1To0p2 << " "
343  << fMVAVar_GammaIso_DR0p2To0p3 << " "
344  << fMVAVar_GammaIso_DR0p3To0p4 << " "
345  << fMVAVar_GammaIso_DR0p4To0p5 << "\n"
346  << "PF Neutral Hadron Iso ( 0.0 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 ): "
347  << fMVAVar_NeutralHadronIso_DR0p0To0p1 << " "
348  << fMVAVar_NeutralHadronIso_DR0p1To0p2 << " "
349  << fMVAVar_NeutralHadronIso_DR0p2To0p3 << " "
350  << fMVAVar_NeutralHadronIso_DR0p3To0p4 << " "
351  << fMVAVar_NeutralHadronIso_DR0p4To0p5 << " "
352  << " ### MVA " << mva << endl;
353  }
354 
355  return mva;
356 
357 }
358 
359 Double_t MuonMVAEstimator::mvaValueIso(HWW& hww, Int_t mu, Double_t rho, MuonEffectiveArea::MuonEffectiveAreaTarget EATarget,
360  std::vector<Int_t> IdentifiedEle, std::vector<Int_t> IdentifiedMu, Bool_t printDebug ) {
361 
362 
363  if (!fisInitialized) {
364  edm::LogError("NotInitialized") << "Error: MuonMVAEstimator not properly initialized.\n";
365  return -9999;
366  }
367 
368 
369  Double_t mvavalue = -9999.;
370  if( hww.mus_trkidx().at(mu)<0 && hww.mus_sta_d0().at(mu)<0) return -9999;
371 
372 
373 /*
374  double rho = 0;
375  if (!(isnan(float(Rho)) || isinf(float(Rho)))) rho = Rho;
376 */
377 
378  // vertex selection
379  Int_t ivtx = 0;
380 
381  Double_t Pt = hww.mus_trkidx().at(mu)>=0 ? hww.mus_trk_p4().at(mu).pt() : hww.mus_sta_p4().at(mu).pt() ;
382  Double_t Eta = hww.mus_trkidx().at(mu)>=0 ? hww.mus_trk_p4().at(mu).eta() : hww.mus_sta_p4().at(mu).eta() ;
383  Bool_t isGlobalMuon = (hww.mus_type().at(mu)) & (1<<1);
384  Bool_t isTrackerMuon = (hww.mus_type().at(mu)) & (1<<2);
385  Double_t Rho = rho;
386  //MuonEffectiveArea::MuonEffectiveAreaTarget EATarget = EATarget;
387 
388  Double_t ChargedIso_DR0p0To0p1 = 0;
389  Double_t ChargedIso_DR0p1To0p2 = 0;
390  Double_t ChargedIso_DR0p2To0p3 = 0;
391  Double_t ChargedIso_DR0p3To0p4 = 0;
392  Double_t ChargedIso_DR0p4To0p5 = 0;
393  Double_t GammaIso_DR0p0To0p1 = 0;
394  Double_t GammaIso_DR0p1To0p2 = 0;
395  Double_t GammaIso_DR0p2To0p3 = 0;
396  Double_t GammaIso_DR0p3To0p4 = 0;
397  Double_t GammaIso_DR0p4To0p5 = 0;
398  Double_t NeutralHadronIso_DR0p0To0p1 = 0;
399  Double_t NeutralHadronIso_DR0p1To0p2 = 0;
400  Double_t NeutralHadronIso_DR0p2To0p3 = 0;
401  Double_t NeutralHadronIso_DR0p3To0p4 = 0;
402  Double_t NeutralHadronIso_DR0p4To0p5 = 0;
403 
404 
405  // Calcluate energy deposit in rings
406  for(unsigned int ipf = 0; ipf < hww.pfcands_p4().size(); ++ipf) {
407 
408 
409  // exclude muon itself
410  if ( hww.pfcands_trkidx().at(ipf)>=0 && hww.mus_trkidx().at(mu)>=0 &&
411  hww.pfcands_trkidx().at(ipf) == hww.mus_trkidx().at(mu) ) continue;
412 
413  //************************************************************
414  // New Isolation Calculations
415  //************************************************************
416  double dr = sqrt( pow( hww.pfcands_p4().at(ipf).eta() - hww.mus_trk_p4().at(mu).eta(), 2)
417  + pow(acos(cos(hww.pfcands_p4().at(ipf).phi() - hww.mus_trk_p4().at(mu).phi() )), 2) );
418 
419  if (dr < 1.0) {
420  Bool_t IsLeptonFootprint = kFALSE;
421 
422  //************************************************************
423  // Lepton Footprint Removal
424  //************************************************************
425 
426 
427  for (unsigned int iele=0; iele<IdentifiedEle.size(); iele++) {
428 
429  int ele = IdentifiedEle[iele];
430  // if pf candidate matches an electron passing ID cuts, then veto it
431  // pfcands_pfelsidx : index of PFElectron
432  // pfels_elsidx : index of els matching PFElectron
433  if( abs(hww.pfcands_particleId().at(ipf))==11 &&
434  hww.pfcands_pfelsidx().at(ipf)>=0 &&
435  hww.pfels_elsidx().at(hww.pfcands_pfelsidx().at(ipf))>=0 &&
436  ele == hww.pfels_elsidx().at(hww.pfcands_pfelsidx().at(ipf)) ) IsLeptonFootprint = true;
437  if( abs(hww.pfcands_particleId().at(ipf))==11 &&
438  hww.pfcands_pfelsidx().at(ipf)>=0 &&
439  hww.pfels_elsidx().at(hww.pfcands_pfelsidx().at(ipf))>=0 &&
440  hww.els_trkidx().at(ele) == hww.pfcands_trkidx().at(ipf) ) IsLeptonFootprint = true;
441 
442  //if pf candidate lies in veto regions of electron passing ID cuts, then veto it
443  double tmpDR = sqrt( pow(hww.pfcands_p4().at(ipf).eta() - hww.els_p4().at(ele).eta(),2)
444  + pow(acos(cos(hww.pfcands_p4().at(ipf).phi() - hww.els_p4().at(ele).phi())),2));
445  if( hww.pfcands_trkidx().at(ipf)>=0 && fabs(hww.els_etaSC().at(ele)) >= 1.479
446  && tmpDR < 0.015) IsLeptonFootprint = kTRUE;
447  if( translatePdgIdToType(hww.pfcands_particleId().at(ipf)) == kgamma && fabs(hww.els_etaSC().at(ele)) >= 1.479
448  && tmpDR < 0.08) IsLeptonFootprint = kTRUE;
449  }
450 
451  for (unsigned int imu=0; imu<IdentifiedMu.size(); imu++) {
452 
453  unsigned int idenmu = IdentifiedMu[imu];
454  //if pf candidate matches an muon passing ID cuts, then veto it
455  if ( hww.pfcands_trkidx().at(ipf)>=0 && hww.mus_trkidx().at(idenmu)>=0 &&
456  hww.pfcands_trkidx().at(ipf) == hww.mus_trkidx().at(idenmu) ) IsLeptonFootprint = true;
457 
458 
459  //if pf candidate lies in veto regions of muon passing ID cuts, then veto it
460  double tmpDR = sqrt( pow(hww.pfcands_p4().at(ipf).eta() - hww.mus_p4().at(idenmu).eta(),2)
461  + pow(acos(cos(hww.pfcands_p4().at(ipf).phi() - hww.mus_p4().at(idenmu).phi())),2));
462  if(hww.pfcands_trkidx().at(ipf)>=0 && tmpDR < 0.01) IsLeptonFootprint = kTRUE;
463  }
464 
465  if( !IsLeptonFootprint ) {
466 
467  if( hww.pfcands_trkidx().at(ipf)>=0) { //Charged
468 
469  // dZ cut
470  if ( (fabs( trks_dz_pv(hww, hww.pfcands_trkidx().at(ipf), ivtx ).first
471  - dzPVmu(hww.mus_vertex_p4().at(mu), hww.mus_trk_p4().at(mu), hww.vtxs_position().at(ivtx)) ) > 0.2) ) continue;
472 
473  // Veto any PFmuon, or PFEle
474  if ( translatePdgIdToType(hww.pfcands_particleId().at(ipf)) == ke
475  || translatePdgIdToType(hww.pfcands_particleId().at(ipf)) == kmu ) continue;
476 
477  // Footprint Veto
478  if (fabs(Eta) > 1.479 && dr < 0.01) continue;
479 
480  if (dr < 0.1) ChargedIso_DR0p0To0p1 += hww.pfcands_p4().at(ipf).pt();
481  if (dr >= 0.1 && dr < 0.2) ChargedIso_DR0p1To0p2 += hww.pfcands_p4().at(ipf).pt();
482  if (dr >= 0.2 && dr < 0.3) ChargedIso_DR0p2To0p3 += hww.pfcands_p4().at(ipf).pt();
483  if (dr >= 0.3 && dr < 0.4) ChargedIso_DR0p3To0p4 += hww.pfcands_p4().at(ipf).pt();
484  if (dr >= 0.4 && dr < 0.5) ChargedIso_DR0p4To0p5 += hww.pfcands_p4().at(ipf).pt();
485  }
486  else if ( translatePdgIdToType(hww.pfcands_particleId().at(ipf)) == kgamma ) { //Gamma
487  if (dr < 0.1) GammaIso_DR0p0To0p1 += hww.pfcands_p4().at(ipf).pt();
488  if (dr >= 0.1 && dr < 0.2) GammaIso_DR0p1To0p2 += hww.pfcands_p4().at(ipf).pt();
489  if (dr >= 0.2 && dr < 0.3) GammaIso_DR0p2To0p3 += hww.pfcands_p4().at(ipf).pt();
490  if (dr >= 0.3 && dr < 0.4) GammaIso_DR0p3To0p4 += hww.pfcands_p4().at(ipf).pt();
491  if (dr >= 0.4 && dr < 0.5) GammaIso_DR0p4To0p5 += hww.pfcands_p4().at(ipf).pt();
492  }
493  else { //NeutralHadron
494  if (dr < 0.1) NeutralHadronIso_DR0p0To0p1 += hww.pfcands_p4().at(ipf).pt();
495  if (dr >= 0.1 && dr < 0.2) NeutralHadronIso_DR0p1To0p2 += hww.pfcands_p4().at(ipf).pt();
496  if (dr >= 0.2 && dr < 0.3) NeutralHadronIso_DR0p2To0p3 += hww.pfcands_p4().at(ipf).pt();
497  if (dr >= 0.3 && dr < 0.4) NeutralHadronIso_DR0p3To0p4 += hww.pfcands_p4().at(ipf).pt();
498  if (dr >= 0.4 && dr < 0.5) NeutralHadronIso_DR0p4To0p5 += hww.pfcands_p4().at(ipf).pt();
499  }
500  } // IsLeptonFootprint
501  } // dR < 1.0
502  } // loop over PF candidates
503 
504 
505  // Get the mva value
506  mvavalue = MuonMVAEstimator::mvaValue_Iso(
507  Pt,
508  Eta,
509  isGlobalMuon,
510  isTrackerMuon,
511  Rho,
512  EATarget,
513  ChargedIso_DR0p0To0p1,
514  ChargedIso_DR0p1To0p2,
515  ChargedIso_DR0p2To0p3,
516  ChargedIso_DR0p3To0p4,
517  ChargedIso_DR0p4To0p5,
518  GammaIso_DR0p0To0p1,
519  GammaIso_DR0p1To0p2,
520  GammaIso_DR0p2To0p3,
521  GammaIso_DR0p3To0p4,
522  GammaIso_DR0p4To0p5,
523  NeutralHadronIso_DR0p0To0p1,
524  NeutralHadronIso_DR0p1To0p2,
525  NeutralHadronIso_DR0p2To0p3,
526  NeutralHadronIso_DR0p3To0p4,
527  NeutralHadronIso_DR0p4To0p5,
528  printDebug);
529 
530  return mvavalue;
531 
532 }
533 
534 
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:501
type
type
Definition: HCALResponse.h:21
i
int i
Definition: DBlmapReader.cc:9
MuonMVAEstimator::fMVAVar_MuNPixelHits
Float_t fMVAVar_MuNPixelHits
Definition: MuonMVAEstimator.h:103
HWW::pfels_elsidx
std::vector< int > & pfels_elsidx()
Definition: HWW.cc:749
MuonMVAEstimator::fMVAVar_MuNValidHits
Float_t fMVAVar_MuNValidHits
Definition: MuonMVAEstimator.h:101
MuonMVAEstimator::fMVAVar_MuCaloCompatibility
Float_t fMVAVar_MuCaloCompatibility
Definition: MuonMVAEstimator.h:110
MuonMVAEstimator::fMVAVar_ChargedIso_DR0p2To0p3
Float_t fMVAVar_ChargedIso_DR0p2To0p3
Definition: MuonMVAEstimator.h:119
MuonMVAEstimator::fTMVAReader
std::vector< TMVA::Reader * > fTMVAReader
Definition: MuonMVAEstimator.h:85
MuonMVAEstimator::fMVAVar_NeutralHadronIso_DR0p1To0p2
Float_t fMVAVar_NeutralHadronIso_DR0p1To0p2
Definition: MuonMVAEstimator.h:128
MuonMVAEstimator::fMVAVar_ChargedIso_DR0p3To0p4
Float_t fMVAVar_ChargedIso_DR0p3To0p4
Definition: MuonMVAEstimator.h:120
rho
Definition: DDAxes.h:10
MuonMVAEstimator::fMVAVar_ChargedIso_DR0p1To0p2
Float_t fMVAVar_ChargedIso_DR0p1To0p2
Definition: MuonMVAEstimator.h:118
MuonMVAEstimator::fMVAVar_GammaIso_DR0p4To0p5
Float_t fMVAVar_GammaIso_DR0p4To0p5
Definition: MuonMVAEstimator.h:126
HWW::mus_trkidx
std::vector< int > & mus_trkidx()
Definition: HWW.cc:535
MuonMVAEstimator::fMVAVar_GammaIso_DR0p3To0p4
Float_t fMVAVar_GammaIso_DR0p3To0p4
Definition: MuonMVAEstimator.h:125
MuonMVAEstimator::fMethodname
std::string fMethodname
Definition: MuonMVAEstimator.h:86
HWW::mus_p4
std::vector< LorentzVector > & mus_p4()
Definition: HWW.cc:367
HWW::pfcands_p4
std::vector< LorentzVector > & pfcands_p4()
Definition: HWW.cc:725
MuonMVAEstimator::fMVAVar_GammaIso_DR0p1To0p2
Float_t fMVAVar_GammaIso_DR0p1To0p2
Definition: MuonMVAEstimator.h:123
Min
T Min(T a, T b)
Definition: MathUtil.h:39
eta
T eta() const
Definition: Basic3DVectorLD.h:177
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
MuonMVAEstimator::fMVAVar_MuSegmentCompatibility
Float_t fMVAVar_MuSegmentCompatibility
Definition: MuonMVAEstimator.h:109
HWW::pfcands_trkidx
std::vector< int > & pfcands_trkidx()
Definition: HWW.cc:729
HWW::pfcands_pfelsidx
std::vector< int > & pfcands_pfelsidx()
Definition: HWW.cc:737
MuonMVAEstimator::GetMVABin
UInt_t GetMVABin(double eta, double pt, Bool_t isGlobal, Bool_t isTrackerMuon) const
Definition: MuonMVAEstimator.cc:223
MuonMVAEstimator::fMVAVar_ChargedIso_DR0p0To0p1
Float_t fMVAVar_ChargedIso_DR0p0To0p1
Definition: MuonMVAEstimator.h:117
HWW::mus_trk_p4
std::vector< LorentzVector > & mus_trk_p4()
Definition: HWW.cc:371
MuonMVAEstimator::mvaValue_Iso
Double_t mvaValue_Iso(Double_t Pt, Double_t Eta, Bool_t isGlobalMuon, Bool_t isTrackerMuon, Double_t Rho, MuonEffectiveArea::MuonEffectiveAreaTarget EATarget, Double_t ChargedIso_DR0p0To0p1, Double_t ChargedIso_DR0p1To0p2, Double_t ChargedIso_DR0p2To0p3, Double_t ChargedIso_DR0p3To0p4, Double_t ChargedIso_DR0p4To0p5, Double_t GammaIso_DR0p0To0p1, Double_t GammaIso_DR0p1To0p2, Double_t GammaIso_DR0p2To0p3, Double_t GammaIso_DR0p3To0p4, Double_t GammaIso_DR0p4To0p5, Double_t NeutralHadronIso_DR0p0To0p1, Double_t NeutralHadronIso_DR0p1To0p2, Double_t NeutralHadronIso_DR0p2To0p3, Double_t NeutralHadronIso_DR0p3To0p4, Double_t NeutralHadronIso_DR0p4To0p5, Bool_t printDebug=kFALSE)
Definition: MuonMVAEstimator.cc:285
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
p4
double p4[4]
Definition: TauolaWrapper.h:92
ParticleType
ParticleType
Definition: MuonMVAEstimator.cc:21
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
HWW::pfcands_particleId
std::vector< int > & pfcands_particleId()
Definition: HWW.cc:733
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
MuonMVAEstimator::mvaValueIso
Double_t mvaValueIso(HWW &, Int_t imu, Double_t rho, MuonEffectiveArea::MuonEffectiveAreaTarget EATarget, std::vector< Int_t > IdentifiedEle, std::vector< Int_t > IdentifiedMu, Bool_t printDebug)
Definition: MuonMVAEstimator.cc:359
HWW::vtxs_position
std::vector< LorentzVector > & vtxs_position()
Definition: HWW.cc:5
MuonEffectiveArea::GetMuonEffectiveArea
static const Double_t GetMuonEffectiveArea(MuonEffectiveAreaType type, Double_t SCEta, MuonEffectiveAreaTarget EffectiveAreaTarget=kMuEAData2011)
Definition: MuonEffectiveArea.h:69
MuonMVAEstimator::fMVAVar_ChargedIso_DR0p4To0p5
Float_t fMVAVar_ChargedIso_DR0p4To0p5
Definition: MuonMVAEstimator.h:121
RPCpg::mu
const int mu
Definition: Constants.h:22
first
bool first
Definition: L1TdeRCT.cc:75
HWW::els_p4
std::vector< LorentzVector > & els_p4()
Definition: HWW.cc:101
HWW::els_etaSC
std::vector< float > & els_etaSC()
Definition: HWW.cc:117
dzPVmu
double dzPVmu(const LorentzVector &vtx, const LorentzVector &p4, const LorentzVector &pv)
Definition: MuonMVAEstimator.cc:14
MuonMVAEstimator::fMVAVar_NeutralHadronIso_DR0p3To0p4
Float_t fMVAVar_NeutralHadronIso_DR0p3To0p4
Definition: MuonMVAEstimator.h:130
MuonMVAEstimator::fMVAVar_NeutralHadronIso_DR0p2To0p3
Float_t fMVAVar_NeutralHadronIso_DR0p2To0p3
Definition: MuonMVAEstimator.h:129
HWW
Definition: HWW.h:12
Max
T Max(T a, T b)
Definition: MathUtil.h:44
MuonMVAEstimator::fMVAVar_MuHadS9Energy
Float_t fMVAVar_MuHadS9Energy
Definition: MuonMVAEstimator.h:113
MuonMVAEstimator::initialize
void initialize(std::string methodName, std::string weightsfile, MuonMVAEstimator::MVAType type)
Definition: MuonMVAEstimator.cc:74
HWW::mus_sta_d0
std::vector< float > & mus_sta_d0()
Definition: HWW.cc:483
MuonMVAEstimator::fMVAVar_MuDeltaRMean
Float_t fMVAVar_MuDeltaRMean
Definition: MuonMVAEstimator.h:138
MuonMVAEstimator::fMVAVar_MuTypeBits
Float_t fMVAVar_MuTypeBits
Definition: MuonMVAEstimator.h:96
HWWFunctions::trks_dz_pv
std::pair< double, double > trks_dz_pv(HWW &, int itrk, int ipv)
Definition: trackSelections.cc:70
MuonMVAEstimator::fMVAVar_GammaIso_DR0p0To0p1
Float_t fMVAVar_GammaIso_DR0p0To0p1
Definition: MuonMVAEstimator.h:122
MuonMVAEstimator::fMVAVar_MuRelIsoPFCharged
Float_t fMVAVar_MuRelIsoPFCharged
Definition: MuonMVAEstimator.h:135
LorentzVector
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< float > > LorentzVector
Definition: analysisEnums.h:9
MuonMVAEstimator::fMVAVar_NeutralHadronIso_DR0p0To0p1
Float_t fMVAVar_NeutralHadronIso_DR0p0To0p1
Definition: MuonMVAEstimator.h:127
MuonMVAEstimator::fMVAVar_NeutralHadronIso_DR0p4To0p5
Float_t fMVAVar_NeutralHadronIso_DR0p4To0p5
Definition: MuonMVAEstimator.h:131
MuonMVAEstimator::fMVAVar_MuGlobalNchi2
Float_t fMVAVar_MuGlobalNchi2
Definition: MuonMVAEstimator.h:100
MuonMVAEstimator::fMVAVar_MuRelIsoPFNeutral
Float_t fMVAVar_MuRelIsoPFNeutral
Definition: MuonMVAEstimator.h:136
MuonMVAEstimator::fMVAVar_GammaIso_DR0p2To0p3
Float_t fMVAVar_GammaIso_DR0p2To0p3
Definition: MuonMVAEstimator.h:124
translatePdgIdToType
int translatePdgIdToType(int pdgid)
Definition: MuonMVAEstimator.cc:35
HWW::mus_vertex_p4
std::vector< LorentzVector > & mus_vertex_p4()
Definition: HWW.cc:375
MuonMVAEstimator::fMVAVar_MuRelIsoPFPhotons
Float_t fMVAVar_MuRelIsoPFPhotons
Definition: MuonMVAEstimator.h:137
HWW::mus_sta_p4
std::vector< LorentzVector > & mus_sta_p4()
Definition: HWW.cc:379
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
MuonMVAEstimator::fMVAVar_MuEmS9Energy
Float_t fMVAVar_MuEmS9Energy
Definition: MuonMVAEstimator.h:114
MuonMVAEstimator::fMVAVar_MuNTrackerHits
Float_t fMVAVar_MuNTrackerHits
Definition: MuonMVAEstimator.h:102
HWW::els_trkidx
std::vector< int > & els_trkidx()
Definition: HWW.cc:341
SQLiteCheck_cfg.printDebug
tuple printDebug
Definition: SQLiteCheck_cfg.py:37
HWW::mus_type
std::vector< int > & mus_type()
Definition: HWW.cc:555
Generated for CMSSW Reference Manual by doxygen 1.8.5