ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MVAMet.cc
Revision: 1.30
Committed: Sat Jan 12 11:49:50 2013 UTC (12 years, 3 months ago) by pharris
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_029c, Mit_029b, Mit_029a
Changes since 1.29: +84 -2 lines
Log Message: 
Updated MVA Met and Jet ID

File Contents

# User Rev Content
1 pharris 1.1 #include "MitPhysics/Utils/interface/MVAMet.h"
2 pharris 1.26 #include "MitPhysics/Utils/interface/MetLeptonTools.h"
3 pharris 1.1 #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.26 #include <utility>
11 pharris 1.1
12     ClassImp(mithep::MVAMet)
13    
14     using namespace mithep;
15    
16     //--------------------------------------------------------------------------------------------------
17     MVAMet::MVAMet() :
18     fRecoilTools(0),
19 pharris 1.4 fPhiMethodName ("PhiCorrection"),
20     fU1MethodName ("U1Correction"),
21     fCovU1MethodName ("CovU1"),
22     fCovU2MethodName ("CovU2"),
23 pharris 1.1 fIsInitialized(kFALSE),
24     fU (0),
25     fUPhi (0),
26     fTKSumEt(0),
27     fTKU (0),
28     fTKUPhi (0),
29     fNPSumEt(0),
30     fNPU (0),
31     fNPUPhi (0),
32     fPUSumEt(0),
33     fPUMet (0),
34     fPUMetPhi(0),
35     fPCSumEt(0),
36     fPCU (0),
37     fPCUPhi (0),
38     fJSPt1 (0),
39     fJSEta1 (0),
40     fJSPhi1 (0),
41     fJSPt2 (0),
42     fJSEta2 (0),
43     fJSPhi2 (0),
44     fNJet (0),
45     fNAllJet(0),
46     fNPV (0),
47     fPhiReader(0),
48 pharris 1.4 fU1Reader(0),
49     fCovU1Reader(0),
50 pharris 1.28 fCovU2Reader(0),
51     fMetLeptonTools(0) { }
52 pharris 1.1 //--------------------------------------------------------------------------------------------------
53     MVAMet::~MVAMet()
54     {
55 pharris 1.4 delete fPhiReader;
56     delete fU1Reader;
57     delete fCovU1Reader;
58     delete fCovU2Reader;
59    
60 pharris 1.1 }
61     //--------------------------------------------------------------------------------------------------
62 pharris 1.2 /*
63 pharris 1.1 void MVAMet::setVariables(TMVA::Reader *iReader,bool iScale) {
64     iReader->AddVariable( "npv" , &fNPV );
65     iReader->AddVariable( "u" , &fU );
66     iReader->AddVariable( "uphi" , &fUPhi );
67     iReader->AddVariable( "chsumet/sumet" , &fTKSumEt );
68     iReader->AddVariable( "tku" , &fTKU );
69     iReader->AddVariable( "tkuphi" , &fTKUPhi );
70     iReader->AddVariable( "nopusumet/sumet" , &fNPSumEt );
71     iReader->AddVariable( "nopuu" , &fNPU );
72     iReader->AddVariable( "nopuuphi" , &fNPUPhi );
73     iReader->AddVariable( "pusumet/sumet" , &fPUSumEt );
74     iReader->AddVariable( "pumet" , &fPUMet );
75     iReader->AddVariable( "pumetphi" , &fPUMetPhi);
76     iReader->AddVariable( "pucsumet/sumet" , &fPCSumEt );
77     iReader->AddVariable( "pucu" , &fPCU );
78     iReader->AddVariable( "pucuphi" , &fPCUPhi );
79     iReader->AddVariable( "jspt_1" , &fJSPt1 );
80     iReader->AddVariable( "jseta_1" , &fJSEta1 );
81     iReader->AddVariable( "jsphi_1" , &fJSPhi1 );
82     iReader->AddVariable( "jspt_2" , &fJSPt2 );
83     iReader->AddVariable( "jseta_2" , &fJSEta2 );
84     iReader->AddVariable( "jsphi_2" , &fJSPhi2 );
85     iReader->AddVariable( "nalljet" , &fNAllJet );
86     iReader->AddVariable( "njet" , &fNJet );
87     if(iScale) iReader->AddVariable( "uphi_mva" , &fUPhiMVA );
88     }
89     //--------------------------------------------------------------------------------------------------
90     void MVAMet::Initialize( TString iU1MethodName,
91     TString iPhiMethodName,
92     TString iJetMVAFile,
93     TString iU1Weights,
94     TString iPhiWeights,
95     MVAMet::MVAType iType) {
96    
97     fIsInitialized = kTRUE;
98     fRecoilTools = new RecoilTools(iJetMVAFile);
99 pharris 1.2
100 pharris 1.1 fU1MethodName = iU1MethodName;
101     fPhiMethodName = iPhiMethodName;
102     fType = iType;
103     fPhiReader = new TMVA::Reader( "!Color:!Silent:Error" );
104     fU1Reader = new TMVA::Reader( "!Color:!Silent:Error" );
105     if (fType == kBaseline) {
106     setVariables(fPhiReader,false);
107     setVariables(fU1Reader ,true);
108     }
109     fPhiReader->BookMVA(fPhiMethodName , iPhiWeights );
110     fU1Reader ->BookMVA(fU1MethodName , iU1Weights );
111 pharris 1.2
112 pharris 1.1 std::cout << "Jet ID MVA Initialization\n";
113     std::cout << "Phi Method Name : " << fPhiMethodName << " , type == " << iType << std::endl;
114     std::cout << "U1 Method Name : " << fU1MethodName << " , type == " << iType << std::endl;
115 pharris 1.2
116    
117 pharris 1.1 }
118 pharris 1.2 */
119     //--------------------------------------------------------------------------------------------------
120 pharris 1.3 void MVAMet::Initialize(TString iJetLowPtFile,
121     TString iJetHighPtFile,
122     TString iJetCutFile,
123     TString iU1Weights,
124     TString iPhiWeights,
125 pharris 1.4 TString iCovU1Weights,
126     TString iCovU2Weights,
127 pharris 1.30 JetIDMVA::MVAType iType) {
128 pharris 1.2
129     fIsInitialized = kTRUE;
130     fType = iType;
131 pharris 1.16 f42 = iU1Weights.Contains("42");
132 pharris 1.30
133     fRecoilTools = new RecoilTools(iJetLowPtFile,iJetHighPtFile,iJetCutFile,f42,fType);
134 pharris 1.16 if(f42) fRecoilTools->fJetIDMVA->fJetPtMin = 1.;
135 pharris 1.2
136 bendavid 1.6 ROOT::Cintex::Cintex::Enable();
137    
138 pharris 1.2 TFile *lPhiForest = new TFile(iPhiWeights,"READ");
139     fPhiReader = (GBRForest*)lPhiForest->Get(fPhiMethodName);
140     lPhiForest->Close();
141    
142     TFile *lU1Forest = new TFile(iU1Weights,"READ");
143     fU1Reader = (GBRForest*)lU1Forest->Get(fU1MethodName);
144     lU1Forest->Close();
145 pharris 1.1
146 pharris 1.4 TFile *lCovU1Forest = new TFile(iCovU1Weights,"READ");
147     fCovU1Reader = (GBRForest*)lCovU1Forest->Get(fCovU1MethodName);
148     lCovU1Forest->Close();
149    
150     TFile *lCovU2Forest = new TFile(iCovU2Weights,"READ");
151     fCovU2Reader = (GBRForest*)lCovU2Forest->Get(fCovU2MethodName);
152     lCovU2Forest->Close();
153 pharris 1.5
154     fCov = new TMatrixD(2,2);
155 pharris 1.7 fPhiVals = new Float_t[23];
156     fU1Vals = new Float_t[25];
157     fCovVals = new Float_t[26];
158 pharris 1.28
159    
160     fMetLeptonTools = new MetLeptonTools();
161 pharris 1.2 }
162     //--------------------------------------------------------------------------------------------------
163 pharris 1.23 Float_t* MVAMet::getVals() {
164     fU1Vals[0] = fSumEt ;
165     fU1Vals[1] = fNPV ;
166     fU1Vals[2] = fU ;
167     fU1Vals[3] = fUPhi ;
168     fU1Vals[4] = fTKSumEt ;
169     fU1Vals[5] = fTKU ;
170     fU1Vals[6] = fTKUPhi ;
171     fU1Vals[7] = fNPSumEt ;
172     fU1Vals[8] = fNPU ;
173     fU1Vals[9] = fNPUPhi ;
174     fU1Vals[10] = fPUSumEt ;
175     fU1Vals[11] = fPUMet ;
176     fU1Vals[12] = fPUMetPhi;
177     fU1Vals[13] = fPCSumEt ;
178     fU1Vals[14] = fPCU ;
179     fU1Vals[15] = fPCUPhi ;
180     fU1Vals[16] = fJSPt1 ;
181     fU1Vals[17] = fJSEta1 ;
182     fU1Vals[18] = fJSPhi1 ;
183     fU1Vals[19] = fJSPt2 ;
184     fU1Vals[20] = fJSEta2 ;
185     fU1Vals[21] = fJSPhi2 ;
186     fU1Vals[22] = fNAllJet ;
187     fU1Vals[23] = fNJet ;
188     fU1Vals[24] = fUPhiMVA ;
189     return fU1Vals;
190     }
191     //--------------------------------------------------------------------------------------------------
192 pharris 1.2 Double_t MVAMet::evaluatePhi() {
193     fPhiVals[0] = fNPV ;
194     fPhiVals[1] = fU ;
195     fPhiVals[2] = fUPhi ;
196     fPhiVals[3] = fTKSumEt ;
197     fPhiVals[4] = fTKU ;
198     fPhiVals[5] = fTKUPhi ;
199     fPhiVals[6] = fNPSumEt ;
200     fPhiVals[7] = fNPU ;
201     fPhiVals[8] = fNPUPhi ;
202     fPhiVals[9] = fPUSumEt ;
203     fPhiVals[10] = fPUMet ;
204     fPhiVals[11] = fPUMetPhi;
205     fPhiVals[12] = fPCSumEt ;
206     fPhiVals[13] = fPCU ;
207     fPhiVals[14] = fPCUPhi ;
208     fPhiVals[15] = fJSPt1 ;
209     fPhiVals[16] = fJSEta1 ;
210     fPhiVals[17] = fJSPhi1 ;
211     fPhiVals[18] = fJSPt2 ;
212     fPhiVals[19] = fJSEta2 ;
213     fPhiVals[20] = fJSPhi2 ;
214     fPhiVals[21] = fNAllJet ;
215     fPhiVals[22] = fNJet ;
216     return fPhiReader->GetResponse(fPhiVals);
217     }
218     //--------------------------------------------------------------------------------------------------
219     Double_t MVAMet::evaluateU1() {
220 pharris 1.5 fU1Vals[0] = fSumEt ;
221     fU1Vals[1] = fNPV ;
222     fU1Vals[2] = fU ;
223     fU1Vals[3] = fUPhi ;
224     fU1Vals[4] = fTKSumEt ;
225     fU1Vals[5] = fTKU ;
226     fU1Vals[6] = fTKUPhi ;
227     fU1Vals[7] = fNPSumEt ;
228     fU1Vals[8] = fNPU ;
229     fU1Vals[9] = fNPUPhi ;
230     fU1Vals[10] = fPUSumEt ;
231     fU1Vals[11] = fPUMet ;
232     fU1Vals[12] = fPUMetPhi;
233     fU1Vals[13] = fPCSumEt ;
234     fU1Vals[14] = fPCU ;
235     fU1Vals[15] = fPCUPhi ;
236     fU1Vals[16] = fJSPt1 ;
237     fU1Vals[17] = fJSEta1 ;
238     fU1Vals[18] = fJSPhi1 ;
239     fU1Vals[19] = fJSPt2 ;
240     fU1Vals[20] = fJSEta2 ;
241     fU1Vals[21] = fJSPhi2 ;
242     fU1Vals[22] = fNAllJet ;
243     fU1Vals[23] = fNJet ;
244     fU1Vals[24] = fUPhiMVA ;
245 pharris 1.2 return fU1Reader->GetResponse(fU1Vals);
246     }
247 pharris 1.1 //--------------------------------------------------------------------------------------------------
248 pharris 1.4 Double_t MVAMet::evaluateCovU1() {
249 pharris 1.5 fCovVals[0] = fSumEt ;
250     fCovVals[1] = fNPV ;
251     fCovVals[2] = fU ;
252     fCovVals[3] = fUPhi ;
253     fCovVals[4] = fTKSumEt ;
254     fCovVals[5] = fTKU ;
255     fCovVals[6] = fTKUPhi ;
256     fCovVals[7] = fNPSumEt ;
257     fCovVals[8] = fNPU ;
258     fCovVals[9] = fNPUPhi ;
259     fCovVals[10] = fPUSumEt ;
260     fCovVals[11] = fPUMet ;
261     fCovVals[12] = fPUMetPhi;
262     fCovVals[13] = fPCSumEt ;
263     fCovVals[14] = fPCU ;
264     fCovVals[15] = fPCUPhi ;
265     fCovVals[16] = fJSPt1 ;
266     fCovVals[17] = fJSEta1 ;
267     fCovVals[18] = fJSPhi1 ;
268     fCovVals[19] = fJSPt2 ;
269     fCovVals[20] = fJSEta2 ;
270     fCovVals[21] = fJSPhi2 ;
271     fCovVals[22] = fNAllJet ;
272     fCovVals[23] = fNJet ;
273     fCovVals[24] = fUPhiMVA ;
274     fCovVals[25] = fUMVA ;
275 pharris 1.28 double lCovU1 = fCovU1Reader->GetResponse(fCovVals);
276     if(!f42) lCovU1 = lCovU1*lCovU1*fUMVA*fUMVA;
277     return lCovU1;
278 pharris 1.4 }
279     //--------------------------------------------------------------------------------------------------
280     Double_t MVAMet::evaluateCovU2() {
281 pharris 1.5 fCovVals[0] = fSumEt ;
282     fCovVals[1] = fNPV ;
283     fCovVals[2] = fU ;
284     fCovVals[3] = fUPhi ;
285     fCovVals[4] = fTKSumEt ;
286     fCovVals[5] = fTKU ;
287     fCovVals[6] = fTKUPhi ;
288     fCovVals[7] = fNPSumEt ;
289     fCovVals[8] = fNPU ;
290     fCovVals[9] = fNPUPhi ;
291     fCovVals[10] = fPUSumEt ;
292     fCovVals[11] = fPUMet ;
293     fCovVals[12] = fPUMetPhi;
294     fCovVals[13] = fPCSumEt ;
295     fCovVals[14] = fPCU ;
296     fCovVals[15] = fPCUPhi ;
297     fCovVals[16] = fJSPt1 ;
298     fCovVals[17] = fJSEta1 ;
299     fCovVals[18] = fJSPhi1 ;
300     fCovVals[19] = fJSPt2 ;
301     fCovVals[20] = fJSEta2 ;
302     fCovVals[21] = fJSPhi2 ;
303     fCovVals[22] = fNAllJet ;
304     fCovVals[23] = fNJet ;
305     fCovVals[24] = fUPhiMVA ;
306     fCovVals[25] = fUMVA ;
307 pharris 1.28 double lCovU2 = fCovU2Reader->GetResponse(fCovVals);
308     if(!f42) lCovU2 = lCovU2*lCovU2*fUMVA*fUMVA;
309     return lCovU2;
310 pharris 1.4 }
311     //--------------------------------------------------------------------------------------------------
312 pharris 1.1 Double_t MVAMet::MVAValue( Bool_t iPhi,
313     Float_t iPFSumEt,
314     Float_t iU ,
315     Float_t iUPhi ,
316     Float_t iTKSumEt,
317     Float_t iTKU ,
318     Float_t iTKUPhi ,
319     Float_t iNPSumEt,
320     Float_t iNPU ,
321     Float_t iNPUPhi ,
322     Float_t iPUSumEt,
323     Float_t iPUMet ,
324     Float_t iPUMetPhi,
325     Float_t iPCSumEt,
326     Float_t iPCU ,
327     Float_t iPCUPhi ,
328     Float_t iJSPt1 ,
329     Float_t iJSEta1 ,
330     Float_t iJSPhi1 ,
331     Float_t iJSPt2 ,
332     Float_t iJSEta2 ,
333     Float_t iJSPhi2 ,
334     Float_t iNAllJet,
335     Float_t iNJet ,
336     Float_t iNPV
337     ){
338    
339     if (!fIsInitialized) {
340     std::cout << "Error: MVA Met not properly initialized.\n";
341     return -9999;
342     }
343 pharris 1.5
344     fSumEt = iPFSumEt;
345 pharris 1.1 fU = iU ;
346     fUPhi = iUPhi ;
347     fTKSumEt = iTKSumEt/iPFSumEt;
348     fTKU = iTKU ;
349     fTKUPhi = iTKUPhi ;
350     fNPSumEt = iNPSumEt/iPFSumEt;
351     fNPU = iNPU ;
352     fNPUPhi = iNPUPhi ;
353     fPUSumEt = iPUSumEt/iPFSumEt;
354     fPUMet = iPUMet ;
355     fPUMetPhi = iPUMetPhi;
356     fPCSumEt = iPCSumEt/iPFSumEt;
357     fPCU = iPCU ;
358     fPCUPhi = iPCUPhi ;
359     fJSPt1 = iJSPt1 ;
360     fJSEta1 = iJSEta1 ;
361     fJSPhi1 = iJSPhi1 ;
362     fJSPt2 = iJSPt2 ;
363     fJSEta2 = iJSEta2 ;
364     fJSPhi2 = iJSPhi2 ;
365     fNAllJet = iNAllJet;
366     fNJet = iNJet ;
367     fNPV = iNPV ;
368 pharris 1.5
369 pharris 1.1 Double_t lMVA = -9999;
370 pharris 1.2 lMVA = evaluatePhi();
371 pharris 1.1 if(!iPhi) fUPhiMVA = iUPhi+lMVA;
372 pharris 1.2 //Not no nice feature of the training
373 pharris 1.4 //fTKSumEt /= iPFSumEt;
374     //fNPSumEt /= iPFSumEt;
375     //fPUSumEt /= iPFSumEt;
376     //fPCSumEt /= iPFSumEt;
377 pharris 1.2 if(!iPhi) lMVA = evaluateU1();
378 pharris 1.1 return lMVA;
379     }
380     //--------------------------------------------------------------------------------------------------
381     //====> Please not that the jet collection must be cleaned => all jets near leptons must be removed
382 pharris 1.23 Met MVAMet::GetMet( Bool_t iPhi,
383     Float_t iPtVis,Float_t iPhiVis,Float_t iSumEtVis,
384     Float_t iPtQ ,Float_t iPhiQ ,Float_t iSumEtQ,
385 pharris 1.1 const PFMet *iMet ,
386 pharris 1.3 const PFCandidateCol *iCands,
387     const Vertex *iVertex,const VertexCol *iVertices,
388 pharris 1.1 const PFJetCol *iJets ,
389     FactorizedJetCorrector *iJetCorrector,
390     const PileupEnergyDensityCol *iPUEnergyDensity,
391     int iNPV,
392     Bool_t printDebug) {
393 pharris 1.14
394     int lNPV = 0;
395     for(unsigned int i0 = 0; i0 < iVertices->GetEntries(); i0++) {
396     const Vertex *lPV = iVertices->At(i0);
397     if(lPV->Ndof() < 4.0) continue;
398     if(fabs(lPV->Z()) > 24.) continue;
399     if(lPV->Position().Rho() > 2.) continue;
400     lNPV++;
401     }
402 pharris 1.1
403 pharris 1.7 Met lPFRec = fRecoilTools->pfRecoil (iPtVis,iPhiVis,iSumEtVis,iCands);
404 pharris 1.23 Met lTKRec = fRecoilTools->trackRecoil(iPtQ ,iPhiQ ,iSumEtQ ,iCands,iVertex);
405     Met lNPRec = fRecoilTools->NoPURecoil (iPtQ ,iPhiQ ,iSumEtQ ,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices);
406 pharris 1.3 Met lPCRec = fRecoilTools->PUCRecoil (iPtVis,iPhiVis,iSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices);
407     Met lPUMet = fRecoilTools->PUMet (iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices);
408 pharris 1.1
409 pharris 1.23
410 pharris 1.1 Double_t lPt0 = 0; const PFJet *lLead = 0;
411     Double_t lPt1 = 0; const PFJet *l2nd = 0;
412     int lNAllJet = 0;
413     int lNJet = 0;
414     for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
415     const PFJet *pJet = iJets->At(i0);
416     Double_t pPt = fRecoilTools->fJetIDMVA->correctedPt(pJet,iJetCorrector,iPUEnergyDensity);
417 pharris 1.24 //if(pPt < 0) continue;
418 pharris 1.1 if(!JetTools::passPFLooseId(pJet)) continue;
419 pharris 1.16 if(f42 && pPt > 1.) lNAllJet++;
420     if(!f42) lNAllJet++;
421 pharris 1.1 if(pPt > 30) lNJet++;
422 pharris 1.24 if(lPt0 < pPt) {lPt1 = lPt0; lPt0 = pPt; l2nd = lLead; lLead = pJet; continue;}
423     if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
424 pharris 1.1 }
425 pharris 1.5
426     fSumEt = lPFRec.SumEt();
427 pharris 1.1 fU = lPFRec.Pt();
428     fUPhi = lPFRec.Phi();
429     fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
430     fTKU = lTKRec.Pt();
431     fTKUPhi = lTKRec.Phi();
432     fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
433     fNPU = lNPRec.Pt();
434     fNPUPhi = lNPRec.Phi();
435     fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
436     fPUMet = lPUMet.Pt();
437     fPUMetPhi = lPUMet.Phi();
438     fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
439     fPCU = lPCRec.Pt() ;
440     fPCUPhi = lPCRec.Phi() ;
441     fJSPt1 = lPt0;
442     fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
443     fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
444     fJSPt2 = lPt1;
445     fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
446     fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
447     fNJet = lNJet ;
448     fNAllJet = lNAllJet;
449 pharris 1.14 fNPV = lNPV ;
450 pharris 1.5 Float_t lMVA = evaluatePhi();
451     if(!iPhi) fUPhiMVA = fUPhi + lMVA;
452 pharris 1.1 //Not no nice feature of teh training
453 pharris 1.4 //fTKSumEt /= lPFRec.SumEt();
454     //fNPSumEt /= lPFRec.SumEt();
455     //fPUSumEt /= lPFRec.SumEt();
456     //fPCSumEt /= lPFRec.SumEt();
457 pharris 1.7 if(!iPhi) lMVA = evaluateU1();//fU1Reader ->EvaluateMVA( fU1MethodName );
458     fUMVA = fU*lMVA;
459    
460     TLorentzVector lUVec(0,0,0,0); lUVec .SetPtEtaPhiM(fUMVA,0,fUPhiMVA,0);
461 pharris 1.1 TLorentzVector lVVec(0,0,0,0); lVVec .SetPtEtaPhiM(iPtVis ,0,iPhiVis ,0);
462     if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
463     lUVec -= lVVec;
464     Met lMet(lUVec.Px(),lUVec.Py());
465 pharris 1.5
466 pharris 1.4 //TMatrixD lCov(2,2);
467     //Covariance matrix perpendicular and parallel to the recoil (sort of correct)
468     double lCovU1 = evaluateCovU1();
469     double lCovU2 = evaluateCovU2();
470    
471     double lCos2 = cos(fUPhiMVA)*cos(fUPhiMVA);
472     double lSin2 = sin(fUPhiMVA)*sin(fUPhiMVA);
473 pharris 1.5
474 pharris 1.4 //Now Compute teh covariance matrix in X and Y
475 pharris 1.5 (*fCov)(0,0) = lCovU1*lCos2+lCovU2*lSin2;
476     (*fCov)(1,0) = -lCovU1*sin(fUPhiMVA)*cos(fUPhiMVA)+lCovU2*sin(fUPhiMVA)*cos(fUPhiMVA);
477     (*fCov)(0,1) = (*fCov)(1,0);
478     (*fCov)(1,1) = lCovU1*lSin2+lCovU2*lCos2;
479 pharris 1.10 TMatrixD lInv(2,2);
480     lInv(0,0) = (*fCov)(0,0); lInv(1,1) = (*fCov)(1,1); lInv(1,0) = (*fCov)(1,0); lInv(0,1) = (*fCov)(0,1);
481     if(lInv.Determinant() != 0) lInv.Invert();
482     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));
483     fUncertainty = sqrt(lCovU1+lCovU2);
484 pharris 1.7
485 pharris 1.1 if (printDebug == kTRUE) {
486 ceballos 1.11 std::cout << "Debug Met MVA: "
487 pharris 1.1 << fU << " : "
488     << fUPhi << " : "
489     << fTKSumEt << " : "
490     << fTKU << " : "
491     << fTKUPhi << " : "
492     << fNPSumEt << " : "
493     << fNPU << " : "
494     << fNPUPhi << " : "
495     << fPUSumEt << " : "
496     << fPUMet << " : "
497     << fPUMetPhi << " : "
498     << fPCUPhi << " : "
499     << fPCSumEt << " : "
500     << fPCU << " : "
501     << fPCUPhi << " : "
502     << fJSPt1 << " : "
503     << fJSEta1 << " : "
504     << fJSPhi1 << " : "
505     << fJSPt2 << " : "
506     << fJSEta2 << " : "
507     << fJSPhi2 << " : "
508     << fNJet << " : "
509     << fNAllJet << " : "
510     << fNPV << " : "
511     << " === : === "
512 ceballos 1.12 << lMet.Pt() << " : "
513     << lMet.Phi() << " : "
514 pharris 1.1 << std::endl;
515     }
516    
517     return lMet;
518     }
519 pharris 1.2 //--------------------------------------------------------------------------------------------------
520     //====> Please not that the jet collection must be cleaned => all jets near leptons must be removed
521     //====> Corrected Jets
522 pharris 1.23 Met MVAMet::GetMet( Bool_t iPhi,
523     Float_t iPtVis,Float_t iPhiVis,Float_t iSumEtVis,
524     Float_t iPtQ ,Float_t iPhiQ ,Float_t iSumEtQ , //Charged components
525 pharris 1.2 const PFMet *iMet ,
526 pharris 1.3 const PFCandidateCol *iCands,
527 pharris 1.7 const Vertex *iVertex ,const VertexCol *iVertices,Double_t iRho,
528 pharris 1.2 const PFJetCol *iJets ,
529     int iNPV,
530     Bool_t printDebug) {
531 pharris 1.14
532     int lNPV = 0;
533     for(unsigned int i0 = 0; i0 < iVertices->GetEntries(); i0++) {
534     const Vertex *lPV = iVertices->At(i0);
535     if(lPV->Ndof() < 4.0) continue;
536     if(fabs(lPV->Z()) > 24.) continue;
537     if(lPV->Position().Rho() > 2.) continue;
538     lNPV++;
539     }
540 pharris 1.7 Met lPFRec = fRecoilTools->pfRecoil (iPtVis,iPhiVis,iSumEtVis,iCands);
541 pharris 1.23 Met lTKRec = fRecoilTools->trackRecoil(iPtQ ,iPhiQ ,iSumEtQ , iCands,iVertex);
542     Met lNPRec = fRecoilTools->NoPURecoil (iPtQ ,iPhiQ ,iSumEtQ ,iJets,iCands,iVertex,iVertices,iRho);
543 pharris 1.7 Met lPCRec = fRecoilTools->PUCRecoil (iPtVis,iPhiVis,iSumEtVis,iJets,iCands,iVertex,iVertices,iRho);
544     Met lPUMet = fRecoilTools->PUMet ( iJets,iCands,iVertex,iVertices,iRho);
545 pharris 1.2
546     Double_t lPt0 = 0; const PFJet *lLead = 0;
547     Double_t lPt1 = 0; const PFJet *l2nd = 0;
548     int lNAllJet = 0;
549     int lNJet = 0;
550     for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
551     const PFJet *pJet = iJets->At(i0);
552     Double_t pPt = pJet->Pt();
553     if(!JetTools::passPFLooseId(pJet)) continue;
554 pharris 1.16 if(f42 && pPt > 1.) lNAllJet++;
555     if(!f42) lNAllJet++;
556 pharris 1.24 if(pPt > 30.) lNJet++;
557     if(lPt0 < pPt) {lPt1 = lPt0; lPt0 = pPt; l2nd = lLead; lLead = pJet; continue;}
558     if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
559 pharris 1.2 }
560    
561 pharris 1.5 fSumEt = lPFRec.SumEt();
562 pharris 1.2 fU = lPFRec.Pt();
563     fUPhi = lPFRec.Phi();
564     fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
565     fTKU = lTKRec.Pt();
566     fTKUPhi = lTKRec.Phi();
567     fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
568     fNPU = lNPRec.Pt();
569     fNPUPhi = lNPRec.Phi();
570     fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
571     fPUMet = lPUMet.Pt();
572     fPUMetPhi = lPUMet.Phi();
573     fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
574     fPCU = lPCRec.Pt() ;
575     fPCUPhi = lPCRec.Phi() ;
576     fJSPt1 = lPt0;
577     fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
578     fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
579     fJSPt2 = lPt1;
580     fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
581     fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
582     fNJet = lNJet ;
583     fNAllJet = lNAllJet;
584 pharris 1.14 fNPV = lNPV ;
585 pharris 1.7
586 pharris 1.2 Float_t lMVA = evaluatePhi();
587    
588     if(!iPhi) fUPhiMVA = fUPhi + lMVA;
589     //Not no nice feature of teh training
590 pharris 1.4 //fTKSumEt /= lPFRec.SumEt();
591     //fNPSumEt /= lPFRec.SumEt();
592     //fPUSumEt /= lPFRec.SumEt();
593     //fPCSumEt /= lPFRec.SumEt();
594 pharris 1.7 if(!iPhi) lMVA = evaluateU1();//fU1Reader ->EvaluateMVA( fU1MethodName );
595     fUMVA = lMVA*fU;
596    
597     TLorentzVector lUVec(0,0,0,0); lUVec .SetPtEtaPhiM(fUMVA ,0,fUPhiMVA,0);
598 pharris 1.2 TLorentzVector lVVec(0,0,0,0); lVVec .SetPtEtaPhiM(iPtVis ,0,iPhiVis ,0);
599     if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
600     lUVec -= lVVec;
601     Met lMet(lUVec.Px(),lUVec.Py());
602 pharris 1.4 //Cov matrix
603     double lCovU1 = evaluateCovU1();
604     double lCovU2 = evaluateCovU2();
605    
606     double lCos2 = cos(fUPhiMVA)*cos(fUPhiMVA);
607     double lSin2 = sin(fUPhiMVA)*sin(fUPhiMVA);
608    
609     //Now Compute teh covariance matrix in X and Y
610 pharris 1.5 (*fCov)(0,0) = lCovU1*lCos2+lCovU2*lSin2;
611     (*fCov)(1,0) = -lCovU1*sin(fUPhiMVA)*cos(fUPhiMVA)+lCovU2*sin(fUPhiMVA)*cos(fUPhiMVA);
612     (*fCov)(0,1) = (*fCov)(1,0);
613     (*fCov)(1,1) = lCovU1*lSin2+lCovU2*lCos2;
614 pharris 1.10 TMatrixD lInv(2,2);
615     lInv(0,0) = (*fCov)(0,0); lInv(1,1) = (*fCov)(1,1); lInv(1,0) = (*fCov)(1,0); lInv(0,1) = (*fCov)(0,1);
616     if(lInv.Determinant() != 0) lInv.Invert();
617     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));
618     fUncertainty = sqrt(lCovU1+lCovU2);
619 pharris 1.2
620     if (printDebug == kTRUE) {
621 ceballos 1.11 std::cout << "Debug Met MVA: "
622 pharris 1.2 << fU << " : "
623     << fUPhi << " : "
624     << fTKSumEt << " : "
625     << fTKU << " : "
626     << fTKUPhi << " : "
627     << fNPSumEt << " : "
628     << fNPU << " : "
629     << fNPUPhi << " : "
630     << fPUSumEt << " : "
631     << fPUMet << " : "
632     << fPUMetPhi << " : "
633     << fPCUPhi << " : "
634     << fPCSumEt << " : "
635     << fPCU << " : "
636     << fPCUPhi << " : "
637     << fJSPt1 << " : "
638     << fJSEta1 << " : "
639     << fJSPhi1 << " : "
640     << fJSPt2 << " : "
641     << fJSEta2 << " : "
642     << fJSPhi2 << " : "
643     << fNJet << " : "
644     << fNAllJet << " : "
645     << fNPV << " : "
646     << " === : === "
647 ceballos 1.12 << lMet.Pt() << " : "
648     << lMet.Phi() << " : "
649 pharris 1.2 << std::endl;
650     }
651    
652     return lMet;
653     }
654     //--------------------------------------------------------------------------------------------------
655 pharris 1.1 //Interms of the two candidates
656 pharris 1.23 Met MVAMet::GetMet( bool iPhi,
657     Float_t iPt1,Float_t iPhi1,Float_t iEta1,Float_t iChargedFrac1,
658     Float_t iPt2,Float_t iPhi2,Float_t iEta2,Float_t iChargedFrac2,
659 pharris 1.1 const PFMet *iMet ,
660 pharris 1.3 const PFCandidateCol *iCands,
661     const Vertex *iVertex,const VertexCol *iVertices,
662 pharris 1.1 const PFJetCol *iJets ,
663     FactorizedJetCorrector *iJetCorrector,
664     const PileupEnergyDensityCol *iPUEnergyDensity,
665     int iNPV,
666     Bool_t printDebug) {
667 pharris 1.14
668     int lNPV = 0;
669     for(unsigned int i0 = 0; i0 < iVertices->GetEntries(); i0++) {
670     const Vertex *lPV = iVertices->At(i0);
671     if(lPV->Ndof() < 4.0) continue;
672     if(fabs(lPV->Z()) > 24.) continue;
673     if(lPV->Position().Rho() > 2.) continue;
674     lNPV++;
675     }
676 pharris 1.23 //Met lMVec1 = fRecoilTools->pfCone(iPhi1,iEta1,iCands);
677     //Met lMVec2 = fRecoilTools->pfCone(iPhi2,iEta2,iCands);
678     //double lChargedFrac1 = 1.;
679     //double lChargedFrac2 = 1.;
680     //if(iId == MVAMet::TauTau){
681     // lChargedFrac1 = chargedFracTau(
682     // }
683    
684     TLorentzVector lVVec1 (0,0,0,0); lVVec1 .SetPtEtaPhiM(iPt1,0 ,iPhi1,0);
685     TLorentzVector lVVec2 (0,0,0,0); lVVec2 .SetPtEtaPhiM(iPt2,0 ,iPhi2,0);
686     TLorentzVector lVVec1Q (0,0,0,0); lVVec1Q .SetPtEtaPhiM(iPt1*iChargedFrac1,0,iPhi1,0);
687     TLorentzVector lVVec2Q (0,0,0,0); lVVec2Q .SetPtEtaPhiM(iPt2*iChargedFrac2,0,iPhi2,0);
688    
689 pharris 1.30 if(iChargedFrac1 < 0 ) {
690     Met lQ1Cone = fRecoilTools->pfCone(iPhi1,iEta1,iCands,iVertex,true);
691     lVVec1Q.SetPtEtaPhiM(lQ1Cone.Pt(),0,lQ1Cone.Phi(),0);
692     }
693     if(iChargedFrac2 < 0 ) {
694     Met lQ2Cone = fRecoilTools->pfCone(iPhi2,iEta2,iCands,iVertex,true);
695     lVVec2Q.SetPtEtaPhiM(lQ2Cone.Pt(),0,lQ2Cone.Phi(),0);
696     }
697 pharris 1.23 //Met lPFRec = fRecoilTools->pfRecoil (iPhi1,iEta1,iPhi2,iEta2, iCands);
698     //Met lTKRec = fRecoilTools->trackRecoil(iPhi1,iEta1,iPhi2,iEta2, iCands,iVertex);
699     //Met lNPRec = fRecoilTools->NoPURecoil (iPhi1,iEta1,iPhi2,iEta2,iJets,iCands,iVertex,iVertices,iJetCorrector,iPUEnergyDensity);
700     //Met lPCRec = fRecoilTools->PUCRecoil (iPhi1,iEta1,iPhi2,iEta2,iJets,iCands,iVertex,iVertices,iJetCorrector,iPUEnergyDensity);
701     //Met lPUMet = fRecoilTools->PUMet (iPhi1,iEta1,iPhi2,iEta2,iJets,iCands,iVertex,iVertices,iJetCorrector,iPUEnergyDensity);
702    
703     //Met lNPMet = fRecoilTools->NoPUMet (iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2);
704     //std::cout << "===> " << lNPMet.Pt() << " -- " << lNPMet.Phi() << std::endl;
705     //std::cout << "===> " << lVVec1.Pt() << " -- " << lVVec1.Phi() << " -- " << iPt1 << " -- " << iPhi1 << std::endl;
706     //std::cout << "===> " << lVVec2.Pt() << " -- " << lVVec2.Phi() << " -- " << iPt2 << " -- " << iPhi2 << std::endl;
707     //std::cout << "===> " << lNPRec.Pt() << " -- " << lNPRec.Phi() << std::endl;
708 ceballos 1.22
709 pharris 1.23 Float_t lSumEtVis = lVVec1.Pt() + lVVec2.Pt();
710 pharris 1.1 lVVec1+=lVVec2;
711 pharris 1.23 Float_t lPtVis = lVVec1.Pt();
712     Float_t lPhiVis = lVVec1.Phi();
713    
714     Float_t lSumEtQ = lVVec1Q.Pt() + lVVec2Q.Pt();
715     lVVec1Q+=lVVec2Q;
716     Float_t lPtVisQ = lVVec1Q.Pt();
717     Float_t lPhiVisQ = lVVec1Q.Phi();
718    
719     Met lPFRec = fRecoilTools->pfRecoil (lPtVis ,lPhiVis, lSumEtVis,iCands);
720     Met lTKRec = fRecoilTools->trackRecoil(lPtVisQ ,lPhiVisQ ,lSumEtQ ,iCands,iVertex);
721     Met lNPRec = fRecoilTools->NoPURecoil (lPtVisQ ,lPhiVisQ ,lSumEtQ ,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2);
722     Met lPCRec = fRecoilTools->PUCRecoil (lPtVis ,lPhiVis ,lSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2);
723     Met lPUMet = fRecoilTools->PUMet ( iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2);
724 pharris 1.1
725 pharris 1.25 //Met lNPRec1 = fRecoilTools->NoPUMet (iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices,0,0,0,0);
726 pharris 1.19 //std::cout << " ====> PF => " << lPFRec.Pt() << " -- " << lPFRec1.Pt() << std::endl;
727     //std::cout << " ====> TK => " << lTKRec.Pt() << " -- " << lTKRec1.Pt() << std::endl;
728     //std::cout << " ====> NP => " << lNPRec.Pt() << " -- " << lNPRec1.Pt() << std::endl;
729     //std::cout << " ====> PC => " << lPCRec.Pt() << " -- " << lPCRec1.Pt() << std::endl;
730    
731 pharris 1.1 Double_t lPt0 = 0; const PFJet *lLead = 0;
732     Double_t lPt1 = 0; const PFJet *l2nd = 0;
733     int lNAllJet = 0;
734     int lNJet = 0;
735     for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
736     const PFJet *pJet = iJets->At(i0);
737     Double_t pPt = fRecoilTools->fJetIDMVA->correctedPt(pJet,iJetCorrector,iPUEnergyDensity);
738     double pDEta1 = pJet->Eta() - iEta1;
739     double pDPhi1 = fabs(pJet->Phi() - iPhi1); if(pDPhi1 > 2.*TMath::Pi()-pDPhi1) pDPhi1 = 2.*TMath::Pi()-pDPhi1;
740     double pDR1 = sqrt(pDEta1*pDEta1 + pDPhi1*pDPhi1);
741     if(pDR1 < 0.5) continue;
742     double pDEta2 = pJet->Eta() - iEta2;
743     double pDPhi2 = fabs(pJet->Phi() - iPhi2); if(pDPhi2 > 2.*TMath::Pi()-pDPhi2) pDPhi2 = 2.*TMath::Pi()-pDPhi2;
744     double pDR2 = sqrt(pDEta2*pDEta2 + pDPhi2*pDPhi2);
745     if(pDR2 < 0.5) continue;
746     if(!JetTools::passPFLooseId(pJet)) continue;
747 pharris 1.16 if(f42 && pPt > 1.) lNAllJet++;
748     if(!f42) lNAllJet++;
749 pharris 1.1 if(pPt > 30.) lNJet++;
750 pharris 1.24 if(lPt0 < pPt) {lPt1 = lPt0; lPt0 = pPt; l2nd = lLead; lLead = pJet; continue;}
751     if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
752 pharris 1.1 }
753    
754 pharris 1.5 fSumEt = lPFRec.SumEt();
755 pharris 1.1 fU = lPFRec.Pt();
756     fUPhi = lPFRec.Phi();
757     fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
758     fTKU = lTKRec.Pt();
759     fTKUPhi = lTKRec.Phi();
760     fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
761     fNPU = lNPRec.Pt();
762     fNPUPhi = lNPRec.Phi();
763     fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
764     fPUMet = lPUMet.Pt();
765     fPUMetPhi= lPUMet.Phi();
766     fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
767     fPCU = lPCRec.Pt() ;
768     fPCUPhi = lPCRec.Phi() ;
769     fJSPt1 = lPt0;
770     fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
771     fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
772     fJSPt2 = lPt1;
773     fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
774     fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
775     fNJet = lNJet ;
776     fNAllJet = lNAllJet;
777 pharris 1.14 fNPV = lNPV ;
778 pharris 1.5
779 pharris 1.2 Float_t lMVA = evaluatePhi();//fPhiReader ->EvaluateMVA( fPhiMethodName );
780     if(!iPhi) fUPhiMVA = fUPhi + lMVA;
781 pharris 1.4 //fTKSumEt /= lPFRec.SumEt();
782     //fNPSumEt /= lPFRec.SumEt();
783     //fPUSumEt /= lPFRec.SumEt();
784     //fPCSumEt /= lPFRec.SumEt();
785 pharris 1.7 if(!iPhi) lMVA = evaluateU1();//fU1Reader ->EvaluateMVA( fU1MethodName );
786     fUMVA = fU*lMVA;
787    
788     TLorentzVector lUVec (0,0,0,0); lUVec .SetPtEtaPhiM(fUMVA ,0,fUPhiMVA,0);
789 pharris 1.2 TLorentzVector lVVec (0,0,0,0); lVVec .SetPtEtaPhiM(lPtVis ,0,lPhiVis ,0);
790     if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
791     lUVec -= lVVec;
792     Met lMet(lUVec.Px(),lUVec.Py());
793 pharris 1.4 //Cov matrix
794     double lCovU1 = evaluateCovU1();
795     double lCovU2 = evaluateCovU2();
796    
797     double lCos2 = cos(fUPhiMVA)*cos(fUPhiMVA);
798     double lSin2 = sin(fUPhiMVA)*sin(fUPhiMVA);
799    
800     //Now Compute teh covariance matrix in X and Y
801 pharris 1.5 (*fCov)(0,0) = lCovU1*lCos2+lCovU2*lSin2;
802     (*fCov)(1,0) = -lCovU1*sin(fUPhiMVA)*cos(fUPhiMVA)+lCovU2*sin(fUPhiMVA)*cos(fUPhiMVA);
803     (*fCov)(0,1) = (*fCov)(1,0);
804     (*fCov)(1,1) = lCovU1*lSin2+lCovU2*lCos2;
805 pharris 1.10 TMatrixD lInv(2,2);
806     lInv(0,0) = (*fCov)(0,0); lInv(1,1) = (*fCov)(1,1); lInv(1,0) = (*fCov)(1,0); lInv(0,1) = (*fCov)(0,1);
807     if(lInv.Determinant() != 0) lInv.Invert();
808     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));
809     fUncertainty = sqrt(lCovU1+lCovU2);
810    
811 pharris 1.2
812     if (printDebug == kTRUE) {
813 ceballos 1.11 std::cout << "Debug Met MVA: "
814 pharris 1.14 << fSumEt << " : "
815 pharris 1.2 << fU << " : "
816     << fUPhi << " : "
817     << fTKSumEt << " : "
818     << fTKU << " : "
819     << fTKUPhi << " : "
820     << fNPSumEt << " : "
821     << fNPU << " : "
822     << fNPUPhi << " : "
823     << fPUSumEt << " : "
824     << fPUMet << " : "
825     << fPUMetPhi << " : "
826     << fPCSumEt << " : "
827     << fPCU << " : "
828     << fPCUPhi << " : "
829     << fJSPt1 << " : "
830     << fJSEta1 << " : "
831     << fJSPhi1 << " : "
832     << fJSPt2 << " : "
833     << fJSEta2 << " : "
834     << fJSPhi2 << " : "
835     << fNJet << " : "
836     << fNAllJet << " : "
837     << fNPV << " : "
838     << " === : === "
839 ceballos 1.12 << lMet.Pt() << " : "
840     << lMet.Phi() << " : "
841 pharris 1.2 << std::endl;
842     }
843    
844     return lMet;
845     }
846     //--------------------------------------------------------------------------------------------------
847     //Interms of the two candidates => corrected Jets
848     Met MVAMet::GetMet( Bool_t iPhi,
849 pharris 1.23 Float_t iPt1,Float_t iPhi1,Float_t iEta1,Float_t iChargedFrac1,
850     Float_t iPt2,Float_t iPhi2,Float_t iEta2,Float_t iChargedFrac2,
851 pharris 1.2 const PFMet *iMet ,
852 pharris 1.3 const PFCandidateCol *iCands,
853 pharris 1.7 const Vertex *iVertex,const VertexCol *iVertices,Double_t iRho,
854 pharris 1.2 const PFJetCol *iJets ,
855     int iNPV,
856     Bool_t printDebug) {
857 pharris 1.14 int lNPV = 0;
858     for(unsigned int i0 = 0; i0 < iVertices->GetEntries(); i0++) {
859     const Vertex *lPV = iVertices->At(i0);
860     if(lPV->Ndof() < 4.0) continue;
861     if(fabs(lPV->Z()) > 24.) continue;
862     if(lPV->Position().Rho() > 2.) continue;
863     lNPV++;
864     }
865 pharris 1.2
866     TLorentzVector lVVec1(0,0,0,0); lVVec1.SetPtEtaPhiM(iPt1,0,iPhi1 ,0);
867     TLorentzVector lVVec2(0,0,0,0); lVVec2.SetPtEtaPhiM(iPt2,0,iPhi2 ,0);
868 pharris 1.23 TLorentzVector lVVec1Q (0,0,0,0); lVVec1Q .SetPtEtaPhiM(iPt1*iChargedFrac1,0,iPhi1,0);
869     TLorentzVector lVVec2Q (0,0,0,0); lVVec2Q .SetPtEtaPhiM(iPt2*iChargedFrac2,0,iPhi2,0);
870    
871 pharris 1.30 if(iChargedFrac1 < 0 ) {
872     Met lQ1Cone = fRecoilTools->pfCone(iPhi1,iEta1,iCands,iVertex,true);
873     lVVec1Q.SetPtEtaPhiM(lQ1Cone.Pt(),0,lQ1Cone.Phi(),0);
874     }
875     if(iChargedFrac2 < 0 ) {
876     Met lQ2Cone = fRecoilTools->pfCone(iPhi2,iEta2,iCands,iVertex,true);
877     lVVec2Q.SetPtEtaPhiM(lQ2Cone.Pt(),0,lQ2Cone.Phi(),0);
878     }
879    
880 pharris 1.23 Float_t lSumEtVis = lVVec1.Pt() + lVVec2.Pt();
881 pharris 1.2 lVVec1+=lVVec2;
882 pharris 1.23 Float_t lPtVis = lVVec1.Pt();
883     Float_t lPhiVis = lVVec1.Phi();
884 pharris 1.19
885 pharris 1.23 Float_t lSumEtQ = lVVec1Q.Pt() + lVVec2Q.Pt();
886     lVVec1Q+=lVVec2Q;
887     Float_t lPtVisQ = lVVec1Q.Pt();
888     Float_t lPhiVisQ = lVVec1Q.Phi();
889    
890     Met lPFRec = fRecoilTools->pfRecoil (lPtVis ,lPhiVis, lSumEtVis,iCands);
891     Met lTKRec = fRecoilTools->trackRecoil(lPtVisQ ,lPhiVisQ ,lSumEtQ ,iCands,iVertex);
892     Met lNPRec = fRecoilTools->NoPURecoil (lPtVisQ ,lPhiVisQ ,lSumEtQ ,iJets,iCands,iVertex,iVertices,1.,iPhi1,iEta1,iPhi2,iEta2);
893     Met lPCRec = fRecoilTools->PUCRecoil (lPtVis ,lPhiVis ,lSumEtVis,iJets,iCands,iVertex,iVertices,1.,iPhi1,iEta1,iPhi2,iEta2);
894     Met lPUMet = fRecoilTools->PUMet ( iJets,iCands,iVertex,iVertices,1.,iPhi1,iEta1,iPhi2,iEta2);
895    
896     //Met lPFRec = fRecoilTools->pfRecoil (iPhi1,iEta1,iPhi2,iEta2, iCands);
897     //Met lTKRec = fRecoilTools->trackRecoil(iPhi1,iEta1,iPhi2,iEta2, iCands,iVertex);
898     //Met lNPRec = fRecoilTools->NoPURecoil (iPhi1,iEta1,iPhi2,iEta2,iJets,iCands,iVertex,iVertices); //Not using rho
899     //Met lPCRec = fRecoilTools->PUCRecoil (iPhi1,iEta1,iPhi2,iEta2,iJets,iCands,iVertex,iVertices);
900     //Met lPUMet = fRecoilTools->PUMet (iPhi1,iEta1,iPhi2,iEta2,iJets,iCands,iVertex,iVertices);
901 pharris 1.2
902     Double_t lPt0 = 0; const PFJet *lLead = 0;
903     Double_t lPt1 = 0; const PFJet *l2nd = 0;
904     int lNAllJet = 0;
905     int lNJet = 0;
906     for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
907     const PFJet *pJet = iJets->At(i0);
908     Double_t pPt = pJet->Pt();
909     double pDEta1 = pJet->Eta() - iEta1;
910     double pDPhi1 = fabs(pJet->Phi() - iPhi1); if(pDPhi1 > 2.*TMath::Pi()-pDPhi1) pDPhi1 = 2.*TMath::Pi()-pDPhi1;
911     double pDR1 = sqrt(pDEta1*pDEta1 + pDPhi1*pDPhi1);
912     if(pDR1 < 0.5) continue;
913     double pDEta2 = pJet->Eta() - iEta2;
914     double pDPhi2 = fabs(pJet->Phi() - iPhi2); if(pDPhi2 > 2.*TMath::Pi()-pDPhi2) pDPhi2 = 2.*TMath::Pi()-pDPhi2;
915     double pDR2 = sqrt(pDEta2*pDEta2 + pDPhi2*pDPhi2);
916     if(pDR2 < 0.5) continue;
917     if(!JetTools::passPFLooseId(pJet)) continue;
918 pharris 1.16 if(f42 && pPt > 1.) lNAllJet++;
919     if(!f42) lNAllJet++;
920 pharris 1.2 if(pPt > 30.) lNJet++;
921 pharris 1.24 if(lPt0 < pPt) {lPt1 = lPt0; lPt0 = pPt; l2nd = lLead; lPt1 = l2nd->Pt(); lLead = pJet; continue;}
922     if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
923 pharris 1.2 }
924 pharris 1.5 fSumEt = lPFRec.SumEt();
925 pharris 1.2 fU = lPFRec.Pt();
926     fUPhi = lPFRec.Phi();
927     fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
928     fTKU = lTKRec.Pt();
929     fTKUPhi = lTKRec.Phi();
930     fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
931     fNPU = lNPRec.Pt();
932     fNPUPhi = lNPRec.Phi();
933     fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
934     fPUMet = lPUMet.Pt();
935     fPUMetPhi= lPUMet.Phi();
936     fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
937     fPCU = lPCRec.Pt() ;
938     fPCUPhi = lPCRec.Phi() ;
939     fJSPt1 = lPt0;
940     fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
941     fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
942     fJSPt2 = lPt1;
943     fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
944     fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
945     fNJet = lNJet ;
946     fNAllJet = lNAllJet;
947 pharris 1.14 fNPV = lNPV ;
948 pharris 1.2
949     Float_t lMVA = evaluatePhi();//fPhiReader ->EvaluateMVA( fPhiMethodName );
950 pharris 1.1
951     if(!iPhi) fUPhiMVA = fUPhi + lMVA;
952 pharris 1.4 //fTKSumEt /= lPFRec.SumEt();
953     //fNPSumEt /= lPFRec.SumEt();
954     //fPUSumEt /= lPFRec.SumEt();
955     //fPCSumEt /= lPFRec.SumEt();
956 pharris 1.7 if(!iPhi) lMVA = evaluateU1();//fU1Reader ->EvaluateMVA( fU1MethodName );
957     fUMVA = lMVA*fU;
958    
959     TLorentzVector lUVec (0,0,0,0); lUVec .SetPtEtaPhiM(fUMVA ,0,fUPhiMVA,0);
960 pharris 1.1 TLorentzVector lVVec (0,0,0,0); lVVec .SetPtEtaPhiM(lPtVis ,0,lPhiVis ,0);
961 pharris 1.7 if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
962 pharris 1.1 lUVec -= lVVec;
963     Met lMet(lUVec.Px(),lUVec.Py());
964 pharris 1.4 //Cov matrix
965     double lCovU1 = evaluateCovU1();
966     double lCovU2 = evaluateCovU2();
967    
968     double lCos2 = cos(fUPhiMVA)*cos(fUPhiMVA);
969     double lSin2 = sin(fUPhiMVA)*sin(fUPhiMVA);
970     //Now Compute teh covariance matrix in X and Y
971 pharris 1.5 (*fCov)(0,0) = lCovU1*lCos2+lCovU2*lSin2;
972     (*fCov)(1,0) = -lCovU1*sin(fUPhiMVA)*cos(fUPhiMVA)+lCovU2*sin(fUPhiMVA)*cos(fUPhiMVA);
973     (*fCov)(0,1) = (*fCov)(1,0);
974     (*fCov)(1,1) = lCovU1*lSin2+lCovU2*lCos2;
975 pharris 1.10 TMatrixD lInv(2,2);
976     lInv(0,0) = (*fCov)(0,0); lInv(1,1) = (*fCov)(1,1); lInv(1,0) = (*fCov)(1,0); lInv(0,1) = (*fCov)(0,1);
977     if(lInv.Determinant() != 0) lInv.Invert();
978     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));
979     fUncertainty = sqrt(lCovU1+lCovU2);
980 pharris 1.7
981 pharris 1.1 if (printDebug == kTRUE) {
982 ceballos 1.11 std::cout << "Debug Met MVA: "
983 pharris 1.1 << fU << " : "
984     << fUPhi << " : "
985     << fTKSumEt << " : "
986     << fTKU << " : "
987     << fTKUPhi << " : "
988     << fNPSumEt << " : "
989     << fNPU << " : "
990     << fNPUPhi << " : "
991     << fPUSumEt << " : "
992     << fPUMet << " : "
993     << fPUMetPhi << " : "
994     << fPCSumEt << " : "
995     << fPCU << " : "
996     << fPCUPhi << " : "
997     << fJSPt1 << " : "
998     << fJSEta1 << " : "
999     << fJSPhi1 << " : "
1000     << fJSPt2 << " : "
1001     << fJSEta2 << " : "
1002     << fJSPhi2 << " : "
1003     << fNJet << " : "
1004     << fNAllJet << " : "
1005     << fNPV << " : "
1006     << " === : === "
1007 ceballos 1.12 << lMet.Pt() << " : "
1008     << lMet.Phi() << " : "
1009 pharris 1.1 << std::endl;
1010     }
1011    
1012     return lMet;
1013     }
1014 pharris 1.26 //--------------------------------------------------------------------------------------------------
1015     Met MVAMet::GetMet(const MuonCol *iMuons,const ElectronCol *iElectrons,const PFTauCol *iTaus,
1016     const PFCandidateCol *iCands,const PFJetCol *iJets,const Vertex *iPV,const VertexCol *iVertices,const PFMetCol *iPFMet,
1017     FactorizedJetCorrector *iJetCorrector,const PileupEnergyDensityCol* iPUEnergyDensity) {
1018     const Vertex *lPV = iPV; if(iPV == 0 && iVertices->GetEntries() > 0) lPV = iVertices->At(0);
1019     FourVectorM lTotVec(0,0,0,0); double lTotSumEt = 0;
1020     FourVectorM lVisVec(0,0,0,0); double lVisSumEt = 0;
1021     std::vector<std::pair<FourVectorM,FourVectorM> > lDecay; std::vector<int> lId;
1022     fNMuons = 0;
1023     for(UInt_t i0 = 0; i0 < iMuons->GetEntries(); i0++) {
1024     const Muon *pMu = iMuons->At(i0);
1025     if(!MetLeptonTools::looseMuId(pMu,iCands,lPV,iVertices)) continue;
1026     std::pair<FourVectorM,FourVectorM> pVec(pMu->Mom(),pMu->Mom());
1027     lDecay .push_back(pVec);
1028     lId .push_back(0);
1029     fNMuons++;
1030     }
1031     fNElectrons=0;
1032     for(UInt_t i0 = 0; i0 < iElectrons->GetEntries(); i0++) {
1033     const Electron *pElectron = iElectrons->At(i0);
1034     if(!MetLeptonTools::looseEleId(pElectron,iPUEnergyDensity,iCands,lPV,iVertices)) continue;
1035     std::pair<FourVectorM,FourVectorM> pVec(pElectron->Mom(),pElectron->Mom());
1036     lDecay .push_back(pVec);
1037     lId .push_back(1);
1038     fNElectrons++;
1039     }
1040     fNTaus = 0;
1041     for(UInt_t i0 = 0; i0 < iTaus->GetEntries(); i0++) {
1042     const PFTau *pTau = iTaus->At(i0);
1043 pharris 1.28 if(!fMetLeptonTools->looseTauId(pTau,iPUEnergyDensity)) continue;
1044 pharris 1.26 FourVectorM pVis(0,0,0,0); pVis.SetCoordinates(pTau->Pt()*MetLeptonTools::vis(pTau),pTau->Eta(),pTau->Phi(),pTau->Mass()*MetLeptonTools::vis(pTau));
1045     std::pair<FourVectorM,FourVectorM> pVec(pTau->Mom(),pVis);
1046     lDecay .push_back(pVec);
1047     lId .push_back(2);
1048     fNTaus++;
1049     }
1050     std::vector<std::pair<FourVectorM,FourVectorM> > lFinalDecay;
1051     std::vector<int> lFinalId;
1052     for(unsigned int i0 = 0; i0 < lDecay.size(); i0++) {
1053     bool pAdd = true;
1054     for(unsigned int i1 = 0; i1 < lDecay.size(); i1++) {
1055     if(i0 == i1) continue;
1056     if(MathUtils::DeltaR(lDecay[i0].first,lDecay[i1].first) < 0.5) pAdd = false;
1057     if(!pAdd && lId[i0] != 2 && lId[i1] == 2) pAdd = true;
1058     if(!pAdd && ((lId[i0] != 2 && lId[i1] != 2) || (lId[i0] == 2 && lId[i1] == 2))
1059     && lDecay[i0].first.pt() > lDecay[i1].first.pt()) pAdd = true;
1060     if(MathUtils::DeltaR(lDecay[i0].first,lDecay[i1].first) < 0.5 && lDecay[i0].first.pt() == lDecay[i1].first.pt()) { pAdd = true;
1061     for(unsigned int i2 = 0; i2 < lFinalDecay.size(); i2++) if(fabs(lFinalDecay[i2].first.pt() - lDecay[i0].first.pt()) < 0.1) pAdd = false;
1062     }
1063     if(!pAdd) break;
1064     //if(!pAdd && lId[i2] == 2 && lId[i1] != 2) pAdd = true;
1065     }
1066     if(pAdd) lFinalDecay.push_back(lDecay[i0]);
1067     if(pAdd) lFinalId .push_back(lId [i0]);
1068     }
1069     for(unsigned int i0 = 0; i0 < lFinalDecay.size(); i0++) {
1070     lTotVec += lFinalDecay[i0].first;
1071     lVisVec += lFinalDecay[i0].second;
1072     lTotSumEt += lFinalDecay[i0].first.pt();
1073     lVisSumEt += lFinalDecay[i0].second.pt();
1074     }
1075     PFJetOArr *lCleanJets = new PFJetOArr();
1076     for(UInt_t i0 = 0; i0 < iJets->GetEntries(); i0++) {
1077     const PFJet *pJet = iJets->At(i0);
1078     bool pClean = false;
1079     for(unsigned int i1 = 0; i1 < lFinalDecay.size(); i1++) {
1080     if(MathUtils::DeltaR(pJet->Mom(),lFinalDecay[i1].first) < 0.5) pClean = true;
1081     }
1082     if(!pClean) lCleanJets->Add(pJet);
1083     }
1084 pharris 1.27 //for(unsigned int i0 = 0; i0 < lFinalDecay.size(); i0++) std::cout << "----> " << lFinalDecay[i0].first.pt() << " -- " << lFinalDecay[i0].first.phi() << " -- " << lFinalDecay[i0].first.eta() << " -- " << lFinalId[i0] << " -- " << MathUtils::DeltaR(lFinalDecay[i0].first,lFinalDecay[0].first)<< std::endl;
1085 pharris 1.26 Met lMVAMet = GetMet( false,
1086     lTotVec.Pt(),lTotVec.Phi(),lTotSumEt,
1087     lVisVec.Pt(),lVisVec.Phi(),lVisSumEt,
1088     iPFMet->At(0),
1089     iCands,lPV,iVertices,
1090     lCleanJets,
1091     iJetCorrector,
1092     iPUEnergyDensity,
1093 pharris 1.27 int(iVertices->GetEntries()));//,true);
1094 pharris 1.26 return lMVAMet;
1095     }
1096 pharris 1.30 Met MVAMet::GetMet(const PhotonCol *iPhotons,
1097     const PFCandidateCol *iCands,const PFJetCol *iJets,const Vertex *iPV,const VertexCol *iVertices,const PFMetCol *iPFMet,
1098     FactorizedJetCorrector *iJetCorrector,const PileupEnergyDensityCol* iPUEnergyDensity) {
1099     const Vertex *lPV = iPV; if(iPV == 0 && iVertices->GetEntries() > 0) lPV = iVertices->At(0);
1100     FourVectorM lTotVec(0,0,0,0); double lTotSumEt = 0;
1101     FourVectorM lVisVec(0,0,0,0); double lVisSumEt = 0;
1102     std::vector<std::pair<FourVectorM,FourVectorM> > lDecay; std::vector<int> lId;
1103     fNPhotons = 0;
1104     for(UInt_t i0 = 0; i0 < iPhotons->GetEntries(); i0++) {
1105     const Photon *pPhoton = iPhotons->At(i0);
1106     if(!MetLeptonTools::loosePhotonId(pPhoton)) continue;
1107     FourVectorM pVis(0,0,0,0); pVis.SetCoordinates(0.,pPhoton->Eta(),pPhoton->Phi(),0);
1108     std::pair<FourVectorM,FourVectorM> pVec(pPhoton->Mom(),pVis);
1109     lDecay .push_back(pVec);
1110     lId .push_back(0);
1111     fNPhotons++;
1112     }
1113     std::vector<std::pair<FourVectorM,FourVectorM> > lFinalDecay;
1114     std::vector<int> lFinalId;
1115     for(unsigned int i0 = 0; i0 < lDecay.size(); i0++) {
1116     bool pAdd = true;
1117     for(unsigned int i1 = 0; i1 < lDecay.size(); i1++) {
1118     if(i0 == i1) continue;
1119     if(MathUtils::DeltaR(lDecay[i0].first,lDecay[i1].first) < 0.5) pAdd = false;
1120     if(!pAdd && lId[i0] != 2 && lId[i1] == 2) pAdd = true;
1121     if(!pAdd && ((lId[i0] != 2 && lId[i1] != 2) || (lId[i0] == 2 && lId[i1] == 2))
1122     && lDecay[i0].first.pt() > lDecay[i1].first.pt()) pAdd = true;
1123     if(MathUtils::DeltaR(lDecay[i0].first,lDecay[i1].first) < 0.5 && lDecay[i0].first.pt() == lDecay[i1].first.pt()) { pAdd = true;
1124     for(unsigned int i2 = 0; i2 < lFinalDecay.size(); i2++) if(fabs(lFinalDecay[i2].first.pt() - lDecay[i0].first.pt()) < 0.1) pAdd = false;
1125     }
1126     if(!pAdd) break;
1127     //if(!pAdd && lId[i2] == 2 && lId[i1] != 2) pAdd = true;
1128     }
1129     if(pAdd) lFinalDecay.push_back(lDecay[i0]);
1130     if(pAdd) lFinalId .push_back(lId [i0]);
1131     }
1132     for(unsigned int i0 = 0; i0 < lFinalDecay.size(); i0++) {
1133     lTotVec += lFinalDecay[i0].first;
1134     lVisVec += lFinalDecay[i0].second;
1135     lTotSumEt += lFinalDecay[i0].first.pt();
1136     lVisSumEt += lFinalDecay[i0].second.pt();
1137     }
1138     PFJetOArr *lCleanJets = new PFJetOArr();
1139     for(UInt_t i0 = 0; i0 < iJets->GetEntries(); i0++) {
1140     const PFJet *pJet = iJets->At(i0);
1141     bool pClean = false;
1142     for(unsigned int i1 = 0; i1 < lFinalDecay.size(); i1++) {
1143     if(MathUtils::DeltaR(pJet->Mom(),lFinalDecay[i1].first) < 0.5) pClean = true;
1144     }
1145     if(!pClean) lCleanJets->Add(pJet);
1146     }
1147     //for(unsigned int i0 = 0; i0 < lFinalDecay.size(); i0++) std::cout << "----> " << lFinalDecay[i0].first.pt() << " -- " << lFinalDecay[i0].first.phi() << " -- " << lFinalDecay[i0].first.eta() << " -- " << lFinalId[i0] << " -- " << MathUtils::DeltaR(lFinalDecay[i0].first,lFinalDecay[0].first)<< std::endl;
1148     Met lMVAMet = GetMet( false,
1149     lTotVec.Pt(),lTotVec.Phi(),lTotSumEt,
1150     lVisVec.Pt(),lVisVec.Phi(),lVisSumEt,
1151     iPFMet->At(0),
1152     iCands,lPV,iVertices,
1153     lCleanJets,
1154     iJetCorrector,
1155     iPUEnergyDensity,
1156     int(iVertices->GetEntries()));//,true);
1157     return lMVAMet;
1158     }
1159    
1160 pharris 1.1