ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/JetIDMVA.cc
Revision: 1.2
Committed: Wed Apr 4 07:57:01 2012 UTC (13 years, 1 month ago) by pharris
Content type: text/plain
Branch: MAIN
Changes since 1.1: +98 -0 lines
Log Message:
Fixed Jet MVA

File Contents

# User Rev Content
1 pharris 1.1 #include "MitPhysics/Utils/interface/JetIDMVA.h"
2     #include "MitPhysics/Utils/interface/JetTools.h"
3     #include "MitAna/DataTree/interface/StableData.h"
4     #include <TFile.h>
5     #include <TRandom3.h>
6     #include "TMVA/Tools.h"
7     #include "TMVA/Reader.h"
8    
9    
10     ClassImp(mithep::JetIDMVA)
11    
12     using namespace mithep;
13    
14     //--------------------------------------------------------------------------------------------------
15     JetIDMVA::JetIDMVA() :
16     fMethodName("JetIDMA"),
17     fIsInitialized(kFALSE),
18     fJetPtMin(10), //We need to lower this
19     fNPV (0),
20     fJPt1 (0),
21     fJEta1 (0),
22     fJPhi1 (0),
23     fJD01 (0),
24     fJDZ1 (0),
25     fJM1 (0),
26     fNPart1 (0),
27     fLPt1 (0),
28     fLEta1 (0),
29     fLPhi1 (0),
30     fSPt1 (0),
31     fSEta1 (0),
32     fSPhi1 (0),
33     fNEPt1 (0),
34     fNEEta1 (0),
35     fNEPhi1 (0),
36     fEMPt1 (0),
37     fEMEta1 (0),
38     fEMPhi1 (0),
39     fChPt1 (0),
40     fChPhi1 (0),
41     fLFr1 (0),
42     fDRLC1 (0),
43     fDRLS1 (0),
44     fDRM1 (0),
45     fDRNE1 (0),
46     fDREM1 (0),
47     fDRCH1 (0)
48     {
49     fReader = 0;
50     }
51     //--------------------------------------------------------------------------------------------------
52     JetIDMVA::~JetIDMVA()
53     {
54     fReader = 0;
55     }
56    
57     //--------------------------------------------------------------------------------------------------
58     void JetIDMVA::Initialize( TString iMethodName,
59     TString iWeights,
60     JetIDMVA::MVAType iType)
61     {
62    
63     fIsInitialized = kTRUE;
64     fMethodName = iMethodName;
65     fType = iType;
66     fReader = 0;
67     fReader = new TMVA::Reader( "!Color:!Silent:Error" );
68     if (fType == kBaseline) {
69     fReader->AddSpectator( "npv" , &fNPV );
70     fReader->AddVariable( "jspt_1" , &fJPt1 );
71     fReader->AddVariable( "jseta_1" , &fJEta1 );
72     fReader->AddVariable( "jsphi_1" , &fJPhi1 );
73     fReader->AddVariable( "jd0_1" , &fJD01 );
74     fReader->AddVariable( "jdZ_1" , &fJDZ1 );
75     fReader->AddVariable( "jm_1" , &fJM1 );
76     fReader->AddVariable( "npart_1" , &fNPart1 );
77     fReader->AddVariable( "lpt_1" , &fLPt1 );
78     fReader->AddVariable( "leta_1" , &fLEta1 );
79     fReader->AddVariable( "lphi_1" , &fLPhi1 );
80     fReader->AddVariable( "spt_1" , &fSPt1 );
81     fReader->AddVariable( "seta_1" , &fSEta1 );
82     fReader->AddVariable( "sphi_1" , &fSPhi1 );
83     fReader->AddVariable( "lnept_1" , &fNEPt1 );
84     fReader->AddVariable( "lneeta_1" , &fNEEta1 );
85     fReader->AddVariable( "lnephi_1" , &fNEPhi1 );
86     fReader->AddVariable( "lempt_1" , &fEMPt1 );
87     fReader->AddVariable( "lemeta_1" , &fEMEta1 );
88     fReader->AddVariable( "lemphi_1" , &fEMPhi1 );
89     fReader->AddVariable( "lchpt_1" , &fChPt1 );
90     fReader->AddVariable( "lchphi_1" , &fChPhi1 );
91     fReader->AddVariable( "lLfr_1" , &fLFr1 );
92     fReader->AddVariable( "drlc_1" , &fDRLC1 );
93     fReader->AddVariable( "drls_1" , &fDRLS1 );
94     fReader->AddVariable( "drm_1" , &fDRM1 );
95     fReader->AddVariable( "drmne_1" , &fDRNE1 );
96     fReader->AddVariable( "drem_1" , &fDREM1 );
97     fReader->AddVariable( "drch_1" , &fDRCH1 );
98     }
99     fReader->BookMVA(fMethodName , iWeights );
100     std::cout << "Jet ID MVA Initialization\n";
101     std::cout << "MethodName : " << fMethodName << " , type == " << fType << std::endl;
102     }
103    
104     //--------------------------------------------------------------------------------------------------
105     Double_t JetIDMVA::MVAValue(
106     Float_t iNPV ,
107     Float_t iJPt1 ,
108     Float_t iJEta1 ,
109     Float_t iJPhi1 ,
110     Float_t iJD01 ,
111     Float_t iJDZ1 ,
112     Float_t iJM1 ,
113     Float_t iNPart1 ,
114     Float_t iLPt1 ,
115     Float_t iLEta1 ,
116     Float_t iLPhi1 ,
117     Float_t iSPt1 ,
118     Float_t iSEta1 ,
119     Float_t iSPhi1 ,
120     Float_t iNEPt1 ,
121     Float_t iNEEta1 ,
122     Float_t iNEPhi1 ,
123     Float_t iEMPt1 ,
124     Float_t iEMEta1 ,
125     Float_t iEMPhi1 ,
126     Float_t iChPt1 ,
127     Float_t iChPhi1 ,
128     Float_t iLFr1 ,
129     Float_t iDRLC1 ,
130     Float_t iDRLS1 ,
131     Float_t iDRM1 ,
132     Float_t iDRNE1 ,
133     Float_t iDREM1 ,
134     Float_t iDRCH1
135     ){
136    
137     if (!fIsInitialized) {
138     std::cout << "Error: JetIDMVA not properly initialized.\n";
139     return -9999;
140     }
141    
142     fNPV = iNPV;
143     fJPt1 = iJPt1;
144     fJEta1 = iJEta1;
145     fJPhi1 = fJPhi1;
146     fJD01 = iJD01;
147     fJDZ1 = iJDZ1;
148     fJM1 = iJM1 ;
149     fNPart1 = iNPart1;
150     fLPt1 = iLPt1;
151     fLEta1 = iLEta1;
152     fLPhi1 = iLPhi1;
153     fSPt1 = iSPt1;
154     fSEta1 = iSEta1;
155     fSPhi1 = iSPhi1;
156     fNEPt1 = iNEPt1;
157     fNEEta1 = iNEEta1;
158     fNEPhi1 = iNEPhi1;
159     fEMPt1 = iEMPt1;
160     fEMEta1 = iEMEta1;
161     fEMPhi1 = iEMPhi1;
162     fChPt1 = iChPt1;
163     fChPhi1 = iChPhi1;
164     fLFr1 = iLFr1;
165     fDRLC1 = iDRLC1;
166     fDRLS1 = iDRLS1;
167     fDRM1 = iDRM1;
168     fDRNE1 = iDRNE1;
169     fDREM1 = iDREM1;
170     fDRCH1 = iDRCH1;
171    
172     Double_t lMVA = -9999;
173     lMVA = fReader->EvaluateMVA( fMethodName );
174    
175     return lMVA;
176     }
177     //--------------------------------------------------------------------------------------------------
178     Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,
179     FactorizedJetCorrector *iJetCorrector,
180     const PileupEnergyDensityCol *iPileupEnergyDensity) {
181     if(!JetTools::passPFLooseId(iJet)) return false;
182     if(correctedPt(iJet,iJetCorrector,iPileupEnergyDensity) < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false; //This line is a bug in the Met training
183     if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
184     double lMVA = MVAValue(iJet,iVertex,iJetCorrector,iPileupEnergyDensity);
185     if(lMVA < -0.8) return false;
186     if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
187     return true;
188     }
189     //--------------------------------------------------------------------------------------------------
190 pharris 1.2 Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex) {
191     if(!JetTools::passPFLooseId(iJet)) return false;
192     if(iJet->Pt() < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false; //This line is a bug in the Met training
193     if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
194     double lMVA = MVAValue(iJet,iVertex);
195     if(lMVA < -0.8) return false;
196     if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
197     return true;
198     }
199     //--------------------------------------------------------------------------------------------------
200 pharris 1.1 Double_t JetIDMVA::MVAValue(const PFJet *iJet,const Vertex *iVertex, //Vertex here is the PV
201     FactorizedJetCorrector *iJetCorrector,
202     const PileupEnergyDensityCol *iPileupEnergyDensity,
203     Bool_t printDebug) {
204    
205     if (!fIsInitialized) {
206     std::cout << "Error: JetIDMVA not properly initialized.\n";
207     return -9999;
208     }
209     if(!JetTools::passPFLooseId(iJet)) return -2.;
210    
211     //set all input variables
212     fJPt1 = correctedPt(iJet,iJetCorrector,iPileupEnergyDensity);
213     fJEta1 = iJet->RawMom().Eta();
214     fJPhi1 = iJet->RawMom().Phi();
215     fJM1 = iJet->Mass();
216    
217     const mithep::PFCandidate *lLead = JetTools::leadCand(iJet,-1);
218     const mithep::PFCandidate *lSecond = JetTools::leadCand(iJet,-1,true);
219     const mithep::PFCandidate *lLeadNeut = JetTools::leadCand(iJet ,5);
220     const mithep::PFCandidate *lLeadEm = JetTools::leadCand(iJet ,4);
221     const mithep::PFCandidate *lLeadCh = JetTools::leadCand(iJet ,1);
222    
223     fJD01 = JetTools::impactParameter(iJet,iVertex);
224     fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
225     fNPart1 = iJet->NPFCands();
226     fLPt1 = lLead ->Pt();
227     fLEta1 = lLead ->Eta();
228     fLPhi1 = lLead ->Phi();
229     fSPt1 = lSecond ->Pt();
230     fSEta1 = lSecond ->Eta();
231     fSPhi1 = lSecond ->Phi();
232     fNEPt1 = lLeadNeut->Pt();
233     fNEEta1 = lLeadNeut->Eta();
234     fNEPhi1 = lLeadNeut->Phi();
235     fEMPt1 = lLeadEm ->Pt();
236     fEMEta1 = lLeadEm ->Eta();
237     fEMPhi1 = lLeadEm ->Phi();
238     fChPt1 = lLeadCh ->Pt();
239     //fChEta1 = lLeadCh ->Eta();
240     fChPhi1 = lLeadCh ->Phi();
241     fLFr1 = lLead->Pt()/iJet->Pt();
242    
243     fDRLC1 = MathUtils::DeltaR(iJet->Mom(),lLead ->Mom());
244     fDRLS1 = MathUtils::DeltaR(iJet->Mom(),lSecond->Mom());
245     fDRM1 = JetTools::dRMean (iJet,-1);
246     fDRNE1 = JetTools::dRMean (iJet, 5);
247     fDREM1 = JetTools::dRMean (iJet, 4);
248     fDRCH1 = JetTools::dRMean (iJet, 1);
249    
250     double lMVA = fReader->EvaluateMVA( fMethodName );
251    
252     if (printDebug == kTRUE) {
253     std::cout << "Debug Jet MVA: "
254     << fNPV << " "
255     << fJPt1 << " "
256     << fJEta1 << " "
257     << fJPhi1 << " "
258     << fJD01 << " "
259     << fJDZ1 << " "
260     << fJM1 << " "
261     << fNPart1 << " "
262     << fLPt1 << " "
263     << fLEta1 << " "
264     << fLPhi1 << " "
265     << fSPt1 << " "
266     << fSEta1 << " "
267     << fSPhi1 << " "
268     << fNEPt1 << " "
269     << fNEEta1 << " "
270     << fNEPhi1 << " "
271     << fEMPt1 << " "
272     << fEMEta1 << " "
273     << fEMPhi1 << " "
274     << fChPt1 << " "
275     << fChPhi1 << " "
276     << fLFr1 << " "
277     << fDRLC1 << " "
278     << fDRLS1 << " "
279     << fDRM1 << " "
280     << fDRNE1 << " "
281     << fDREM1 << " "
282     << fDRCH1 << " "
283     << " === : === "
284     << lMVA << " "
285     << std::endl;
286     }
287    
288     return lMVA;
289     }
290 pharris 1.2 Double_t JetIDMVA::MVAValue(const PFJet *iJet,const Vertex *iVertex, //Vertex here is the PV
291     Bool_t printDebug) {
292    
293     if (!fIsInitialized) {
294     std::cout << "Error: JetIDMVA not properly initialized.\n";
295     return -9999;
296     }
297     if(!JetTools::passPFLooseId(iJet)) return -2.;
298    
299     //set all input variables
300     fJPt1 = iJet->Pt();
301     fJEta1 = iJet->RawMom().Eta();
302     fJPhi1 = iJet->RawMom().Phi();
303     fJM1 = iJet->Mass();
304    
305     const mithep::PFCandidate *lLead = JetTools::leadCand(iJet,-1);
306     const mithep::PFCandidate *lSecond = JetTools::leadCand(iJet,-1,true);
307     const mithep::PFCandidate *lLeadNeut = JetTools::leadCand(iJet ,5);
308     const mithep::PFCandidate *lLeadEm = JetTools::leadCand(iJet ,4);
309     const mithep::PFCandidate *lLeadCh = JetTools::leadCand(iJet ,1);
310    
311     fJD01 = JetTools::impactParameter(iJet,iVertex);
312     fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
313     fNPart1 = iJet->NPFCands();
314     fLPt1 = lLead ->Pt();
315     fLEta1 = lLead ->Eta();
316     fLPhi1 = lLead ->Phi();
317     fSPt1 = lSecond ->Pt();
318     fSEta1 = lSecond ->Eta();
319     fSPhi1 = lSecond ->Phi();
320     fNEPt1 = lLeadNeut->Pt();
321     fNEEta1 = lLeadNeut->Eta();
322     fNEPhi1 = lLeadNeut->Phi();
323     fEMPt1 = lLeadEm ->Pt();
324     fEMEta1 = lLeadEm ->Eta();
325     fEMPhi1 = lLeadEm ->Phi();
326     fChPt1 = lLeadCh ->Pt();
327     //fChEta1 = lLeadCh ->Eta();
328     fChPhi1 = lLeadCh ->Phi();
329     fLFr1 = lLead->Pt()/iJet->RawMom().Pt();
330    
331     fDRLC1 = MathUtils::DeltaR(iJet->Mom(),lLead ->Mom());
332     fDRLS1 = MathUtils::DeltaR(iJet->Mom(),lSecond->Mom());
333     fDRM1 = JetTools::dRMean (iJet,-1);
334     fDRNE1 = JetTools::dRMean (iJet, 5);
335     fDREM1 = JetTools::dRMean (iJet, 4);
336     fDRCH1 = JetTools::dRMean (iJet, 1);
337    
338     double lMVA = fReader->EvaluateMVA( fMethodName );
339    
340     if (printDebug == kTRUE) {
341     std::cout << "Debug Jet MVA: "
342     << fNPV << " "
343     << fJPt1 << " "
344     << fJEta1 << " "
345     << fJPhi1 << " "
346     << fJD01 << " "
347     << fJDZ1 << " "
348     << fJM1 << " "
349     << fNPart1 << " "
350     << fLPt1 << " "
351     << fLEta1 << " "
352     << fLPhi1 << " "
353     << fSPt1 << " "
354     << fSEta1 << " "
355     << fSPhi1 << " "
356     << fNEPt1 << " "
357     << fNEEta1 << " "
358     << fNEPhi1 << " "
359     << fEMPt1 << " "
360     << fEMEta1 << " "
361     << fEMPhi1 << " "
362     << fChPt1 << " "
363     << fChPhi1 << " "
364     << fLFr1 << " "
365     << fDRLC1 << " "
366     << fDRLS1 << " "
367     << fDRM1 << " "
368     << fDRNE1 << " "
369     << fDREM1 << " "
370     << fDRCH1 << " "
371     << " === : === "
372     << lMVA << " "
373     << std::endl;
374     }
375    
376     return lMVA;
377     }
378 pharris 1.1 Double_t JetIDMVA::correctedPt(const PFJet *iJet, FactorizedJetCorrector *iJetCorrector,const PileupEnergyDensityCol *iPUEnergyDensity) {
379     const FourVectorM rawMom = iJet->RawMom();
380     iJetCorrector->setJetEta(rawMom.Eta());
381     iJetCorrector->setJetPt (rawMom.Pt());
382     iJetCorrector->setJetPhi(rawMom.Phi());
383     iJetCorrector->setJetE (rawMom.E());
384     iJetCorrector->setRho (iPUEnergyDensity->At(0)->RhoHighEta());
385     iJetCorrector->setJetA (iJet->JetArea());
386     iJetCorrector->setJetEMF(-99.0);
387     Double_t correction = iJetCorrector->getCorrection();
388     return rawMom.Pt()*correction;
389     }