ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MVAMet.cc
Revision: 1.13
Committed: Mon Apr 30 17:39:37 2012 UTC (13 years ago) by pharris
Content type: text/plain
Branch: MAIN
Changes since 1.12: +0 -1 lines
Log Message:
Fixed MVA Met bug

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 bendavid 1.6 #include "CondFormats/EgammaObjects/interface/GBRForest.h"
8     #include "Cintex/Cintex.h"
9 pharris 1.1
10     ClassImp(mithep::MVAMet)
11    
12     using namespace mithep;
13    
14     //--------------------------------------------------------------------------------------------------
15     MVAMet::MVAMet() :
16     fRecoilTools(0),
17 pharris 1.4 fPhiMethodName ("PhiCorrection"),
18     fU1MethodName ("U1Correction"),
19     fCovU1MethodName ("CovU1"),
20     fCovU2MethodName ("CovU2"),
21 pharris 1.1 fIsInitialized(kFALSE),
22     fU (0),
23     fUPhi (0),
24     fTKSumEt(0),
25     fTKU (0),
26     fTKUPhi (0),
27     fNPSumEt(0),
28     fNPU (0),
29     fNPUPhi (0),
30     fPUSumEt(0),
31     fPUMet (0),
32     fPUMetPhi(0),
33     fPCSumEt(0),
34     fPCU (0),
35     fPCUPhi (0),
36     fJSPt1 (0),
37     fJSEta1 (0),
38     fJSPhi1 (0),
39     fJSPt2 (0),
40     fJSEta2 (0),
41     fJSPhi2 (0),
42     fNJet (0),
43     fNAllJet(0),
44     fNPV (0),
45     fPhiReader(0),
46 pharris 1.4 fU1Reader(0),
47     fCovU1Reader(0),
48     fCovU2Reader(0) { }
49 pharris 1.1 //--------------------------------------------------------------------------------------------------
50     MVAMet::~MVAMet()
51     {
52 pharris 1.4 delete fPhiReader;
53     delete fU1Reader;
54     delete fCovU1Reader;
55     delete fCovU2Reader;
56    
57 pharris 1.1 }
58     //--------------------------------------------------------------------------------------------------
59 pharris 1.2 /*
60 pharris 1.1 void MVAMet::setVariables(TMVA::Reader *iReader,bool iScale) {
61     iReader->AddVariable( "npv" , &fNPV );
62     iReader->AddVariable( "u" , &fU );
63     iReader->AddVariable( "uphi" , &fUPhi );
64     iReader->AddVariable( "chsumet/sumet" , &fTKSumEt );
65     iReader->AddVariable( "tku" , &fTKU );
66     iReader->AddVariable( "tkuphi" , &fTKUPhi );
67     iReader->AddVariable( "nopusumet/sumet" , &fNPSumEt );
68     iReader->AddVariable( "nopuu" , &fNPU );
69     iReader->AddVariable( "nopuuphi" , &fNPUPhi );
70     iReader->AddVariable( "pusumet/sumet" , &fPUSumEt );
71     iReader->AddVariable( "pumet" , &fPUMet );
72     iReader->AddVariable( "pumetphi" , &fPUMetPhi);
73     iReader->AddVariable( "pucsumet/sumet" , &fPCSumEt );
74     iReader->AddVariable( "pucu" , &fPCU );
75     iReader->AddVariable( "pucuphi" , &fPCUPhi );
76     iReader->AddVariable( "jspt_1" , &fJSPt1 );
77     iReader->AddVariable( "jseta_1" , &fJSEta1 );
78     iReader->AddVariable( "jsphi_1" , &fJSPhi1 );
79     iReader->AddVariable( "jspt_2" , &fJSPt2 );
80     iReader->AddVariable( "jseta_2" , &fJSEta2 );
81     iReader->AddVariable( "jsphi_2" , &fJSPhi2 );
82     iReader->AddVariable( "nalljet" , &fNAllJet );
83     iReader->AddVariable( "njet" , &fNJet );
84     if(iScale) iReader->AddVariable( "uphi_mva" , &fUPhiMVA );
85     }
86     //--------------------------------------------------------------------------------------------------
87     void MVAMet::Initialize( TString iU1MethodName,
88     TString iPhiMethodName,
89     TString iJetMVAFile,
90     TString iU1Weights,
91     TString iPhiWeights,
92     MVAMet::MVAType iType) {
93    
94     fIsInitialized = kTRUE;
95     fRecoilTools = new RecoilTools(iJetMVAFile);
96 pharris 1.2
97 pharris 1.1 fU1MethodName = iU1MethodName;
98     fPhiMethodName = iPhiMethodName;
99     fType = iType;
100     fPhiReader = new TMVA::Reader( "!Color:!Silent:Error" );
101     fU1Reader = new TMVA::Reader( "!Color:!Silent:Error" );
102     if (fType == kBaseline) {
103     setVariables(fPhiReader,false);
104     setVariables(fU1Reader ,true);
105     }
106     fPhiReader->BookMVA(fPhiMethodName , iPhiWeights );
107     fU1Reader ->BookMVA(fU1MethodName , iU1Weights );
108 pharris 1.2
109 pharris 1.1 std::cout << "Jet ID MVA Initialization\n";
110     std::cout << "Phi Method Name : " << fPhiMethodName << " , type == " << iType << std::endl;
111     std::cout << "U1 Method Name : " << fU1MethodName << " , type == " << iType << std::endl;
112 pharris 1.2
113    
114 pharris 1.1 }
115 pharris 1.2 */
116     //--------------------------------------------------------------------------------------------------
117 pharris 1.3 void MVAMet::Initialize(TString iJetLowPtFile,
118     TString iJetHighPtFile,
119     TString iJetCutFile,
120     TString iU1Weights,
121     TString iPhiWeights,
122 pharris 1.4 TString iCovU1Weights,
123     TString iCovU2Weights,
124 pharris 1.3 MVAMet::MVAType iType) {
125 pharris 1.2
126     fIsInitialized = kTRUE;
127 pharris 1.3 fRecoilTools = new RecoilTools(iJetLowPtFile,iJetHighPtFile,iJetCutFile);
128    
129 pharris 1.2 fType = iType;
130    
131 bendavid 1.6 ROOT::Cintex::Cintex::Enable();
132    
133 pharris 1.2 TFile *lPhiForest = new TFile(iPhiWeights,"READ");
134     fPhiReader = (GBRForest*)lPhiForest->Get(fPhiMethodName);
135     lPhiForest->Close();
136    
137     TFile *lU1Forest = new TFile(iU1Weights,"READ");
138     fU1Reader = (GBRForest*)lU1Forest->Get(fU1MethodName);
139     lU1Forest->Close();
140 pharris 1.1
141 pharris 1.4 TFile *lCovU1Forest = new TFile(iCovU1Weights,"READ");
142     fCovU1Reader = (GBRForest*)lCovU1Forest->Get(fCovU1MethodName);
143     lCovU1Forest->Close();
144    
145     TFile *lCovU2Forest = new TFile(iCovU2Weights,"READ");
146     fCovU2Reader = (GBRForest*)lCovU2Forest->Get(fCovU2MethodName);
147     lCovU2Forest->Close();
148 pharris 1.5
149     fCov = new TMatrixD(2,2);
150 pharris 1.7 fPhiVals = new Float_t[23];
151     fU1Vals = new Float_t[25];
152     fCovVals = new Float_t[26];
153 pharris 1.2 }
154     //--------------------------------------------------------------------------------------------------
155     Double_t MVAMet::evaluatePhi() {
156     fPhiVals[0] = fNPV ;
157     fPhiVals[1] = fU ;
158     fPhiVals[2] = fUPhi ;
159     fPhiVals[3] = fTKSumEt ;
160     fPhiVals[4] = fTKU ;
161     fPhiVals[5] = fTKUPhi ;
162     fPhiVals[6] = fNPSumEt ;
163     fPhiVals[7] = fNPU ;
164     fPhiVals[8] = fNPUPhi ;
165     fPhiVals[9] = fPUSumEt ;
166     fPhiVals[10] = fPUMet ;
167     fPhiVals[11] = fPUMetPhi;
168     fPhiVals[12] = fPCSumEt ;
169     fPhiVals[13] = fPCU ;
170     fPhiVals[14] = fPCUPhi ;
171     fPhiVals[15] = fJSPt1 ;
172     fPhiVals[16] = fJSEta1 ;
173     fPhiVals[17] = fJSPhi1 ;
174     fPhiVals[18] = fJSPt2 ;
175     fPhiVals[19] = fJSEta2 ;
176     fPhiVals[20] = fJSPhi2 ;
177     fPhiVals[21] = fNAllJet ;
178     fPhiVals[22] = fNJet ;
179     return fPhiReader->GetResponse(fPhiVals);
180     }
181     //--------------------------------------------------------------------------------------------------
182     Double_t MVAMet::evaluateU1() {
183 pharris 1.5 fU1Vals[0] = fSumEt ;
184     fU1Vals[1] = fNPV ;
185     fU1Vals[2] = fU ;
186     fU1Vals[3] = fUPhi ;
187     fU1Vals[4] = fTKSumEt ;
188     fU1Vals[5] = fTKU ;
189     fU1Vals[6] = fTKUPhi ;
190     fU1Vals[7] = fNPSumEt ;
191     fU1Vals[8] = fNPU ;
192     fU1Vals[9] = fNPUPhi ;
193     fU1Vals[10] = fPUSumEt ;
194     fU1Vals[11] = fPUMet ;
195     fU1Vals[12] = fPUMetPhi;
196     fU1Vals[13] = fPCSumEt ;
197     fU1Vals[14] = fPCU ;
198     fU1Vals[15] = fPCUPhi ;
199     fU1Vals[16] = fJSPt1 ;
200     fU1Vals[17] = fJSEta1 ;
201     fU1Vals[18] = fJSPhi1 ;
202     fU1Vals[19] = fJSPt2 ;
203     fU1Vals[20] = fJSEta2 ;
204     fU1Vals[21] = fJSPhi2 ;
205     fU1Vals[22] = fNAllJet ;
206     fU1Vals[23] = fNJet ;
207     fU1Vals[24] = fUPhiMVA ;
208 pharris 1.2 return fU1Reader->GetResponse(fU1Vals);
209     }
210 pharris 1.1 //--------------------------------------------------------------------------------------------------
211 pharris 1.4 Double_t MVAMet::evaluateCovU1() {
212 pharris 1.5 fCovVals[0] = fSumEt ;
213     fCovVals[1] = fNPV ;
214     fCovVals[2] = fU ;
215     fCovVals[3] = fUPhi ;
216     fCovVals[4] = fTKSumEt ;
217     fCovVals[5] = fTKU ;
218     fCovVals[6] = fTKUPhi ;
219     fCovVals[7] = fNPSumEt ;
220     fCovVals[8] = fNPU ;
221     fCovVals[9] = fNPUPhi ;
222     fCovVals[10] = fPUSumEt ;
223     fCovVals[11] = fPUMet ;
224     fCovVals[12] = fPUMetPhi;
225     fCovVals[13] = fPCSumEt ;
226     fCovVals[14] = fPCU ;
227     fCovVals[15] = fPCUPhi ;
228     fCovVals[16] = fJSPt1 ;
229     fCovVals[17] = fJSEta1 ;
230     fCovVals[18] = fJSPhi1 ;
231     fCovVals[19] = fJSPt2 ;
232     fCovVals[20] = fJSEta2 ;
233     fCovVals[21] = fJSPhi2 ;
234     fCovVals[22] = fNAllJet ;
235     fCovVals[23] = fNJet ;
236     fCovVals[24] = fUPhiMVA ;
237     fCovVals[25] = fUMVA ;
238 pharris 1.4 return fCovU1Reader->GetResponse(fCovVals);
239     }
240     //--------------------------------------------------------------------------------------------------
241     Double_t MVAMet::evaluateCovU2() {
242 pharris 1.5 fCovVals[0] = fSumEt ;
243     fCovVals[1] = fNPV ;
244     fCovVals[2] = fU ;
245     fCovVals[3] = fUPhi ;
246     fCovVals[4] = fTKSumEt ;
247     fCovVals[5] = fTKU ;
248     fCovVals[6] = fTKUPhi ;
249     fCovVals[7] = fNPSumEt ;
250     fCovVals[8] = fNPU ;
251     fCovVals[9] = fNPUPhi ;
252     fCovVals[10] = fPUSumEt ;
253     fCovVals[11] = fPUMet ;
254     fCovVals[12] = fPUMetPhi;
255     fCovVals[13] = fPCSumEt ;
256     fCovVals[14] = fPCU ;
257     fCovVals[15] = fPCUPhi ;
258     fCovVals[16] = fJSPt1 ;
259     fCovVals[17] = fJSEta1 ;
260     fCovVals[18] = fJSPhi1 ;
261     fCovVals[19] = fJSPt2 ;
262     fCovVals[20] = fJSEta2 ;
263     fCovVals[21] = fJSPhi2 ;
264     fCovVals[22] = fNAllJet ;
265     fCovVals[23] = fNJet ;
266     fCovVals[24] = fUPhiMVA ;
267     fCovVals[25] = fUMVA ;
268 pharris 1.4 return fCovU2Reader->GetResponse(fCovVals);
269     }
270     //--------------------------------------------------------------------------------------------------
271 pharris 1.1 Double_t MVAMet::MVAValue( Bool_t iPhi,
272     Float_t iPFSumEt,
273     Float_t iU ,
274     Float_t iUPhi ,
275     Float_t iTKSumEt,
276     Float_t iTKU ,
277     Float_t iTKUPhi ,
278     Float_t iNPSumEt,
279     Float_t iNPU ,
280     Float_t iNPUPhi ,
281     Float_t iPUSumEt,
282     Float_t iPUMet ,
283     Float_t iPUMetPhi,
284     Float_t iPCSumEt,
285     Float_t iPCU ,
286     Float_t iPCUPhi ,
287     Float_t iJSPt1 ,
288     Float_t iJSEta1 ,
289     Float_t iJSPhi1 ,
290     Float_t iJSPt2 ,
291     Float_t iJSEta2 ,
292     Float_t iJSPhi2 ,
293     Float_t iNAllJet,
294     Float_t iNJet ,
295     Float_t iNPV
296     ){
297    
298     if (!fIsInitialized) {
299     std::cout << "Error: MVA Met not properly initialized.\n";
300     return -9999;
301     }
302 pharris 1.5
303     fSumEt = iPFSumEt;
304 pharris 1.1 fU = iU ;
305     fUPhi = iUPhi ;
306     fTKSumEt = iTKSumEt/iPFSumEt;
307     fTKU = iTKU ;
308     fTKUPhi = iTKUPhi ;
309     fNPSumEt = iNPSumEt/iPFSumEt;
310     fNPU = iNPU ;
311     fNPUPhi = iNPUPhi ;
312     fPUSumEt = iPUSumEt/iPFSumEt;
313     fPUMet = iPUMet ;
314     fPUMetPhi = iPUMetPhi;
315     fPCSumEt = iPCSumEt/iPFSumEt;
316     fPCU = iPCU ;
317     fPCUPhi = iPCUPhi ;
318     fJSPt1 = iJSPt1 ;
319     fJSEta1 = iJSEta1 ;
320     fJSPhi1 = iJSPhi1 ;
321     fJSPt2 = iJSPt2 ;
322     fJSEta2 = iJSEta2 ;
323     fJSPhi2 = iJSPhi2 ;
324     fNAllJet = iNAllJet;
325     fNJet = iNJet ;
326     fNPV = iNPV ;
327 pharris 1.5
328 pharris 1.1 Double_t lMVA = -9999;
329 pharris 1.2 lMVA = evaluatePhi();
330 pharris 1.1 if(!iPhi) fUPhiMVA = iUPhi+lMVA;
331 pharris 1.2 //Not no nice feature of the training
332 pharris 1.4 //fTKSumEt /= iPFSumEt;
333     //fNPSumEt /= iPFSumEt;
334     //fPUSumEt /= iPFSumEt;
335     //fPCSumEt /= iPFSumEt;
336 pharris 1.2 if(!iPhi) lMVA = evaluateU1();
337 pharris 1.1 return lMVA;
338     }
339     //--------------------------------------------------------------------------------------------------
340     //====> Please not that the jet collection must be cleaned => all jets near leptons must be removed
341     Met MVAMet::GetMet( Bool_t iPhi,Float_t iPtVis,Float_t iPhiVis,Float_t iSumEtVis,
342     const PFMet *iMet ,
343 pharris 1.3 const PFCandidateCol *iCands,
344     const Vertex *iVertex,const VertexCol *iVertices,
345 pharris 1.1 const PFJetCol *iJets ,
346     FactorizedJetCorrector *iJetCorrector,
347     const PileupEnergyDensityCol *iPUEnergyDensity,
348     int iNPV,
349     Bool_t printDebug) {
350    
351 pharris 1.7 Met lPFRec = fRecoilTools->pfRecoil (iPtVis,iPhiVis,iSumEtVis,iCands);
352 pharris 1.1 Met lTKRec = fRecoilTools->trackRecoil(iPtVis,iPhiVis,iSumEtVis,iCands,iVertex);
353 pharris 1.3 Met lNPRec = fRecoilTools->NoPURecoil (iPtVis,iPhiVis,iSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices);
354     Met lPCRec = fRecoilTools->PUCRecoil (iPtVis,iPhiVis,iSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices);
355     Met lPUMet = fRecoilTools->PUMet (iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices);
356 pharris 1.1
357     Double_t lPt0 = 0; const PFJet *lLead = 0;
358     Double_t lPt1 = 0; const PFJet *l2nd = 0;
359     int lNAllJet = 0;
360     int lNJet = 0;
361     for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
362     const PFJet *pJet = iJets->At(i0);
363     Double_t pPt = fRecoilTools->fJetIDMVA->correctedPt(pJet,iJetCorrector,iPUEnergyDensity);
364     if(!JetTools::passPFLooseId(pJet)) continue;
365     lNAllJet++;
366     if(pPt > 30) lNJet++;
367     if(lPt0 < pPt) {lPt0 = pPt; lLead = pJet; continue;}
368     if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
369     }
370 pharris 1.5
371     fSumEt = lPFRec.SumEt();
372 pharris 1.1 fU = lPFRec.Pt();
373     fUPhi = lPFRec.Phi();
374     fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
375     fTKU = lTKRec.Pt();
376     fTKUPhi = lTKRec.Phi();
377     fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
378     fNPU = lNPRec.Pt();
379     fNPUPhi = lNPRec.Phi();
380     fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
381     fPUMet = lPUMet.Pt();
382     fPUMetPhi = lPUMet.Phi();
383     fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
384     fPCU = lPCRec.Pt() ;
385     fPCUPhi = lPCRec.Phi() ;
386     fJSPt1 = lPt0;
387     fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
388     fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
389     fJSPt2 = lPt1;
390     fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
391     fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
392     fNJet = lNJet ;
393     fNAllJet = lNAllJet;
394     fNPV = iNPV ;
395 pharris 1.5 Float_t lMVA = evaluatePhi();
396     if(!iPhi) fUPhiMVA = fUPhi + lMVA;
397 pharris 1.1 //Not no nice feature of teh training
398 pharris 1.4 //fTKSumEt /= lPFRec.SumEt();
399     //fNPSumEt /= lPFRec.SumEt();
400     //fPUSumEt /= lPFRec.SumEt();
401     //fPCSumEt /= lPFRec.SumEt();
402 pharris 1.7 if(!iPhi) lMVA = evaluateU1();//fU1Reader ->EvaluateMVA( fU1MethodName );
403     fUMVA = fU*lMVA;
404    
405     TLorentzVector lUVec(0,0,0,0); lUVec .SetPtEtaPhiM(fUMVA,0,fUPhiMVA,0);
406 pharris 1.1 TLorentzVector lVVec(0,0,0,0); lVVec .SetPtEtaPhiM(iPtVis ,0,iPhiVis ,0);
407     if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
408     lUVec -= lVVec;
409     Met lMet(lUVec.Px(),lUVec.Py());
410 pharris 1.5
411 pharris 1.4 //TMatrixD lCov(2,2);
412     //Covariance matrix perpendicular and parallel to the recoil (sort of correct)
413     double lCovU1 = evaluateCovU1();
414     double lCovU2 = evaluateCovU2();
415    
416     double lCos2 = cos(fUPhiMVA)*cos(fUPhiMVA);
417     double lSin2 = sin(fUPhiMVA)*sin(fUPhiMVA);
418 pharris 1.5
419 pharris 1.4 //Now Compute teh covariance matrix in X and Y
420 pharris 1.5 (*fCov)(0,0) = lCovU1*lCos2+lCovU2*lSin2;
421     (*fCov)(1,0) = -lCovU1*sin(fUPhiMVA)*cos(fUPhiMVA)+lCovU2*sin(fUPhiMVA)*cos(fUPhiMVA);
422     (*fCov)(0,1) = (*fCov)(1,0);
423     (*fCov)(1,1) = lCovU1*lSin2+lCovU2*lCos2;
424 pharris 1.10 TMatrixD lInv(2,2);
425     lInv(0,0) = (*fCov)(0,0); lInv(1,1) = (*fCov)(1,1); lInv(1,0) = (*fCov)(1,0); lInv(0,1) = (*fCov)(0,1);
426     if(lInv.Determinant() != 0) lInv.Invert();
427     fSignificance = TMath::Sqrt(lUVec.Px()*lUVec.Px()*(lInv)(0,0) + 2.*lUVec.Px()*lUVec.Py()*(lInv)(1,0) + lUVec.Py()*lUVec.Py()*(lInv)(1,1));
428     fUncertainty = sqrt(lCovU1+lCovU2);
429 pharris 1.7
430 pharris 1.1 if (printDebug == kTRUE) {
431 ceballos 1.11 std::cout << "Debug Met MVA: "
432 pharris 1.1 << fU << " : "
433     << fUPhi << " : "
434     << fTKSumEt << " : "
435     << fTKU << " : "
436     << fTKUPhi << " : "
437     << fNPSumEt << " : "
438     << fNPU << " : "
439     << fNPUPhi << " : "
440     << fPUSumEt << " : "
441     << fPUMet << " : "
442     << fPUMetPhi << " : "
443     << fPCUPhi << " : "
444     << fPCSumEt << " : "
445     << fPCU << " : "
446     << fPCUPhi << " : "
447     << fJSPt1 << " : "
448     << fJSEta1 << " : "
449     << fJSPhi1 << " : "
450     << fJSPt2 << " : "
451     << fJSEta2 << " : "
452     << fJSPhi2 << " : "
453     << fNJet << " : "
454     << fNAllJet << " : "
455     << fNPV << " : "
456     << " === : === "
457 ceballos 1.12 << lMet.Pt() << " : "
458     << lMet.Phi() << " : "
459 pharris 1.1 << std::endl;
460     }
461    
462     return lMet;
463     }
464 pharris 1.2 //--------------------------------------------------------------------------------------------------
465     //====> Please not that the jet collection must be cleaned => all jets near leptons must be removed
466     //====> Corrected Jets
467     Met MVAMet::GetMet( Bool_t iPhi,Float_t iPtVis,Float_t iPhiVis,Float_t iSumEtVis,
468     const PFMet *iMet ,
469 pharris 1.3 const PFCandidateCol *iCands,
470 pharris 1.7 const Vertex *iVertex ,const VertexCol *iVertices,Double_t iRho,
471 pharris 1.2 const PFJetCol *iJets ,
472     int iNPV,
473     Bool_t printDebug) {
474    
475 pharris 1.7 Met lPFRec = fRecoilTools->pfRecoil (iPtVis,iPhiVis,iSumEtVis,iCands);
476 pharris 1.3 Met lTKRec = fRecoilTools->trackRecoil(iPtVis,iPhiVis,iSumEtVis, iCands,iVertex);
477 pharris 1.7 Met lNPRec = fRecoilTools->NoPURecoil (iPtVis,iPhiVis,iSumEtVis,iJets,iCands,iVertex,iVertices,iRho);
478     Met lPCRec = fRecoilTools->PUCRecoil (iPtVis,iPhiVis,iSumEtVis,iJets,iCands,iVertex,iVertices,iRho);
479     Met lPUMet = fRecoilTools->PUMet ( iJets,iCands,iVertex,iVertices,iRho);
480 pharris 1.2
481     Double_t lPt0 = 0; const PFJet *lLead = 0;
482     Double_t lPt1 = 0; const PFJet *l2nd = 0;
483     int lNAllJet = 0;
484     int lNJet = 0;
485     for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
486     const PFJet *pJet = iJets->At(i0);
487     Double_t pPt = pJet->Pt();
488     if(!JetTools::passPFLooseId(pJet)) continue;
489     lNAllJet++;
490     if(pPt > 30) lNJet++;
491     if(lPt0 < pPt) {lPt0 = pPt; lLead = pJet; continue;}
492     if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
493     }
494    
495 pharris 1.5 fSumEt = lPFRec.SumEt();
496 pharris 1.2 fU = lPFRec.Pt();
497     fUPhi = lPFRec.Phi();
498     fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
499     fTKU = lTKRec.Pt();
500     fTKUPhi = lTKRec.Phi();
501     fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
502     fNPU = lNPRec.Pt();
503     fNPUPhi = lNPRec.Phi();
504     fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
505     fPUMet = lPUMet.Pt();
506     fPUMetPhi = lPUMet.Phi();
507     fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
508     fPCU = lPCRec.Pt() ;
509     fPCUPhi = lPCRec.Phi() ;
510     fJSPt1 = lPt0;
511     fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
512     fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
513     fJSPt2 = lPt1;
514     fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
515     fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
516     fNJet = lNJet ;
517     fNAllJet = lNAllJet;
518     fNPV = iNPV ;
519 pharris 1.7
520 pharris 1.2 Float_t lMVA = evaluatePhi();
521    
522     if(!iPhi) fUPhiMVA = fUPhi + lMVA;
523     //Not no nice feature of teh training
524 pharris 1.4 //fTKSumEt /= lPFRec.SumEt();
525     //fNPSumEt /= lPFRec.SumEt();
526     //fPUSumEt /= lPFRec.SumEt();
527     //fPCSumEt /= lPFRec.SumEt();
528 pharris 1.7 if(!iPhi) lMVA = evaluateU1();//fU1Reader ->EvaluateMVA( fU1MethodName );
529     fUMVA = lMVA*fU;
530    
531     TLorentzVector lUVec(0,0,0,0); lUVec .SetPtEtaPhiM(fUMVA ,0,fUPhiMVA,0);
532 pharris 1.2 TLorentzVector lVVec(0,0,0,0); lVVec .SetPtEtaPhiM(iPtVis ,0,iPhiVis ,0);
533     if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
534     lUVec -= lVVec;
535     Met lMet(lUVec.Px(),lUVec.Py());
536 pharris 1.4 //Cov matrix
537     double lCovU1 = evaluateCovU1();
538     double lCovU2 = evaluateCovU2();
539    
540     double lCos2 = cos(fUPhiMVA)*cos(fUPhiMVA);
541     double lSin2 = sin(fUPhiMVA)*sin(fUPhiMVA);
542    
543     //Now Compute teh covariance matrix in X and Y
544 pharris 1.5 (*fCov)(0,0) = lCovU1*lCos2+lCovU2*lSin2;
545     (*fCov)(1,0) = -lCovU1*sin(fUPhiMVA)*cos(fUPhiMVA)+lCovU2*sin(fUPhiMVA)*cos(fUPhiMVA);
546     (*fCov)(0,1) = (*fCov)(1,0);
547     (*fCov)(1,1) = lCovU1*lSin2+lCovU2*lCos2;
548 pharris 1.10 TMatrixD lInv(2,2);
549     lInv(0,0) = (*fCov)(0,0); lInv(1,1) = (*fCov)(1,1); lInv(1,0) = (*fCov)(1,0); lInv(0,1) = (*fCov)(0,1);
550     if(lInv.Determinant() != 0) lInv.Invert();
551     fSignificance = TMath::Sqrt(lUVec.Px()*lUVec.Px()*(lInv)(0,0) + 2.*lUVec.Px()*lUVec.Py()*(lInv)(1,0) + lUVec.Py()*lUVec.Py()*(lInv)(1,1));
552     fUncertainty = sqrt(lCovU1+lCovU2);
553 pharris 1.2
554     if (printDebug == kTRUE) {
555 ceballos 1.11 std::cout << "Debug Met MVA: "
556 pharris 1.2 << fU << " : "
557     << fUPhi << " : "
558     << fTKSumEt << " : "
559     << fTKU << " : "
560     << fTKUPhi << " : "
561     << fNPSumEt << " : "
562     << fNPU << " : "
563     << fNPUPhi << " : "
564     << fPUSumEt << " : "
565     << fPUMet << " : "
566     << fPUMetPhi << " : "
567     << fPCUPhi << " : "
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 ceballos 1.12 << lMet.Pt() << " : "
582     << lMet.Phi() << " : "
583 pharris 1.2 << std::endl;
584     }
585    
586     return lMet;
587     }
588     //--------------------------------------------------------------------------------------------------
589 pharris 1.1 //Interms of the two candidates
590     Met MVAMet::GetMet( Bool_t iPhi,
591     Float_t iPt1,Float_t iPhi1,Float_t iEta1,
592     Float_t iPt2,Float_t iPhi2,Float_t iEta2,
593     const PFMet *iMet ,
594 pharris 1.3 const PFCandidateCol *iCands,
595     const Vertex *iVertex,const VertexCol *iVertices,
596 pharris 1.1 const PFJetCol *iJets ,
597     FactorizedJetCorrector *iJetCorrector,
598     const PileupEnergyDensityCol *iPUEnergyDensity,
599     int iNPV,
600     Bool_t printDebug) {
601    
602     TLorentzVector lVVec1(0,0,0,0); lVVec1.SetPtEtaPhiM(iPt1,0,iPhi1 ,0);
603     TLorentzVector lVVec2(0,0,0,0); lVVec2.SetPtEtaPhiM(iPt2,0,iPhi2 ,0);
604     lVVec1+=lVVec2;
605     Float_t lPtVis = lVVec1.Pt();
606     Float_t lPhiVis = lVVec1.Phi();
607     Float_t lSumEtVis = iPt1 + iPt2;
608 pharris 1.7 Met lPFRec = fRecoilTools->pfRecoil (lPtVis,lPhiVis,lSumEtVis,iCands);
609 pharris 1.1 Met lTKRec = fRecoilTools->trackRecoil(lPtVis,lPhiVis,lSumEtVis,iCands,iVertex);
610 pharris 1.3 Met lNPRec = fRecoilTools->NoPURecoil (lPtVis,lPhiVis,lSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices,
611     iPhi1,iEta1,iPhi2,iEta2);
612     Met lPCRec = fRecoilTools->PUCRecoil (lPtVis,lPhiVis,lSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices,
613     iPhi1,iEta1,iPhi2,iEta2);
614     Met lPUMet = fRecoilTools->PUMet ( iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices,
615     iPhi1,iEta1,iPhi2,iEta2);
616 pharris 1.1
617     Double_t lPt0 = 0; const PFJet *lLead = 0;
618     Double_t lPt1 = 0; const PFJet *l2nd = 0;
619     int lNAllJet = 0;
620     int lNJet = 0;
621     for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
622     const PFJet *pJet = iJets->At(i0);
623     Double_t pPt = fRecoilTools->fJetIDMVA->correctedPt(pJet,iJetCorrector,iPUEnergyDensity);
624     double pDEta1 = pJet->Eta() - iEta1;
625     double pDPhi1 = fabs(pJet->Phi() - iPhi1); if(pDPhi1 > 2.*TMath::Pi()-pDPhi1) pDPhi1 = 2.*TMath::Pi()-pDPhi1;
626     double pDR1 = sqrt(pDEta1*pDEta1 + pDPhi1*pDPhi1);
627     if(pDR1 < 0.5) continue;
628     double pDEta2 = pJet->Eta() - iEta2;
629     double pDPhi2 = fabs(pJet->Phi() - iPhi2); if(pDPhi2 > 2.*TMath::Pi()-pDPhi2) pDPhi2 = 2.*TMath::Pi()-pDPhi2;
630     double pDR2 = sqrt(pDEta2*pDEta2 + pDPhi2*pDPhi2);
631     if(pDR2 < 0.5) continue;
632     if(!JetTools::passPFLooseId(pJet)) continue;
633     lNAllJet++;
634     if(pPt > 30.) lNJet++;
635     if(lPt0 < pPt) {lPt0 = pPt; lLead = pJet; continue;}
636     if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
637     }
638    
639 pharris 1.5 fSumEt = lPFRec.SumEt();
640 pharris 1.1 fU = lPFRec.Pt();
641     fUPhi = lPFRec.Phi();
642     fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
643     fTKU = lTKRec.Pt();
644     fTKUPhi = lTKRec.Phi();
645     fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
646     fNPU = lNPRec.Pt();
647     fNPUPhi = lNPRec.Phi();
648     fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
649     fPUMet = lPUMet.Pt();
650     fPUMetPhi= lPUMet.Phi();
651     fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
652     fPCU = lPCRec.Pt() ;
653     fPCUPhi = lPCRec.Phi() ;
654     fJSPt1 = lPt0;
655     fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
656     fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
657     fJSPt2 = lPt1;
658     fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
659     fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
660     fNJet = lNJet ;
661     fNAllJet = lNAllJet;
662     fNPV = iNPV ;
663 pharris 1.5
664 pharris 1.2 Float_t lMVA = evaluatePhi();//fPhiReader ->EvaluateMVA( fPhiMethodName );
665     if(!iPhi) fUPhiMVA = fUPhi + lMVA;
666 pharris 1.4 //fTKSumEt /= lPFRec.SumEt();
667     //fNPSumEt /= lPFRec.SumEt();
668     //fPUSumEt /= lPFRec.SumEt();
669     //fPCSumEt /= lPFRec.SumEt();
670 pharris 1.7 if(!iPhi) lMVA = evaluateU1();//fU1Reader ->EvaluateMVA( fU1MethodName );
671     fUMVA = fU*lMVA;
672    
673     TLorentzVector lUVec (0,0,0,0); lUVec .SetPtEtaPhiM(fUMVA ,0,fUPhiMVA,0);
674 pharris 1.2 TLorentzVector lVVec (0,0,0,0); lVVec .SetPtEtaPhiM(lPtVis ,0,lPhiVis ,0);
675     if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
676     lUVec -= lVVec;
677     Met lMet(lUVec.Px(),lUVec.Py());
678 pharris 1.4 //Cov matrix
679     double lCovU1 = evaluateCovU1();
680     double lCovU2 = evaluateCovU2();
681    
682     double lCos2 = cos(fUPhiMVA)*cos(fUPhiMVA);
683     double lSin2 = sin(fUPhiMVA)*sin(fUPhiMVA);
684    
685     //Now Compute teh covariance matrix in X and Y
686 pharris 1.5 (*fCov)(0,0) = lCovU1*lCos2+lCovU2*lSin2;
687     (*fCov)(1,0) = -lCovU1*sin(fUPhiMVA)*cos(fUPhiMVA)+lCovU2*sin(fUPhiMVA)*cos(fUPhiMVA);
688     (*fCov)(0,1) = (*fCov)(1,0);
689     (*fCov)(1,1) = lCovU1*lSin2+lCovU2*lCos2;
690 pharris 1.10 TMatrixD lInv(2,2);
691     lInv(0,0) = (*fCov)(0,0); lInv(1,1) = (*fCov)(1,1); lInv(1,0) = (*fCov)(1,0); lInv(0,1) = (*fCov)(0,1);
692     if(lInv.Determinant() != 0) lInv.Invert();
693     fSignificance = TMath::Sqrt(lUVec.Px()*lUVec.Px()*(lInv)(0,0) + 2.*lUVec.Px()*lUVec.Py()*(lInv)(1,0) + lUVec.Py()*lUVec.Py()*(lInv)(1,1));
694     fUncertainty = sqrt(lCovU1+lCovU2);
695    
696 pharris 1.2
697     if (printDebug == kTRUE) {
698 ceballos 1.11 std::cout << "Debug Met MVA: "
699 pharris 1.2 << fU << " : "
700     << fUPhi << " : "
701     << fTKSumEt << " : "
702     << fTKU << " : "
703     << fTKUPhi << " : "
704     << fNPSumEt << " : "
705     << fNPU << " : "
706     << fNPUPhi << " : "
707     << fPUSumEt << " : "
708     << fPUMet << " : "
709     << fPUMetPhi << " : "
710     << fPCSumEt << " : "
711     << fPCU << " : "
712     << fPCUPhi << " : "
713     << fJSPt1 << " : "
714     << fJSEta1 << " : "
715     << fJSPhi1 << " : "
716     << fJSPt2 << " : "
717     << fJSEta2 << " : "
718     << fJSPhi2 << " : "
719     << fNJet << " : "
720     << fNAllJet << " : "
721     << fNPV << " : "
722     << " === : === "
723 ceballos 1.12 << lMet.Pt() << " : "
724     << lMet.Phi() << " : "
725 pharris 1.2 << std::endl;
726     }
727    
728     return lMet;
729     }
730     //--------------------------------------------------------------------------------------------------
731     //Interms of the two candidates => corrected Jets
732     Met MVAMet::GetMet( Bool_t iPhi,
733     Float_t iPt1,Float_t iPhi1,Float_t iEta1,
734     Float_t iPt2,Float_t iPhi2,Float_t iEta2,
735     const PFMet *iMet ,
736 pharris 1.3 const PFCandidateCol *iCands,
737 pharris 1.7 const Vertex *iVertex,const VertexCol *iVertices,Double_t iRho,
738 pharris 1.2 const PFJetCol *iJets ,
739     int iNPV,
740     Bool_t printDebug) {
741    
742     TLorentzVector lVVec1(0,0,0,0); lVVec1.SetPtEtaPhiM(iPt1,0,iPhi1 ,0);
743     TLorentzVector lVVec2(0,0,0,0); lVVec2.SetPtEtaPhiM(iPt2,0,iPhi2 ,0);
744     lVVec1+=lVVec2;
745     Float_t lPtVis = lVVec1.Pt();
746     Float_t lPhiVis = lVVec1.Phi();
747     Float_t lSumEtVis = iPt1 + iPt2;
748 pharris 1.7 Met lPFRec = fRecoilTools->pfRecoil (lPtVis,lPhiVis,lSumEtVis,iCands);
749 pharris 1.2 Met lTKRec = fRecoilTools->trackRecoil(lPtVis,lPhiVis,lSumEtVis,iCands,iVertex);
750 pharris 1.3 Met lNPRec = fRecoilTools->NoPURecoil (lPtVis,lPhiVis,lSumEtVis,iJets,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2);
751     Met lPCRec = fRecoilTools->PUCRecoil (lPtVis,lPhiVis,lSumEtVis,iJets,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2);
752     Met lPUMet = fRecoilTools->PUMet ( iJets,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2);
753 pharris 1.2
754     Double_t lPt0 = 0; const PFJet *lLead = 0;
755     Double_t lPt1 = 0; const PFJet *l2nd = 0;
756     int lNAllJet = 0;
757     int lNJet = 0;
758     for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
759     const PFJet *pJet = iJets->At(i0);
760     Double_t pPt = pJet->Pt();
761     double pDEta1 = pJet->Eta() - iEta1;
762     double pDPhi1 = fabs(pJet->Phi() - iPhi1); if(pDPhi1 > 2.*TMath::Pi()-pDPhi1) pDPhi1 = 2.*TMath::Pi()-pDPhi1;
763     double pDR1 = sqrt(pDEta1*pDEta1 + pDPhi1*pDPhi1);
764     if(pDR1 < 0.5) continue;
765     double pDEta2 = pJet->Eta() - iEta2;
766     double pDPhi2 = fabs(pJet->Phi() - iPhi2); if(pDPhi2 > 2.*TMath::Pi()-pDPhi2) pDPhi2 = 2.*TMath::Pi()-pDPhi2;
767     double pDR2 = sqrt(pDEta2*pDEta2 + pDPhi2*pDPhi2);
768     if(pDR2 < 0.5) continue;
769     if(!JetTools::passPFLooseId(pJet)) continue;
770     lNAllJet++;
771     if(pPt > 30.) lNJet++;
772     if(lPt0 < pPt) {lPt0 = pPt; lLead = pJet; continue;}
773     if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
774     }
775 pharris 1.5 fSumEt = lPFRec.SumEt();
776 pharris 1.2 fU = lPFRec.Pt();
777     fUPhi = lPFRec.Phi();
778     fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
779     fTKU = lTKRec.Pt();
780     fTKUPhi = lTKRec.Phi();
781     fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
782     fNPU = lNPRec.Pt();
783     fNPUPhi = lNPRec.Phi();
784     fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
785     fPUMet = lPUMet.Pt();
786     fPUMetPhi= lPUMet.Phi();
787     fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
788     fPCU = lPCRec.Pt() ;
789     fPCUPhi = lPCRec.Phi() ;
790     fJSPt1 = lPt0;
791     fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
792     fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
793     fJSPt2 = lPt1;
794     fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
795     fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
796     fNJet = lNJet ;
797     fNAllJet = lNAllJet;
798     fNPV = iNPV ;
799    
800     Float_t lMVA = evaluatePhi();//fPhiReader ->EvaluateMVA( fPhiMethodName );
801 pharris 1.1
802     if(!iPhi) fUPhiMVA = fUPhi + lMVA;
803 pharris 1.4 //fTKSumEt /= lPFRec.SumEt();
804     //fNPSumEt /= lPFRec.SumEt();
805     //fPUSumEt /= lPFRec.SumEt();
806     //fPCSumEt /= lPFRec.SumEt();
807 pharris 1.7 if(!iPhi) lMVA = evaluateU1();//fU1Reader ->EvaluateMVA( fU1MethodName );
808     fUMVA = lMVA*fU;
809    
810     TLorentzVector lUVec (0,0,0,0); lUVec .SetPtEtaPhiM(fUMVA ,0,fUPhiMVA,0);
811 pharris 1.1 TLorentzVector lVVec (0,0,0,0); lVVec .SetPtEtaPhiM(lPtVis ,0,lPhiVis ,0);
812 pharris 1.7 if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
813 pharris 1.1 lUVec -= lVVec;
814     Met lMet(lUVec.Px(),lUVec.Py());
815 pharris 1.4 //Cov matrix
816     double lCovU1 = evaluateCovU1();
817     double lCovU2 = evaluateCovU2();
818    
819     double lCos2 = cos(fUPhiMVA)*cos(fUPhiMVA);
820     double lSin2 = sin(fUPhiMVA)*sin(fUPhiMVA);
821     //Now Compute teh covariance matrix in X and Y
822 pharris 1.5 (*fCov)(0,0) = lCovU1*lCos2+lCovU2*lSin2;
823     (*fCov)(1,0) = -lCovU1*sin(fUPhiMVA)*cos(fUPhiMVA)+lCovU2*sin(fUPhiMVA)*cos(fUPhiMVA);
824     (*fCov)(0,1) = (*fCov)(1,0);
825     (*fCov)(1,1) = lCovU1*lSin2+lCovU2*lCos2;
826 pharris 1.10 TMatrixD lInv(2,2);
827     lInv(0,0) = (*fCov)(0,0); lInv(1,1) = (*fCov)(1,1); lInv(1,0) = (*fCov)(1,0); lInv(0,1) = (*fCov)(0,1);
828     if(lInv.Determinant() != 0) lInv.Invert();
829     fSignificance = TMath::Sqrt(lUVec.Px()*lUVec.Px()*(lInv)(0,0) + 2.*lUVec.Px()*lUVec.Py()*(lInv)(1,0) + lUVec.Py()*lUVec.Py()*(lInv)(1,1));
830     fUncertainty = sqrt(lCovU1+lCovU2);
831 pharris 1.7
832 pharris 1.1 if (printDebug == kTRUE) {
833 ceballos 1.11 std::cout << "Debug Met MVA: "
834 pharris 1.1 << fU << " : "
835     << fUPhi << " : "
836     << fTKSumEt << " : "
837     << fTKU << " : "
838     << fTKUPhi << " : "
839     << fNPSumEt << " : "
840     << fNPU << " : "
841     << fNPUPhi << " : "
842     << fPUSumEt << " : "
843     << fPUMet << " : "
844     << fPUMetPhi << " : "
845     << fPCSumEt << " : "
846     << fPCU << " : "
847     << fPCUPhi << " : "
848     << fJSPt1 << " : "
849     << fJSEta1 << " : "
850     << fJSPhi1 << " : "
851     << fJSPt2 << " : "
852     << fJSEta2 << " : "
853     << fJSPhi2 << " : "
854     << fNJet << " : "
855     << fNAllJet << " : "
856     << fNPV << " : "
857     << " === : === "
858 ceballos 1.12 << lMet.Pt() << " : "
859     << lMet.Phi() << " : "
860 pharris 1.1 << std::endl;
861     }
862    
863     return lMet;
864     }
865