ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MVAMet.cc
Revision: 1.16
Committed: Thu May 3 14:01:48 2012 UTC (13 years ago) by pharris
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_028, Mit_027, Mit_027a
Changes since 1.15: +11 -5 lines
Log Message:
Add 42

File Contents

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