ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/JetIDMVA.cc
Revision: 1.16
Committed: Tue Sep 25 15:39:15 2012 UTC (12 years, 7 months ago) by pharris
Content type: text/plain
Branch: MAIN
Changes since 1.15: +105 -2 lines
Log Message:
Updated MVA Met tool

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