ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/JetIDMVA.cc
Revision: 1.20
Committed: Fri Apr 5 13:14:28 2013 UTC (12 years, 1 month ago) by arapyan
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_029c, Mit_029b, Mit_029a
Changes since 1.19: +0 -1 lines
Log Message:
turn of the debug mode

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.1 if (printDebug == kTRUE) {
426     std::cout << "Debug Jet MVA: "
427 pharris 1.4 << fNVtx << " "
428     << fJPt1 << " "
429     << fJEta1 << " "
430     << fJPhi1 << " "
431     << fJD01 << " "
432     << fJDZ1 << " "
433     << fBeta << " "
434     << fBetaStar << " "
435     << fNCharged << " "
436     << fNNeutrals << " "
437     << fDRMean << " "
438 pharris 1.17 << fPtD << " "
439 pharris 1.4 << fFrac01 << " "
440     << fFrac02 << " "
441     << fFrac03 << " "
442     << fFrac04 << " "
443 pharris 1.11 << fFrac05 << " "
444 pharris 1.17 << fDR2Mean
445 pharris 1.1 << " === : === "
446     << lMVA << " "
447 pharris 1.19 << " -----> raw pt " << iJet->Pt() << std::endl;
448 pharris 1.1 }
449 pharris 1.19 //std::cout << " === " << iJet->Pt() << " -- " << iJet->Eta() << " -- " << fPtD << " -- " << lMVA << std::endl;
450 pharris 1.1 return lMVA;
451     }
452 pharris 1.4 Double_t JetIDMVA::MVAValue(const PFJet *iJet,const Vertex *iVertex, const VertexCol *iVertices,//Vertex here is the PV
453 pharris 1.2 Bool_t printDebug) {
454    
455     if (!fIsInitialized) {
456     std::cout << "Error: JetIDMVA not properly initialized.\n";
457     return -9999;
458     }
459     if(!JetTools::passPFLooseId(iJet)) return -2.;
460    
461     //set all input variables
462 pharris 1.4 fNVtx = iVertices->GetEntries();
463 pharris 1.2 fJPt1 = iJet->Pt();
464     fJEta1 = iJet->RawMom().Eta();
465     fJPhi1 = iJet->RawMom().Phi();
466 pharris 1.4 fJD01 = JetTools::impactParameter(iJet,iVertex);
467     fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
468     fBeta = JetTools::Beta(iJet,iVertex,fDZCut);
469     fBetaStar = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
470     fNCharged = iJet->ChargedMultiplicity();
471     fNNeutrals = iJet->NeutralMultiplicity();
472 pharris 1.17 fPtD = JetTools::W(iJet,-1,0);
473 pharris 1.11 fDRMean = JetTools::dRMean (iJet,-1);
474     fDR2Mean = JetTools::dR2Mean(iJet,-1);
475     fFrac01 = JetTools::frac (iJet,0.1,0. ,-1);
476     fFrac02 = JetTools::frac (iJet,0.2,0.1,-1);
477     fFrac03 = JetTools::frac (iJet,0.3,0.2,-1);
478     fFrac04 = JetTools::frac (iJet,0.4,0.3,-1);
479     fFrac05 = JetTools::frac (iJet,0.5,0.4,-1);
480 pharris 1.4
481     double lMVA = 0;
482 pharris 1.11 if(fJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName );
483     if(fJPt1 > 10) lMVA = fReader ->EvaluateMVA( fHighPtMethodName );
484 pharris 1.5
485 pharris 1.2 if (printDebug == kTRUE) {
486     std::cout << "Debug Jet MVA: "
487 pharris 1.4 << fNVtx << " "
488     << fJPt1 << " "
489     << fJEta1 << " "
490     << fJPhi1 << " "
491     << fJD01 << " "
492     << fJDZ1 << " "
493     << fBeta << " "
494     << fBetaStar << " "
495     << fNCharged << " "
496     << fNNeutrals << " "
497     << fDRMean << " "
498 pharris 1.17 << fPtD << " "
499 pharris 1.4 << fFrac01 << " "
500     << fFrac02 << " "
501     << fFrac03 << " "
502     << fFrac04 << " "
503 pharris 1.11 << fFrac05 << " "
504 pharris 1.17 << fDR2Mean
505 pharris 1.2 << " === : === "
506     << lMVA << " "
507     << std::endl;
508     }
509    
510     return lMVA;
511     }
512 pharris 1.16 //--------------------------------------------------------------------------------------------------
513     Double_t* JetIDMVA::QGValue(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices, //Vertex here is the PV
514     FactorizedJetCorrector *iJetCorrector,
515     const PileupEnergyDensityCol *iPileupEnergyDensity,
516     Bool_t printDebug) {
517    
518     Double_t *lId = new double[3];
519     lId[0] = -2;
520     lId[1] = -2;
521     lId[2] = -2;
522     if (!fIsInitialized) {
523     std::cout << "Error: JetIDMVA not properly initialized.\n";
524     return lId;
525     }
526     if(!JetTools::passPFLooseId(iJet)) return lId;
527    
528     fJPt1 = correctedPt(iJet,iJetCorrector,iPileupEnergyDensity);
529     if(fJPt1 < 20) return lId;
530    
531     //set all input variables
532     fNVtx = iVertices->GetEntries();
533     fJEta1 = iJet->RawMom().Eta();
534     fJPhi1 = iJet->RawMom().Phi();
535     fJD01 = JetTools::impactParameter(iJet,iVertex);
536     fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
537     fBeta = JetTools::Beta(iJet,iVertex,fDZCut);
538     fBetaStar = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
539     fNCharged = iJet->ChargedMultiplicity();
540     fNNeutrals = iJet->NeutralMultiplicity();
541     fNParticles = fNCharged+fNNeutrals;
542     fPtD = JetTools::W(iJet,-1,0);
543    
544     fDRMean = JetTools::dRMean(iJet,-1);
545     fDR2Mean = JetTools::dR2Mean(iJet,-1);
546     fFrac01 = JetTools::frac (iJet,0.1,0. ,-1);
547     fFrac02 = JetTools::frac (iJet,0.2,0.1,-1);
548     fFrac03 = JetTools::frac (iJet,0.3,0.2,-1);
549     fFrac04 = JetTools::frac (iJet,0.4,0.3,-1);
550     fFrac05 = JetTools::frac (iJet,0.5,0.4,-1);
551    
552     double lMVA = 0;
553     lId[0] = fReader->EvaluateMulticlass( fHighPtMethodName )[0];
554     lId[1] = fReader->EvaluateMulticlass( fHighPtMethodName )[1];
555     lId[2] = fReader->EvaluateMulticlass( fHighPtMethodName )[2];
556     if (printDebug == kTRUE) {
557     std::cout << "Debug Jet MVA: "
558     << fNVtx << " "
559     << fJPt1 << " "
560     << fJEta1 << " "
561     << fJPhi1 << " "
562     << fJD01 << " "
563     << fJDZ1 << " "
564     << fBeta << " "
565     << fBetaStar << " "
566     << fNCharged << " "
567     << fNNeutrals << " "
568     << fDRMean << " "
569     << fFrac01 << " "
570     << fFrac02 << " "
571     << fFrac03 << " "
572     << fFrac04 << " "
573     << fFrac05 << " "
574     << fDRMean
575     << " === : === "
576     << lMVA << " "
577     << std::endl;
578     }
579     return lId;
580     }
581 ceballos 1.10 Double_t JetIDMVA::correctedPt(const PFJet *iJet, FactorizedJetCorrector *iJetCorrector,
582     const PileupEnergyDensityCol *iPUEnergyDensity,
583 pharris 1.16 RhoUtilities::RhoType type,int iId) {
584 ceballos 1.10 Double_t Rho = 0.0;
585     switch(type) {
586     case RhoUtilities::MIT_RHO_VORONOI_HIGH_ETA:
587     Rho = iPUEnergyDensity->At(0)->Rho();
588     break;
589     case RhoUtilities::MIT_RHO_VORONOI_LOW_ETA:
590     Rho = iPUEnergyDensity->At(0)->RhoLowEta();
591     break;
592     case RhoUtilities::MIT_RHO_RANDOM_HIGH_ETA:
593     Rho = iPUEnergyDensity->At(0)->RhoRandom();
594     break;
595     case RhoUtilities::MIT_RHO_RANDOM_LOW_ETA:
596     Rho = iPUEnergyDensity->At(0)->RhoRandomLowEta();
597     break;
598     case RhoUtilities::CMS_RHO_RHOKT6PFJETS:
599     Rho = iPUEnergyDensity->At(0)->RhoKt6PFJets();
600     break;
601     default:
602     // use the old default
603     Rho = iPUEnergyDensity->At(0)->Rho();
604     break;
605     }
606    
607 pharris 1.1 const FourVectorM rawMom = iJet->RawMom();
608     iJetCorrector->setJetEta(rawMom.Eta());
609     iJetCorrector->setJetPt (rawMom.Pt());
610     iJetCorrector->setJetPhi(rawMom.Phi());
611     iJetCorrector->setJetE (rawMom.E());
612 ceballos 1.10 iJetCorrector->setRho (Rho);
613 pharris 1.1 iJetCorrector->setJetA (iJet->JetArea());
614     iJetCorrector->setJetEMF(-99.0);
615 pharris 1.16 Double_t correction = 1.;
616     if(iId < 0 || iId == 100) correction = iJetCorrector->getCorrection();
617     std::vector<float> lCorrections; if(iId != -1 && iId != 100) lCorrections = iJetCorrector->getSubCorrections();
618     if(iId > -1 && iId < int(lCorrections.size())) correction = lCorrections[iId];
619     if(iId == 100) {
620     iJetCorrector->setJetEta(rawMom.Eta());
621     iJetCorrector->setJetPt (rawMom.Pt());
622     iJetCorrector->setJetPhi(rawMom.Phi());
623     iJetCorrector->setJetE (rawMom.E());
624     iJetCorrector->setRho (Rho);
625     iJetCorrector->setJetA (iJet->JetArea());
626     iJetCorrector->setJetEMF(-99.0);
627     lCorrections = iJetCorrector->getSubCorrections();
628     double correction2 = 1;
629     correction2 *= lCorrections[0];
630     correction = correction-correction2;
631     }
632    
633 pharris 1.1 return rawMom.Pt()*correction;
634     }