ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MVAMet.cc
Revision: 1.3
Committed: Sat Apr 7 09:36:33 2012 UTC (13 years, 1 month ago) by pharris
Content type: text/plain
Branch: MAIN
Changes since 1.2: +32 -26 lines
Log Message:
Update of Met regresion

File Contents

# User Rev Content
1 pharris 1.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 pharris 1.2 //#include "TMVA/Tools.h"
8     //#include "TMVA/Reader.h"
9 pharris 1.1
10    
11     ClassImp(mithep::MVAMet)
12    
13     using namespace mithep;
14    
15     //--------------------------------------------------------------------------------------------------
16     MVAMet::MVAMet() :
17     fRecoilTools(0),
18 pharris 1.2 fPhiMethodName("PhiCorrection"),
19     fU1MethodName ("U1Correction"),
20 pharris 1.1 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 pharris 1.2 /*
55 pharris 1.1 void MVAMet::setVariables(TMVA::Reader *iReader,bool iScale) {
56     iReader->AddVariable( "npv" , &fNPV );
57     iReader->AddVariable( "u" , &fU );
58     iReader->AddVariable( "uphi" , &fUPhi );
59     iReader->AddVariable( "chsumet/sumet" , &fTKSumEt );
60     iReader->AddVariable( "tku" , &fTKU );
61     iReader->AddVariable( "tkuphi" , &fTKUPhi );
62     iReader->AddVariable( "nopusumet/sumet" , &fNPSumEt );
63     iReader->AddVariable( "nopuu" , &fNPU );
64     iReader->AddVariable( "nopuuphi" , &fNPUPhi );
65     iReader->AddVariable( "pusumet/sumet" , &fPUSumEt );
66     iReader->AddVariable( "pumet" , &fPUMet );
67     iReader->AddVariable( "pumetphi" , &fPUMetPhi);
68     iReader->AddVariable( "pucsumet/sumet" , &fPCSumEt );
69     iReader->AddVariable( "pucu" , &fPCU );
70     iReader->AddVariable( "pucuphi" , &fPCUPhi );
71     iReader->AddVariable( "jspt_1" , &fJSPt1 );
72     iReader->AddVariable( "jseta_1" , &fJSEta1 );
73     iReader->AddVariable( "jsphi_1" , &fJSPhi1 );
74     iReader->AddVariable( "jspt_2" , &fJSPt2 );
75     iReader->AddVariable( "jseta_2" , &fJSEta2 );
76     iReader->AddVariable( "jsphi_2" , &fJSPhi2 );
77     iReader->AddVariable( "nalljet" , &fNAllJet );
78     iReader->AddVariable( "njet" , &fNJet );
79     if(iScale) iReader->AddVariable( "uphi_mva" , &fUPhiMVA );
80     }
81     //--------------------------------------------------------------------------------------------------
82     void MVAMet::Initialize( TString iU1MethodName,
83     TString iPhiMethodName,
84     TString iJetMVAFile,
85     TString iU1Weights,
86     TString iPhiWeights,
87     MVAMet::MVAType iType) {
88    
89     fIsInitialized = kTRUE;
90     fRecoilTools = new RecoilTools(iJetMVAFile);
91 pharris 1.2
92 pharris 1.1 fU1MethodName = iU1MethodName;
93     fPhiMethodName = iPhiMethodName;
94     fType = iType;
95     fPhiReader = new TMVA::Reader( "!Color:!Silent:Error" );
96     fU1Reader = new TMVA::Reader( "!Color:!Silent:Error" );
97     if (fType == kBaseline) {
98     setVariables(fPhiReader,false);
99     setVariables(fU1Reader ,true);
100     }
101     fPhiReader->BookMVA(fPhiMethodName , iPhiWeights );
102     fU1Reader ->BookMVA(fU1MethodName , iU1Weights );
103 pharris 1.2
104 pharris 1.1 std::cout << "Jet ID MVA Initialization\n";
105     std::cout << "Phi Method Name : " << fPhiMethodName << " , type == " << iType << std::endl;
106     std::cout << "U1 Method Name : " << fU1MethodName << " , type == " << iType << std::endl;
107 pharris 1.2
108    
109 pharris 1.1 }
110 pharris 1.2 */
111     //--------------------------------------------------------------------------------------------------
112 pharris 1.3 void MVAMet::Initialize(TString iJetLowPtFile,
113     TString iJetHighPtFile,
114     TString iJetCutFile,
115     TString iU1Weights,
116     TString iPhiWeights,
117     MVAMet::MVAType iType) {
118 pharris 1.2
119     fIsInitialized = kTRUE;
120 pharris 1.3 fRecoilTools = new RecoilTools(iJetLowPtFile,iJetHighPtFile,iJetCutFile);
121    
122 pharris 1.2 fType = iType;
123    
124     TFile *lPhiForest = new TFile(iPhiWeights,"READ");
125     fPhiReader = (GBRForest*)lPhiForest->Get(fPhiMethodName);
126     lPhiForest->Close();
127    
128     TFile *lU1Forest = new TFile(iU1Weights,"READ");
129     fU1Reader = (GBRForest*)lU1Forest->Get(fU1MethodName);
130     lU1Forest->Close();
131 pharris 1.1
132 pharris 1.2 fPhiVals = new Float_t[23];
133     fU1Vals = new Float_t[24];
134     }
135     //--------------------------------------------------------------------------------------------------
136     Double_t MVAMet::evaluatePhi() {
137     fPhiVals[0] = fNPV ;
138     fPhiVals[1] = fU ;
139     fPhiVals[2] = fUPhi ;
140     fPhiVals[3] = fTKSumEt ;
141     fPhiVals[4] = fTKU ;
142     fPhiVals[5] = fTKUPhi ;
143     fPhiVals[6] = fNPSumEt ;
144     fPhiVals[7] = fNPU ;
145     fPhiVals[8] = fNPUPhi ;
146     fPhiVals[9] = fPUSumEt ;
147     fPhiVals[10] = fPUMet ;
148     fPhiVals[11] = fPUMetPhi;
149     fPhiVals[12] = fPCSumEt ;
150     fPhiVals[13] = fPCU ;
151     fPhiVals[14] = fPCUPhi ;
152     fPhiVals[15] = fJSPt1 ;
153     fPhiVals[16] = fJSEta1 ;
154     fPhiVals[17] = fJSPhi1 ;
155     fPhiVals[18] = fJSPt2 ;
156     fPhiVals[19] = fJSEta2 ;
157     fPhiVals[20] = fJSPhi2 ;
158     fPhiVals[21] = fNAllJet ;
159     fPhiVals[22] = fNJet ;
160     return fPhiReader->GetResponse(fPhiVals);
161     }
162     //--------------------------------------------------------------------------------------------------
163     Double_t MVAMet::evaluateU1() {
164     fU1Vals[0] = fNPV ;
165     fU1Vals[1] = fU ;
166     fU1Vals[2] = fUPhi ;
167     fU1Vals[3] = fTKSumEt ;
168     fU1Vals[4] = fTKU ;
169     fU1Vals[5] = fTKUPhi ;
170     fU1Vals[6] = fNPSumEt ;
171     fU1Vals[7] = fNPU ;
172     fU1Vals[8] = fNPUPhi ;
173     fU1Vals[9] = fPUSumEt ;
174     fU1Vals[10] = fPUMet ;
175     fU1Vals[11] = fPUMetPhi;
176     fU1Vals[12] = fPCSumEt ;
177     fU1Vals[13] = fPCU ;
178     fU1Vals[14] = fPCUPhi ;
179     fU1Vals[15] = fJSPt1 ;
180     fU1Vals[16] = fJSEta1 ;
181     fU1Vals[17] = fJSPhi1 ;
182     fU1Vals[18] = fJSPt2 ;
183     fU1Vals[19] = fJSEta2 ;
184     fU1Vals[20] = fJSPhi2 ;
185     fU1Vals[21] = fNAllJet ;
186     fU1Vals[22] = fNJet ;
187     fU1Vals[23] = fUPhiMVA ;
188     return fU1Reader->GetResponse(fU1Vals);
189     }
190 pharris 1.1 //--------------------------------------------------------------------------------------------------
191     Double_t MVAMet::MVAValue( Bool_t iPhi,
192     Float_t iPFSumEt,
193     Float_t iU ,
194     Float_t iUPhi ,
195     Float_t iTKSumEt,
196     Float_t iTKU ,
197     Float_t iTKUPhi ,
198     Float_t iNPSumEt,
199     Float_t iNPU ,
200     Float_t iNPUPhi ,
201     Float_t iPUSumEt,
202     Float_t iPUMet ,
203     Float_t iPUMetPhi,
204     Float_t iPCSumEt,
205     Float_t iPCU ,
206     Float_t iPCUPhi ,
207     Float_t iJSPt1 ,
208     Float_t iJSEta1 ,
209     Float_t iJSPhi1 ,
210     Float_t iJSPt2 ,
211     Float_t iJSEta2 ,
212     Float_t iJSPhi2 ,
213     Float_t iNAllJet,
214     Float_t iNJet ,
215     Float_t iNPV
216     ){
217    
218     if (!fIsInitialized) {
219     std::cout << "Error: MVA Met not properly initialized.\n";
220     return -9999;
221     }
222 pharris 1.2
223 pharris 1.1 fU = iU ;
224     fUPhi = iUPhi ;
225     fTKSumEt = iTKSumEt/iPFSumEt;
226     fTKU = iTKU ;
227     fTKUPhi = iTKUPhi ;
228     fNPSumEt = iNPSumEt/iPFSumEt;
229     fNPU = iNPU ;
230     fNPUPhi = iNPUPhi ;
231     fPUSumEt = iPUSumEt/iPFSumEt;
232     fPUMet = iPUMet ;
233     fPUMetPhi = iPUMetPhi;
234     fPCSumEt = iPCSumEt/iPFSumEt;
235     fPCU = iPCU ;
236     fPCUPhi = iPCUPhi ;
237     fJSPt1 = iJSPt1 ;
238     fJSEta1 = iJSEta1 ;
239     fJSPhi1 = iJSPhi1 ;
240     fJSPt2 = iJSPt2 ;
241     fJSEta2 = iJSEta2 ;
242     fJSPhi2 = iJSPhi2 ;
243     fNAllJet = iNAllJet;
244     fNJet = iNJet ;
245     fNPV = iNPV ;
246    
247     Double_t lMVA = -9999;
248 pharris 1.2 lMVA = evaluatePhi();
249 pharris 1.1 if(!iPhi) fUPhiMVA = iUPhi+lMVA;
250 pharris 1.2 //Not no nice feature of the training
251 pharris 1.1 fTKSumEt /= iPFSumEt;
252     fNPSumEt /= iPFSumEt;
253     fPUSumEt /= iPFSumEt;
254     fPCSumEt /= iPFSumEt;
255 pharris 1.2 if(!iPhi) lMVA = evaluateU1();
256 pharris 1.1 return lMVA;
257     }
258     //--------------------------------------------------------------------------------------------------
259     //====> Please not that the jet collection must be cleaned => all jets near leptons must be removed
260     Met MVAMet::GetMet( Bool_t iPhi,Float_t iPtVis,Float_t iPhiVis,Float_t iSumEtVis,
261     const PFMet *iMet ,
262 pharris 1.3 const PFCandidateCol *iCands,
263     const Vertex *iVertex,const VertexCol *iVertices,
264 pharris 1.1 const PFJetCol *iJets ,
265     FactorizedJetCorrector *iJetCorrector,
266     const PileupEnergyDensityCol *iPUEnergyDensity,
267     int iNPV,
268     Bool_t printDebug) {
269    
270     Met lPFRec = fRecoilTools->pfRecoil (iPtVis,iPhiVis,iSumEtVis,iMet);
271     Met lTKRec = fRecoilTools->trackRecoil(iPtVis,iPhiVis,iSumEtVis,iCands,iVertex);
272 pharris 1.3 Met lNPRec = fRecoilTools->NoPURecoil (iPtVis,iPhiVis,iSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices);
273     Met lPCRec = fRecoilTools->PUCRecoil (iPtVis,iPhiVis,iSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices);
274     Met lPUMet = fRecoilTools->PUMet (iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices);
275 pharris 1.1
276     Double_t lPt0 = 0; const PFJet *lLead = 0;
277     Double_t lPt1 = 0; const PFJet *l2nd = 0;
278     int lNAllJet = 0;
279     int lNJet = 0;
280     for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
281     const PFJet *pJet = iJets->At(i0);
282     Double_t pPt = fRecoilTools->fJetIDMVA->correctedPt(pJet,iJetCorrector,iPUEnergyDensity);
283     if(!JetTools::passPFLooseId(pJet)) continue;
284     lNAllJet++;
285     if(pPt > 30) lNJet++;
286     if(lPt0 < pPt) {lPt0 = pPt; lLead = pJet; continue;}
287     if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
288     }
289    
290     fU = lPFRec.Pt();
291     fUPhi = lPFRec.Phi();
292     fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
293     fTKU = lTKRec.Pt();
294     fTKUPhi = lTKRec.Phi();
295     fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
296     fNPU = lNPRec.Pt();
297     fNPUPhi = lNPRec.Phi();
298     fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
299     fPUMet = lPUMet.Pt();
300     fPUMetPhi = lPUMet.Phi();
301     fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
302     fPCU = lPCRec.Pt() ;
303     fPCUPhi = lPCRec.Phi() ;
304     fJSPt1 = lPt0;
305     fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
306     fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
307     fJSPt2 = lPt1;
308     fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
309     fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
310     fNJet = lNJet ;
311     fNAllJet = lNAllJet;
312     fNPV = iNPV ;
313    
314 pharris 1.2 Float_t lMVA = evaluatePhi();
315 pharris 1.1
316     if(!iPhi) fUPhiMVA = fUPhi + lMVA;
317     //Not no nice feature of teh training
318     fTKSumEt /= lPFRec.SumEt();
319     fNPSumEt /= lPFRec.SumEt();
320     fPUSumEt /= lPFRec.SumEt();
321     fPCSumEt /= lPFRec.SumEt();
322 pharris 1.2 if(!iPhi) lMVA = evaluateU1();//fU1Reader ->EvaluateMVA( fU1MethodName );
323 pharris 1.1
324     TLorentzVector lUVec(0,0,0,0); lUVec .SetPtEtaPhiM(fU*lMVA,0,fUPhiMVA,0);
325     TLorentzVector lVVec(0,0,0,0); lVVec .SetPtEtaPhiM(iPtVis ,0,iPhiVis ,0);
326     if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
327     lUVec -= lVVec;
328     Met lMet(lUVec.Px(),lUVec.Py());
329    
330     if (printDebug == kTRUE) {
331     std::cout << "Debug Jet MVA: "
332     << fU << " : "
333     << fUPhi << " : "
334     << fTKSumEt << " : "
335     << fTKU << " : "
336     << fTKUPhi << " : "
337     << fNPSumEt << " : "
338     << fNPU << " : "
339     << fNPUPhi << " : "
340     << fPUSumEt << " : "
341     << fPUMet << " : "
342     << fPUMetPhi << " : "
343     << fPCUPhi << " : "
344     << fPCSumEt << " : "
345     << fPCU << " : "
346     << fPCUPhi << " : "
347     << fJSPt1 << " : "
348     << fJSEta1 << " : "
349     << fJSPhi1 << " : "
350     << fJSPt2 << " : "
351     << fJSEta2 << " : "
352     << fJSPhi2 << " : "
353     << fNJet << " : "
354     << fNAllJet << " : "
355     << fNPV << " : "
356     << " === : === "
357     << std::endl;
358     }
359    
360     return lMet;
361     }
362 pharris 1.2 //--------------------------------------------------------------------------------------------------
363     //====> Please not that the jet collection must be cleaned => all jets near leptons must be removed
364     //====> Corrected Jets
365     Met MVAMet::GetMet( Bool_t iPhi,Float_t iPtVis,Float_t iPhiVis,Float_t iSumEtVis,
366     const PFMet *iMet ,
367 pharris 1.3 const PFCandidateCol *iCands,
368     const Vertex *iVertex ,const VertexCol *iVertices,
369 pharris 1.2 const PFJetCol *iJets ,
370     int iNPV,
371     Bool_t printDebug) {
372    
373     Met lPFRec = fRecoilTools->pfRecoil (iPtVis,iPhiVis,iSumEtVis,iMet);
374 pharris 1.3 Met lTKRec = fRecoilTools->trackRecoil(iPtVis,iPhiVis,iSumEtVis, iCands,iVertex);
375     Met lNPRec = fRecoilTools->NoPURecoil (iPtVis,iPhiVis,iSumEtVis,iJets,iCands,iVertex,iVertices);
376     Met lPCRec = fRecoilTools->PUCRecoil (iPtVis,iPhiVis,iSumEtVis,iJets,iCands,iVertex,iVertices);
377     Met lPUMet = fRecoilTools->PUMet ( iJets,iCands,iVertex,iVertices);
378 pharris 1.2
379     Double_t lPt0 = 0; const PFJet *lLead = 0;
380     Double_t lPt1 = 0; const PFJet *l2nd = 0;
381     int lNAllJet = 0;
382     int lNJet = 0;
383     for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
384     const PFJet *pJet = iJets->At(i0);
385     Double_t pPt = pJet->Pt();
386     if(!JetTools::passPFLooseId(pJet)) continue;
387     lNAllJet++;
388     if(pPt > 30) lNJet++;
389     if(lPt0 < pPt) {lPt0 = pPt; lLead = pJet; continue;}
390     if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
391     }
392    
393     fU = lPFRec.Pt();
394     fUPhi = lPFRec.Phi();
395     fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
396     fTKU = lTKRec.Pt();
397     fTKUPhi = lTKRec.Phi();
398     fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
399     fNPU = lNPRec.Pt();
400     fNPUPhi = lNPRec.Phi();
401     fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
402     fPUMet = lPUMet.Pt();
403     fPUMetPhi = lPUMet.Phi();
404     fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
405     fPCU = lPCRec.Pt() ;
406     fPCUPhi = lPCRec.Phi() ;
407     fJSPt1 = lPt0;
408     fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
409     fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
410     fJSPt2 = lPt1;
411     fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
412     fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
413     fNJet = lNJet ;
414     fNAllJet = lNAllJet;
415     fNPV = iNPV ;
416    
417     Float_t lMVA = evaluatePhi();
418    
419     if(!iPhi) fUPhiMVA = fUPhi + lMVA;
420     //Not no nice feature of teh training
421     fTKSumEt /= lPFRec.SumEt();
422     fNPSumEt /= lPFRec.SumEt();
423     fPUSumEt /= lPFRec.SumEt();
424     fPCSumEt /= lPFRec.SumEt();
425     if(!iPhi) lMVA = evaluateU1();//fU1Reader ->EvaluateMVA( fU1MethodName );
426    
427     TLorentzVector lUVec(0,0,0,0); lUVec .SetPtEtaPhiM(fU*lMVA,0,fUPhiMVA,0);
428     TLorentzVector lVVec(0,0,0,0); lVVec .SetPtEtaPhiM(iPtVis ,0,iPhiVis ,0);
429     if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
430     lUVec -= lVVec;
431     Met lMet(lUVec.Px(),lUVec.Py());
432    
433     if (printDebug == kTRUE) {
434     std::cout << "Debug Jet MVA: "
435     << fU << " : "
436     << fUPhi << " : "
437     << fTKSumEt << " : "
438     << fTKU << " : "
439     << fTKUPhi << " : "
440     << fNPSumEt << " : "
441     << fNPU << " : "
442     << fNPUPhi << " : "
443     << fPUSumEt << " : "
444     << fPUMet << " : "
445     << fPUMetPhi << " : "
446     << fPCUPhi << " : "
447     << fPCSumEt << " : "
448     << fPCU << " : "
449     << fPCUPhi << " : "
450     << fJSPt1 << " : "
451     << fJSEta1 << " : "
452     << fJSPhi1 << " : "
453     << fJSPt2 << " : "
454     << fJSEta2 << " : "
455     << fJSPhi2 << " : "
456     << fNJet << " : "
457     << fNAllJet << " : "
458     << fNPV << " : "
459     << " === : === "
460     << std::endl;
461     }
462    
463     return lMet;
464     }
465     //--------------------------------------------------------------------------------------------------
466 pharris 1.1 //Interms of the two candidates
467     Met MVAMet::GetMet( Bool_t iPhi,
468     Float_t iPt1,Float_t iPhi1,Float_t iEta1,
469     Float_t iPt2,Float_t iPhi2,Float_t iEta2,
470     const PFMet *iMet ,
471 pharris 1.3 const PFCandidateCol *iCands,
472     const Vertex *iVertex,const VertexCol *iVertices,
473 pharris 1.1 const PFJetCol *iJets ,
474     FactorizedJetCorrector *iJetCorrector,
475     const PileupEnergyDensityCol *iPUEnergyDensity,
476     int iNPV,
477     Bool_t printDebug) {
478    
479     TLorentzVector lVVec1(0,0,0,0); lVVec1.SetPtEtaPhiM(iPt1,0,iPhi1 ,0);
480     TLorentzVector lVVec2(0,0,0,0); lVVec2.SetPtEtaPhiM(iPt2,0,iPhi2 ,0);
481     lVVec1+=lVVec2;
482     Float_t lPtVis = lVVec1.Pt();
483     Float_t lPhiVis = lVVec1.Phi();
484     Float_t lSumEtVis = iPt1 + iPt2;
485     Met lPFRec = fRecoilTools->pfRecoil (lPtVis,lPhiVis,lSumEtVis,iMet);
486     Met lTKRec = fRecoilTools->trackRecoil(lPtVis,lPhiVis,lSumEtVis,iCands,iVertex);
487 pharris 1.3 Met lNPRec = fRecoilTools->NoPURecoil (lPtVis,lPhiVis,lSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices,
488     iPhi1,iEta1,iPhi2,iEta2);
489     Met lPCRec = fRecoilTools->PUCRecoil (lPtVis,lPhiVis,lSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices,
490     iPhi1,iEta1,iPhi2,iEta2);
491     Met lPUMet = fRecoilTools->PUMet ( iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices,
492     iPhi1,iEta1,iPhi2,iEta2);
493 pharris 1.1
494     Double_t lPt0 = 0; const PFJet *lLead = 0;
495     Double_t lPt1 = 0; const PFJet *l2nd = 0;
496     int lNAllJet = 0;
497     int lNJet = 0;
498     for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
499     const PFJet *pJet = iJets->At(i0);
500     Double_t pPt = fRecoilTools->fJetIDMVA->correctedPt(pJet,iJetCorrector,iPUEnergyDensity);
501     double pDEta1 = pJet->Eta() - iEta1;
502     double pDPhi1 = fabs(pJet->Phi() - iPhi1); if(pDPhi1 > 2.*TMath::Pi()-pDPhi1) pDPhi1 = 2.*TMath::Pi()-pDPhi1;
503     double pDR1 = sqrt(pDEta1*pDEta1 + pDPhi1*pDPhi1);
504     if(pDR1 < 0.5) continue;
505     double pDEta2 = pJet->Eta() - iEta2;
506     double pDPhi2 = fabs(pJet->Phi() - iPhi2); if(pDPhi2 > 2.*TMath::Pi()-pDPhi2) pDPhi2 = 2.*TMath::Pi()-pDPhi2;
507     double pDR2 = sqrt(pDEta2*pDEta2 + pDPhi2*pDPhi2);
508     if(pDR2 < 0.5) continue;
509     if(!JetTools::passPFLooseId(pJet)) continue;
510     lNAllJet++;
511     if(pPt > 30.) lNJet++;
512     if(lPt0 < pPt) {lPt0 = pPt; lLead = pJet; continue;}
513     if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
514     }
515    
516     fU = lPFRec.Pt();
517     fUPhi = lPFRec.Phi();
518     fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
519     fTKU = lTKRec.Pt();
520     fTKUPhi = lTKRec.Phi();
521     fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
522     fNPU = lNPRec.Pt();
523     fNPUPhi = lNPRec.Phi();
524     fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
525     fPUMet = lPUMet.Pt();
526     fPUMetPhi= lPUMet.Phi();
527     fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
528     fPCU = lPCRec.Pt() ;
529     fPCUPhi = lPCRec.Phi() ;
530     fJSPt1 = lPt0;
531     fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
532     fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
533     fJSPt2 = lPt1;
534     fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
535     fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
536     fNJet = lNJet ;
537     fNAllJet = lNAllJet;
538     fNPV = iNPV ;
539    
540 pharris 1.2 Float_t lMVA = evaluatePhi();//fPhiReader ->EvaluateMVA( fPhiMethodName );
541    
542     if(!iPhi) fUPhiMVA = fUPhi + lMVA;
543     fTKSumEt /= lPFRec.SumEt();
544     fNPSumEt /= lPFRec.SumEt();
545     fPUSumEt /= lPFRec.SumEt();
546     fPCSumEt /= lPFRec.SumEt();
547     if(!iPhi) lMVA = evaluateU1();//fU1Reader ->EvaluateMVA( fU1MethodName );
548    
549     TLorentzVector lUVec (0,0,0,0); lUVec .SetPtEtaPhiM(fU*lMVA,0,fUPhiMVA,0);
550     TLorentzVector lVVec (0,0,0,0); lVVec .SetPtEtaPhiM(lPtVis ,0,lPhiVis ,0);
551     if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
552     lUVec -= lVVec;
553     Met lMet(lUVec.Px(),lUVec.Py());
554    
555     if (printDebug == kTRUE) {
556     std::cout << "Debug Jet MVA: "
557     << fU << " : "
558     << fUPhi << " : "
559     << fTKSumEt << " : "
560     << fTKU << " : "
561     << fTKUPhi << " : "
562     << fNPSumEt << " : "
563     << fNPU << " : "
564     << fNPUPhi << " : "
565     << fPUSumEt << " : "
566     << fPUMet << " : "
567     << fPUMetPhi << " : "
568     << fPCSumEt << " : "
569     << fPCU << " : "
570     << fPCUPhi << " : "
571     << fJSPt1 << " : "
572     << fJSEta1 << " : "
573     << fJSPhi1 << " : "
574     << fJSPt2 << " : "
575     << fJSEta2 << " : "
576     << fJSPhi2 << " : "
577     << fNJet << " : "
578     << fNAllJet << " : "
579     << fNPV << " : "
580     << " === : === "
581     << std::endl;
582     }
583    
584     return lMet;
585     }
586     //--------------------------------------------------------------------------------------------------
587     //Interms of the two candidates => corrected Jets
588     Met MVAMet::GetMet( Bool_t iPhi,
589     Float_t iPt1,Float_t iPhi1,Float_t iEta1,
590     Float_t iPt2,Float_t iPhi2,Float_t iEta2,
591     const PFMet *iMet ,
592 pharris 1.3 const PFCandidateCol *iCands,
593     const Vertex *iVertex,const VertexCol *iVertices,
594 pharris 1.2 const PFJetCol *iJets ,
595     int iNPV,
596     Bool_t printDebug) {
597    
598     TLorentzVector lVVec1(0,0,0,0); lVVec1.SetPtEtaPhiM(iPt1,0,iPhi1 ,0);
599     TLorentzVector lVVec2(0,0,0,0); lVVec2.SetPtEtaPhiM(iPt2,0,iPhi2 ,0);
600     lVVec1+=lVVec2;
601     Float_t lPtVis = lVVec1.Pt();
602     Float_t lPhiVis = lVVec1.Phi();
603     Float_t lSumEtVis = iPt1 + iPt2;
604     Met lPFRec = fRecoilTools->pfRecoil (lPtVis,lPhiVis,lSumEtVis,iMet);
605     Met lTKRec = fRecoilTools->trackRecoil(lPtVis,lPhiVis,lSumEtVis,iCands,iVertex);
606 pharris 1.3 Met lNPRec = fRecoilTools->NoPURecoil (lPtVis,lPhiVis,lSumEtVis,iJets,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2);
607     Met lPCRec = fRecoilTools->PUCRecoil (lPtVis,lPhiVis,lSumEtVis,iJets,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2);
608     Met lPUMet = fRecoilTools->PUMet ( iJets,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2);
609 pharris 1.2
610     Double_t lPt0 = 0; const PFJet *lLead = 0;
611     Double_t lPt1 = 0; const PFJet *l2nd = 0;
612     int lNAllJet = 0;
613     int lNJet = 0;
614     for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
615     const PFJet *pJet = iJets->At(i0);
616     Double_t pPt = pJet->Pt();
617     if( pJet->TrackCountingHighEffBJetTagsDisc() == -100 && pPt < 10.) continue;
618     double pDEta1 = pJet->Eta() - iEta1;
619     double pDPhi1 = fabs(pJet->Phi() - iPhi1); if(pDPhi1 > 2.*TMath::Pi()-pDPhi1) pDPhi1 = 2.*TMath::Pi()-pDPhi1;
620     double pDR1 = sqrt(pDEta1*pDEta1 + pDPhi1*pDPhi1);
621     if(pDR1 < 0.5) continue;
622     double pDEta2 = pJet->Eta() - iEta2;
623     double pDPhi2 = fabs(pJet->Phi() - iPhi2); if(pDPhi2 > 2.*TMath::Pi()-pDPhi2) pDPhi2 = 2.*TMath::Pi()-pDPhi2;
624     double pDR2 = sqrt(pDEta2*pDEta2 + pDPhi2*pDPhi2);
625     if(pDR2 < 0.5) continue;
626     if(!JetTools::passPFLooseId(pJet)) continue;
627     lNAllJet++;
628     if(pPt > 30.) lNJet++;
629     if(lPt0 < pPt) {lPt0 = pPt; lLead = pJet; continue;}
630     if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
631     }
632    
633     fU = lPFRec.Pt();
634     fUPhi = lPFRec.Phi();
635     fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
636     fTKU = lTKRec.Pt();
637     fTKUPhi = lTKRec.Phi();
638     fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
639     fNPU = lNPRec.Pt();
640     fNPUPhi = lNPRec.Phi();
641     fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
642     fPUMet = lPUMet.Pt();
643     fPUMetPhi= lPUMet.Phi();
644     fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
645     fPCU = lPCRec.Pt() ;
646     fPCUPhi = lPCRec.Phi() ;
647     fJSPt1 = lPt0;
648     fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
649     fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
650     fJSPt2 = lPt1;
651     fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
652     fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
653     fNJet = lNJet ;
654     fNAllJet = lNAllJet;
655     fNPV = iNPV ;
656    
657     Float_t lMVA = evaluatePhi();//fPhiReader ->EvaluateMVA( fPhiMethodName );
658 pharris 1.1
659     if(!iPhi) fUPhiMVA = fUPhi + lMVA;
660     fTKSumEt /= lPFRec.SumEt();
661     fNPSumEt /= lPFRec.SumEt();
662     fPUSumEt /= lPFRec.SumEt();
663     fPCSumEt /= lPFRec.SumEt();
664 pharris 1.2 if(!iPhi) lMVA = evaluateU1();//fU1Reader ->EvaluateMVA( fU1MethodName );
665 pharris 1.1
666     TLorentzVector lUVec (0,0,0,0); lUVec .SetPtEtaPhiM(fU*lMVA,0,fUPhiMVA,0);
667     TLorentzVector lVVec (0,0,0,0); lVVec .SetPtEtaPhiM(lPtVis ,0,lPhiVis ,0);
668     if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
669     lUVec -= lVVec;
670     Met lMet(lUVec.Px(),lUVec.Py());
671    
672     if (printDebug == kTRUE) {
673     std::cout << "Debug Jet MVA: "
674     << fU << " : "
675     << fUPhi << " : "
676     << fTKSumEt << " : "
677     << fTKU << " : "
678     << fTKUPhi << " : "
679     << fNPSumEt << " : "
680     << fNPU << " : "
681     << fNPUPhi << " : "
682     << fPUSumEt << " : "
683     << fPUMet << " : "
684     << fPUMetPhi << " : "
685     << fPCSumEt << " : "
686     << fPCU << " : "
687     << fPCUPhi << " : "
688     << fJSPt1 << " : "
689     << fJSEta1 << " : "
690     << fJSPhi1 << " : "
691     << fJSPt2 << " : "
692     << fJSEta2 << " : "
693     << fJSPhi2 << " : "
694     << fNJet << " : "
695     << fNAllJet << " : "
696     << fNPV << " : "
697     << " === : === "
698     << std::endl;
699     }
700    
701     return lMet;
702     }
703