ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/JetIDMVA.cc
Revision: 1.19
Committed: Sun Mar 31 20:27:52 2013 UTC (12 years, 1 month ago) by pharris
Content type: text/plain
Branch: MAIN
Changes since 1.18: +5 -5 lines
Log Message:
updated jet id mva

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     TString iCutFileName)
57 pharris 1.1 {
58    
59     fIsInitialized = kTRUE;
60     fType = iType;
61 pharris 1.4 fCutType = iCutType;
62 bendavid 1.14
63     std::string lCutId = "JetIdParams";
64 pharris 1.17 if(fType == k42) lCutId = "PuJetIdOptMVA_wp";
65     if(fType == k52) lCutId = "full_5x_wp";
66     if(fType == kCut) lCutId = "PuJetIdCutBased_wp";
67     if(fType == k53) lCutId = "full_53x_wp";
68     if(fType == k53MET) lCutId = "met_53x_wp";
69    
70     //Load Cut Matrix
71 bendavid 1.14 edm::ParameterSet lDefConfig = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>("JetIdParams");
72     edm::ParameterSet lConfig = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>(lCutId);
73     std::string lCutType = "Tight";
74     if(fCutType == kMedium) lCutType = "Medium";
75     if(fCutType == kLoose ) lCutType = "Loose";
76     if(fCutType == kMET ) lCutType = "MET";
77     if(fType != kCut) {
78     std::string lLowPtCut = "MET";
79     std::vector<double> lPt010 = lDefConfig.getParameter<std::vector<double> >(("Pt010_" +lLowPtCut).c_str());
80     std::vector<double> lPt1020 = lConfig.getParameter<std::vector<double> >(("Pt1020_"+lCutType).c_str());
81     std::vector<double> lPt2030 = lConfig.getParameter<std::vector<double> >(("Pt2030_"+lCutType).c_str());
82     std::vector<double> lPt3050 = lConfig.getParameter<std::vector<double> >(("Pt3050_"+lCutType).c_str());
83     for(int i0 = 0; i0 < 4; i0++) fMVACut[0][i0] = lPt010 [i0];
84     for(int i0 = 0; i0 < 4; i0++) fMVACut[1][i0] = lPt1020[i0];
85     for(int i0 = 0; i0 < 4; i0++) fMVACut[2][i0] = lPt2030[i0];
86     for(int i0 = 0; i0 < 4; i0++) fMVACut[3][i0] = lPt3050[i0];
87     }
88     if(fType == kCut) {
89     for(int i0 = 0; i0 < 2; i0++) {
90     std::string lFullCutType = lCutType;
91     if(i0 == 0) lFullCutType = "BetaStar"+ lCutType;
92     if(i0 == 1) lFullCutType = "RMS" + lCutType;
93     std::vector<double> pt010 = lConfig.getParameter<std::vector<double> >(("Pt010_" +lFullCutType).c_str());
94     std::vector<double> pt1020 = lConfig.getParameter<std::vector<double> >(("Pt1020_"+lFullCutType).c_str());
95     std::vector<double> pt2030 = lConfig.getParameter<std::vector<double> >(("Pt2030_"+lFullCutType).c_str());
96     std::vector<double> pt3050 = lConfig.getParameter<std::vector<double> >(("Pt3050_"+lFullCutType).c_str());
97     if(i0 == 0) {
98     for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[0][i2] = pt010 [i2];
99     for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[1][i2] = pt1020[i2];
100     for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[2][i2] = pt2030[i2];
101     for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[3][i2] = pt3050[i2];
102     }
103     if(i0 == 1) {
104     for(int i2 = 0; i2 < 4; i2++) fRMSCut[0][i2] = pt010 [i2];
105     for(int i2 = 0; i2 < 4; i2++) fRMSCut[1][i2] = pt1020[i2];
106     for(int i2 = 0; i2 < 4; i2++) fRMSCut[2][i2] = pt2030[i2];
107     for(int i2 = 0; i2 < 4; i2++) fRMSCut[3][i2] = pt3050[i2];
108     }
109     }
110     return;
111     }
112    
113    
114 pharris 1.11 fLowPtReader = 0;
115     fLowPtReader = new TMVA::Reader( "!Color:!Silent:Error" );
116     fLowPtReader->AddVariable( "nvtx" , &fNVtx );
117 pharris 1.17 if(fType != k53) fLowPtReader->AddVariable( "jetPt" , &fJPt1 );
118     if(fType != k53) fLowPtReader->AddVariable( "jetEta" , &fJEta1 );
119     if(fType != k53) fLowPtReader->AddVariable( "jetPhi" , &fJPhi1 );
120 pharris 1.11 fLowPtReader->AddVariable( "dZ" , &fJDZ1 );
121 pharris 1.18 if(fType != k53 && fType != k53MET) fLowPtReader->AddVariable( "d0" , &fJD01 );
122 pharris 1.11 fLowPtReader->AddVariable( "beta" , &fBeta );
123     fLowPtReader->AddVariable( "betaStar" , &fBetaStar );
124     fLowPtReader->AddVariable( "nCharged" , &fNCharged );
125     fLowPtReader->AddVariable( "nNeutrals", &fNNeutrals );
126 pharris 1.17 if(fType != k53 && fType != k53MET) fLowPtReader->AddVariable( "dRMean" , &fDRMean );
127     if(fType == k53 || fType == k53MET) fLowPtReader->AddVariable( "dR2Mean" , &fDR2Mean );
128     if(fType == k53 || fType == k53MET) fLowPtReader->AddVariable( "ptD" , &fPtD );
129 pharris 1.11 fLowPtReader->AddVariable( "frac01" , &fFrac01 );
130     fLowPtReader->AddVariable( "frac02" , &fFrac02 );
131     fLowPtReader->AddVariable( "frac03" , &fFrac03 );
132     fLowPtReader->AddVariable( "frac04" , &fFrac04 );
133     fLowPtReader->AddVariable( "frac05" , &fFrac05 );
134 pharris 1.18 if(fType == k53) fLowPtReader->AddSpectator( "jetPt" , &fJPt1 );
135     if(fType == k53) fLowPtReader->AddSpectator( "jetEta" , &fJEta1 );
136     if(fType == k53) fLowPtReader->AddSpectator( "jetPhi" , &fJPhi1 );
137 pharris 1.11
138 pharris 1.1 fReader = 0;
139     fReader = new TMVA::Reader( "!Color:!Silent:Error" );
140     if (fType == kBaseline) {
141 pharris 1.4 fReader->AddVariable( "nvtx" , &fNVtx );
142     fReader->AddVariable( "jetPt" , &fJPt1 );
143     fReader->AddVariable( "jetEta" , &fJEta1 );
144     fReader->AddVariable( "jetPhi" , &fJPhi1 );
145     fReader->AddVariable( "dZ" , &fJDZ1 );
146     fReader->AddVariable( "d0" , &fJD01 );
147     fReader->AddVariable( "beta" , &fBeta );
148     fReader->AddVariable( "betaStar" , &fBetaStar );
149     fReader->AddVariable( "nCharged" , &fNCharged );
150     fReader->AddVariable( "nNeutrals", &fNNeutrals );
151     fReader->AddVariable( "dRMean" , &fDRMean );
152     fReader->AddVariable( "frac01" , &fFrac01 );
153     fReader->AddVariable( "frac02" , &fFrac02 );
154     fReader->AddVariable( "frac03" , &fFrac03 );
155     fReader->AddVariable( "frac04" , &fFrac04 );
156     fReader->AddVariable( "frac05" , &fFrac05 );
157 pharris 1.1 }
158 pharris 1.11 if (fType == k42) {
159     fReader->AddVariable( "frac01" , &fFrac01 );
160     fReader->AddVariable( "frac02" , &fFrac02 );
161     fReader->AddVariable( "frac03" , &fFrac03 );
162     fReader->AddVariable( "frac04" , &fFrac04 );
163     fReader->AddVariable( "frac05" , &fFrac05 );
164     fReader->AddVariable( "nvtx" , &fNVtx );
165     fReader->AddVariable( "nNeutrals", &fNNeutrals );
166     fReader->AddVariable( "beta" , &fBeta );
167     fReader->AddVariable( "betaStar" , &fBetaStar );
168     fReader->AddVariable( "dZ" , &fJDZ1 );
169     fReader->AddVariable( "nCharged" , &fNCharged );
170     fReader->AddSpectator( "jetPt" , &fJPt1 );
171     fReader->AddSpectator( "jetEta" , &fJEta1 );
172     }
173     if (fType == k52) {
174     fReader->AddVariable( "frac01" , &fFrac01 );
175     fReader->AddVariable( "frac02" , &fFrac02 );
176     fReader->AddVariable( "frac03" , &fFrac03 );
177     fReader->AddVariable( "frac04" , &fFrac04 );
178     fReader->AddVariable( "frac05" , &fFrac05 );
179     fReader->AddVariable( "dR2Mean" , &fDR2Mean );
180     fReader->AddVariable( "nvtx" , &fNVtx );
181     fReader->AddVariable( "nNeutrals", &fNNeutrals );
182     fReader->AddVariable( "beta" , &fBeta );
183     fReader->AddVariable( "betaStar" , &fBetaStar );
184     fReader->AddVariable( "dZ" , &fJDZ1 );
185     fReader->AddVariable( "nCharged" , &fNCharged );
186     fReader->AddSpectator( "jetPt" , &fJPt1 );
187     fReader->AddSpectator( "jetEta" , &fJEta1 );
188     }
189 pharris 1.17 if (fType == k53) {
190     fReader->AddVariable( "nvtx" , &fNVtx );
191     fReader->AddVariable( "dZ" , &fJDZ1 );
192     fReader->AddVariable( "beta" , &fBeta );
193     fReader->AddVariable( "betaStar" , &fBetaStar );
194     fReader->AddVariable( "nCharged" , &fNCharged );
195     fReader->AddVariable( "nNeutrals", &fNNeutrals );
196 pharris 1.19 fReader->AddVariable( "dR2Mean" , &fDR2Mean );
197 pharris 1.17 fReader->AddVariable( "ptD" , &fPtD );
198     fReader->AddVariable( "frac01" , &fFrac01 );
199     fReader->AddVariable( "frac02" , &fFrac02 );
200     fReader->AddVariable( "frac03" , &fFrac03 );
201     fReader->AddVariable( "frac04" , &fFrac04 );
202     fReader->AddVariable( "frac05" , &fFrac05 );
203     fReader->AddSpectator( "jetPt" , &fJPt1 );
204     fReader->AddSpectator( "jetEta" , &fJEta1 );
205     fReader->AddSpectator( "jetPhi" , &fJPhi1 );
206     }
207     if (fType == k53MET) {
208     fReader->AddVariable( "nvtx" , &fNVtx );
209     fReader->AddVariable( "jetPt" , &fJPt1 );
210     fReader->AddVariable( "jetEta" , &fJEta1 );
211     fReader->AddVariable( "jetPhi" , &fJPhi1 );
212     fReader->AddVariable( "dZ" , &fJDZ1 );
213     fReader->AddVariable( "beta" , &fBeta );
214     fReader->AddVariable( "betaStar" , &fBetaStar );
215     fReader->AddVariable( "nCharged" , &fNCharged );
216     fReader->AddVariable( "nNeutrals", &fNNeutrals );
217 pharris 1.19 fReader->AddVariable( "dR2Mean" , &fDR2Mean );
218 pharris 1.17 fReader->AddVariable( "ptD" , &fPtD );
219     fReader->AddVariable( "frac01" , &fFrac01 );
220     fReader->AddVariable( "frac02" , &fFrac02 );
221     fReader->AddVariable( "frac03" , &fFrac03 );
222     fReader->AddVariable( "frac04" , &fFrac04 );
223     fReader->AddVariable( "frac05" , &fFrac05 );
224     }
225 pharris 1.16 if (fType == kQGP) {
226     fReader->AddVariable( "nvtx" , &fNVtx );
227     fReader->AddVariable( "jetPt" , &fJPt1 );
228     fReader->AddVariable( "jetEta" , &fJEta1 );
229     fReader->AddVariable( "jetPhi" , &fJPhi1 );
230     fReader->AddVariable( "beta" , &fBeta );
231     fReader->AddVariable( "betaStar" , &fBetaStar );
232     fReader->AddVariable( "nParticles", &fNNeutrals );
233     fReader->AddVariable( "nCharged" , &fNCharged );
234     fReader->AddVariable( "dRMean" , &fDRMean );
235     fReader->AddVariable( "ptD" , &fPtD );
236     fReader->AddVariable( "frac01" , &fFrac01 );
237     fReader->AddVariable( "frac02" , &fFrac02 );
238     fReader->AddVariable( "frac03" , &fFrac03 );
239     fReader->AddVariable( "frac04" , &fFrac04 );
240     fReader->AddVariable( "frac05" , &fFrac05 );
241     }
242 pharris 1.11
243     fLowPtReader->BookMVA(fLowPtMethodName , iLowPtWeights );
244 pharris 1.4 fReader->BookMVA(fHighPtMethodName , iHighPtWeights );
245 pharris 1.1 std::cout << "Jet ID MVA Initialization\n";
246 pharris 1.4 std::cout << "MethodName : " << fLowPtMethodName << " , type == " << fType << std::endl;
247    
248 pharris 1.1 }
249    
250     //--------------------------------------------------------------------------------------------------
251     Double_t JetIDMVA::MVAValue(
252     Float_t iNPV ,
253     Float_t iJPt1 ,
254     Float_t iJEta1 ,
255     Float_t iJPhi1 ,
256     Float_t iJD01 ,
257     Float_t iJDZ1 ,
258 pharris 1.4 Float_t iBeta ,
259     Float_t iBetaStar,
260     Float_t iNCharged,
261     Float_t iNNeutrals,
262     Float_t iDRMean ,
263     Float_t iFrac01 ,
264     Float_t iFrac02 ,
265     Float_t iFrac03 ,
266     Float_t iFrac04 ,
267 pharris 1.11 Float_t iFrac05 ,
268 pharris 1.17 Float_t iDR2Mean ,
269     Float_t iPtD
270 pharris 1.1 ){
271    
272 pharris 1.4 if(!fIsInitialized) {
273 pharris 1.1 std::cout << "Error: JetIDMVA not properly initialized.\n";
274     return -9999;
275     }
276    
277 pharris 1.4 fNVtx = iNPV;
278     fJPt1 = iJPt1;
279     fJEta1 = iJEta1;
280     fJPhi1 = fJPhi1;
281     fJD01 = iJD01;
282     fJDZ1 = iJDZ1;
283     fBeta = iBeta;
284     fBetaStar = iBetaStar;
285     fNCharged = iNCharged;
286     fNNeutrals = iNNeutrals;
287     fDRMean = iDRMean;
288 pharris 1.17 fPtD = iPtD;
289 pharris 1.4 fFrac01 = iFrac01;
290     fFrac02 = iFrac02;
291     fFrac03 = iFrac03;
292     fFrac04 = iFrac04;
293     fFrac05 = iFrac05;
294 pharris 1.11 fDR2Mean = iDR2Mean;
295 pharris 1.1
296     Double_t lMVA = -9999;
297 pharris 1.11 if(iJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName );
298 pharris 1.4 if(iJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
299 pharris 1.1 return lMVA;
300     }
301     //--------------------------------------------------------------------------------------------------
302 pharris 1.15 Bool_t JetIDMVA::passPt(const PFJet *iJet, FactorizedJetCorrector *iJetCorrector,
303     const PileupEnergyDensityCol *iPileupEnergyDensity,
304     RhoUtilities::RhoType type) {
305     if(iJetCorrector == 0) return (iJet->Pt() > fJetPtMin);
306     double lPt = correctedPt(iJet, iJetCorrector,iPileupEnergyDensity,type);
307     if(lPt < fJetPtMin) return false;
308     return true;
309     }
310     //--------------------------------------------------------------------------------------------------
311 pharris 1.4 Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices,
312 pharris 1.1 FactorizedJetCorrector *iJetCorrector,
313 ceballos 1.10 const PileupEnergyDensityCol *iPileupEnergyDensity,
314     RhoUtilities::RhoType type) {
315 pharris 1.4
316 pharris 1.1 if(!JetTools::passPFLooseId(iJet)) return false;
317 pharris 1.7 if(fabs(iJet->Eta()) > 4.99) return false;
318 pharris 1.5
319 pharris 1.4 double lMVA = MVAValue (iJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity);
320 ceballos 1.10 double lPt = correctedPt(iJet, iJetCorrector,iPileupEnergyDensity,type);
321 pharris 1.8 if(lPt < fJetPtMin) return false;
322 pharris 1.7
323 pharris 1.4 int lPtId = 0;
324 pharris 1.7 if(lPt > 10 && lPt < 20) lPtId = 1;
325     if(lPt > 20 && lPt < 30) lPtId = 2;
326 pharris 1.4 if(lPt > 30 ) lPtId = 3;
327    
328     int lEtaId = 0;
329     if(fabs(iJet->Eta()) > 2.5 && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
330     if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
331     if(fabs(iJet->Eta()) > 3.0 && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
332    
333     double lMVACut = fMVACut[lPtId][lEtaId];
334     if(lMVA < lMVACut) return false;
335 pharris 1.1 return true;
336     }
337     //--------------------------------------------------------------------------------------------------
338 pharris 1.12 Bool_t JetIDMVA::passCut(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices) {
339     if(!JetTools::passPFLooseId(iJet)) return false;
340     if(iJet->Pt() < fJetPtMin) return false;
341     if(fabs(iJet->Eta()) > 4.99) return false;
342 bendavid 1.13 //if(fType == kCut) passCut(iJet,iVertex,iVertices);
343 pharris 1.12
344     double lPt = iJet->Pt();
345     int lPtId = 0;
346     if(lPt > 10 && lPt < 20) lPtId = 1;
347     if(lPt > 20 && lPt < 30) lPtId = 2;
348     if(lPt > 30 ) lPtId = 3;
349    
350     int lEtaId = 0;
351     if(fabs(iJet->Eta()) > 2.5 && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
352     if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
353     if(fabs(iJet->Eta()) > 3.0 && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
354     float betaStarModified = JetTools::betaStarClassic(iJet,iVertex,iVertices)/log(iVertices ->GetEntries()-0.64);
355     float dR2Mean = JetTools::dR2Mean(iJet,-1);
356    
357     if(betaStarModified < fBetaStarCut[lPtId][lEtaId] &&
358     dR2Mean < fRMSCut [lPtId][lEtaId]
359     ) return true;
360    
361     return false;
362     }
363     //--------------------------------------------------------------------------------------------------
364 pharris 1.4 Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices) {
365 pharris 1.2 if(!JetTools::passPFLooseId(iJet)) return false;
366 pharris 1.5 if(iJet->Pt() < fJetPtMin) return false;
367 pharris 1.7 if(fabs(iJet->Eta()) > 4.99) return false;
368 bendavid 1.13 if(fType == kCut) return passCut(iJet,iVertex,iVertices);
369 pharris 1.4 double lMVA = MVAValue(iJet,iVertex,iVertices);
370    
371     int lPtId = 0;
372     if(iJet->Pt() > 10 && iJet->Pt() < 20) lPtId = 1;
373     if(iJet->Pt() > 20 && iJet->Pt() < 30) lPtId = 2;
374     if(iJet->Pt() > 30 ) lPtId = 3;
375    
376     int lEtaId = 0;
377     if(fabs(iJet->Eta()) > 2.5 && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
378     if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
379     if(fabs(iJet->Eta()) > 3.0 && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
380    
381     double lMVACut = fMVACut[lPtId][lEtaId];
382     if(lMVA < lMVACut) return false;
383 pharris 1.2 return true;
384 pharris 1.4 //if(lMVA < -0.8) return false;
385     //if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
386     //if(iJet->Pt() < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false; //This line is a bug in the Met training
387     //if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
388 pharris 1.2 }
389     //--------------------------------------------------------------------------------------------------
390 pharris 1.4 Double_t JetIDMVA::MVAValue(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices, //Vertex here is the PV
391 pharris 1.1 FactorizedJetCorrector *iJetCorrector,
392     const PileupEnergyDensityCol *iPileupEnergyDensity,
393     Bool_t printDebug) {
394    
395     if (!fIsInitialized) {
396     std::cout << "Error: JetIDMVA not properly initialized.\n";
397     return -9999;
398     }
399     if(!JetTools::passPFLooseId(iJet)) return -2.;
400    
401     //set all input variables
402 pharris 1.4 fNVtx = iVertices->GetEntries();
403 pharris 1.1 fJPt1 = correctedPt(iJet,iJetCorrector,iPileupEnergyDensity);
404     fJEta1 = iJet->RawMom().Eta();
405     fJPhi1 = iJet->RawMom().Phi();
406 pharris 1.4 fJD01 = JetTools::impactParameter(iJet,iVertex);
407     fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
408     fBeta = JetTools::Beta(iJet,iVertex,fDZCut);
409     fBetaStar = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
410     fNCharged = iJet->ChargedMultiplicity();
411     fNNeutrals = iJet->NeutralMultiplicity();
412    
413     fDRMean = JetTools::dRMean(iJet,-1);
414 pharris 1.11 fDR2Mean = JetTools::dR2Mean(iJet,-1);
415 pharris 1.17 fPtD = JetTools::W(iJet,-1,0);
416 pharris 1.4 fFrac01 = JetTools::frac (iJet,0.1,0. ,-1);
417     fFrac02 = JetTools::frac (iJet,0.2,0.1,-1);
418     fFrac03 = JetTools::frac (iJet,0.3,0.2,-1);
419     fFrac04 = JetTools::frac (iJet,0.4,0.3,-1);
420     fFrac05 = JetTools::frac (iJet,0.5,0.4,-1);
421    
422     double lMVA = 0;
423 pharris 1.11 if(fJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName );
424 pharris 1.4 if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
425 pharris 1.19 printDebug = kTRUE;
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     Rho = iPUEnergyDensity->At(0)->Rho();
605     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     }