ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MVAMet.cc
Revision: 1.7
Committed: Tue Apr 24 21:28:01 2012 UTC (13 years ago) by pharris
Content type: text/plain
Branch: MAIN
Changes since 1.6: +40 -33 lines
Log Message:
Final debug of MVA Met

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     //--------------------------------------------------------------------------------------------------
242     Double_t MVAMet::evaluateCovU2() {
243 pharris 1.5 fCovVals[0] = fSumEt ;
244     fCovVals[1] = fNPV ;
245     fCovVals[2] = fU ;
246     fCovVals[3] = fUPhi ;
247     fCovVals[4] = fTKSumEt ;
248     fCovVals[5] = fTKU ;
249     fCovVals[6] = fTKUPhi ;
250     fCovVals[7] = fNPSumEt ;
251     fCovVals[8] = fNPU ;
252     fCovVals[9] = fNPUPhi ;
253     fCovVals[10] = fPUSumEt ;
254     fCovVals[11] = fPUMet ;
255     fCovVals[12] = fPUMetPhi;
256     fCovVals[13] = fPCSumEt ;
257     fCovVals[14] = fPCU ;
258     fCovVals[15] = fPCUPhi ;
259     fCovVals[16] = fJSPt1 ;
260     fCovVals[17] = fJSEta1 ;
261     fCovVals[18] = fJSPhi1 ;
262     fCovVals[19] = fJSPt2 ;
263     fCovVals[20] = fJSEta2 ;
264     fCovVals[21] = fJSPhi2 ;
265     fCovVals[22] = fNAllJet ;
266     fCovVals[23] = fNJet ;
267     fCovVals[24] = fUPhiMVA ;
268     fCovVals[25] = fUMVA ;
269 pharris 1.4 return fCovU2Reader->GetResponse(fCovVals);
270     }
271     //--------------------------------------------------------------------------------------------------
272 pharris 1.1 Double_t MVAMet::MVAValue( Bool_t iPhi,
273     Float_t iPFSumEt,
274     Float_t iU ,
275     Float_t iUPhi ,
276     Float_t iTKSumEt,
277     Float_t iTKU ,
278     Float_t iTKUPhi ,
279     Float_t iNPSumEt,
280     Float_t iNPU ,
281     Float_t iNPUPhi ,
282     Float_t iPUSumEt,
283     Float_t iPUMet ,
284     Float_t iPUMetPhi,
285     Float_t iPCSumEt,
286     Float_t iPCU ,
287     Float_t iPCUPhi ,
288     Float_t iJSPt1 ,
289     Float_t iJSEta1 ,
290     Float_t iJSPhi1 ,
291     Float_t iJSPt2 ,
292     Float_t iJSEta2 ,
293     Float_t iJSPhi2 ,
294     Float_t iNAllJet,
295     Float_t iNJet ,
296     Float_t iNPV
297     ){
298    
299     if (!fIsInitialized) {
300     std::cout << "Error: MVA Met not properly initialized.\n";
301     return -9999;
302     }
303 pharris 1.5
304     fSumEt = iPFSumEt;
305 pharris 1.1 fU = iU ;
306     fUPhi = iUPhi ;
307     fTKSumEt = iTKSumEt/iPFSumEt;
308     fTKU = iTKU ;
309     fTKUPhi = iTKUPhi ;
310     fNPSumEt = iNPSumEt/iPFSumEt;
311     fNPU = iNPU ;
312     fNPUPhi = iNPUPhi ;
313     fPUSumEt = iPUSumEt/iPFSumEt;
314     fPUMet = iPUMet ;
315     fPUMetPhi = iPUMetPhi;
316     fPCSumEt = iPCSumEt/iPFSumEt;
317     fPCU = iPCU ;
318     fPCUPhi = iPCUPhi ;
319     fJSPt1 = iJSPt1 ;
320     fJSEta1 = iJSEta1 ;
321     fJSPhi1 = iJSPhi1 ;
322     fJSPt2 = iJSPt2 ;
323     fJSEta2 = iJSEta2 ;
324     fJSPhi2 = iJSPhi2 ;
325     fNAllJet = iNAllJet;
326     fNJet = iNJet ;
327     fNPV = iNPV ;
328 pharris 1.5
329 pharris 1.1 Double_t lMVA = -9999;
330 pharris 1.2 lMVA = evaluatePhi();
331 pharris 1.1 if(!iPhi) fUPhiMVA = iUPhi+lMVA;
332 pharris 1.2 //Not no nice feature of the training
333 pharris 1.4 //fTKSumEt /= iPFSumEt;
334     //fNPSumEt /= iPFSumEt;
335     //fPUSumEt /= iPFSumEt;
336     //fPCSumEt /= iPFSumEt;
337 pharris 1.2 if(!iPhi) lMVA = evaluateU1();
338 pharris 1.1 return lMVA;
339     }
340     //--------------------------------------------------------------------------------------------------
341     //====> Please not that the jet collection must be cleaned => all jets near leptons must be removed
342     Met MVAMet::GetMet( Bool_t iPhi,Float_t iPtVis,Float_t iPhiVis,Float_t iSumEtVis,
343     const PFMet *iMet ,
344 pharris 1.3 const PFCandidateCol *iCands,
345     const Vertex *iVertex,const VertexCol *iVertices,
346 pharris 1.1 const PFJetCol *iJets ,
347     FactorizedJetCorrector *iJetCorrector,
348     const PileupEnergyDensityCol *iPUEnergyDensity,
349     int iNPV,
350     Bool_t printDebug) {
351    
352 pharris 1.7 Met lPFRec = fRecoilTools->pfRecoil (iPtVis,iPhiVis,iSumEtVis,iCands);
353 pharris 1.1 Met lTKRec = fRecoilTools->trackRecoil(iPtVis,iPhiVis,iSumEtVis,iCands,iVertex);
354 pharris 1.3 Met lNPRec = fRecoilTools->NoPURecoil (iPtVis,iPhiVis,iSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices);
355     Met lPCRec = fRecoilTools->PUCRecoil (iPtVis,iPhiVis,iSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices);
356     Met lPUMet = fRecoilTools->PUMet (iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices);
357 pharris 1.1
358     Double_t lPt0 = 0; const PFJet *lLead = 0;
359     Double_t lPt1 = 0; const PFJet *l2nd = 0;
360     int lNAllJet = 0;
361     int lNJet = 0;
362     for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
363     const PFJet *pJet = iJets->At(i0);
364     Double_t pPt = fRecoilTools->fJetIDMVA->correctedPt(pJet,iJetCorrector,iPUEnergyDensity);
365     if(!JetTools::passPFLooseId(pJet)) continue;
366     lNAllJet++;
367     if(pPt > 30) lNJet++;
368     if(lPt0 < pPt) {lPt0 = pPt; lLead = pJet; continue;}
369     if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
370     }
371 pharris 1.5
372     fSumEt = lPFRec.SumEt();
373 pharris 1.1 fU = lPFRec.Pt();
374     fUPhi = lPFRec.Phi();
375     fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
376     fTKU = lTKRec.Pt();
377     fTKUPhi = lTKRec.Phi();
378     fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
379     fNPU = lNPRec.Pt();
380     fNPUPhi = lNPRec.Phi();
381     fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
382     fPUMet = lPUMet.Pt();
383     fPUMetPhi = lPUMet.Phi();
384     fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
385     fPCU = lPCRec.Pt() ;
386     fPCUPhi = lPCRec.Phi() ;
387     fJSPt1 = lPt0;
388     fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
389     fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
390     fJSPt2 = lPt1;
391     fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
392     fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
393     fNJet = lNJet ;
394     fNAllJet = lNAllJet;
395     fNPV = iNPV ;
396 pharris 1.5 Float_t lMVA = evaluatePhi();
397     if(!iPhi) fUPhiMVA = fUPhi + lMVA;
398 pharris 1.1 //Not no nice feature of teh training
399 pharris 1.4 //fTKSumEt /= lPFRec.SumEt();
400     //fNPSumEt /= lPFRec.SumEt();
401     //fPUSumEt /= lPFRec.SumEt();
402     //fPCSumEt /= lPFRec.SumEt();
403 pharris 1.7 if(!iPhi) lMVA = evaluateU1();//fU1Reader ->EvaluateMVA( fU1MethodName );
404     fUMVA = fU*lMVA;
405    
406     TLorentzVector lUVec(0,0,0,0); lUVec .SetPtEtaPhiM(fUMVA,0,fUPhiMVA,0);
407 pharris 1.1 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 pharris 1.7 fUncertainty = lUVec.Px()*lUVec.Px()*(*fCov)(0,0) + 2.*lUVec.Px()*lUVec.Py()*(*fCov)(1,0) + lUVec.Py()*lUVec.Py()*(*fCov)(1,1);
427     fSignificance = lMet.Pt()/fUncertainty;
428    
429 pharris 1.1 if (printDebug == kTRUE) {
430     std::cout << "Debug Jet MVA: "
431     << fU << " : "
432     << fUPhi << " : "
433     << fTKSumEt << " : "
434     << fTKU << " : "
435     << fTKUPhi << " : "
436     << fNPSumEt << " : "
437     << fNPU << " : "
438     << fNPUPhi << " : "
439     << fPUSumEt << " : "
440     << fPUMet << " : "
441     << fPUMetPhi << " : "
442     << fPCUPhi << " : "
443     << fPCSumEt << " : "
444     << fPCU << " : "
445     << fPCUPhi << " : "
446     << fJSPt1 << " : "
447     << fJSEta1 << " : "
448     << fJSPhi1 << " : "
449     << fJSPt2 << " : "
450     << fJSEta2 << " : "
451     << fJSPhi2 << " : "
452     << fNJet << " : "
453     << fNAllJet << " : "
454     << fNPV << " : "
455     << " === : === "
456     << std::endl;
457     }
458    
459     return lMet;
460     }
461 pharris 1.2 //--------------------------------------------------------------------------------------------------
462     //====> Please not that the jet collection must be cleaned => all jets near leptons must be removed
463     //====> Corrected Jets
464     Met MVAMet::GetMet( Bool_t iPhi,Float_t iPtVis,Float_t iPhiVis,Float_t iSumEtVis,
465     const PFMet *iMet ,
466 pharris 1.3 const PFCandidateCol *iCands,
467 pharris 1.7 const Vertex *iVertex ,const VertexCol *iVertices,Double_t iRho,
468 pharris 1.2 const PFJetCol *iJets ,
469     int iNPV,
470     Bool_t printDebug) {
471    
472 pharris 1.7 Met lPFRec = fRecoilTools->pfRecoil (iPtVis,iPhiVis,iSumEtVis,iCands);
473 pharris 1.3 Met lTKRec = fRecoilTools->trackRecoil(iPtVis,iPhiVis,iSumEtVis, iCands,iVertex);
474 pharris 1.7 Met lNPRec = fRecoilTools->NoPURecoil (iPtVis,iPhiVis,iSumEtVis,iJets,iCands,iVertex,iVertices,iRho);
475     Met lPCRec = fRecoilTools->PUCRecoil (iPtVis,iPhiVis,iSumEtVis,iJets,iCands,iVertex,iVertices,iRho);
476     Met lPUMet = fRecoilTools->PUMet ( iJets,iCands,iVertex,iVertices,iRho);
477 pharris 1.2
478     Double_t lPt0 = 0; const PFJet *lLead = 0;
479     Double_t lPt1 = 0; const PFJet *l2nd = 0;
480     int lNAllJet = 0;
481     int lNJet = 0;
482     for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
483     const PFJet *pJet = iJets->At(i0);
484     Double_t pPt = pJet->Pt();
485     if(!JetTools::passPFLooseId(pJet)) continue;
486     lNAllJet++;
487     if(pPt > 30) lNJet++;
488     if(lPt0 < pPt) {lPt0 = pPt; lLead = pJet; continue;}
489     if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
490     }
491    
492 pharris 1.5 fSumEt = lPFRec.SumEt();
493 pharris 1.2 fU = lPFRec.Pt();
494     fUPhi = lPFRec.Phi();
495     fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
496     fTKU = lTKRec.Pt();
497     fTKUPhi = lTKRec.Phi();
498     fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
499     fNPU = lNPRec.Pt();
500     fNPUPhi = lNPRec.Phi();
501     fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
502     fPUMet = lPUMet.Pt();
503     fPUMetPhi = lPUMet.Phi();
504     fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
505     fPCU = lPCRec.Pt() ;
506     fPCUPhi = lPCRec.Phi() ;
507     fJSPt1 = lPt0;
508     fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
509     fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
510     fJSPt2 = lPt1;
511     fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
512     fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
513     fNJet = lNJet ;
514     fNAllJet = lNAllJet;
515     fNPV = iNPV ;
516 pharris 1.7
517 pharris 1.2 Float_t lMVA = evaluatePhi();
518    
519     if(!iPhi) fUPhiMVA = fUPhi + lMVA;
520     //Not no nice feature of teh training
521 pharris 1.4 //fTKSumEt /= lPFRec.SumEt();
522     //fNPSumEt /= lPFRec.SumEt();
523     //fPUSumEt /= lPFRec.SumEt();
524     //fPCSumEt /= lPFRec.SumEt();
525 pharris 1.7 if(!iPhi) lMVA = evaluateU1();//fU1Reader ->EvaluateMVA( fU1MethodName );
526     fUMVA = lMVA*fU;
527    
528     TLorentzVector lUVec(0,0,0,0); lUVec .SetPtEtaPhiM(fUMVA ,0,fUPhiMVA,0);
529 pharris 1.2 TLorentzVector lVVec(0,0,0,0); lVVec .SetPtEtaPhiM(iPtVis ,0,iPhiVis ,0);
530     if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
531     lUVec -= lVVec;
532     Met lMet(lUVec.Px(),lUVec.Py());
533 pharris 1.4 //Cov matrix
534     double lCovU1 = evaluateCovU1();
535     double lCovU2 = evaluateCovU2();
536    
537     double lCos2 = cos(fUPhiMVA)*cos(fUPhiMVA);
538     double lSin2 = sin(fUPhiMVA)*sin(fUPhiMVA);
539    
540     //Now Compute teh covariance matrix in X and Y
541 pharris 1.5 (*fCov)(0,0) = lCovU1*lCos2+lCovU2*lSin2;
542     (*fCov)(1,0) = -lCovU1*sin(fUPhiMVA)*cos(fUPhiMVA)+lCovU2*sin(fUPhiMVA)*cos(fUPhiMVA);
543     (*fCov)(0,1) = (*fCov)(1,0);
544     (*fCov)(1,1) = lCovU1*lSin2+lCovU2*lCos2;
545 pharris 1.7 fUncertainty = lUVec.Px()*lUVec.Px()*(*fCov)(0,0) + 2.*lUVec.Px()*lUVec.Py()*(*fCov)(1,0) + lUVec.Py()*lUVec.Py()*(*fCov)(1,1);
546     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.7 fUncertainty = lUVec.Px()*lUVec.Px()*(*fCov)(0,0) + 2.*lUVec.Px()*lUVec.Py()*(*fCov)(1,0) + lUVec.Py()*lUVec.Py()*(*fCov)(1,1);
683     fSignificance = lMet.Pt()/fUncertainty;
684 pharris 1.2
685     if (printDebug == kTRUE) {
686     std::cout << "Debug Jet MVA: "
687     << fU << " : "
688     << fUPhi << " : "
689     << fTKSumEt << " : "
690     << fTKU << " : "
691     << fTKUPhi << " : "
692     << fNPSumEt << " : "
693     << fNPU << " : "
694     << fNPUPhi << " : "
695     << fPUSumEt << " : "
696     << fPUMet << " : "
697     << fPUMetPhi << " : "
698     << fPCSumEt << " : "
699     << fPCU << " : "
700     << fPCUPhi << " : "
701     << fJSPt1 << " : "
702     << fJSEta1 << " : "
703     << fJSPhi1 << " : "
704     << fJSPt2 << " : "
705     << fJSEta2 << " : "
706     << fJSPhi2 << " : "
707     << fNJet << " : "
708     << fNAllJet << " : "
709     << fNPV << " : "
710     << " === : === "
711     << std::endl;
712     }
713    
714     return lMet;
715     }
716     //--------------------------------------------------------------------------------------------------
717     //Interms of the two candidates => corrected Jets
718     Met MVAMet::GetMet( Bool_t iPhi,
719     Float_t iPt1,Float_t iPhi1,Float_t iEta1,
720     Float_t iPt2,Float_t iPhi2,Float_t iEta2,
721     const PFMet *iMet ,
722 pharris 1.3 const PFCandidateCol *iCands,
723 pharris 1.7 const Vertex *iVertex,const VertexCol *iVertices,Double_t iRho,
724 pharris 1.2 const PFJetCol *iJets ,
725     int iNPV,
726     Bool_t printDebug) {
727    
728     TLorentzVector lVVec1(0,0,0,0); lVVec1.SetPtEtaPhiM(iPt1,0,iPhi1 ,0);
729     TLorentzVector lVVec2(0,0,0,0); lVVec2.SetPtEtaPhiM(iPt2,0,iPhi2 ,0);
730     lVVec1+=lVVec2;
731     Float_t lPtVis = lVVec1.Pt();
732     Float_t lPhiVis = lVVec1.Phi();
733     Float_t lSumEtVis = iPt1 + iPt2;
734 pharris 1.7 Met lPFRec = fRecoilTools->pfRecoil (lPtVis,lPhiVis,lSumEtVis,iCands);
735 pharris 1.2 Met lTKRec = fRecoilTools->trackRecoil(lPtVis,lPhiVis,lSumEtVis,iCands,iVertex);
736 pharris 1.3 Met lNPRec = fRecoilTools->NoPURecoil (lPtVis,lPhiVis,lSumEtVis,iJets,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2);
737     Met lPCRec = fRecoilTools->PUCRecoil (lPtVis,lPhiVis,lSumEtVis,iJets,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2);
738     Met lPUMet = fRecoilTools->PUMet ( iJets,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2);
739 pharris 1.2
740     Double_t lPt0 = 0; const PFJet *lLead = 0;
741     Double_t lPt1 = 0; const PFJet *l2nd = 0;
742     int lNAllJet = 0;
743     int lNJet = 0;
744     for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
745     const PFJet *pJet = iJets->At(i0);
746     Double_t pPt = pJet->Pt();
747     if( pJet->TrackCountingHighEffBJetTagsDisc() == -100 && pPt < 10.) continue;
748     double pDEta1 = pJet->Eta() - iEta1;
749     double pDPhi1 = fabs(pJet->Phi() - iPhi1); if(pDPhi1 > 2.*TMath::Pi()-pDPhi1) pDPhi1 = 2.*TMath::Pi()-pDPhi1;
750     double pDR1 = sqrt(pDEta1*pDEta1 + pDPhi1*pDPhi1);
751     if(pDR1 < 0.5) continue;
752     double pDEta2 = pJet->Eta() - iEta2;
753     double pDPhi2 = fabs(pJet->Phi() - iPhi2); if(pDPhi2 > 2.*TMath::Pi()-pDPhi2) pDPhi2 = 2.*TMath::Pi()-pDPhi2;
754     double pDR2 = sqrt(pDEta2*pDEta2 + pDPhi2*pDPhi2);
755     if(pDR2 < 0.5) continue;
756     if(!JetTools::passPFLooseId(pJet)) continue;
757     lNAllJet++;
758     if(pPt > 30.) lNJet++;
759     if(lPt0 < pPt) {lPt0 = pPt; lLead = pJet; continue;}
760     if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
761     }
762 pharris 1.5 fSumEt = lPFRec.SumEt();
763 pharris 1.2 fU = lPFRec.Pt();
764     fUPhi = lPFRec.Phi();
765     fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
766     fTKU = lTKRec.Pt();
767     fTKUPhi = lTKRec.Phi();
768     fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
769     fNPU = lNPRec.Pt();
770     fNPUPhi = lNPRec.Phi();
771     fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
772     fPUMet = lPUMet.Pt();
773     fPUMetPhi= lPUMet.Phi();
774     fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
775     fPCU = lPCRec.Pt() ;
776     fPCUPhi = lPCRec.Phi() ;
777     fJSPt1 = lPt0;
778     fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
779     fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
780     fJSPt2 = lPt1;
781     fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
782     fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
783     fNJet = lNJet ;
784     fNAllJet = lNAllJet;
785     fNPV = iNPV ;
786    
787     Float_t lMVA = evaluatePhi();//fPhiReader ->EvaluateMVA( fPhiMethodName );
788 pharris 1.1
789     if(!iPhi) fUPhiMVA = fUPhi + lMVA;
790 pharris 1.4 //fTKSumEt /= lPFRec.SumEt();
791     //fNPSumEt /= lPFRec.SumEt();
792     //fPUSumEt /= lPFRec.SumEt();
793     //fPCSumEt /= lPFRec.SumEt();
794 pharris 1.7 if(!iPhi) lMVA = evaluateU1();//fU1Reader ->EvaluateMVA( fU1MethodName );
795     fUMVA = lMVA*fU;
796    
797     TLorentzVector lUVec (0,0,0,0); lUVec .SetPtEtaPhiM(fUMVA ,0,fUPhiMVA,0);
798 pharris 1.1 TLorentzVector lVVec (0,0,0,0); lVVec .SetPtEtaPhiM(lPtVis ,0,lPhiVis ,0);
799 pharris 1.7 if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
800 pharris 1.1 lUVec -= lVVec;
801     Met lMet(lUVec.Px(),lUVec.Py());
802 pharris 1.4 //Cov matrix
803     double lCovU1 = evaluateCovU1();
804     double lCovU2 = evaluateCovU2();
805    
806     double lCos2 = cos(fUPhiMVA)*cos(fUPhiMVA);
807     double lSin2 = sin(fUPhiMVA)*sin(fUPhiMVA);
808     //Now Compute teh covariance matrix in X and Y
809 pharris 1.5 (*fCov)(0,0) = lCovU1*lCos2+lCovU2*lSin2;
810     (*fCov)(1,0) = -lCovU1*sin(fUPhiMVA)*cos(fUPhiMVA)+lCovU2*sin(fUPhiMVA)*cos(fUPhiMVA);
811     (*fCov)(0,1) = (*fCov)(1,0);
812     (*fCov)(1,1) = lCovU1*lSin2+lCovU2*lCos2;
813 pharris 1.7 fUncertainty = lUVec.Px()*lUVec.Px()*(*fCov)(0,0) + 2.*lUVec.Px()*lUVec.Py()*(*fCov)(1,0) + lUVec.Py()*lUVec.Py()*(*fCov)(1,1);
814     fSignificance = lMet.Pt()/fUncertainty;
815    
816 pharris 1.1 if (printDebug == kTRUE) {
817     std::cout << "Debug Jet MVA: "
818     << fU << " : "
819     << fUPhi << " : "
820     << fTKSumEt << " : "
821     << fTKU << " : "
822     << fTKUPhi << " : "
823     << fNPSumEt << " : "
824     << fNPU << " : "
825     << fNPUPhi << " : "
826     << fPUSumEt << " : "
827     << fPUMet << " : "
828     << fPUMetPhi << " : "
829     << fPCSumEt << " : "
830     << fPCU << " : "
831     << fPCUPhi << " : "
832     << fJSPt1 << " : "
833     << fJSEta1 << " : "
834     << fJSPhi1 << " : "
835     << fJSPt2 << " : "
836     << fJSEta2 << " : "
837     << fJSPhi2 << " : "
838     << fNJet << " : "
839     << fNAllJet << " : "
840     << fNPV << " : "
841     << " === : === "
842     << std::endl;
843     }
844    
845     return lMet;
846     }
847