ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MVAMet.cc
Revision: 1.1
Committed: Wed Mar 21 18:56:26 2012 UTC (13 years, 1 month ago) by pharris
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_025e, Mit_025d
Log Message:
Adding MET regression

File Contents

# Content
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