ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/JetIDMVA.cc
Revision: 1.21
Committed: Tue Jul 16 22:29:52 2013 UTC (11 years, 9 months ago) by pharris
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.20: +4 -3 lines
Log Message:
synch for 42X

File Contents

# User Rev Content
1 pharris 1.4 #include "FWCore/PythonParameterSet/interface/MakeParameterSets.h"
2     #include "FWCore/ParameterSet/interface/ParameterSet.h"
3 pharris 1.1 #include "MitPhysics/Utils/interface/JetIDMVA.h"
4     #include "MitPhysics/Utils/interface/JetTools.h"
5     #include "MitAna/DataTree/interface/StableData.h"
6     #include <TFile.h>
7     #include <TRandom3.h>
8     #include "TMVA/Tools.h"
9     #include "TMVA/Reader.h"
10    
11    
12     ClassImp(mithep::JetIDMVA)
13    
14     using namespace mithep;
15    
16     //--------------------------------------------------------------------------------------------------
17     JetIDMVA::JetIDMVA() :
18 pharris 1.8 fJetPtMin(0.) , //We need to lower this
19 pharris 1.4 fDZCut (0.2),
20     fLowPtMethodName ("JetIDMVALowPt" ),
21     fHighPtMethodName("JetIDMVAHighPt"),
22 pharris 1.1 fIsInitialized(kFALSE),
23 pharris 1.4 fNVtx (0),
24     fJPt1 (0),
25     fJEta1 (0),
26     fJPhi1 (0),
27     fJD01 (0),
28     fJDZ1 (0),
29     fBeta (0),
30     fBetaStar (0),
31     fNCharged (0),
32     fNNeutrals(0),
33     fDRMean (0),
34 pharris 1.17 fPtD (0),
35 pharris 1.4 fFrac01 (0),
36     fFrac02 (0),
37     fFrac03 (0),
38     fFrac04 (0),
39 pharris 1.11 fFrac05 (0),
40     fDR2Mean (0)
41 pharris 1.1 {
42     fReader = 0;
43     }
44     //--------------------------------------------------------------------------------------------------
45 pharris 1.4 JetIDMVA::~JetIDMVA() {
46    
47 pharris 1.11 delete fReader;
48     delete fLowPtReader;
49 pharris 1.1 }
50    
51     //--------------------------------------------------------------------------------------------------
52 pharris 1.4 void JetIDMVA::Initialize( JetIDMVA::CutType iCutType,
53     TString iLowPtWeights,
54     TString iHighPtWeights,
55     JetIDMVA::MVAType iType,
56 pharris 1.21 TString iCutFileName,bool i42)
57 pharris 1.1 {
58    
59     fIsInitialized = kTRUE;
60     fType = iType;
61 pharris 1.4 fCutType = iCutType;
62 pharris 1.21 f42 = i42;
63    
64 bendavid 1.14 std::string lCutId = "JetIdParams";
65 pharris 1.17 if(fType == k42) lCutId = "PuJetIdOptMVA_wp";
66     if(fType == k52) lCutId = "full_5x_wp";
67     if(fType == kCut) lCutId = "PuJetIdCutBased_wp";
68     if(fType == k53) lCutId = "full_53x_wp";
69     if(fType == k53MET) lCutId = "met_53x_wp";
70    
71     //Load Cut Matrix
72 bendavid 1.14 edm::ParameterSet lDefConfig = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>("JetIdParams");
73     edm::ParameterSet lConfig = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>(lCutId);
74     std::string lCutType = "Tight";
75     if(fCutType == kMedium) lCutType = "Medium";
76     if(fCutType == kLoose ) lCutType = "Loose";
77     if(fCutType == kMET ) lCutType = "MET";
78     if(fType != kCut) {
79     std::string lLowPtCut = "MET";
80     std::vector<double> lPt010 = lDefConfig.getParameter<std::vector<double> >(("Pt010_" +lLowPtCut).c_str());
81     std::vector<double> lPt1020 = lConfig.getParameter<std::vector<double> >(("Pt1020_"+lCutType).c_str());
82     std::vector<double> lPt2030 = lConfig.getParameter<std::vector<double> >(("Pt2030_"+lCutType).c_str());
83     std::vector<double> lPt3050 = lConfig.getParameter<std::vector<double> >(("Pt3050_"+lCutType).c_str());
84     for(int i0 = 0; i0 < 4; i0++) fMVACut[0][i0] = lPt010 [i0];
85     for(int i0 = 0; i0 < 4; i0++) fMVACut[1][i0] = lPt1020[i0];
86     for(int i0 = 0; i0 < 4; i0++) fMVACut[2][i0] = lPt2030[i0];
87     for(int i0 = 0; i0 < 4; i0++) fMVACut[3][i0] = lPt3050[i0];
88     }
89     if(fType == kCut) {
90     for(int i0 = 0; i0 < 2; i0++) {
91     std::string lFullCutType = lCutType;
92     if(i0 == 0) lFullCutType = "BetaStar"+ lCutType;
93     if(i0 == 1) lFullCutType = "RMS" + lCutType;
94     std::vector<double> pt010 = lConfig.getParameter<std::vector<double> >(("Pt010_" +lFullCutType).c_str());
95     std::vector<double> pt1020 = lConfig.getParameter<std::vector<double> >(("Pt1020_"+lFullCutType).c_str());
96     std::vector<double> pt2030 = lConfig.getParameter<std::vector<double> >(("Pt2030_"+lFullCutType).c_str());
97     std::vector<double> pt3050 = lConfig.getParameter<std::vector<double> >(("Pt3050_"+lFullCutType).c_str());
98     if(i0 == 0) {
99     for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[0][i2] = pt010 [i2];
100     for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[1][i2] = pt1020[i2];
101     for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[2][i2] = pt2030[i2];
102     for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[3][i2] = pt3050[i2];
103     }
104     if(i0 == 1) {
105     for(int i2 = 0; i2 < 4; i2++) fRMSCut[0][i2] = pt010 [i2];
106     for(int i2 = 0; i2 < 4; i2++) fRMSCut[1][i2] = pt1020[i2];
107     for(int i2 = 0; i2 < 4; i2++) fRMSCut[2][i2] = pt2030[i2];
108     for(int i2 = 0; i2 < 4; i2++) fRMSCut[3][i2] = pt3050[i2];
109     }
110     }
111     return;
112     }
113    
114    
115 pharris 1.11 fLowPtReader = 0;
116     fLowPtReader = new TMVA::Reader( "!Color:!Silent:Error" );
117     fLowPtReader->AddVariable( "nvtx" , &fNVtx );
118 pharris 1.17 if(fType != k53) fLowPtReader->AddVariable( "jetPt" , &fJPt1 );
119     if(fType != k53) fLowPtReader->AddVariable( "jetEta" , &fJEta1 );
120     if(fType != k53) fLowPtReader->AddVariable( "jetPhi" , &fJPhi1 );
121 pharris 1.11 fLowPtReader->AddVariable( "dZ" , &fJDZ1 );
122 pharris 1.18 if(fType != k53 && fType != k53MET) fLowPtReader->AddVariable( "d0" , &fJD01 );
123 pharris 1.11 fLowPtReader->AddVariable( "beta" , &fBeta );
124     fLowPtReader->AddVariable( "betaStar" , &fBetaStar );
125     fLowPtReader->AddVariable( "nCharged" , &fNCharged );
126     fLowPtReader->AddVariable( "nNeutrals", &fNNeutrals );
127 pharris 1.17 if(fType != k53 && fType != k53MET) fLowPtReader->AddVariable( "dRMean" , &fDRMean );
128     if(fType == k53 || fType == k53MET) fLowPtReader->AddVariable( "dR2Mean" , &fDR2Mean );
129     if(fType == k53 || fType == k53MET) fLowPtReader->AddVariable( "ptD" , &fPtD );
130 pharris 1.11 fLowPtReader->AddVariable( "frac01" , &fFrac01 );
131     fLowPtReader->AddVariable( "frac02" , &fFrac02 );
132     fLowPtReader->AddVariable( "frac03" , &fFrac03 );
133     fLowPtReader->AddVariable( "frac04" , &fFrac04 );
134     fLowPtReader->AddVariable( "frac05" , &fFrac05 );
135 pharris 1.18 if(fType == k53) fLowPtReader->AddSpectator( "jetPt" , &fJPt1 );
136     if(fType == k53) fLowPtReader->AddSpectator( "jetEta" , &fJEta1 );
137     if(fType == k53) fLowPtReader->AddSpectator( "jetPhi" , &fJPhi1 );
138 pharris 1.11
139 pharris 1.1 fReader = 0;
140     fReader = new TMVA::Reader( "!Color:!Silent:Error" );
141     if (fType == kBaseline) {
142 pharris 1.4 fReader->AddVariable( "nvtx" , &fNVtx );
143     fReader->AddVariable( "jetPt" , &fJPt1 );
144     fReader->AddVariable( "jetEta" , &fJEta1 );
145     fReader->AddVariable( "jetPhi" , &fJPhi1 );
146     fReader->AddVariable( "dZ" , &fJDZ1 );
147     fReader->AddVariable( "d0" , &fJD01 );
148     fReader->AddVariable( "beta" , &fBeta );
149     fReader->AddVariable( "betaStar" , &fBetaStar );
150     fReader->AddVariable( "nCharged" , &fNCharged );
151     fReader->AddVariable( "nNeutrals", &fNNeutrals );
152     fReader->AddVariable( "dRMean" , &fDRMean );
153     fReader->AddVariable( "frac01" , &fFrac01 );
154     fReader->AddVariable( "frac02" , &fFrac02 );
155     fReader->AddVariable( "frac03" , &fFrac03 );
156     fReader->AddVariable( "frac04" , &fFrac04 );
157     fReader->AddVariable( "frac05" , &fFrac05 );
158 pharris 1.1 }
159 pharris 1.11 if (fType == k42) {
160     fReader->AddVariable( "frac01" , &fFrac01 );
161     fReader->AddVariable( "frac02" , &fFrac02 );
162     fReader->AddVariable( "frac03" , &fFrac03 );
163     fReader->AddVariable( "frac04" , &fFrac04 );
164     fReader->AddVariable( "frac05" , &fFrac05 );
165     fReader->AddVariable( "nvtx" , &fNVtx );
166     fReader->AddVariable( "nNeutrals", &fNNeutrals );
167     fReader->AddVariable( "beta" , &fBeta );
168     fReader->AddVariable( "betaStar" , &fBetaStar );
169     fReader->AddVariable( "dZ" , &fJDZ1 );
170     fReader->AddVariable( "nCharged" , &fNCharged );
171     fReader->AddSpectator( "jetPt" , &fJPt1 );
172     fReader->AddSpectator( "jetEta" , &fJEta1 );
173     }
174     if (fType == k52) {
175     fReader->AddVariable( "frac01" , &fFrac01 );
176     fReader->AddVariable( "frac02" , &fFrac02 );
177     fReader->AddVariable( "frac03" , &fFrac03 );
178     fReader->AddVariable( "frac04" , &fFrac04 );
179     fReader->AddVariable( "frac05" , &fFrac05 );
180     fReader->AddVariable( "dR2Mean" , &fDR2Mean );
181     fReader->AddVariable( "nvtx" , &fNVtx );
182     fReader->AddVariable( "nNeutrals", &fNNeutrals );
183     fReader->AddVariable( "beta" , &fBeta );
184     fReader->AddVariable( "betaStar" , &fBetaStar );
185     fReader->AddVariable( "dZ" , &fJDZ1 );
186     fReader->AddVariable( "nCharged" , &fNCharged );
187     fReader->AddSpectator( "jetPt" , &fJPt1 );
188     fReader->AddSpectator( "jetEta" , &fJEta1 );
189     }
190 pharris 1.17 if (fType == k53) {
191     fReader->AddVariable( "nvtx" , &fNVtx );
192     fReader->AddVariable( "dZ" , &fJDZ1 );
193     fReader->AddVariable( "beta" , &fBeta );
194     fReader->AddVariable( "betaStar" , &fBetaStar );
195     fReader->AddVariable( "nCharged" , &fNCharged );
196     fReader->AddVariable( "nNeutrals", &fNNeutrals );
197 pharris 1.19 fReader->AddVariable( "dR2Mean" , &fDR2Mean );
198 pharris 1.17 fReader->AddVariable( "ptD" , &fPtD );
199     fReader->AddVariable( "frac01" , &fFrac01 );
200     fReader->AddVariable( "frac02" , &fFrac02 );
201     fReader->AddVariable( "frac03" , &fFrac03 );
202     fReader->AddVariable( "frac04" , &fFrac04 );
203     fReader->AddVariable( "frac05" , &fFrac05 );
204     fReader->AddSpectator( "jetPt" , &fJPt1 );
205     fReader->AddSpectator( "jetEta" , &fJEta1 );
206     fReader->AddSpectator( "jetPhi" , &fJPhi1 );
207     }
208     if (fType == k53MET) {
209     fReader->AddVariable( "nvtx" , &fNVtx );
210     fReader->AddVariable( "jetPt" , &fJPt1 );
211     fReader->AddVariable( "jetEta" , &fJEta1 );
212     fReader->AddVariable( "jetPhi" , &fJPhi1 );
213     fReader->AddVariable( "dZ" , &fJDZ1 );
214     fReader->AddVariable( "beta" , &fBeta );
215     fReader->AddVariable( "betaStar" , &fBetaStar );
216     fReader->AddVariable( "nCharged" , &fNCharged );
217     fReader->AddVariable( "nNeutrals", &fNNeutrals );
218 pharris 1.19 fReader->AddVariable( "dR2Mean" , &fDR2Mean );
219 pharris 1.17 fReader->AddVariable( "ptD" , &fPtD );
220     fReader->AddVariable( "frac01" , &fFrac01 );
221     fReader->AddVariable( "frac02" , &fFrac02 );
222     fReader->AddVariable( "frac03" , &fFrac03 );
223     fReader->AddVariable( "frac04" , &fFrac04 );
224     fReader->AddVariable( "frac05" , &fFrac05 );
225     }
226 pharris 1.16 if (fType == kQGP) {
227     fReader->AddVariable( "nvtx" , &fNVtx );
228     fReader->AddVariable( "jetPt" , &fJPt1 );
229     fReader->AddVariable( "jetEta" , &fJEta1 );
230     fReader->AddVariable( "jetPhi" , &fJPhi1 );
231     fReader->AddVariable( "beta" , &fBeta );
232     fReader->AddVariable( "betaStar" , &fBetaStar );
233     fReader->AddVariable( "nParticles", &fNNeutrals );
234     fReader->AddVariable( "nCharged" , &fNCharged );
235     fReader->AddVariable( "dRMean" , &fDRMean );
236     fReader->AddVariable( "ptD" , &fPtD );
237     fReader->AddVariable( "frac01" , &fFrac01 );
238     fReader->AddVariable( "frac02" , &fFrac02 );
239     fReader->AddVariable( "frac03" , &fFrac03 );
240     fReader->AddVariable( "frac04" , &fFrac04 );
241     fReader->AddVariable( "frac05" , &fFrac05 );
242     }
243 pharris 1.11
244     fLowPtReader->BookMVA(fLowPtMethodName , iLowPtWeights );
245 pharris 1.4 fReader->BookMVA(fHighPtMethodName , iHighPtWeights );
246 pharris 1.1 std::cout << "Jet ID MVA Initialization\n";
247 pharris 1.4 std::cout << "MethodName : " << fLowPtMethodName << " , type == " << fType << std::endl;
248    
249 pharris 1.1 }
250    
251     //--------------------------------------------------------------------------------------------------
252     Double_t JetIDMVA::MVAValue(
253     Float_t iNPV ,
254     Float_t iJPt1 ,
255     Float_t iJEta1 ,
256     Float_t iJPhi1 ,
257     Float_t iJD01 ,
258     Float_t iJDZ1 ,
259 pharris 1.4 Float_t iBeta ,
260     Float_t iBetaStar,
261     Float_t iNCharged,
262     Float_t iNNeutrals,
263     Float_t iDRMean ,
264     Float_t iFrac01 ,
265     Float_t iFrac02 ,
266     Float_t iFrac03 ,
267     Float_t iFrac04 ,
268 pharris 1.11 Float_t iFrac05 ,
269 pharris 1.17 Float_t iDR2Mean ,
270     Float_t iPtD
271 pharris 1.1 ){
272    
273 pharris 1.4 if(!fIsInitialized) {
274 pharris 1.1 std::cout << "Error: JetIDMVA not properly initialized.\n";
275     return -9999;
276     }
277    
278 pharris 1.4 fNVtx = iNPV;
279     fJPt1 = iJPt1;
280     fJEta1 = iJEta1;
281     fJPhi1 = fJPhi1;
282     fJD01 = iJD01;
283     fJDZ1 = iJDZ1;
284     fBeta = iBeta;
285     fBetaStar = iBetaStar;
286     fNCharged = iNCharged;
287     fNNeutrals = iNNeutrals;
288     fDRMean = iDRMean;
289 pharris 1.17 fPtD = iPtD;
290 pharris 1.4 fFrac01 = iFrac01;
291     fFrac02 = iFrac02;
292     fFrac03 = iFrac03;
293     fFrac04 = iFrac04;
294     fFrac05 = iFrac05;
295 pharris 1.11 fDR2Mean = iDR2Mean;
296 pharris 1.1
297     Double_t lMVA = -9999;
298 pharris 1.11 if(iJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName );
299 pharris 1.4 if(iJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
300 pharris 1.1 return lMVA;
301     }
302     //--------------------------------------------------------------------------------------------------
303 pharris 1.15 Bool_t JetIDMVA::passPt(const PFJet *iJet, FactorizedJetCorrector *iJetCorrector,
304     const PileupEnergyDensityCol *iPileupEnergyDensity,
305     RhoUtilities::RhoType type) {
306     if(iJetCorrector == 0) return (iJet->Pt() > fJetPtMin);
307     double lPt = correctedPt(iJet, iJetCorrector,iPileupEnergyDensity,type);
308     if(lPt < fJetPtMin) return false;
309     return true;
310     }
311     //--------------------------------------------------------------------------------------------------
312 pharris 1.4 Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices,
313 pharris 1.1 FactorizedJetCorrector *iJetCorrector,
314 ceballos 1.10 const PileupEnergyDensityCol *iPileupEnergyDensity,
315     RhoUtilities::RhoType type) {
316 pharris 1.4
317 pharris 1.1 if(!JetTools::passPFLooseId(iJet)) return false;
318 pharris 1.7 if(fabs(iJet->Eta()) > 4.99) return false;
319 pharris 1.5
320 pharris 1.4 double lMVA = MVAValue (iJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity);
321 ceballos 1.10 double lPt = correctedPt(iJet, iJetCorrector,iPileupEnergyDensity,type);
322 pharris 1.8 if(lPt < fJetPtMin) return false;
323 pharris 1.7
324 pharris 1.4 int lPtId = 0;
325 pharris 1.7 if(lPt > 10 && lPt < 20) lPtId = 1;
326     if(lPt > 20 && lPt < 30) lPtId = 2;
327 pharris 1.4 if(lPt > 30 ) lPtId = 3;
328    
329     int lEtaId = 0;
330     if(fabs(iJet->Eta()) > 2.5 && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
331     if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
332     if(fabs(iJet->Eta()) > 3.0 && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
333    
334     double lMVACut = fMVACut[lPtId][lEtaId];
335     if(lMVA < lMVACut) return false;
336 pharris 1.1 return true;
337     }
338     //--------------------------------------------------------------------------------------------------
339 pharris 1.12 Bool_t JetIDMVA::passCut(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices) {
340     if(!JetTools::passPFLooseId(iJet)) return false;
341     if(iJet->Pt() < fJetPtMin) return false;
342     if(fabs(iJet->Eta()) > 4.99) return false;
343 bendavid 1.13 //if(fType == kCut) passCut(iJet,iVertex,iVertices);
344 pharris 1.12
345     double lPt = iJet->Pt();
346     int lPtId = 0;
347     if(lPt > 10 && lPt < 20) lPtId = 1;
348     if(lPt > 20 && lPt < 30) lPtId = 2;
349     if(lPt > 30 ) lPtId = 3;
350    
351     int lEtaId = 0;
352     if(fabs(iJet->Eta()) > 2.5 && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
353     if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
354     if(fabs(iJet->Eta()) > 3.0 && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
355     float betaStarModified = JetTools::betaStarClassic(iJet,iVertex,iVertices)/log(iVertices ->GetEntries()-0.64);
356     float dR2Mean = JetTools::dR2Mean(iJet,-1);
357    
358     if(betaStarModified < fBetaStarCut[lPtId][lEtaId] &&
359     dR2Mean < fRMSCut [lPtId][lEtaId]
360     ) return true;
361    
362     return false;
363     }
364     //--------------------------------------------------------------------------------------------------
365 pharris 1.4 Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices) {
366 pharris 1.2 if(!JetTools::passPFLooseId(iJet)) return false;
367 pharris 1.5 if(iJet->Pt() < fJetPtMin) return false;
368 pharris 1.7 if(fabs(iJet->Eta()) > 4.99) return false;
369 bendavid 1.13 if(fType == kCut) return passCut(iJet,iVertex,iVertices);
370 pharris 1.4 double lMVA = MVAValue(iJet,iVertex,iVertices);
371    
372     int lPtId = 0;
373     if(iJet->Pt() > 10 && iJet->Pt() < 20) lPtId = 1;
374     if(iJet->Pt() > 20 && iJet->Pt() < 30) lPtId = 2;
375     if(iJet->Pt() > 30 ) lPtId = 3;
376    
377     int lEtaId = 0;
378     if(fabs(iJet->Eta()) > 2.5 && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
379     if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
380     if(fabs(iJet->Eta()) > 3.0 && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
381    
382     double lMVACut = fMVACut[lPtId][lEtaId];
383     if(lMVA < lMVACut) return false;
384 pharris 1.2 return true;
385 pharris 1.4 //if(lMVA < -0.8) return false;
386     //if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
387     //if(iJet->Pt() < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false; //This line is a bug in the Met training
388     //if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
389 pharris 1.2 }
390     //--------------------------------------------------------------------------------------------------
391 pharris 1.4 Double_t JetIDMVA::MVAValue(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices, //Vertex here is the PV
392 pharris 1.1 FactorizedJetCorrector *iJetCorrector,
393     const PileupEnergyDensityCol *iPileupEnergyDensity,
394     Bool_t printDebug) {
395    
396     if (!fIsInitialized) {
397     std::cout << "Error: JetIDMVA not properly initialized.\n";
398     return -9999;
399     }
400     if(!JetTools::passPFLooseId(iJet)) return -2.;
401    
402     //set all input variables
403 pharris 1.4 fNVtx = iVertices->GetEntries();
404 pharris 1.1 fJPt1 = correctedPt(iJet,iJetCorrector,iPileupEnergyDensity);
405     fJEta1 = iJet->RawMom().Eta();
406     fJPhi1 = iJet->RawMom().Phi();
407 pharris 1.4 fJD01 = JetTools::impactParameter(iJet,iVertex);
408     fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
409     fBeta = JetTools::Beta(iJet,iVertex,fDZCut);
410     fBetaStar = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
411     fNCharged = iJet->ChargedMultiplicity();
412     fNNeutrals = iJet->NeutralMultiplicity();
413    
414     fDRMean = JetTools::dRMean(iJet,-1);
415 pharris 1.11 fDR2Mean = JetTools::dR2Mean(iJet,-1);
416 pharris 1.17 fPtD = JetTools::W(iJet,-1,0);
417 pharris 1.4 fFrac01 = JetTools::frac (iJet,0.1,0. ,-1);
418     fFrac02 = JetTools::frac (iJet,0.2,0.1,-1);
419     fFrac03 = JetTools::frac (iJet,0.3,0.2,-1);
420     fFrac04 = JetTools::frac (iJet,0.4,0.3,-1);
421     fFrac05 = JetTools::frac (iJet,0.5,0.4,-1);
422    
423     double lMVA = 0;
424 pharris 1.11 if(fJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName );
425 pharris 1.4 if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
426 pharris 1.1 if (printDebug == kTRUE) {
427     std::cout << "Debug Jet MVA: "
428 pharris 1.4 << fNVtx << " "
429     << fJPt1 << " "
430     << fJEta1 << " "
431     << fJPhi1 << " "
432     << fJD01 << " "
433     << fJDZ1 << " "
434     << fBeta << " "
435     << fBetaStar << " "
436     << fNCharged << " "
437     << fNNeutrals << " "
438     << fDRMean << " "
439 pharris 1.17 << fPtD << " "
440 pharris 1.4 << fFrac01 << " "
441     << fFrac02 << " "
442     << fFrac03 << " "
443     << fFrac04 << " "
444 pharris 1.11 << fFrac05 << " "
445 pharris 1.17 << fDR2Mean
446 pharris 1.1 << " === : === "
447     << lMVA << " "
448 pharris 1.19 << " -----> raw pt " << iJet->Pt() << std::endl;
449 pharris 1.1 }
450 pharris 1.19 //std::cout << " === " << iJet->Pt() << " -- " << iJet->Eta() << " -- " << fPtD << " -- " << lMVA << std::endl;
451 pharris 1.1 return lMVA;
452     }
453 pharris 1.4 Double_t JetIDMVA::MVAValue(const PFJet *iJet,const Vertex *iVertex, const VertexCol *iVertices,//Vertex here is the PV
454 pharris 1.2 Bool_t printDebug) {
455    
456     if (!fIsInitialized) {
457     std::cout << "Error: JetIDMVA not properly initialized.\n";
458     return -9999;
459     }
460     if(!JetTools::passPFLooseId(iJet)) return -2.;
461    
462     //set all input variables
463 pharris 1.4 fNVtx = iVertices->GetEntries();
464 pharris 1.2 fJPt1 = iJet->Pt();
465     fJEta1 = iJet->RawMom().Eta();
466     fJPhi1 = iJet->RawMom().Phi();
467 pharris 1.4 fJD01 = JetTools::impactParameter(iJet,iVertex);
468     fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
469     fBeta = JetTools::Beta(iJet,iVertex,fDZCut);
470     fBetaStar = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
471     fNCharged = iJet->ChargedMultiplicity();
472     fNNeutrals = iJet->NeutralMultiplicity();
473 pharris 1.17 fPtD = JetTools::W(iJet,-1,0);
474 pharris 1.11 fDRMean = JetTools::dRMean (iJet,-1);
475     fDR2Mean = JetTools::dR2Mean(iJet,-1);
476     fFrac01 = JetTools::frac (iJet,0.1,0. ,-1);
477     fFrac02 = JetTools::frac (iJet,0.2,0.1,-1);
478     fFrac03 = JetTools::frac (iJet,0.3,0.2,-1);
479     fFrac04 = JetTools::frac (iJet,0.4,0.3,-1);
480     fFrac05 = JetTools::frac (iJet,0.5,0.4,-1);
481 pharris 1.4
482     double lMVA = 0;
483 pharris 1.11 if(fJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName );
484     if(fJPt1 > 10) lMVA = fReader ->EvaluateMVA( fHighPtMethodName );
485 pharris 1.5
486 pharris 1.2 if (printDebug == kTRUE) {
487     std::cout << "Debug Jet MVA: "
488 pharris 1.4 << fNVtx << " "
489     << fJPt1 << " "
490     << fJEta1 << " "
491     << fJPhi1 << " "
492     << fJD01 << " "
493     << fJDZ1 << " "
494     << fBeta << " "
495     << fBetaStar << " "
496     << fNCharged << " "
497     << fNNeutrals << " "
498     << fDRMean << " "
499 pharris 1.17 << fPtD << " "
500 pharris 1.4 << fFrac01 << " "
501     << fFrac02 << " "
502     << fFrac03 << " "
503     << fFrac04 << " "
504 pharris 1.11 << fFrac05 << " "
505 pharris 1.17 << fDR2Mean
506 pharris 1.2 << " === : === "
507     << lMVA << " "
508     << std::endl;
509     }
510    
511     return lMVA;
512     }
513 pharris 1.16 //--------------------------------------------------------------------------------------------------
514     Double_t* JetIDMVA::QGValue(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices, //Vertex here is the PV
515     FactorizedJetCorrector *iJetCorrector,
516     const PileupEnergyDensityCol *iPileupEnergyDensity,
517     Bool_t printDebug) {
518    
519     Double_t *lId = new double[3];
520     lId[0] = -2;
521     lId[1] = -2;
522     lId[2] = -2;
523     if (!fIsInitialized) {
524     std::cout << "Error: JetIDMVA not properly initialized.\n";
525     return lId;
526     }
527     if(!JetTools::passPFLooseId(iJet)) return lId;
528    
529     fJPt1 = correctedPt(iJet,iJetCorrector,iPileupEnergyDensity);
530     if(fJPt1 < 20) return lId;
531    
532     //set all input variables
533     fNVtx = iVertices->GetEntries();
534     fJEta1 = iJet->RawMom().Eta();
535     fJPhi1 = iJet->RawMom().Phi();
536     fJD01 = JetTools::impactParameter(iJet,iVertex);
537     fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
538     fBeta = JetTools::Beta(iJet,iVertex,fDZCut);
539     fBetaStar = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
540     fNCharged = iJet->ChargedMultiplicity();
541     fNNeutrals = iJet->NeutralMultiplicity();
542     fNParticles = fNCharged+fNNeutrals;
543     fPtD = JetTools::W(iJet,-1,0);
544    
545     fDRMean = JetTools::dRMean(iJet,-1);
546     fDR2Mean = JetTools::dR2Mean(iJet,-1);
547     fFrac01 = JetTools::frac (iJet,0.1,0. ,-1);
548     fFrac02 = JetTools::frac (iJet,0.2,0.1,-1);
549     fFrac03 = JetTools::frac (iJet,0.3,0.2,-1);
550     fFrac04 = JetTools::frac (iJet,0.4,0.3,-1);
551     fFrac05 = JetTools::frac (iJet,0.5,0.4,-1);
552    
553     double lMVA = 0;
554     lId[0] = fReader->EvaluateMulticlass( fHighPtMethodName )[0];
555     lId[1] = fReader->EvaluateMulticlass( fHighPtMethodName )[1];
556     lId[2] = fReader->EvaluateMulticlass( fHighPtMethodName )[2];
557     if (printDebug == kTRUE) {
558     std::cout << "Debug Jet MVA: "
559     << fNVtx << " "
560     << fJPt1 << " "
561     << fJEta1 << " "
562     << fJPhi1 << " "
563     << fJD01 << " "
564     << fJDZ1 << " "
565     << fBeta << " "
566     << fBetaStar << " "
567     << fNCharged << " "
568     << fNNeutrals << " "
569     << fDRMean << " "
570     << fFrac01 << " "
571     << fFrac02 << " "
572     << fFrac03 << " "
573     << fFrac04 << " "
574     << fFrac05 << " "
575     << fDRMean
576     << " === : === "
577     << lMVA << " "
578     << std::endl;
579     }
580     return lId;
581     }
582 ceballos 1.10 Double_t JetIDMVA::correctedPt(const PFJet *iJet, FactorizedJetCorrector *iJetCorrector,
583     const PileupEnergyDensityCol *iPUEnergyDensity,
584 pharris 1.16 RhoUtilities::RhoType type,int iId) {
585 ceballos 1.10 Double_t Rho = 0.0;
586     switch(type) {
587     case RhoUtilities::MIT_RHO_VORONOI_HIGH_ETA:
588     Rho = iPUEnergyDensity->At(0)->Rho();
589     break;
590     case RhoUtilities::MIT_RHO_VORONOI_LOW_ETA:
591     Rho = iPUEnergyDensity->At(0)->RhoLowEta();
592     break;
593     case RhoUtilities::MIT_RHO_RANDOM_HIGH_ETA:
594     Rho = iPUEnergyDensity->At(0)->RhoRandom();
595     break;
596     case RhoUtilities::MIT_RHO_RANDOM_LOW_ETA:
597     Rho = iPUEnergyDensity->At(0)->RhoRandomLowEta();
598     break;
599     case RhoUtilities::CMS_RHO_RHOKT6PFJETS:
600     Rho = iPUEnergyDensity->At(0)->RhoKt6PFJets();
601     break;
602     default:
603     // use the old default
604 pharris 1.21 f42 ? Rho = iPUEnergyDensity->At(0)->RhoRandom() : Rho = iPUEnergyDensity->At(0)->Rho();
605 ceballos 1.10 break;
606     }
607    
608 pharris 1.1 const FourVectorM rawMom = iJet->RawMom();
609     iJetCorrector->setJetEta(rawMom.Eta());
610     iJetCorrector->setJetPt (rawMom.Pt());
611     iJetCorrector->setJetPhi(rawMom.Phi());
612     iJetCorrector->setJetE (rawMom.E());
613 ceballos 1.10 iJetCorrector->setRho (Rho);
614 pharris 1.1 iJetCorrector->setJetA (iJet->JetArea());
615     iJetCorrector->setJetEMF(-99.0);
616 pharris 1.16 Double_t correction = 1.;
617     if(iId < 0 || iId == 100) correction = iJetCorrector->getCorrection();
618     std::vector<float> lCorrections; if(iId != -1 && iId != 100) lCorrections = iJetCorrector->getSubCorrections();
619     if(iId > -1 && iId < int(lCorrections.size())) correction = lCorrections[iId];
620     if(iId == 100) {
621     iJetCorrector->setJetEta(rawMom.Eta());
622     iJetCorrector->setJetPt (rawMom.Pt());
623     iJetCorrector->setJetPhi(rawMom.Phi());
624     iJetCorrector->setJetE (rawMom.E());
625     iJetCorrector->setRho (Rho);
626     iJetCorrector->setJetA (iJet->JetArea());
627     iJetCorrector->setJetEMF(-99.0);
628     lCorrections = iJetCorrector->getSubCorrections();
629     double correction2 = 1;
630     correction2 *= lCorrections[0];
631     correction = correction-correction2;
632     }
633    
634 pharris 1.1 return rawMom.Pt()*correction;
635     }