ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MVAMet.cc
Revision: 1.29
Committed: Mon Oct 1 17:37:48 2012 UTC (12 years, 7 months ago) by pharris
Content type: text/plain
Branch: MAIN
Changes since 1.28: +1 -1 lines
Log Message:
fixed 42X bug

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