ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MVAMet.cc
Revision: 1.2
Committed: Wed Apr 4 14:43:12 2012 UTC (13 years, 1 month ago) by pharris
Content type: text/plain
Branch: MAIN
Changes since 1.1: +314 -14 lines
Log Message:
Added Regression with corrected jets and updated gbr forest

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