ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MVAMet.cc
Revision: 1.5
Committed: Mon Apr 23 19:45:16 2012 UTC (13 years ago) by pharris
Content type: text/plain
Branch: MAIN
Changes since 1.4: +119 -113 lines
Log Message:
Small fixes

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