ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MVAMet.cc
Revision: 1.22
Committed: Sat Jun 16 10:14:39 2012 UTC (12 years, 10 months ago) by ceballos
Content type: text/plain
Branch: MAIN
Changes since 1.21: +4 -6 lines
Log Message:
fixing MVAMet again

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