ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/JetIDMVA.cc
(Generate patch)

Comparing UserCode/MitPhysics/Utils/src/JetIDMVA.cc (file contents):
Revision 1.3 by pharris, Wed Apr 4 09:51:31 2012 UTC vs.
Revision 1.20 by arapyan, Fri Apr 5 13:14:28 2013 UTC

# Line 1 | Line 1
1 + #include "FWCore/PythonParameterSet/interface/MakeParameterSets.h"
2 + #include "FWCore/ParameterSet/interface/ParameterSet.h"
3   #include "MitPhysics/Utils/interface/JetIDMVA.h"
4   #include "MitPhysics/Utils/interface/JetTools.h"
5   #include "MitAna/DataTree/interface/StableData.h"
# Line 13 | Line 15 | using namespace mithep;
15  
16   //--------------------------------------------------------------------------------------------------
17   JetIDMVA::JetIDMVA() :
18 <  fJetPtMin(10), //We need to lower this
19 <  fMethodName("JetIDMA"),
18 >  fJetPtMin(0.)  , //We need to lower this
19 >  fDZCut   (0.2),
20 >  fLowPtMethodName ("JetIDMVALowPt" ),
21 >  fHighPtMethodName("JetIDMVAHighPt"),
22    fIsInitialized(kFALSE),
23 <  fNPV    (0),
24 <  fJPt1   (0),
25 <  fJEta1  (0),
26 <  fJPhi1  (0),
27 <  fJD01   (0),
28 <  fJDZ1   (0),
29 <  fJM1    (0),
30 <  fNPart1 (0),
31 <  fLPt1   (0),
32 <  fLEta1  (0),
33 <  fLPhi1  (0),
34 <  fSPt1   (0),
35 <  fSEta1  (0),
36 <  fSPhi1  (0),
37 <  fNEPt1  (0),
38 <  fNEEta1 (0),
39 <  fNEPhi1 (0),
40 <  fEMPt1  (0),
37 <  fEMEta1 (0),
38 <  fEMPhi1 (0),
39 <  fChPt1  (0),
40 <  fChPhi1 (0),
41 <  fLFr1   (0),
42 <  fDRLC1  (0),
43 <  fDRLS1  (0),
44 <  fDRM1   (0),
45 <  fDRNE1  (0),
46 <  fDREM1  (0),
47 <  fDRCH1  (0)
23 >  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 >  fPtD      (0),
35 >  fFrac01   (0),
36 >  fFrac02   (0),
37 >  fFrac03   (0),
38 >  fFrac04   (0),
39 >  fFrac05   (0),
40 >  fDR2Mean  (0)
41   {    
42    fReader = 0;
43   }
44   //--------------------------------------------------------------------------------------------------
45 < JetIDMVA::~JetIDMVA()
46 < {
47 <  fReader = 0;
45 > JetIDMVA::~JetIDMVA() {
46 >
47 >  delete fReader;
48 >  delete fLowPtReader;
49   }
50  
51   //--------------------------------------------------------------------------------------------------
52 < void JetIDMVA::Initialize( TString iMethodName,
53 <                           TString iWeights,
54 <                            JetIDMVA::MVAType iType)
52 > void JetIDMVA::Initialize( JetIDMVA::CutType iCutType,
53 >                           TString iLowPtWeights,
54 >                           TString iHighPtWeights,
55 >                           JetIDMVA::MVAType iType,
56 >                           TString iCutFileName)
57   {
58    
59    fIsInitialized = kTRUE;
64  fMethodName    = iMethodName;
60    fType          = iType;
61 +  fCutType       = iCutType;
62 +  
63 +  std::string lCutId = "JetIdParams";
64 +  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 +  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 +  fLowPtReader   = 0;
115 +  fLowPtReader   = new TMVA::Reader( "!Color:!Silent:Error" );  
116 +  fLowPtReader->AddVariable( "nvtx"     , &fNVtx      );
117 +  if(fType != k53) fLowPtReader->AddVariable( "jetPt"    , &fJPt1      );  
118 +  if(fType != k53) fLowPtReader->AddVariable( "jetEta"   , &fJEta1     );
119 +  if(fType != k53) fLowPtReader->AddVariable( "jetPhi"   , &fJPhi1     );            
120 +  fLowPtReader->AddVariable( "dZ"       , &fJDZ1      );
121 +  if(fType != k53 && fType != k53MET) fLowPtReader->AddVariable( "d0"       , &fJD01      );
122 +  fLowPtReader->AddVariable( "beta"     , &fBeta      );
123 +  fLowPtReader->AddVariable( "betaStar" , &fBetaStar  );
124 +  fLowPtReader->AddVariable( "nCharged" , &fNCharged  );
125 +  fLowPtReader->AddVariable( "nNeutrals", &fNNeutrals );
126 +  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 +  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 +  if(fType == k53) fLowPtReader->AddSpectator( "jetPt"   , &fJPt1      );  
135 +  if(fType == k53) fLowPtReader->AddSpectator( "jetEta"  , &fJEta1     );
136 +  if(fType == k53) fLowPtReader->AddSpectator( "jetPhi"  , &fJPhi1     );  
137 +  
138    fReader        = 0;
139    fReader        = new TMVA::Reader( "!Color:!Silent:Error" );  
140    if (fType == kBaseline) {
141 <    fReader->AddSpectator( "npv"     , &fNPV    );
142 <    fReader->AddVariable( "jspt_1"   , &fJPt1   );
143 <    fReader->AddVariable( "jseta_1"  , &fJEta1  );
144 <    fReader->AddVariable( "jsphi_1"  , &fJPhi1  );
145 <    fReader->AddVariable( "jd0_1"    , &fJD01   );
146 <    fReader->AddVariable( "jdZ_1"    , &fJDZ1   );
147 <    fReader->AddVariable( "jm_1"     , &fJM1    );
148 <    fReader->AddVariable( "npart_1"  , &fNPart1 );
149 <    fReader->AddVariable( "lpt_1"    , &fLPt1   );
150 <    fReader->AddVariable( "leta_1"   , &fLEta1  );
151 <    fReader->AddVariable( "lphi_1"   , &fLPhi1  );
152 <    fReader->AddVariable( "spt_1"    , &fSPt1   );
153 <    fReader->AddVariable( "seta_1"   , &fSEta1  );
154 <    fReader->AddVariable( "sphi_1"   , &fSPhi1  );
155 <    fReader->AddVariable( "lnept_1"  , &fNEPt1  );
156 <    fReader->AddVariable( "lneeta_1" , &fNEEta1 );
85 <    fReader->AddVariable( "lnephi_1" , &fNEPhi1 );
86 <    fReader->AddVariable( "lempt_1"  , &fEMPt1  );
87 <    fReader->AddVariable( "lemeta_1" , &fEMEta1 );
88 <    fReader->AddVariable( "lemphi_1" , &fEMPhi1 );
89 <    fReader->AddVariable( "lchpt_1"  , &fChPt1  );
90 <    fReader->AddVariable( "lchphi_1" , &fChPhi1 );
91 <    fReader->AddVariable( "lLfr_1"   , &fLFr1   );
92 <    fReader->AddVariable( "drlc_1"   , &fDRLC1  );
93 <    fReader->AddVariable( "drls_1"   , &fDRLS1  );
94 <    fReader->AddVariable( "drm_1"    , &fDRM1   );
95 <    fReader->AddVariable( "drmne_1"  , &fDRNE1  );
96 <    fReader->AddVariable( "drem_1"   , &fDREM1  );
97 <    fReader->AddVariable( "drch_1"   , &fDRCH1  );
141 >    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    }
158 <  fReader->BookMVA(fMethodName , iWeights );
158 >  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 >  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 >    fReader->AddVariable( "dR2Mean"  , &fDR2Mean   );
197 >    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 >    fReader->AddVariable( "dR2Mean"  , &fDR2Mean   );
218 >    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 >  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 >
243 >  fLowPtReader->BookMVA(fLowPtMethodName  , iLowPtWeights );
244 >  fReader->BookMVA(fHighPtMethodName , iHighPtWeights );
245    std::cout << "Jet ID MVA Initialization\n";
246 <  std::cout << "MethodName : " << fMethodName << " , type == " << fType << std::endl;
246 >  std::cout << "MethodName : " << fLowPtMethodName << " , type == " << fType << std::endl;
247 >
248   }
249  
250   //--------------------------------------------------------------------------------------------------
# Line 109 | Line 255 | Double_t JetIDMVA::MVAValue(
255                              Float_t iJPhi1  ,
256                              Float_t iJD01   ,
257                              Float_t iJDZ1   ,
258 <                            Float_t iJM1    ,
259 <                            Float_t iNPart1 ,
260 <                            Float_t iLPt1   ,
261 <                            Float_t iLEta1  ,
262 <                            Float_t iLPhi1  ,
263 <                            Float_t iSPt1   ,
264 <                            Float_t iSEta1  ,
265 <                            Float_t iSPhi1  ,
266 <                            Float_t iNEPt1  ,
267 <                            Float_t iNEEta1 ,
268 <                            Float_t iNEPhi1 ,
269 <                            Float_t iEMPt1  ,
124 <                            Float_t iEMEta1 ,
125 <                            Float_t iEMPhi1 ,
126 <                            Float_t iChPt1  ,
127 <                            Float_t iChPhi1 ,
128 <                            Float_t iLFr1   ,
129 <                            Float_t iDRLC1  ,
130 <                            Float_t iDRLS1  ,
131 <                            Float_t iDRM1   ,
132 <                            Float_t iDRNE1 ,
133 <                            Float_t iDREM1  ,
134 <                            Float_t iDRCH1  
258 >                            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 >                            Float_t iFrac05  ,
268 >                            Float_t iDR2Mean ,
269 >                            Float_t iPtD
270                              ){
271    
272 <  if (!fIsInitialized) {
272 >  if(!fIsInitialized) {
273      std::cout << "Error: JetIDMVA not properly initialized.\n";
274      return -9999;
275    }
276    
277 <  fNPV    = iNPV;
278 <  fJPt1   = iJPt1;
279 <  fJEta1  = iJEta1;
280 <  fJPhi1  = fJPhi1;
281 <  fJD01   = iJD01;
282 <  fJDZ1   = iJDZ1;
283 <  fJM1    = iJM1 ;
284 <  fNPart1 = iNPart1;
285 <  fLPt1   = iLPt1;
286 <  fLEta1  = iLEta1;
287 <  fLPhi1  = iLPhi1;
288 <  fSPt1   = iSPt1;
289 <  fSEta1  = iSEta1;
290 <  fSPhi1  = iSPhi1;
291 <  fNEPt1  = iNEPt1;
292 <  fNEEta1 = iNEEta1;
293 <  fNEPhi1 = iNEPhi1;
294 <  fEMPt1  = iEMPt1;
160 <  fEMEta1 = iEMEta1;
161 <  fEMPhi1 = iEMPhi1;
162 <  fChPt1  = iChPt1;
163 <  fChPhi1 = iChPhi1;
164 <  fLFr1   = iLFr1;
165 <  fDRLC1  = iDRLC1;
166 <  fDRLS1  = iDRLS1;
167 <  fDRM1   = iDRM1;
168 <  fDRNE1  = iDRNE1;
169 <  fDREM1  = iDREM1;
170 <  fDRCH1  = iDRCH1;  
277 >  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 >  fPtD       = iPtD;
289 >  fFrac01    = iFrac01;
290 >  fFrac02    = iFrac02;
291 >  fFrac03    = iFrac03;
292 >  fFrac04    = iFrac04;
293 >  fFrac05    = iFrac05;
294 >  fDR2Mean   = iDR2Mean;
295  
296    Double_t lMVA = -9999;  
297 <  lMVA = fReader->EvaluateMVA( fMethodName );
298 <
297 >  if(iJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName  );
298 >  if(iJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
299    return lMVA;
300   }
301   //--------------------------------------------------------------------------------------------------
302 < Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,
302 > 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 > Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices,
312                        FactorizedJetCorrector *iJetCorrector,
313 <                      const PileupEnergyDensityCol *iPileupEnergyDensity) {
313 >                      const PileupEnergyDensityCol *iPileupEnergyDensity,
314 >                      RhoUtilities::RhoType type) {
315 >  
316    if(!JetTools::passPFLooseId(iJet))                 return false;
317 <  if(correctedPt(iJet,iJetCorrector,iPileupEnergyDensity) < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false; //This line is a bug in the Met training
318 <  if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
319 <  double lMVA = MVAValue(iJet,iVertex,iJetCorrector,iPileupEnergyDensity);
320 <  if(lMVA < -0.8)                            return false;
321 <  if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
317 >  if(fabs(iJet->Eta()) > 4.99)                       return false;
318 >  
319 >  double lMVA = MVAValue   (iJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity);
320 >  double lPt  = correctedPt(iJet,                  iJetCorrector,iPileupEnergyDensity,type);
321 >  if(lPt < fJetPtMin)                         return false;
322 >  
323 >  int lPtId = 0;
324 >  if(lPt > 10 && lPt < 20) lPtId = 1;
325 >  if(lPt > 20 && lPt < 30) lPtId = 2;
326 >  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    return true;
336   }
337   //--------------------------------------------------------------------------------------------------
338 < Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex) {
338 > 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 && 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 <  double lMVA = MVAValue(iJet,iVertex);
343 <  if(lMVA < -0.8)                            return false;
344 <  if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
340 >  if(iJet->Pt()        < fJetPtMin) return false;
341 >  if(fabs(iJet->Eta()) > 4.99)      return false;
342 >  //if(fType == kCut) passCut(iJet,iVertex,iVertices);
343 >
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 > Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices) {
365 >  if(!JetTools::passPFLooseId(iJet))                 return false;
366 >  if(iJet->Pt()        < fJetPtMin) return false;
367 >  if(fabs(iJet->Eta()) > 4.99)      return false;
368 >  if(fType == kCut) return passCut(iJet,iVertex,iVertices);
369 >  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    return true;
384 +  //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   }
389   //--------------------------------------------------------------------------------------------------
390 < Double_t JetIDMVA::MVAValue(const PFJet *iJet,const Vertex *iVertex, //Vertex here is the PV
390 > Double_t JetIDMVA::MVAValue(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices, //Vertex here is the PV
391                              FactorizedJetCorrector *iJetCorrector,
392                              const PileupEnergyDensityCol *iPileupEnergyDensity,
393                              Bool_t printDebug) {
# Line 209 | Line 399 | Double_t JetIDMVA::MVAValue(const PFJet
399    if(!JetTools::passPFLooseId(iJet)) return -2.;
400  
401    //set all input variables
402 +  fNVtx      = iVertices->GetEntries();
403    fJPt1      = correctedPt(iJet,iJetCorrector,iPileupEnergyDensity);
404    fJEta1     = iJet->RawMom().Eta();
405    fJPhi1     = iJet->RawMom().Phi();
406 <  fJM1       = iJet->Mass();
407 <
408 <  const mithep::PFCandidate *lLead     = JetTools::leadCand(iJet,-1);
409 <  const mithep::PFCandidate *lSecond   = JetTools::leadCand(iJet,-1,true);
410 <  const mithep::PFCandidate *lLeadNeut = JetTools::leadCand(iJet ,5);
411 <  const mithep::PFCandidate *lLeadEm   = JetTools::leadCand(iJet ,4);
412 <  const mithep::PFCandidate *lLeadCh   = JetTools::leadCand(iJet ,1);
413 <
414 <  fJD01         = JetTools::impactParameter(iJet,iVertex);  
415 <  fJDZ1         = JetTools::impactParameter(iJet,iVertex,true);
416 <  fNPart1       = iJet->NPFCands();
417 <  fLPt1         = lLead    ->Pt();
418 <  fLEta1        = lLead    ->Eta();
419 <  fLPhi1        = lLead    ->Phi();
420 <  fSPt1         = lSecond  ->Pt();
421 <  fSEta1        = lSecond  ->Eta();
422 <  fSPhi1        = lSecond  ->Phi();
423 <  fNEPt1        = lLeadNeut->Pt();
424 <  fNEEta1       = lLeadNeut->Eta();
234 <  fNEPhi1       = lLeadNeut->Phi();
235 <  fEMPt1        = lLeadEm  ->Pt();
236 <  fEMEta1       = lLeadEm  ->Eta();
237 <  fEMPhi1       = lLeadEm  ->Phi();
238 <  fChPt1        = lLeadCh  ->Pt();
239 <  //fChEta1       = lLeadCh  ->Eta();
240 <  fChPhi1       = lLeadCh  ->Phi();
241 <  fLFr1         = lLead->Pt()/iJet->Pt();
242 <
243 <  fDRLC1        = MathUtils::DeltaR(iJet->Mom(),lLead  ->Mom());
244 <  fDRLS1        = MathUtils::DeltaR(iJet->Mom(),lSecond->Mom());
245 <  fDRM1         = JetTools::dRMean (iJet,-1);
246 <  fDRNE1        = JetTools::dRMean (iJet, 5);
247 <  fDREM1        = JetTools::dRMean (iJet, 4);
248 <  fDRCH1        = JetTools::dRMean (iJet, 1);
249 <
250 <  double lMVA = fReader->EvaluateMVA( fMethodName );
251 <  
406 >  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 >  fDR2Mean   = JetTools::dR2Mean(iJet,-1);
415 >  fPtD       = JetTools::W(iJet,-1,0);  
416 >  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 >  if(fJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName  );
424 >  if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
425    if (printDebug == kTRUE) {
426      std::cout << "Debug Jet MVA: "
427 <              << fNPV    << " "
428 <              << fJPt1   << " "
429 <              << fJEta1  << " "
430 <              << fJPhi1  << " "
431 <              << fJD01   << " "
432 <              << fJDZ1   << " "
433 <              << fJM1    << " "
434 <              << fNPart1 << " "
435 <              << fLPt1   << " "
436 <              << fLEta1  << " "
437 <              << fLPhi1  << " "
438 <              << fSPt1   << " "
439 <              << fSEta1  << " "
440 <              << fSPhi1  << " "
441 <              << fNEPt1  << " "
442 <              << fNEEta1 << " "
443 <              << fNEPhi1 << " "
444 <              << fEMPt1  << " "
272 <              << fEMEta1 << " "
273 <              << fEMPhi1 << " "
274 <              << fChPt1  << " "
275 <              << fChPhi1 << " "
276 <              << fLFr1   << " "
277 <              << fDRLC1  << " "
278 <              << fDRLS1  << " "
279 <              << fDRM1   << " "
280 <              << fDRNE1 << " "
281 <              << fDREM1  << " "
282 <              << fDRCH1  << " "
427 >              << fNVtx      << " "
428 >              << fJPt1      << " "
429 >              << fJEta1     << " "
430 >              << fJPhi1     << " "
431 >              << fJD01      << " "
432 >              << fJDZ1      << " "
433 >              << fBeta      << " "
434 >              << fBetaStar  << " "
435 >              << fNCharged  << " "
436 >              << fNNeutrals << " "
437 >              << fDRMean    << " "
438 >              << fPtD       << " "
439 >              << fFrac01    << " "
440 >              << fFrac02    << " "
441 >              << fFrac03    << " "
442 >              << fFrac04    << " "
443 >              << fFrac05    << " "
444 >              << fDR2Mean    
445                << " === : === "
446                << lMVA << " "    
447 <              << std::endl;
447 >              << " -----> raw pt " << iJet->Pt() << std::endl;
448    }
449 <
449 >  //std::cout << " === " << iJet->Pt() << " -- " << iJet->Eta() << " -- "  << fPtD << " -- " << lMVA << std::endl;
450    return lMVA;
451   }
452 < Double_t JetIDMVA::MVAValue(const PFJet *iJet,const Vertex *iVertex, //Vertex here is the PV
452 > Double_t JetIDMVA::MVAValue(const PFJet *iJet,const Vertex *iVertex, const VertexCol *iVertices,//Vertex here is the PV
453                              Bool_t printDebug) {
454    
455    if (!fIsInitialized) {
# Line 297 | Line 459 | Double_t JetIDMVA::MVAValue(const PFJet
459    if(!JetTools::passPFLooseId(iJet)) return -2.;
460  
461    //set all input variables
462 +  fNVtx      = iVertices->GetEntries();
463    fJPt1      = iJet->Pt();
464    fJEta1     = iJet->RawMom().Eta();
465    fJPhi1     = iJet->RawMom().Phi();
466 <  fJM1       = iJet->Mass();
467 <
468 <  const mithep::PFCandidate *lLead     = JetTools::leadCand(iJet,-1);
469 <  const mithep::PFCandidate *lSecond   = JetTools::leadCand(iJet,-1,true);
470 <  const mithep::PFCandidate *lLeadNeut = JetTools::leadCand(iJet ,5);
471 <  const mithep::PFCandidate *lLeadEm   = JetTools::leadCand(iJet ,4);
472 <  const mithep::PFCandidate *lLeadCh   = JetTools::leadCand(iJet ,1);
473 <
474 <  fJD01         = JetTools::impactParameter(iJet,iVertex);  
475 <  fJDZ1         = JetTools::impactParameter(iJet,iVertex,true);
476 <  fNPart1       = iJet->NPFCands();
477 <  fLPt1         = lLead    ->Pt();
478 <  fLEta1        = lLead    ->Eta();
479 <  fLPhi1        = lLead    ->Phi();
480 <  fSPt1         = lSecond  ->Pt();
481 <  fSEta1        = lSecond  ->Eta();
482 <  fSPhi1        = lSecond  ->Phi();
483 <  fNEPt1        = lLeadNeut->Pt();
484 <  fNEEta1       = lLeadNeut->Eta();
322 <  fNEPhi1       = lLeadNeut->Phi();
323 <  fEMPt1        = lLeadEm  ->Pt();
324 <  fEMEta1       = lLeadEm  ->Eta();
325 <  fEMPhi1       = lLeadEm  ->Phi();
326 <  fChPt1        = lLeadCh  ->Pt();
327 <  //fChEta1       = lLeadCh  ->Eta();
328 <  fChPhi1       = lLeadCh  ->Phi();
329 <  fLFr1         = lLead->Pt()/iJet->RawMom().Pt();
330 <
331 <  fDRLC1        = MathUtils::DeltaR(iJet->Mom(),lLead  ->Mom());
332 <  fDRLS1        = MathUtils::DeltaR(iJet->Mom(),lSecond->Mom());
333 <  fDRM1         = JetTools::dRMean (iJet,-1);
334 <  fDRNE1        = JetTools::dRMean (iJet, 5);
335 <  fDREM1        = JetTools::dRMean (iJet, 4);
336 <  fDRCH1        = JetTools::dRMean (iJet, 1);
337 <
338 <  double lMVA = fReader->EvaluateMVA( fMethodName );
339 <  
466 >  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 >  fPtD       = JetTools::W(iJet,-1,0);  
473 >  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 >
481 >  double lMVA = 0;
482 >  if(fJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName  );
483 >  if(fJPt1 > 10) lMVA = fReader     ->EvaluateMVA( fHighPtMethodName );
484 >  
485    if (printDebug == kTRUE) {
486      std::cout << "Debug Jet MVA: "
487 <              << fNPV    << " "
488 <              << fJPt1   << " "
489 <              << fJEta1  << " "
490 <              << fJPhi1  << " "
491 <              << fJD01   << " "
492 <              << fJDZ1   << " "
493 <              << fJM1    << " "
494 <              << fNPart1 << " "
495 <              << fLPt1   << " "
496 <              << fLEta1  << " "
497 <              << fLPhi1  << " "
498 <              << fSPt1   << " "
499 <              << fSEta1  << " "
500 <              << fSPhi1  << " "
501 <              << fNEPt1  << " "
502 <              << fNEEta1 << " "
503 <              << fNEPhi1 << " "
504 <              << fEMPt1  << " "
360 <              << fEMEta1 << " "
361 <              << fEMPhi1 << " "
362 <              << fChPt1  << " "
363 <              << fChPhi1 << " "
364 <              << fLFr1   << " "
365 <              << fDRLC1  << " "
366 <              << fDRLS1  << " "
367 <              << fDRM1   << " "
368 <              << fDRNE1 << " "
369 <              << fDREM1  << " "
370 <              << fDRCH1  << " "
487 >              << fNVtx      << " "
488 >              << fJPt1      << " "
489 >              << fJEta1     << " "
490 >              << fJPhi1     << " "
491 >              << fJD01      << " "
492 >              << fJDZ1      << " "
493 >              << fBeta      << " "
494 >              << fBetaStar  << " "
495 >              << fNCharged  << " "
496 >              << fNNeutrals << " "
497 >              << fDRMean    << " "
498 >              << fPtD       << " "
499 >              << fFrac01    << " "
500 >              << fFrac02    << " "
501 >              << fFrac03    << " "
502 >              << fFrac04    << " "
503 >              << fFrac05    << " "
504 >              << fDR2Mean    
505                << " === : === "
506                << lMVA << " "    
507                << std::endl;
# Line 375 | Line 509 | Double_t JetIDMVA::MVAValue(const PFJet
509  
510    return lMVA;
511   }
512 < Double_t JetIDMVA::correctedPt(const PFJet *iJet, FactorizedJetCorrector *iJetCorrector,const PileupEnergyDensityCol *iPUEnergyDensity) {
512 > //--------------------------------------------------------------------------------------------------
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 > Double_t JetIDMVA::correctedPt(const PFJet *iJet, FactorizedJetCorrector *iJetCorrector,
582 >                               const PileupEnergyDensityCol *iPUEnergyDensity,
583 >                               RhoUtilities::RhoType type,int iId) {
584 >  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    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 <  iJetCorrector->setRho   (iPUEnergyDensity->At(0)->RhoHighEta());
612 >  iJetCorrector->setRho   (Rho);
613    iJetCorrector->setJetA  (iJet->JetArea());
614    iJetCorrector->setJetEMF(-99.0);    
615 <  Double_t correction = iJetCorrector->getCorrection();
615 >  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    return rawMom.Pt()*correction;
634   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines