1 |
#include "MitPhysics/Utils/interface/MVAMet.h"
|
2 |
#include "MitPhysics/Utils/interface/JetTools.h"
|
3 |
#include "MitPhysics/Utils/interface/RecoilTools.h"
|
4 |
#include "MitAna/DataTree/interface/StableData.h"
|
5 |
#include <TFile.h>
|
6 |
#include <TRandom3.h>
|
7 |
#include "TMVA/Tools.h"
|
8 |
#include "TMVA/Reader.h"
|
9 |
|
10 |
|
11 |
ClassImp(mithep::MVAMet)
|
12 |
|
13 |
using namespace mithep;
|
14 |
|
15 |
//--------------------------------------------------------------------------------------------------
|
16 |
MVAMet::MVAMet() :
|
17 |
fRecoilTools(0),
|
18 |
fPhiMethodName("PhiMVAMet"),
|
19 |
fU1MethodName ("U1MVAMet"),
|
20 |
fIsInitialized(kFALSE),
|
21 |
fU (0),
|
22 |
fUPhi (0),
|
23 |
fTKSumEt(0),
|
24 |
fTKU (0),
|
25 |
fTKUPhi (0),
|
26 |
fNPSumEt(0),
|
27 |
fNPU (0),
|
28 |
fNPUPhi (0),
|
29 |
fPUSumEt(0),
|
30 |
fPUMet (0),
|
31 |
fPUMetPhi(0),
|
32 |
fPCSumEt(0),
|
33 |
fPCU (0),
|
34 |
fPCUPhi (0),
|
35 |
fJSPt1 (0),
|
36 |
fJSEta1 (0),
|
37 |
fJSPhi1 (0),
|
38 |
fJSPt2 (0),
|
39 |
fJSEta2 (0),
|
40 |
fJSPhi2 (0),
|
41 |
fNJet (0),
|
42 |
fNAllJet(0),
|
43 |
fNPV (0),
|
44 |
fPhiReader(0),
|
45 |
fU1Reader(0)
|
46 |
{ }
|
47 |
//--------------------------------------------------------------------------------------------------
|
48 |
MVAMet::~MVAMet()
|
49 |
{
|
50 |
fPhiReader = 0;
|
51 |
fU1Reader = 0;
|
52 |
}
|
53 |
//--------------------------------------------------------------------------------------------------
|
54 |
void MVAMet::setVariables(TMVA::Reader *iReader,bool iScale) {
|
55 |
iReader->AddVariable( "npv" , &fNPV );
|
56 |
iReader->AddVariable( "u" , &fU );
|
57 |
iReader->AddVariable( "uphi" , &fUPhi );
|
58 |
iReader->AddVariable( "chsumet/sumet" , &fTKSumEt );
|
59 |
iReader->AddVariable( "tku" , &fTKU );
|
60 |
iReader->AddVariable( "tkuphi" , &fTKUPhi );
|
61 |
iReader->AddVariable( "nopusumet/sumet" , &fNPSumEt );
|
62 |
iReader->AddVariable( "nopuu" , &fNPU );
|
63 |
iReader->AddVariable( "nopuuphi" , &fNPUPhi );
|
64 |
iReader->AddVariable( "pusumet/sumet" , &fPUSumEt );
|
65 |
iReader->AddVariable( "pumet" , &fPUMet );
|
66 |
iReader->AddVariable( "pumetphi" , &fPUMetPhi);
|
67 |
iReader->AddVariable( "pucsumet/sumet" , &fPCSumEt );
|
68 |
iReader->AddVariable( "pucu" , &fPCU );
|
69 |
iReader->AddVariable( "pucuphi" , &fPCUPhi );
|
70 |
iReader->AddVariable( "jspt_1" , &fJSPt1 );
|
71 |
iReader->AddVariable( "jseta_1" , &fJSEta1 );
|
72 |
iReader->AddVariable( "jsphi_1" , &fJSPhi1 );
|
73 |
iReader->AddVariable( "jspt_2" , &fJSPt2 );
|
74 |
iReader->AddVariable( "jseta_2" , &fJSEta2 );
|
75 |
iReader->AddVariable( "jsphi_2" , &fJSPhi2 );
|
76 |
iReader->AddVariable( "nalljet" , &fNAllJet );
|
77 |
iReader->AddVariable( "njet" , &fNJet );
|
78 |
if(iScale) iReader->AddVariable( "uphi_mva" , &fUPhiMVA );
|
79 |
}
|
80 |
//--------------------------------------------------------------------------------------------------
|
81 |
void MVAMet::Initialize( TString iU1MethodName,
|
82 |
TString iPhiMethodName,
|
83 |
TString iJetMVAFile,
|
84 |
TString iU1Weights,
|
85 |
TString iPhiWeights,
|
86 |
MVAMet::MVAType iType) {
|
87 |
|
88 |
fIsInitialized = kTRUE;
|
89 |
fRecoilTools = new RecoilTools(iJetMVAFile);
|
90 |
|
91 |
fU1MethodName = iU1MethodName;
|
92 |
fPhiMethodName = iPhiMethodName;
|
93 |
fType = iType;
|
94 |
fPhiReader = new TMVA::Reader( "!Color:!Silent:Error" );
|
95 |
fU1Reader = new TMVA::Reader( "!Color:!Silent:Error" );
|
96 |
if (fType == kBaseline) {
|
97 |
setVariables(fPhiReader,false);
|
98 |
setVariables(fU1Reader ,true);
|
99 |
}
|
100 |
fPhiReader->BookMVA(fPhiMethodName , iPhiWeights );
|
101 |
fU1Reader ->BookMVA(fU1MethodName , iU1Weights );
|
102 |
std::cout << "Jet ID MVA Initialization\n";
|
103 |
std::cout << "Phi Method Name : " << fPhiMethodName << " , type == " << iType << std::endl;
|
104 |
std::cout << "U1 Method Name : " << fU1MethodName << " , type == " << iType << std::endl;
|
105 |
}
|
106 |
|
107 |
//--------------------------------------------------------------------------------------------------
|
108 |
Double_t MVAMet::MVAValue( Bool_t iPhi,
|
109 |
Float_t iPFSumEt,
|
110 |
Float_t iU ,
|
111 |
Float_t iUPhi ,
|
112 |
Float_t iTKSumEt,
|
113 |
Float_t iTKU ,
|
114 |
Float_t iTKUPhi ,
|
115 |
Float_t iNPSumEt,
|
116 |
Float_t iNPU ,
|
117 |
Float_t iNPUPhi ,
|
118 |
Float_t iPUSumEt,
|
119 |
Float_t iPUMet ,
|
120 |
Float_t iPUMetPhi,
|
121 |
Float_t iPCSumEt,
|
122 |
Float_t iPCU ,
|
123 |
Float_t iPCUPhi ,
|
124 |
Float_t iJSPt1 ,
|
125 |
Float_t iJSEta1 ,
|
126 |
Float_t iJSPhi1 ,
|
127 |
Float_t iJSPt2 ,
|
128 |
Float_t iJSEta2 ,
|
129 |
Float_t iJSPhi2 ,
|
130 |
Float_t iNAllJet,
|
131 |
Float_t iNJet ,
|
132 |
Float_t iNPV
|
133 |
){
|
134 |
|
135 |
if (!fIsInitialized) {
|
136 |
std::cout << "Error: MVA Met not properly initialized.\n";
|
137 |
return -9999;
|
138 |
}
|
139 |
|
140 |
fU = iU ;
|
141 |
fUPhi = iUPhi ;
|
142 |
fTKSumEt = iTKSumEt/iPFSumEt;
|
143 |
fTKU = iTKU ;
|
144 |
fTKUPhi = iTKUPhi ;
|
145 |
fNPSumEt = iNPSumEt/iPFSumEt;
|
146 |
fNPU = iNPU ;
|
147 |
fNPUPhi = iNPUPhi ;
|
148 |
fPUSumEt = iPUSumEt/iPFSumEt;
|
149 |
fPUMet = iPUMet ;
|
150 |
fPUMetPhi = iPUMetPhi;
|
151 |
fPCSumEt = iPCSumEt/iPFSumEt;
|
152 |
fPCU = iPCU ;
|
153 |
fPCUPhi = iPCUPhi ;
|
154 |
fJSPt1 = iJSPt1 ;
|
155 |
fJSEta1 = iJSEta1 ;
|
156 |
fJSPhi1 = iJSPhi1 ;
|
157 |
fJSPt2 = iJSPt2 ;
|
158 |
fJSEta2 = iJSEta2 ;
|
159 |
fJSPhi2 = iJSPhi2 ;
|
160 |
fNAllJet = iNAllJet;
|
161 |
fNJet = iNJet ;
|
162 |
fNPV = iNPV ;
|
163 |
|
164 |
Double_t lMVA = -9999;
|
165 |
lMVA = fPhiReader->EvaluateMVA( fPhiMethodName );
|
166 |
if(!iPhi) fUPhiMVA = iUPhi+lMVA;
|
167 |
//Not no nice feature of teh training
|
168 |
fTKSumEt /= iPFSumEt;
|
169 |
fNPSumEt /= iPFSumEt;
|
170 |
fPUSumEt /= iPFSumEt;
|
171 |
fPCSumEt /= iPFSumEt;
|
172 |
if(!iPhi) lMVA = fU1Reader ->EvaluateMVA( fU1MethodName );
|
173 |
|
174 |
return lMVA;
|
175 |
}
|
176 |
//--------------------------------------------------------------------------------------------------
|
177 |
//====> Please not that the jet collection must be cleaned => all jets near leptons must be removed
|
178 |
Met MVAMet::GetMet( Bool_t iPhi,Float_t iPtVis,Float_t iPhiVis,Float_t iSumEtVis,
|
179 |
const PFMet *iMet ,
|
180 |
const PFCandidateCol *iCands,const Vertex *iVertex,
|
181 |
const PFJetCol *iJets ,
|
182 |
FactorizedJetCorrector *iJetCorrector,
|
183 |
const PileupEnergyDensityCol *iPUEnergyDensity,
|
184 |
int iNPV,
|
185 |
Bool_t printDebug) {
|
186 |
|
187 |
Met lPFRec = fRecoilTools->pfRecoil (iPtVis,iPhiVis,iSumEtVis,iMet);
|
188 |
Met lTKRec = fRecoilTools->trackRecoil(iPtVis,iPhiVis,iSumEtVis,iCands,iVertex);
|
189 |
Met lNPRec = fRecoilTools->NoPURecoil (iPtVis,iPhiVis,iSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex);
|
190 |
Met lPCRec = fRecoilTools->PUCRecoil (iPtVis,iPhiVis,iSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex);
|
191 |
Met lPUMet = fRecoilTools->PUMet (iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex);
|
192 |
|
193 |
Double_t lPt0 = 0; const PFJet *lLead = 0;
|
194 |
Double_t lPt1 = 0; const PFJet *l2nd = 0;
|
195 |
int lNAllJet = 0;
|
196 |
int lNJet = 0;
|
197 |
for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
|
198 |
const PFJet *pJet = iJets->At(i0);
|
199 |
Double_t pPt = fRecoilTools->fJetIDMVA->correctedPt(pJet,iJetCorrector,iPUEnergyDensity);
|
200 |
if( pJet->TrackCountingHighEffBJetTagsDisc() == -100 && pPt < 10.) continue;
|
201 |
if(!JetTools::passPFLooseId(pJet)) continue;
|
202 |
lNAllJet++;
|
203 |
if(pPt > 30) lNJet++;
|
204 |
if(lPt0 < pPt) {lPt0 = pPt; lLead = pJet; continue;}
|
205 |
if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
|
206 |
}
|
207 |
|
208 |
fU = lPFRec.Pt();
|
209 |
fUPhi = lPFRec.Phi();
|
210 |
fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
|
211 |
fTKU = lTKRec.Pt();
|
212 |
fTKUPhi = lTKRec.Phi();
|
213 |
fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
|
214 |
fNPU = lNPRec.Pt();
|
215 |
fNPUPhi = lNPRec.Phi();
|
216 |
fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
|
217 |
fPUMet = lPUMet.Pt();
|
218 |
fPUMetPhi = lPUMet.Phi();
|
219 |
fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
|
220 |
fPCU = lPCRec.Pt() ;
|
221 |
fPCUPhi = lPCRec.Phi() ;
|
222 |
fJSPt1 = lPt0;
|
223 |
fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
|
224 |
fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
|
225 |
fJSPt2 = lPt1;
|
226 |
fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
|
227 |
fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
|
228 |
fNJet = lNJet ;
|
229 |
fNAllJet = lNAllJet;
|
230 |
fNPV = iNPV ;
|
231 |
|
232 |
Float_t lMVA = fPhiReader ->EvaluateMVA( fPhiMethodName );
|
233 |
|
234 |
if(!iPhi) fUPhiMVA = fUPhi + lMVA;
|
235 |
//Not no nice feature of teh training
|
236 |
fTKSumEt /= lPFRec.SumEt();
|
237 |
fNPSumEt /= lPFRec.SumEt();
|
238 |
fPUSumEt /= lPFRec.SumEt();
|
239 |
fPCSumEt /= lPFRec.SumEt();
|
240 |
if(!iPhi) lMVA = fU1Reader ->EvaluateMVA( fU1MethodName );
|
241 |
|
242 |
TLorentzVector lUVec(0,0,0,0); lUVec .SetPtEtaPhiM(fU*lMVA,0,fUPhiMVA,0);
|
243 |
TLorentzVector lVVec(0,0,0,0); lVVec .SetPtEtaPhiM(iPtVis ,0,iPhiVis ,0);
|
244 |
if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
|
245 |
lUVec -= lVVec;
|
246 |
Met lMet(lUVec.Px(),lUVec.Py());
|
247 |
|
248 |
if (printDebug == kTRUE) {
|
249 |
std::cout << "Debug Jet MVA: "
|
250 |
<< fU << " : "
|
251 |
<< fUPhi << " : "
|
252 |
<< fTKSumEt << " : "
|
253 |
<< fTKU << " : "
|
254 |
<< fTKUPhi << " : "
|
255 |
<< fNPSumEt << " : "
|
256 |
<< fNPU << " : "
|
257 |
<< fNPUPhi << " : "
|
258 |
<< fPUSumEt << " : "
|
259 |
<< fPUMet << " : "
|
260 |
<< fPUMetPhi << " : "
|
261 |
<< fPCUPhi << " : "
|
262 |
<< fPCSumEt << " : "
|
263 |
<< fPCU << " : "
|
264 |
<< fPCUPhi << " : "
|
265 |
<< fJSPt1 << " : "
|
266 |
<< fJSEta1 << " : "
|
267 |
<< fJSPhi1 << " : "
|
268 |
<< fJSPt2 << " : "
|
269 |
<< fJSEta2 << " : "
|
270 |
<< fJSPhi2 << " : "
|
271 |
<< fNJet << " : "
|
272 |
<< fNAllJet << " : "
|
273 |
<< fNPV << " : "
|
274 |
<< " === : === "
|
275 |
<< std::endl;
|
276 |
}
|
277 |
|
278 |
return lMet;
|
279 |
}
|
280 |
//Interms of the two candidates
|
281 |
Met MVAMet::GetMet( Bool_t iPhi,
|
282 |
Float_t iPt1,Float_t iPhi1,Float_t iEta1,
|
283 |
Float_t iPt2,Float_t iPhi2,Float_t iEta2,
|
284 |
const PFMet *iMet ,
|
285 |
const PFCandidateCol *iCands,const Vertex *iVertex,
|
286 |
const PFJetCol *iJets ,
|
287 |
FactorizedJetCorrector *iJetCorrector,
|
288 |
const PileupEnergyDensityCol *iPUEnergyDensity,
|
289 |
int iNPV,
|
290 |
Bool_t printDebug) {
|
291 |
|
292 |
TLorentzVector lVVec1(0,0,0,0); lVVec1.SetPtEtaPhiM(iPt1,0,iPhi1 ,0);
|
293 |
TLorentzVector lVVec2(0,0,0,0); lVVec2.SetPtEtaPhiM(iPt2,0,iPhi2 ,0);
|
294 |
lVVec1+=lVVec2;
|
295 |
Float_t lPtVis = lVVec1.Pt();
|
296 |
Float_t lPhiVis = lVVec1.Phi();
|
297 |
Float_t lSumEtVis = iPt1 + iPt2;
|
298 |
Met lPFRec = fRecoilTools->pfRecoil (lPtVis,lPhiVis,lSumEtVis,iMet);
|
299 |
Met lTKRec = fRecoilTools->trackRecoil(lPtVis,lPhiVis,lSumEtVis,iCands,iVertex);
|
300 |
Met lNPRec = fRecoilTools->NoPURecoil (lPtVis,lPhiVis,lSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iPhi1,iEta1,iPhi2,iEta2);
|
301 |
Met lPCRec = fRecoilTools->PUCRecoil (lPtVis,lPhiVis,lSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iPhi1,iEta1,iPhi2,iEta2);
|
302 |
Met lPUMet = fRecoilTools->PUMet ( iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iPhi1,iEta1,iPhi2,iEta2);
|
303 |
|
304 |
Double_t lPt0 = 0; const PFJet *lLead = 0;
|
305 |
Double_t lPt1 = 0; const PFJet *l2nd = 0;
|
306 |
int lNAllJet = 0;
|
307 |
int lNJet = 0;
|
308 |
for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
|
309 |
const PFJet *pJet = iJets->At(i0);
|
310 |
Double_t pPt = fRecoilTools->fJetIDMVA->correctedPt(pJet,iJetCorrector,iPUEnergyDensity);
|
311 |
if( pJet->TrackCountingHighEffBJetTagsDisc() == -100 && pPt < 10.) continue;
|
312 |
double pDEta1 = pJet->Eta() - iEta1;
|
313 |
double pDPhi1 = fabs(pJet->Phi() - iPhi1); if(pDPhi1 > 2.*TMath::Pi()-pDPhi1) pDPhi1 = 2.*TMath::Pi()-pDPhi1;
|
314 |
double pDR1 = sqrt(pDEta1*pDEta1 + pDPhi1*pDPhi1);
|
315 |
if(pDR1 < 0.5) continue;
|
316 |
double pDEta2 = pJet->Eta() - iEta2;
|
317 |
double pDPhi2 = fabs(pJet->Phi() - iPhi2); if(pDPhi2 > 2.*TMath::Pi()-pDPhi2) pDPhi2 = 2.*TMath::Pi()-pDPhi2;
|
318 |
double pDR2 = sqrt(pDEta2*pDEta2 + pDPhi2*pDPhi2);
|
319 |
if(pDR2 < 0.5) continue;
|
320 |
if(!JetTools::passPFLooseId(pJet)) continue;
|
321 |
lNAllJet++;
|
322 |
if(pPt > 30.) lNJet++;
|
323 |
if(lPt0 < pPt) {lPt0 = pPt; lLead = pJet; continue;}
|
324 |
if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
|
325 |
}
|
326 |
|
327 |
fU = lPFRec.Pt();
|
328 |
fUPhi = lPFRec.Phi();
|
329 |
fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
|
330 |
fTKU = lTKRec.Pt();
|
331 |
fTKUPhi = lTKRec.Phi();
|
332 |
fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
|
333 |
fNPU = lNPRec.Pt();
|
334 |
fNPUPhi = lNPRec.Phi();
|
335 |
fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
|
336 |
fPUMet = lPUMet.Pt();
|
337 |
fPUMetPhi= lPUMet.Phi();
|
338 |
fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
|
339 |
fPCU = lPCRec.Pt() ;
|
340 |
fPCUPhi = lPCRec.Phi() ;
|
341 |
fJSPt1 = lPt0;
|
342 |
fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
|
343 |
fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
|
344 |
fJSPt2 = lPt1;
|
345 |
fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
|
346 |
fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
|
347 |
fNJet = lNJet ;
|
348 |
fNAllJet = lNAllJet;
|
349 |
fNPV = iNPV ;
|
350 |
|
351 |
Float_t lMVA = fPhiReader ->EvaluateMVA( fPhiMethodName );
|
352 |
|
353 |
if(!iPhi) fUPhiMVA = fUPhi + lMVA;
|
354 |
fTKSumEt /= lPFRec.SumEt();
|
355 |
fNPSumEt /= lPFRec.SumEt();
|
356 |
fPUSumEt /= lPFRec.SumEt();
|
357 |
fPCSumEt /= lPFRec.SumEt();
|
358 |
if(!iPhi) lMVA = fU1Reader ->EvaluateMVA( fU1MethodName );
|
359 |
|
360 |
TLorentzVector lUVec (0,0,0,0); lUVec .SetPtEtaPhiM(fU*lMVA,0,fUPhiMVA,0);
|
361 |
TLorentzVector lVVec (0,0,0,0); lVVec .SetPtEtaPhiM(lPtVis ,0,lPhiVis ,0);
|
362 |
if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
|
363 |
lUVec -= lVVec;
|
364 |
Met lMet(lUVec.Px(),lUVec.Py());
|
365 |
|
366 |
if (printDebug == kTRUE) {
|
367 |
std::cout << "Debug Jet MVA: "
|
368 |
<< fU << " : "
|
369 |
<< fUPhi << " : "
|
370 |
<< fTKSumEt << " : "
|
371 |
<< fTKU << " : "
|
372 |
<< fTKUPhi << " : "
|
373 |
<< fNPSumEt << " : "
|
374 |
<< fNPU << " : "
|
375 |
<< fNPUPhi << " : "
|
376 |
<< fPUSumEt << " : "
|
377 |
<< fPUMet << " : "
|
378 |
<< fPUMetPhi << " : "
|
379 |
<< fPCSumEt << " : "
|
380 |
<< fPCU << " : "
|
381 |
<< fPCUPhi << " : "
|
382 |
<< fJSPt1 << " : "
|
383 |
<< fJSEta1 << " : "
|
384 |
<< fJSPhi1 << " : "
|
385 |
<< fJSPt2 << " : "
|
386 |
<< fJSEta2 << " : "
|
387 |
<< fJSPhi2 << " : "
|
388 |
<< fNJet << " : "
|
389 |
<< fNAllJet << " : "
|
390 |
<< fNPV << " : "
|
391 |
<< " === : === "
|
392 |
<< std::endl;
|
393 |
}
|
394 |
|
395 |
return lMet;
|
396 |
}
|
397 |
|