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.1 by pharris, Wed Mar 21 18:56:25 2012 UTC vs.
Revision 1.15 by pharris, Tue Jun 12 09:43:57 2012 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 <  fMethodName("JetIDMA"),
18 >  fJetPtMin(0.)  , //We need to lower this
19 >  fDZCut   (0.2),
20 >  fLowPtMethodName ("JetIDMVALowPt" ),
21 >  fHighPtMethodName("JetIDMVAHighPt"),
22    fIsInitialized(kFALSE),
23 <  fJetPtMin(10), //We need to lower this
24 <  fNPV    (0),
25 <  fJPt1   (0),
26 <  fJEta1  (0),
27 <  fJPhi1  (0),
28 <  fJD01   (0),
29 <  fJDZ1   (0),
30 <  fJM1    (0),
31 <  fNPart1 (0),
32 <  fLPt1   (0),
33 <  fLEta1  (0),
34 <  fLPhi1  (0),
35 <  fSPt1   (0),
36 <  fSEta1  (0),
37 <  fSPhi1  (0),
38 <  fNEPt1  (0),
39 <  fNEEta1 (0),
35 <  fNEPhi1 (0),
36 <  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 >  fFrac01   (0),
35 >  fFrac02   (0),
36 >  fFrac03   (0),
37 >  fFrac04   (0),
38 >  fFrac05   (0),
39 >  fDR2Mean  (0)
40   {    
41    fReader = 0;
42   }
43   //--------------------------------------------------------------------------------------------------
44 < JetIDMVA::~JetIDMVA()
45 < {
46 <  fReader = 0;
44 > JetIDMVA::~JetIDMVA() {
45 >
46 >  delete fReader;
47 >  delete fLowPtReader;
48   }
49  
50   //--------------------------------------------------------------------------------------------------
51 < void JetIDMVA::Initialize( TString iMethodName,
52 <                           TString iWeights,
53 <                            JetIDMVA::MVAType iType)
51 > void JetIDMVA::Initialize( JetIDMVA::CutType iCutType,
52 >                           TString iLowPtWeights,
53 >                           TString iHighPtWeights,
54 >                           JetIDMVA::MVAType iType,
55 >                           TString iCutFileName)
56   {
57    
58    fIsInitialized = kTRUE;
64  fMethodName    = iMethodName;
59    fType          = iType;
60 +  fCutType       = iCutType;
61 +  
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 +  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    fReader        = 0;
130    fReader        = new TMVA::Reader( "!Color:!Silent:Error" );  
131    if (fType == kBaseline) {
132 <    fReader->AddSpectator( "npv"     , &fNPV    );
133 <    fReader->AddVariable( "jspt_1"   , &fJPt1   );
134 <    fReader->AddVariable( "jseta_1"  , &fJEta1  );
135 <    fReader->AddVariable( "jsphi_1"  , &fJPhi1  );
136 <    fReader->AddVariable( "jd0_1"    , &fJD01   );
137 <    fReader->AddVariable( "jdZ_1"    , &fJDZ1   );
138 <    fReader->AddVariable( "jm_1"     , &fJM1    );
139 <    fReader->AddVariable( "npart_1"  , &fNPart1 );
140 <    fReader->AddVariable( "lpt_1"    , &fLPt1   );
141 <    fReader->AddVariable( "leta_1"   , &fLEta1  );
142 <    fReader->AddVariable( "lphi_1"   , &fLPhi1  );
143 <    fReader->AddVariable( "spt_1"    , &fSPt1   );
144 <    fReader->AddVariable( "seta_1"   , &fSEta1  );
145 <    fReader->AddVariable( "sphi_1"   , &fSPhi1  );
146 <    fReader->AddVariable( "lnept_1"  , &fNEPt1  );
147 <    fReader->AddVariable( "lneeta_1" , &fNEEta1 );
148 <    fReader->AddVariable( "lnephi_1" , &fNEPhi1 );
149 <    fReader->AddVariable( "lempt_1"  , &fEMPt1  );
150 <    fReader->AddVariable( "lemeta_1" , &fEMEta1 );
151 <    fReader->AddVariable( "lemphi_1" , &fEMPhi1 );
152 <    fReader->AddVariable( "lchpt_1"  , &fChPt1  );
153 <    fReader->AddVariable( "lchphi_1" , &fChPhi1 );
154 <    fReader->AddVariable( "lLfr_1"   , &fLFr1   );
155 <    fReader->AddVariable( "drlc_1"   , &fDRLC1  );
156 <    fReader->AddVariable( "drls_1"   , &fDRLS1  );
157 <    fReader->AddVariable( "drm_1"    , &fDRM1   );
158 <    fReader->AddVariable( "drmne_1"  , &fDRNE1  );
159 <    fReader->AddVariable( "drem_1"   , &fDREM1  );
160 <    fReader->AddVariable( "drch_1"   , &fDRCH1  );
132 >    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 >  }
149 >  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 <  fReader->BookMVA(fMethodName , iWeights );
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 >
181 >  fLowPtReader->BookMVA(fLowPtMethodName  , iLowPtWeights );
182 >  fReader->BookMVA(fHighPtMethodName , iHighPtWeights );
183    std::cout << "Jet ID MVA Initialization\n";
184 <  std::cout << "MethodName : " << fMethodName << " , type == " << fType << std::endl;
184 >  std::cout << "MethodName : " << fLowPtMethodName << " , type == " << fType << std::endl;
185 >
186   }
187  
188   //--------------------------------------------------------------------------------------------------
# Line 109 | Line 193 | Double_t JetIDMVA::MVAValue(
193                              Float_t iJPhi1  ,
194                              Float_t iJD01   ,
195                              Float_t iJDZ1   ,
196 <                            Float_t iJM1    ,
197 <                            Float_t iNPart1 ,
198 <                            Float_t iLPt1   ,
199 <                            Float_t iLEta1  ,
200 <                            Float_t iLPhi1  ,
201 <                            Float_t iSPt1   ,
202 <                            Float_t iSEta1  ,
203 <                            Float_t iSPhi1  ,
204 <                            Float_t iNEPt1  ,
205 <                            Float_t iNEEta1 ,
206 <                            Float_t iNEPhi1 ,
123 <                            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  
196 >                            Float_t iBeta   ,
197 >                            Float_t iBetaStar,
198 >                            Float_t iNCharged,
199 >                            Float_t iNNeutrals,
200 >                            Float_t iDRMean  ,
201 >                            Float_t iFrac01  ,
202 >                            Float_t iFrac02  ,
203 >                            Float_t iFrac03  ,
204 >                            Float_t iFrac04  ,
205 >                            Float_t iFrac05  ,
206 >                            Float_t iDR2Mean  
207                              ){
208    
209 <  if (!fIsInitialized) {
209 >  if(!fIsInitialized) {
210      std::cout << "Error: JetIDMVA not properly initialized.\n";
211      return -9999;
212    }
213    
214 <  fNPV    = iNPV;
215 <  fJPt1   = iJPt1;
216 <  fJEta1  = iJEta1;
217 <  fJPhi1  = fJPhi1;
218 <  fJD01   = iJD01;
219 <  fJDZ1   = iJDZ1;
220 <  fJM1    = iJM1 ;
221 <  fNPart1 = iNPart1;
222 <  fLPt1   = iLPt1;
223 <  fLEta1  = iLEta1;
224 <  fLPhi1  = iLPhi1;
225 <  fSPt1   = iSPt1;
226 <  fSEta1  = iSEta1;
227 <  fSPhi1  = iSPhi1;
228 <  fNEPt1  = iNEPt1;
229 <  fNEEta1 = iNEEta1;
230 <  fNEPhi1 = iNEPhi1;
159 <  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;  
214 >  fNVtx      = iNPV;
215 >  fJPt1      = iJPt1;
216 >  fJEta1     = iJEta1;
217 >  fJPhi1     = fJPhi1;
218 >  fJD01      = iJD01;
219 >  fJDZ1      = iJDZ1;
220 >  fBeta      = iBeta;
221 >  fBetaStar  = iBetaStar;
222 >  fNCharged  = iNCharged;
223 >  fNNeutrals = iNNeutrals;
224 >  fDRMean    = iDRMean;
225 >  fFrac01    = iFrac01;
226 >  fFrac02    = iFrac02;
227 >  fFrac03    = iFrac03;
228 >  fFrac04    = iFrac04;
229 >  fFrac05    = iFrac05;
230 >  fDR2Mean   = iDR2Mean;
231  
232    Double_t lMVA = -9999;  
233 <  lMVA = fReader->EvaluateMVA( fMethodName );
233 >  if(iJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName  );
234 >  if(iJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
235  
236    return lMVA;
237   }
238   //--------------------------------------------------------------------------------------------------
239 < Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,
239 > Bool_t JetIDMVA::passPt(const PFJet *iJet, FactorizedJetCorrector *iJetCorrector,
240 >                        const PileupEnergyDensityCol *iPileupEnergyDensity,
241 >                        RhoUtilities::RhoType type) {
242 >  if(iJetCorrector == 0) return (iJet->Pt() > fJetPtMin);
243 >  double lPt  = correctedPt(iJet,                  iJetCorrector,iPileupEnergyDensity,type);
244 >  if(lPt < fJetPtMin)                         return false;
245 >  return true;
246 > }
247 > //--------------------------------------------------------------------------------------------------
248 > Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices,
249                        FactorizedJetCorrector *iJetCorrector,
250 <                      const PileupEnergyDensityCol *iPileupEnergyDensity) {
250 >                      const PileupEnergyDensityCol *iPileupEnergyDensity,
251 >                      RhoUtilities::RhoType type) {
252 >  
253 >  if(!JetTools::passPFLooseId(iJet))                 return false;
254 >  if(fabs(iJet->Eta()) > 4.99)                       return false;
255 >  
256 >  double lMVA = MVAValue   (iJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity);
257 >  double lPt  = correctedPt(iJet,                  iJetCorrector,iPileupEnergyDensity,type);
258 >  if(lPt < fJetPtMin)                         return false;
259 >  
260 >  int lPtId = 0;
261 >  if(lPt > 10 && lPt < 20) lPtId = 1;
262 >  if(lPt > 20 && lPt < 30) lPtId = 2;
263 >  if(lPt > 30                   ) lPtId = 3;
264 >  
265 >  int lEtaId = 0;
266 >  if(fabs(iJet->Eta()) > 2.5  && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
267 >  if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
268 >  if(fabs(iJet->Eta()) > 3.0  && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
269 >  
270 >  double lMVACut = fMVACut[lPtId][lEtaId];
271 >  if(lMVA < lMVACut) return false;
272 >  return true;
273 > }
274 > //--------------------------------------------------------------------------------------------------
275 > Bool_t JetIDMVA::passCut(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices) {
276 >  if(!JetTools::passPFLooseId(iJet))                 return false;
277 >  if(iJet->Pt()        < fJetPtMin) return false;
278 >  if(fabs(iJet->Eta()) > 4.99)      return false;
279 >  //if(fType == kCut) passCut(iJet,iVertex,iVertices);
280 >
281 >  double lPt = iJet->Pt();
282 >  int lPtId = 0;
283 >  if(lPt > 10 && lPt < 20) lPtId = 1;
284 >  if(lPt > 20 && lPt < 30) lPtId = 2;
285 >  if(lPt > 30                   ) lPtId = 3;
286 >  
287 >  int lEtaId = 0;
288 >  if(fabs(iJet->Eta()) > 2.5  && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
289 >  if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
290 >  if(fabs(iJet->Eta()) > 3.0  && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
291 >  float betaStarModified = JetTools::betaStarClassic(iJet,iVertex,iVertices)/log(iVertices ->GetEntries()-0.64);
292 >  float dR2Mean          = JetTools::dR2Mean(iJet,-1);
293 >  
294 >  if(betaStarModified < fBetaStarCut[lPtId][lEtaId] &&
295 >     dR2Mean          < fRMSCut     [lPtId][lEtaId]
296 >     ) return true;
297 >  
298 >  return false;
299 > }
300 > //--------------------------------------------------------------------------------------------------
301 > Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices) {
302    if(!JetTools::passPFLooseId(iJet))                 return false;
303 <  if(correctedPt(iJet,iJetCorrector,iPileupEnergyDensity) < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false; //This line is a bug in the Met training
304 <  if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
305 <  double lMVA = MVAValue(iJet,iVertex,iJetCorrector,iPileupEnergyDensity);
306 <  if(lMVA < -0.8)                            return false;
307 <  if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
303 >  if(iJet->Pt()        < fJetPtMin) return false;
304 >  if(fabs(iJet->Eta()) > 4.99)      return false;
305 >  if(fType == kCut) return passCut(iJet,iVertex,iVertices);
306 >  double lMVA = MVAValue(iJet,iVertex,iVertices);
307 >  
308 >  int lPtId = 0;
309 >  if(iJet->Pt() > 10 && iJet->Pt() < 20) lPtId = 1;
310 >  if(iJet->Pt() > 20 && iJet->Pt() < 30) lPtId = 2;
311 >  if(iJet->Pt() > 30                   ) lPtId = 3;
312 >  
313 >  int lEtaId = 0;
314 >  if(fabs(iJet->Eta()) > 2.5  && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
315 >  if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
316 >  if(fabs(iJet->Eta()) > 3.0  && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
317 >  
318 >  double lMVACut = fMVACut[lPtId][lEtaId];
319 >  if(lMVA < lMVACut) return false;
320    return true;
321 +  //if(lMVA < -0.8)                            return false;
322 +  //if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
323 +  //if(iJet->Pt() < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false; //This line is a bug in the Met training
324 +  //if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
325   }
326   //--------------------------------------------------------------------------------------------------
327 < Double_t JetIDMVA::MVAValue(const PFJet *iJet,const Vertex *iVertex, //Vertex here is the PV
327 > Double_t JetIDMVA::MVAValue(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices, //Vertex here is the PV
328                              FactorizedJetCorrector *iJetCorrector,
329                              const PileupEnergyDensityCol *iPileupEnergyDensity,
330                              Bool_t printDebug) {
# Line 199 | Line 336 | Double_t JetIDMVA::MVAValue(const PFJet
336    if(!JetTools::passPFLooseId(iJet)) return -2.;
337  
338    //set all input variables
339 +  fNVtx      = iVertices->GetEntries();
340    fJPt1      = correctedPt(iJet,iJetCorrector,iPileupEnergyDensity);
341    fJEta1     = iJet->RawMom().Eta();
342    fJPhi1     = iJet->RawMom().Phi();
343 <  fJM1       = iJet->Mass();
343 >  fJD01      = JetTools::impactParameter(iJet,iVertex);  
344 >  fJDZ1      = JetTools::impactParameter(iJet,iVertex,true);
345 >  fBeta      = JetTools::Beta(iJet,iVertex,fDZCut);
346 >  fBetaStar  = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
347 >  fNCharged  = iJet->ChargedMultiplicity();
348 >  fNNeutrals = iJet->NeutralMultiplicity();
349 >
350 >  fDRMean    = JetTools::dRMean(iJet,-1);
351 >  fDR2Mean   = JetTools::dR2Mean(iJet,-1);
352 >  fFrac01    = JetTools::frac  (iJet,0.1,0. ,-1);
353 >  fFrac02    = JetTools::frac  (iJet,0.2,0.1,-1);
354 >  fFrac03    = JetTools::frac  (iJet,0.3,0.2,-1);
355 >  fFrac04    = JetTools::frac  (iJet,0.4,0.3,-1);
356 >  fFrac05    = JetTools::frac  (iJet,0.5,0.4,-1);
357 >
358 >  double lMVA = 0;
359 >  if(fJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName  );
360 >  if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
361 >  if (printDebug == kTRUE) {
362 >    std::cout << "Debug Jet MVA: "
363 >              << fNVtx      << " "
364 >              << fJPt1      << " "
365 >              << fJEta1     << " "
366 >              << fJPhi1     << " "
367 >              << fJD01      << " "
368 >              << fJDZ1      << " "
369 >              << fBeta      << " "
370 >              << fBetaStar  << " "
371 >              << fNCharged  << " "
372 >              << fNNeutrals << " "
373 >              << fDRMean    << " "
374 >              << fFrac01    << " "
375 >              << fFrac02    << " "
376 >              << fFrac03    << " "
377 >              << fFrac04    << " "
378 >              << fFrac05    << " "
379 >              << fDRMean    
380 >              << " === : === "
381 >              << lMVA << " "    
382 >              << std::endl;
383 >  }
384  
385 <  const mithep::PFCandidate *lLead     = JetTools::leadCand(iJet,-1);
386 <  const mithep::PFCandidate *lSecond   = JetTools::leadCand(iJet,-1,true);
387 <  const mithep::PFCandidate *lLeadNeut = JetTools::leadCand(iJet ,5);
388 <  const mithep::PFCandidate *lLeadEm   = JetTools::leadCand(iJet ,4);
389 <  const mithep::PFCandidate *lLeadCh   = JetTools::leadCand(iJet ,1);
390 <
391 <  fJD01         = JetTools::impactParameter(iJet,iVertex);  
392 <  fJDZ1         = JetTools::impactParameter(iJet,iVertex,true);
393 <  fNPart1       = iJet->NPFCands();
394 <  fLPt1         = lLead    ->Pt();
217 <  fLEta1        = lLead    ->Eta();
218 <  fLPhi1        = lLead    ->Phi();
219 <  fSPt1         = lSecond  ->Pt();
220 <  fSEta1        = lSecond  ->Eta();
221 <  fSPhi1        = lSecond  ->Phi();
222 <  fNEPt1        = lLeadNeut->Pt();
223 <  fNEEta1       = lLeadNeut->Eta();
224 <  fNEPhi1       = lLeadNeut->Phi();
225 <  fEMPt1        = lLeadEm  ->Pt();
226 <  fEMEta1       = lLeadEm  ->Eta();
227 <  fEMPhi1       = lLeadEm  ->Phi();
228 <  fChPt1        = lLeadCh  ->Pt();
229 <  //fChEta1       = lLeadCh  ->Eta();
230 <  fChPhi1       = lLeadCh  ->Phi();
231 <  fLFr1         = lLead->Pt()/iJet->Pt();
232 <
233 <  fDRLC1        = MathUtils::DeltaR(iJet->Mom(),lLead  ->Mom());
234 <  fDRLS1        = MathUtils::DeltaR(iJet->Mom(),lSecond->Mom());
235 <  fDRM1         = JetTools::dRMean (iJet,-1);
236 <  fDRNE1        = JetTools::dRMean (iJet, 5);
237 <  fDREM1        = JetTools::dRMean (iJet, 4);
238 <  fDRCH1        = JetTools::dRMean (iJet, 1);
385 >  return lMVA;
386 > }
387 > Double_t JetIDMVA::MVAValue(const PFJet *iJet,const Vertex *iVertex, const VertexCol *iVertices,//Vertex here is the PV
388 >                            Bool_t printDebug) {
389 >  
390 >  if (!fIsInitialized) {
391 >    std::cout << "Error: JetIDMVA not properly initialized.\n";
392 >    return -9999;
393 >  }
394 >  if(!JetTools::passPFLooseId(iJet)) return -2.;
395  
396 <  double lMVA = fReader->EvaluateMVA( fMethodName );
397 <  
396 >  //set all input variables
397 >  fNVtx      = iVertices->GetEntries();
398 >  fJPt1      = iJet->Pt();
399 >  fJEta1     = iJet->RawMom().Eta();
400 >  fJPhi1     = iJet->RawMom().Phi();
401 >  fJD01      = JetTools::impactParameter(iJet,iVertex);  
402 >  fJDZ1      = JetTools::impactParameter(iJet,iVertex,true);
403 >  fBeta      = JetTools::Beta(iJet,iVertex,fDZCut);
404 >  fBetaStar  = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
405 >  fNCharged  = iJet->ChargedMultiplicity();
406 >  fNNeutrals = iJet->NeutralMultiplicity();
407 >
408 >  fDRMean    = JetTools::dRMean (iJet,-1);
409 >  fDR2Mean   = JetTools::dR2Mean(iJet,-1);
410 >  fFrac01    = JetTools::frac   (iJet,0.1,0. ,-1);
411 >  fFrac02    = JetTools::frac   (iJet,0.2,0.1,-1);
412 >  fFrac03    = JetTools::frac   (iJet,0.3,0.2,-1);
413 >  fFrac04    = JetTools::frac   (iJet,0.4,0.3,-1);
414 >  fFrac05    = JetTools::frac   (iJet,0.5,0.4,-1);
415 >
416 >  double lMVA = 0;
417 >  if(fJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName  );
418 >  if(fJPt1 > 10) lMVA = fReader     ->EvaluateMVA( fHighPtMethodName );
419 >  
420    if (printDebug == kTRUE) {
421      std::cout << "Debug Jet MVA: "
422 <              << fNPV    << " "
423 <              << fJPt1   << " "
424 <              << fJEta1  << " "
425 <              << fJPhi1  << " "
426 <              << fJD01   << " "
427 <              << fJDZ1   << " "
428 <              << fJM1    << " "
429 <              << fNPart1 << " "
430 <              << fLPt1   << " "
431 <              << fLEta1  << " "
432 <              << fLPhi1  << " "
433 <              << fSPt1   << " "
434 <              << fSEta1  << " "
435 <              << fSPhi1  << " "
436 <              << fNEPt1  << " "
437 <              << fNEEta1 << " "
438 <              << fNEPhi1 << " "
261 <              << fEMPt1  << " "
262 <              << fEMEta1 << " "
263 <              << fEMPhi1 << " "
264 <              << fChPt1  << " "
265 <              << fChPhi1 << " "
266 <              << fLFr1   << " "
267 <              << fDRLC1  << " "
268 <              << fDRLS1  << " "
269 <              << fDRM1   << " "
270 <              << fDRNE1 << " "
271 <              << fDREM1  << " "
272 <              << fDRCH1  << " "
422 >              << fNVtx      << " "
423 >              << fJPt1      << " "
424 >              << fJEta1     << " "
425 >              << fJPhi1     << " "
426 >              << fJD01      << " "
427 >              << fJDZ1      << " "
428 >              << fBeta      << " "
429 >              << fBetaStar  << " "
430 >              << fNCharged  << " "
431 >              << fNNeutrals << " "
432 >              << fDRMean    << " "
433 >              << fFrac01    << " "
434 >              << fFrac02    << " "
435 >              << fFrac03    << " "
436 >              << fFrac04    << " "
437 >              << fFrac05    << " "
438 >              << fDRMean    
439                << " === : === "
440                << lMVA << " "    
441                << std::endl;
# Line 277 | Line 443 | Double_t JetIDMVA::MVAValue(const PFJet
443  
444    return lMVA;
445   }
446 < Double_t JetIDMVA::correctedPt(const PFJet *iJet, FactorizedJetCorrector *iJetCorrector,const PileupEnergyDensityCol *iPUEnergyDensity) {
446 > Double_t JetIDMVA::correctedPt(const PFJet *iJet, FactorizedJetCorrector *iJetCorrector,
447 >                               const PileupEnergyDensityCol *iPUEnergyDensity,
448 >                               RhoUtilities::RhoType type) {
449 >  Double_t Rho = 0.0;
450 >  switch(type) {
451 >  case RhoUtilities::MIT_RHO_VORONOI_HIGH_ETA:
452 >    Rho = iPUEnergyDensity->At(0)->Rho();
453 >    break;
454 >  case RhoUtilities::MIT_RHO_VORONOI_LOW_ETA:
455 >    Rho = iPUEnergyDensity->At(0)->RhoLowEta();
456 >    break;
457 >  case RhoUtilities::MIT_RHO_RANDOM_HIGH_ETA:
458 >    Rho = iPUEnergyDensity->At(0)->RhoRandom();
459 >    break;
460 >  case RhoUtilities::MIT_RHO_RANDOM_LOW_ETA:
461 >    Rho = iPUEnergyDensity->At(0)->RhoRandomLowEta();
462 >    break;
463 >  case RhoUtilities::CMS_RHO_RHOKT6PFJETS:
464 >    Rho = iPUEnergyDensity->At(0)->RhoKt6PFJets();
465 >    break;
466 >  default:
467 >    // use the old default
468 >    Rho = iPUEnergyDensity->At(0)->Rho();
469 >    break;
470 >  }
471 >    
472    const FourVectorM rawMom = iJet->RawMom();
473    iJetCorrector->setJetEta(rawMom.Eta());
474    iJetCorrector->setJetPt (rawMom.Pt());
475    iJetCorrector->setJetPhi(rawMom.Phi());
476    iJetCorrector->setJetE  (rawMom.E());
477 <  iJetCorrector->setRho   (iPUEnergyDensity->At(0)->RhoHighEta());
477 >  iJetCorrector->setRho   (Rho);
478    iJetCorrector->setJetA  (iJet->JetArea());
479    iJetCorrector->setJetEMF(-99.0);    
480    Double_t correction = iJetCorrector->getCorrection();

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines