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.4 |
fJetPtMin(0) , //We need to lower this
|
19 |
|
|
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 |
|
|
fFrac05 (0)
|
39 |
pharris |
1.1 |
{
|
40 |
|
|
fReader = 0;
|
41 |
|
|
}
|
42 |
|
|
//--------------------------------------------------------------------------------------------------
|
43 |
pharris |
1.4 |
JetIDMVA::~JetIDMVA() {
|
44 |
|
|
|
45 |
pharris |
1.1 |
fReader = 0;
|
46 |
|
|
}
|
47 |
|
|
|
48 |
|
|
//--------------------------------------------------------------------------------------------------
|
49 |
pharris |
1.4 |
void JetIDMVA::Initialize( JetIDMVA::CutType iCutType,
|
50 |
|
|
TString iLowPtWeights,
|
51 |
|
|
TString iHighPtWeights,
|
52 |
|
|
JetIDMVA::MVAType iType,
|
53 |
|
|
TString iCutFileName)
|
54 |
pharris |
1.1 |
{
|
55 |
|
|
|
56 |
|
|
fIsInitialized = kTRUE;
|
57 |
|
|
fType = iType;
|
58 |
pharris |
1.4 |
fCutType = iCutType;
|
59 |
pharris |
1.1 |
fReader = 0;
|
60 |
|
|
fReader = new TMVA::Reader( "!Color:!Silent:Error" );
|
61 |
|
|
if (fType == kBaseline) {
|
62 |
pharris |
1.4 |
fReader->AddVariable( "nvtx" , &fNVtx );
|
63 |
|
|
fReader->AddVariable( "jetPt" , &fJPt1 );
|
64 |
|
|
fReader->AddVariable( "jetEta" , &fJEta1 );
|
65 |
|
|
fReader->AddVariable( "jetPhi" , &fJPhi1 );
|
66 |
|
|
fReader->AddVariable( "dZ" , &fJDZ1 );
|
67 |
|
|
fReader->AddVariable( "d0" , &fJD01 );
|
68 |
|
|
fReader->AddVariable( "beta" , &fBeta );
|
69 |
|
|
fReader->AddVariable( "betaStar" , &fBetaStar );
|
70 |
|
|
fReader->AddVariable( "nCharged" , &fNCharged );
|
71 |
|
|
fReader->AddVariable( "nNeutrals", &fNNeutrals );
|
72 |
|
|
fReader->AddVariable( "dRMean" , &fDRMean );
|
73 |
|
|
fReader->AddVariable( "frac01" , &fFrac01 );
|
74 |
|
|
fReader->AddVariable( "frac02" , &fFrac02 );
|
75 |
|
|
fReader->AddVariable( "frac03" , &fFrac03 );
|
76 |
|
|
fReader->AddVariable( "frac04" , &fFrac04 );
|
77 |
|
|
fReader->AddVariable( "frac05" , &fFrac05 );
|
78 |
pharris |
1.1 |
}
|
79 |
pharris |
1.4 |
fReader->BookMVA(fLowPtMethodName , iLowPtWeights );
|
80 |
|
|
fReader->BookMVA(fHighPtMethodName , iHighPtWeights );
|
81 |
pharris |
1.1 |
std::cout << "Jet ID MVA Initialization\n";
|
82 |
pharris |
1.4 |
std::cout << "MethodName : " << fLowPtMethodName << " , type == " << fType << std::endl;
|
83 |
|
|
|
84 |
|
|
//Load Cut Matrix
|
85 |
|
|
edm::ParameterSet lConfig = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>("JetIdParams");
|
86 |
|
|
std::string lCutType = "Tight";
|
87 |
|
|
if(fCutType == kMedium) lCutType = "Medium";
|
88 |
|
|
if(fCutType == kLoose ) lCutType = "Loose";
|
89 |
|
|
std::vector<double> lPt010 = lConfig.getParameter<std::vector<double> >(("Pt010_" +lCutType).c_str());
|
90 |
|
|
std::vector<double> lPt1020 = lConfig.getParameter<std::vector<double> >(("Pt1020_"+lCutType).c_str());
|
91 |
|
|
std::vector<double> lPt2030 = lConfig.getParameter<std::vector<double> >(("Pt2030_"+lCutType).c_str());
|
92 |
|
|
std::vector<double> lPt3050 = lConfig.getParameter<std::vector<double> >(("Pt3050_"+lCutType).c_str());
|
93 |
|
|
for(int i0 = 0; i0 < 4; i0++) fMVACut[0][i0] = lPt010 [i0];
|
94 |
|
|
for(int i0 = 0; i0 < 4; i0++) fMVACut[1][i0] = lPt1020[i0];
|
95 |
|
|
for(int i0 = 0; i0 < 4; i0++) fMVACut[2][i0] = lPt2030[i0];
|
96 |
|
|
for(int i0 = 0; i0 < 4; i0++) fMVACut[3][i0] = lPt3050[i0];
|
97 |
pharris |
1.1 |
}
|
98 |
|
|
|
99 |
|
|
//--------------------------------------------------------------------------------------------------
|
100 |
|
|
Double_t JetIDMVA::MVAValue(
|
101 |
|
|
Float_t iNPV ,
|
102 |
|
|
Float_t iJPt1 ,
|
103 |
|
|
Float_t iJEta1 ,
|
104 |
|
|
Float_t iJPhi1 ,
|
105 |
|
|
Float_t iJD01 ,
|
106 |
|
|
Float_t iJDZ1 ,
|
107 |
pharris |
1.4 |
Float_t iBeta ,
|
108 |
|
|
Float_t iBetaStar,
|
109 |
|
|
Float_t iNCharged,
|
110 |
|
|
Float_t iNNeutrals,
|
111 |
|
|
Float_t iDRMean ,
|
112 |
|
|
Float_t iFrac01 ,
|
113 |
|
|
Float_t iFrac02 ,
|
114 |
|
|
Float_t iFrac03 ,
|
115 |
|
|
Float_t iFrac04 ,
|
116 |
|
|
Float_t iFrac05
|
117 |
pharris |
1.1 |
){
|
118 |
|
|
|
119 |
pharris |
1.4 |
if(!fIsInitialized) {
|
120 |
pharris |
1.1 |
std::cout << "Error: JetIDMVA not properly initialized.\n";
|
121 |
|
|
return -9999;
|
122 |
|
|
}
|
123 |
|
|
|
124 |
pharris |
1.4 |
fNVtx = iNPV;
|
125 |
|
|
fJPt1 = iJPt1;
|
126 |
|
|
fJEta1 = iJEta1;
|
127 |
|
|
fJPhi1 = fJPhi1;
|
128 |
|
|
fJD01 = iJD01;
|
129 |
|
|
fJDZ1 = iJDZ1;
|
130 |
|
|
fBeta = iBeta;
|
131 |
|
|
fBetaStar = iBetaStar;
|
132 |
|
|
fNCharged = iNCharged;
|
133 |
|
|
fNNeutrals = iNNeutrals;
|
134 |
|
|
fDRMean = iDRMean;
|
135 |
|
|
fFrac01 = iFrac01;
|
136 |
|
|
fFrac02 = iFrac02;
|
137 |
|
|
fFrac03 = iFrac03;
|
138 |
|
|
fFrac04 = iFrac04;
|
139 |
|
|
fFrac05 = iFrac05;
|
140 |
pharris |
1.1 |
|
141 |
|
|
Double_t lMVA = -9999;
|
142 |
pharris |
1.4 |
if(iJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName );
|
143 |
|
|
if(iJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
|
144 |
pharris |
1.1 |
|
145 |
|
|
return lMVA;
|
146 |
|
|
}
|
147 |
|
|
//--------------------------------------------------------------------------------------------------
|
148 |
pharris |
1.4 |
Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices,
|
149 |
pharris |
1.1 |
FactorizedJetCorrector *iJetCorrector,
|
150 |
|
|
const PileupEnergyDensityCol *iPileupEnergyDensity) {
|
151 |
pharris |
1.4 |
|
152 |
pharris |
1.1 |
if(!JetTools::passPFLooseId(iJet)) return false;
|
153 |
pharris |
1.4 |
if(iJet->Pt() < fJetPtMin) return false;
|
154 |
|
|
if(iJet->Pt() > 50) return true; //==> we can raise this
|
155 |
|
|
if(fabs(iJet->Eta()) > 4.99) return false;//Castor
|
156 |
|
|
|
157 |
|
|
double lMVA = MVAValue (iJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity);
|
158 |
|
|
double lPt = correctedPt(iJet, iJetCorrector,iPileupEnergyDensity);
|
159 |
|
|
|
160 |
|
|
int lPtId = 0;
|
161 |
|
|
if(lPt > 10 && iJet->Pt() < 20) lPtId = 1;
|
162 |
|
|
if(lPt > 20 && iJet->Pt() < 30) lPtId = 2;
|
163 |
|
|
if(lPt > 30 ) lPtId = 3;
|
164 |
|
|
|
165 |
|
|
int lEtaId = 0;
|
166 |
|
|
if(fabs(iJet->Eta()) > 2.5 && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
|
167 |
|
|
if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
|
168 |
|
|
if(fabs(iJet->Eta()) > 3.0 && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
|
169 |
|
|
|
170 |
|
|
double lMVACut = fMVACut[lPtId][lEtaId];
|
171 |
|
|
if(lMVA < lMVACut) return false;
|
172 |
pharris |
1.1 |
return true;
|
173 |
pharris |
1.4 |
//if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
|
174 |
|
|
//if(correctedPt(iJet,iJetCorrector,iPileupEnergyDensity) < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false;
|
175 |
|
|
//This line is a bug in the Met training
|
176 |
|
|
//if(lMVA < -0.8) return false;
|
177 |
|
|
//if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
|
178 |
pharris |
1.1 |
}
|
179 |
|
|
//--------------------------------------------------------------------------------------------------
|
180 |
pharris |
1.4 |
Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices) {
|
181 |
pharris |
1.2 |
if(!JetTools::passPFLooseId(iJet)) return false;
|
182 |
pharris |
1.4 |
if(iJet->Pt() < fJetPtMin) return false;
|
183 |
|
|
if(iJet->Pt() > 50) return true; //==> we can raise this
|
184 |
|
|
if(fabs(iJet->Eta()) > 4.99) return false;//Castor
|
185 |
|
|
|
186 |
|
|
double lMVA = MVAValue(iJet,iVertex,iVertices);
|
187 |
|
|
|
188 |
|
|
int lPtId = 0;
|
189 |
|
|
if(iJet->Pt() > 10 && iJet->Pt() < 20) lPtId = 1;
|
190 |
|
|
if(iJet->Pt() > 20 && iJet->Pt() < 30) lPtId = 2;
|
191 |
|
|
if(iJet->Pt() > 30 ) lPtId = 3;
|
192 |
|
|
|
193 |
|
|
int lEtaId = 0;
|
194 |
|
|
if(fabs(iJet->Eta()) > 2.5 && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
|
195 |
|
|
if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
|
196 |
|
|
if(fabs(iJet->Eta()) > 3.0 && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
|
197 |
|
|
|
198 |
|
|
double lMVACut = fMVACut[lPtId][lEtaId];
|
199 |
|
|
if(lMVA < lMVACut) return false;
|
200 |
pharris |
1.2 |
return true;
|
201 |
pharris |
1.4 |
//if(lMVA < -0.8) return false;
|
202 |
|
|
//if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
|
203 |
|
|
//if(iJet->Pt() < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false; //This line is a bug in the Met training
|
204 |
|
|
//if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
|
205 |
pharris |
1.2 |
}
|
206 |
|
|
//--------------------------------------------------------------------------------------------------
|
207 |
pharris |
1.4 |
Double_t JetIDMVA::MVAValue(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices, //Vertex here is the PV
|
208 |
pharris |
1.1 |
FactorizedJetCorrector *iJetCorrector,
|
209 |
|
|
const PileupEnergyDensityCol *iPileupEnergyDensity,
|
210 |
|
|
Bool_t printDebug) {
|
211 |
|
|
|
212 |
|
|
if (!fIsInitialized) {
|
213 |
|
|
std::cout << "Error: JetIDMVA not properly initialized.\n";
|
214 |
|
|
return -9999;
|
215 |
|
|
}
|
216 |
|
|
if(!JetTools::passPFLooseId(iJet)) return -2.;
|
217 |
|
|
|
218 |
|
|
//set all input variables
|
219 |
pharris |
1.4 |
fNVtx = iVertices->GetEntries();
|
220 |
pharris |
1.1 |
fJPt1 = correctedPt(iJet,iJetCorrector,iPileupEnergyDensity);
|
221 |
|
|
fJEta1 = iJet->RawMom().Eta();
|
222 |
|
|
fJPhi1 = iJet->RawMom().Phi();
|
223 |
pharris |
1.4 |
fJD01 = JetTools::impactParameter(iJet,iVertex);
|
224 |
|
|
fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
|
225 |
|
|
fBeta = JetTools::Beta(iJet,iVertex,fDZCut);
|
226 |
|
|
fBetaStar = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
|
227 |
|
|
fNCharged = iJet->ChargedMultiplicity();
|
228 |
|
|
fNNeutrals = iJet->NeutralMultiplicity();
|
229 |
|
|
|
230 |
|
|
fDRMean = JetTools::dRMean(iJet,-1);
|
231 |
|
|
fFrac01 = JetTools::frac (iJet,0.1,0. ,-1);
|
232 |
|
|
fFrac02 = JetTools::frac (iJet,0.2,0.1,-1);
|
233 |
|
|
fFrac03 = JetTools::frac (iJet,0.3,0.2,-1);
|
234 |
|
|
fFrac04 = JetTools::frac (iJet,0.4,0.3,-1);
|
235 |
|
|
fFrac05 = JetTools::frac (iJet,0.5,0.4,-1);
|
236 |
|
|
|
237 |
|
|
double lMVA = 0;
|
238 |
|
|
if(fJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName );
|
239 |
|
|
if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
|
240 |
pharris |
1.1 |
|
241 |
|
|
if (printDebug == kTRUE) {
|
242 |
|
|
std::cout << "Debug Jet MVA: "
|
243 |
pharris |
1.4 |
<< fNVtx << " "
|
244 |
|
|
<< fJPt1 << " "
|
245 |
|
|
<< fJEta1 << " "
|
246 |
|
|
<< fJPhi1 << " "
|
247 |
|
|
<< fJD01 << " "
|
248 |
|
|
<< fJDZ1 << " "
|
249 |
|
|
<< fBeta << " "
|
250 |
|
|
<< fBetaStar << " "
|
251 |
|
|
<< fNCharged << " "
|
252 |
|
|
<< fNNeutrals << " "
|
253 |
|
|
<< fDRMean << " "
|
254 |
|
|
<< fFrac01 << " "
|
255 |
|
|
<< fFrac02 << " "
|
256 |
|
|
<< fFrac03 << " "
|
257 |
|
|
<< fFrac04 << " "
|
258 |
|
|
<< fFrac05
|
259 |
pharris |
1.1 |
<< " === : === "
|
260 |
|
|
<< lMVA << " "
|
261 |
|
|
<< std::endl;
|
262 |
|
|
}
|
263 |
|
|
|
264 |
|
|
return lMVA;
|
265 |
|
|
}
|
266 |
pharris |
1.4 |
Double_t JetIDMVA::MVAValue(const PFJet *iJet,const Vertex *iVertex, const VertexCol *iVertices,//Vertex here is the PV
|
267 |
pharris |
1.2 |
Bool_t printDebug) {
|
268 |
|
|
|
269 |
|
|
if (!fIsInitialized) {
|
270 |
|
|
std::cout << "Error: JetIDMVA not properly initialized.\n";
|
271 |
|
|
return -9999;
|
272 |
|
|
}
|
273 |
|
|
if(!JetTools::passPFLooseId(iJet)) return -2.;
|
274 |
|
|
|
275 |
|
|
//set all input variables
|
276 |
pharris |
1.4 |
fNVtx = iVertices->GetEntries();
|
277 |
pharris |
1.2 |
fJPt1 = iJet->Pt();
|
278 |
|
|
fJEta1 = iJet->RawMom().Eta();
|
279 |
|
|
fJPhi1 = iJet->RawMom().Phi();
|
280 |
pharris |
1.4 |
fJD01 = JetTools::impactParameter(iJet,iVertex);
|
281 |
|
|
fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
|
282 |
|
|
fBeta = JetTools::Beta(iJet,iVertex,fDZCut);
|
283 |
|
|
fBetaStar = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
|
284 |
|
|
fNCharged = iJet->ChargedMultiplicity();
|
285 |
|
|
fNNeutrals = iJet->NeutralMultiplicity();
|
286 |
|
|
|
287 |
|
|
fDRMean = JetTools::dRMean(iJet,-1);
|
288 |
|
|
fFrac01 = JetTools::frac (iJet,0.1,0. ,-1);
|
289 |
|
|
fFrac02 = JetTools::frac (iJet,0.2,0.1,-1);
|
290 |
|
|
fFrac03 = JetTools::frac (iJet,0.3,0.2,-1);
|
291 |
|
|
fFrac04 = JetTools::frac (iJet,0.4,0.3,-1);
|
292 |
|
|
fFrac05 = JetTools::frac (iJet,0.5,0.4,-1);
|
293 |
|
|
|
294 |
|
|
double lMVA = 0;
|
295 |
|
|
if(fJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName );
|
296 |
|
|
if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
|
297 |
pharris |
1.2 |
|
298 |
|
|
if (printDebug == kTRUE) {
|
299 |
|
|
std::cout << "Debug Jet MVA: "
|
300 |
pharris |
1.4 |
<< fNVtx << " "
|
301 |
|
|
<< fJPt1 << " "
|
302 |
|
|
<< fJEta1 << " "
|
303 |
|
|
<< fJPhi1 << " "
|
304 |
|
|
<< fJD01 << " "
|
305 |
|
|
<< fJDZ1 << " "
|
306 |
|
|
<< fBeta << " "
|
307 |
|
|
<< fBetaStar << " "
|
308 |
|
|
<< fNCharged << " "
|
309 |
|
|
<< fNNeutrals << " "
|
310 |
|
|
<< fDRMean << " "
|
311 |
|
|
<< fFrac01 << " "
|
312 |
|
|
<< fFrac02 << " "
|
313 |
|
|
<< fFrac03 << " "
|
314 |
|
|
<< fFrac04 << " "
|
315 |
|
|
<< fFrac05
|
316 |
pharris |
1.2 |
<< " === : === "
|
317 |
|
|
<< lMVA << " "
|
318 |
|
|
<< std::endl;
|
319 |
|
|
}
|
320 |
|
|
|
321 |
|
|
return lMVA;
|
322 |
|
|
}
|
323 |
pharris |
1.1 |
Double_t JetIDMVA::correctedPt(const PFJet *iJet, FactorizedJetCorrector *iJetCorrector,const PileupEnergyDensityCol *iPUEnergyDensity) {
|
324 |
|
|
const FourVectorM rawMom = iJet->RawMom();
|
325 |
|
|
iJetCorrector->setJetEta(rawMom.Eta());
|
326 |
|
|
iJetCorrector->setJetPt (rawMom.Pt());
|
327 |
|
|
iJetCorrector->setJetPhi(rawMom.Phi());
|
328 |
|
|
iJetCorrector->setJetE (rawMom.E());
|
329 |
|
|
iJetCorrector->setRho (iPUEnergyDensity->At(0)->RhoHighEta());
|
330 |
|
|
iJetCorrector->setJetA (iJet->JetArea());
|
331 |
|
|
iJetCorrector->setJetEMF(-99.0);
|
332 |
|
|
Double_t correction = iJetCorrector->getCorrection();
|
333 |
|
|
return rawMom.Pt()*correction;
|
334 |
|
|
}
|