ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MVAMet.cc
Revision: 1.26
Committed: Thu Aug 9 13:15:43 2012 UTC (12 years, 9 months ago) by pharris
Content type: text/plain
Branch: MAIN
Changes since 1.25: +86 -1 lines
Log Message:
Added New MVA Met

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