ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MVAMet.cc
Revision: 1.4
Committed: Mon Apr 23 15:38:00 2012 UTC (13 years ago) by pharris
Content type: text/plain
Branch: MAIN
Changes since 1.3: +156 -26 lines
Log Message:
52 MET update w/Covariance matrix

File Contents

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