ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MVAMet.cc
Revision: 1.9
Committed: Tue Apr 24 23:56:34 2012 UTC (13 years ago) by pharris
Content type: text/plain
Branch: MAIN
Changes since 1.8: +7 -5 lines
Log Message:
*** empty log message ***

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