ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MVAMet.cc
Revision: 1.11
Committed: Mon Apr 30 11:45:19 2012 UTC (13 years ago) by ceballos
Content type: text/plain
Branch: MAIN
Changes since 1.10: +4 -4 lines
Log Message:
small typo

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.10 TMatrixD lInv(2,2);
425     lInv(0,0) = (*fCov)(0,0); lInv(1,1) = (*fCov)(1,1); lInv(1,0) = (*fCov)(1,0); lInv(0,1) = (*fCov)(0,1);
426     if(lInv.Determinant() != 0) lInv.Invert();
427     fSignificance = TMath::Sqrt(lUVec.Px()*lUVec.Px()*(lInv)(0,0) + 2.*lUVec.Px()*lUVec.Py()*(lInv)(1,0) + lUVec.Py()*lUVec.Py()*(lInv)(1,1));
428     fUncertainty = sqrt(lCovU1+lCovU2);
429 pharris 1.7
430 pharris 1.1 if (printDebug == kTRUE) {
431 ceballos 1.11 std::cout << "Debug Met MVA: "
432 pharris 1.1 << fU << " : "
433     << fUPhi << " : "
434     << fTKSumEt << " : "
435     << fTKU << " : "
436     << fTKUPhi << " : "
437     << fNPSumEt << " : "
438     << fNPU << " : "
439     << fNPUPhi << " : "
440     << fPUSumEt << " : "
441     << fPUMet << " : "
442     << fPUMetPhi << " : "
443     << fPCUPhi << " : "
444     << fPCSumEt << " : "
445     << fPCU << " : "
446     << fPCUPhi << " : "
447     << fJSPt1 << " : "
448     << fJSEta1 << " : "
449     << fJSPhi1 << " : "
450     << fJSPt2 << " : "
451     << fJSEta2 << " : "
452     << fJSPhi2 << " : "
453     << fNJet << " : "
454     << fNAllJet << " : "
455     << fNPV << " : "
456     << " === : === "
457     << std::endl;
458     }
459    
460     return lMet;
461     }
462 pharris 1.2 //--------------------------------------------------------------------------------------------------
463     //====> Please not that the jet collection must be cleaned => all jets near leptons must be removed
464     //====> Corrected Jets
465     Met MVAMet::GetMet( Bool_t iPhi,Float_t iPtVis,Float_t iPhiVis,Float_t iSumEtVis,
466     const PFMet *iMet ,
467 pharris 1.3 const PFCandidateCol *iCands,
468 pharris 1.7 const Vertex *iVertex ,const VertexCol *iVertices,Double_t iRho,
469 pharris 1.2 const PFJetCol *iJets ,
470     int iNPV,
471     Bool_t printDebug) {
472    
473 pharris 1.7 Met lPFRec = fRecoilTools->pfRecoil (iPtVis,iPhiVis,iSumEtVis,iCands);
474 pharris 1.3 Met lTKRec = fRecoilTools->trackRecoil(iPtVis,iPhiVis,iSumEtVis, iCands,iVertex);
475 pharris 1.7 Met lNPRec = fRecoilTools->NoPURecoil (iPtVis,iPhiVis,iSumEtVis,iJets,iCands,iVertex,iVertices,iRho);
476     Met lPCRec = fRecoilTools->PUCRecoil (iPtVis,iPhiVis,iSumEtVis,iJets,iCands,iVertex,iVertices,iRho);
477     Met lPUMet = fRecoilTools->PUMet ( iJets,iCands,iVertex,iVertices,iRho);
478 pharris 1.2
479     Double_t lPt0 = 0; const PFJet *lLead = 0;
480     Double_t lPt1 = 0; const PFJet *l2nd = 0;
481     int lNAllJet = 0;
482     int lNJet = 0;
483     for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
484     const PFJet *pJet = iJets->At(i0);
485     Double_t pPt = pJet->Pt();
486     if(!JetTools::passPFLooseId(pJet)) continue;
487     lNAllJet++;
488     if(pPt > 30) lNJet++;
489     if(lPt0 < pPt) {lPt0 = pPt; lLead = pJet; continue;}
490     if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
491     }
492    
493 pharris 1.5 fSumEt = lPFRec.SumEt();
494 pharris 1.2 fU = lPFRec.Pt();
495     fUPhi = lPFRec.Phi();
496     fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
497     fTKU = lTKRec.Pt();
498     fTKUPhi = lTKRec.Phi();
499     fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
500     fNPU = lNPRec.Pt();
501     fNPUPhi = lNPRec.Phi();
502     fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
503     fPUMet = lPUMet.Pt();
504     fPUMetPhi = lPUMet.Phi();
505     fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
506     fPCU = lPCRec.Pt() ;
507     fPCUPhi = lPCRec.Phi() ;
508     fJSPt1 = lPt0;
509     fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
510     fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
511     fJSPt2 = lPt1;
512     fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
513     fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
514     fNJet = lNJet ;
515     fNAllJet = lNAllJet;
516     fNPV = iNPV ;
517 pharris 1.7
518 pharris 1.2 Float_t lMVA = evaluatePhi();
519    
520     if(!iPhi) fUPhiMVA = fUPhi + lMVA;
521     //Not no nice feature of teh training
522 pharris 1.4 //fTKSumEt /= lPFRec.SumEt();
523     //fNPSumEt /= lPFRec.SumEt();
524     //fPUSumEt /= lPFRec.SumEt();
525     //fPCSumEt /= lPFRec.SumEt();
526 pharris 1.7 if(!iPhi) lMVA = evaluateU1();//fU1Reader ->EvaluateMVA( fU1MethodName );
527     fUMVA = lMVA*fU;
528    
529     TLorentzVector lUVec(0,0,0,0); lUVec .SetPtEtaPhiM(fUMVA ,0,fUPhiMVA,0);
530 pharris 1.2 TLorentzVector lVVec(0,0,0,0); lVVec .SetPtEtaPhiM(iPtVis ,0,iPhiVis ,0);
531     if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
532     lUVec -= lVVec;
533     Met lMet(lUVec.Px(),lUVec.Py());
534 pharris 1.4 //Cov matrix
535     double lCovU1 = evaluateCovU1();
536     double lCovU2 = evaluateCovU2();
537    
538     double lCos2 = cos(fUPhiMVA)*cos(fUPhiMVA);
539     double lSin2 = sin(fUPhiMVA)*sin(fUPhiMVA);
540    
541     //Now Compute teh covariance matrix in X and Y
542 pharris 1.5 (*fCov)(0,0) = lCovU1*lCos2+lCovU2*lSin2;
543     (*fCov)(1,0) = -lCovU1*sin(fUPhiMVA)*cos(fUPhiMVA)+lCovU2*sin(fUPhiMVA)*cos(fUPhiMVA);
544     (*fCov)(0,1) = (*fCov)(1,0);
545     (*fCov)(1,1) = lCovU1*lSin2+lCovU2*lCos2;
546 pharris 1.10 TMatrixD lInv(2,2);
547     lInv(0,0) = (*fCov)(0,0); lInv(1,1) = (*fCov)(1,1); lInv(1,0) = (*fCov)(1,0); lInv(0,1) = (*fCov)(0,1);
548     if(lInv.Determinant() != 0) lInv.Invert();
549     fSignificance = TMath::Sqrt(lUVec.Px()*lUVec.Px()*(lInv)(0,0) + 2.*lUVec.Px()*lUVec.Py()*(lInv)(1,0) + lUVec.Py()*lUVec.Py()*(lInv)(1,1));
550     fUncertainty = sqrt(lCovU1+lCovU2);
551 pharris 1.2
552     if (printDebug == kTRUE) {
553 ceballos 1.11 std::cout << "Debug Met MVA: "
554 pharris 1.2 << fU << " : "
555     << fUPhi << " : "
556     << fTKSumEt << " : "
557     << fTKU << " : "
558     << fTKUPhi << " : "
559     << fNPSumEt << " : "
560     << fNPU << " : "
561     << fNPUPhi << " : "
562     << fPUSumEt << " : "
563     << fPUMet << " : "
564     << fPUMetPhi << " : "
565     << fPCUPhi << " : "
566     << fPCSumEt << " : "
567     << fPCU << " : "
568     << fPCUPhi << " : "
569     << fJSPt1 << " : "
570     << fJSEta1 << " : "
571     << fJSPhi1 << " : "
572     << fJSPt2 << " : "
573     << fJSEta2 << " : "
574     << fJSPhi2 << " : "
575     << fNJet << " : "
576     << fNAllJet << " : "
577     << fNPV << " : "
578     << " === : === "
579     << std::endl;
580     }
581    
582     return lMet;
583     }
584     //--------------------------------------------------------------------------------------------------
585 pharris 1.1 //Interms of the two candidates
586     Met MVAMet::GetMet( Bool_t iPhi,
587     Float_t iPt1,Float_t iPhi1,Float_t iEta1,
588     Float_t iPt2,Float_t iPhi2,Float_t iEta2,
589     const PFMet *iMet ,
590 pharris 1.3 const PFCandidateCol *iCands,
591     const Vertex *iVertex,const VertexCol *iVertices,
592 pharris 1.1 const PFJetCol *iJets ,
593     FactorizedJetCorrector *iJetCorrector,
594     const PileupEnergyDensityCol *iPUEnergyDensity,
595     int iNPV,
596     Bool_t printDebug) {
597    
598     TLorentzVector lVVec1(0,0,0,0); lVVec1.SetPtEtaPhiM(iPt1,0,iPhi1 ,0);
599     TLorentzVector lVVec2(0,0,0,0); lVVec2.SetPtEtaPhiM(iPt2,0,iPhi2 ,0);
600     lVVec1+=lVVec2;
601     Float_t lPtVis = lVVec1.Pt();
602     Float_t lPhiVis = lVVec1.Phi();
603     Float_t lSumEtVis = iPt1 + iPt2;
604 pharris 1.7 Met lPFRec = fRecoilTools->pfRecoil (lPtVis,lPhiVis,lSumEtVis,iCands);
605 pharris 1.1 Met lTKRec = fRecoilTools->trackRecoil(lPtVis,lPhiVis,lSumEtVis,iCands,iVertex);
606 pharris 1.3 Met lNPRec = fRecoilTools->NoPURecoil (lPtVis,lPhiVis,lSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices,
607     iPhi1,iEta1,iPhi2,iEta2);
608     Met lPCRec = fRecoilTools->PUCRecoil (lPtVis,lPhiVis,lSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices,
609     iPhi1,iEta1,iPhi2,iEta2);
610     Met lPUMet = fRecoilTools->PUMet ( iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices,
611     iPhi1,iEta1,iPhi2,iEta2);
612 pharris 1.1
613     Double_t lPt0 = 0; const PFJet *lLead = 0;
614     Double_t lPt1 = 0; const PFJet *l2nd = 0;
615     int lNAllJet = 0;
616     int lNJet = 0;
617     for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
618     const PFJet *pJet = iJets->At(i0);
619     Double_t pPt = fRecoilTools->fJetIDMVA->correctedPt(pJet,iJetCorrector,iPUEnergyDensity);
620     double pDEta1 = pJet->Eta() - iEta1;
621     double pDPhi1 = fabs(pJet->Phi() - iPhi1); if(pDPhi1 > 2.*TMath::Pi()-pDPhi1) pDPhi1 = 2.*TMath::Pi()-pDPhi1;
622     double pDR1 = sqrt(pDEta1*pDEta1 + pDPhi1*pDPhi1);
623     if(pDR1 < 0.5) continue;
624     double pDEta2 = pJet->Eta() - iEta2;
625     double pDPhi2 = fabs(pJet->Phi() - iPhi2); if(pDPhi2 > 2.*TMath::Pi()-pDPhi2) pDPhi2 = 2.*TMath::Pi()-pDPhi2;
626     double pDR2 = sqrt(pDEta2*pDEta2 + pDPhi2*pDPhi2);
627     if(pDR2 < 0.5) continue;
628     if(!JetTools::passPFLooseId(pJet)) continue;
629     lNAllJet++;
630     if(pPt > 30.) lNJet++;
631     if(lPt0 < pPt) {lPt0 = pPt; lLead = pJet; continue;}
632     if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
633     }
634    
635 pharris 1.5 fSumEt = lPFRec.SumEt();
636 pharris 1.1 fU = lPFRec.Pt();
637     fUPhi = lPFRec.Phi();
638     fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
639     fTKU = lTKRec.Pt();
640     fTKUPhi = lTKRec.Phi();
641     fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
642     fNPU = lNPRec.Pt();
643     fNPUPhi = lNPRec.Phi();
644     fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
645     fPUMet = lPUMet.Pt();
646     fPUMetPhi= lPUMet.Phi();
647     fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
648     fPCU = lPCRec.Pt() ;
649     fPCUPhi = lPCRec.Phi() ;
650     fJSPt1 = lPt0;
651     fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
652     fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
653     fJSPt2 = lPt1;
654     fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
655     fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
656     fNJet = lNJet ;
657     fNAllJet = lNAllJet;
658     fNPV = iNPV ;
659 pharris 1.5
660 pharris 1.2 Float_t lMVA = evaluatePhi();//fPhiReader ->EvaluateMVA( fPhiMethodName );
661     if(!iPhi) fUPhiMVA = fUPhi + lMVA;
662 pharris 1.4 //fTKSumEt /= lPFRec.SumEt();
663     //fNPSumEt /= lPFRec.SumEt();
664     //fPUSumEt /= lPFRec.SumEt();
665     //fPCSumEt /= lPFRec.SumEt();
666 pharris 1.7 if(!iPhi) lMVA = evaluateU1();//fU1Reader ->EvaluateMVA( fU1MethodName );
667     fUMVA = fU*lMVA;
668    
669     TLorentzVector lUVec (0,0,0,0); lUVec .SetPtEtaPhiM(fUMVA ,0,fUPhiMVA,0);
670 pharris 1.2 TLorentzVector lVVec (0,0,0,0); lVVec .SetPtEtaPhiM(lPtVis ,0,lPhiVis ,0);
671     if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
672     lUVec -= lVVec;
673     Met lMet(lUVec.Px(),lUVec.Py());
674 pharris 1.4 //Cov matrix
675     double lCovU1 = evaluateCovU1();
676     double lCovU2 = evaluateCovU2();
677    
678     double lCos2 = cos(fUPhiMVA)*cos(fUPhiMVA);
679     double lSin2 = sin(fUPhiMVA)*sin(fUPhiMVA);
680    
681     //Now Compute teh covariance matrix in X and Y
682 pharris 1.5 (*fCov)(0,0) = lCovU1*lCos2+lCovU2*lSin2;
683     (*fCov)(1,0) = -lCovU1*sin(fUPhiMVA)*cos(fUPhiMVA)+lCovU2*sin(fUPhiMVA)*cos(fUPhiMVA);
684     (*fCov)(0,1) = (*fCov)(1,0);
685     (*fCov)(1,1) = lCovU1*lSin2+lCovU2*lCos2;
686 pharris 1.10 TMatrixD lInv(2,2);
687     lInv(0,0) = (*fCov)(0,0); lInv(1,1) = (*fCov)(1,1); lInv(1,0) = (*fCov)(1,0); lInv(0,1) = (*fCov)(0,1);
688     if(lInv.Determinant() != 0) lInv.Invert();
689     fSignificance = TMath::Sqrt(lUVec.Px()*lUVec.Px()*(lInv)(0,0) + 2.*lUVec.Px()*lUVec.Py()*(lInv)(1,0) + lUVec.Py()*lUVec.Py()*(lInv)(1,1));
690     fUncertainty = sqrt(lCovU1+lCovU2);
691    
692 pharris 1.2
693     if (printDebug == kTRUE) {
694 ceballos 1.11 std::cout << "Debug Met MVA: "
695 pharris 1.2 << fU << " : "
696     << fUPhi << " : "
697     << fTKSumEt << " : "
698     << fTKU << " : "
699     << fTKUPhi << " : "
700     << fNPSumEt << " : "
701     << fNPU << " : "
702     << fNPUPhi << " : "
703     << fPUSumEt << " : "
704     << fPUMet << " : "
705     << fPUMetPhi << " : "
706     << fPCSumEt << " : "
707     << fPCU << " : "
708     << fPCUPhi << " : "
709     << fJSPt1 << " : "
710     << fJSEta1 << " : "
711     << fJSPhi1 << " : "
712     << fJSPt2 << " : "
713     << fJSEta2 << " : "
714     << fJSPhi2 << " : "
715     << fNJet << " : "
716     << fNAllJet << " : "
717     << fNPV << " : "
718     << " === : === "
719     << std::endl;
720     }
721    
722     return lMet;
723     }
724     //--------------------------------------------------------------------------------------------------
725     //Interms of the two candidates => corrected Jets
726     Met MVAMet::GetMet( Bool_t iPhi,
727     Float_t iPt1,Float_t iPhi1,Float_t iEta1,
728     Float_t iPt2,Float_t iPhi2,Float_t iEta2,
729     const PFMet *iMet ,
730 pharris 1.3 const PFCandidateCol *iCands,
731 pharris 1.7 const Vertex *iVertex,const VertexCol *iVertices,Double_t iRho,
732 pharris 1.2 const PFJetCol *iJets ,
733     int iNPV,
734     Bool_t printDebug) {
735    
736     TLorentzVector lVVec1(0,0,0,0); lVVec1.SetPtEtaPhiM(iPt1,0,iPhi1 ,0);
737     TLorentzVector lVVec2(0,0,0,0); lVVec2.SetPtEtaPhiM(iPt2,0,iPhi2 ,0);
738     lVVec1+=lVVec2;
739     Float_t lPtVis = lVVec1.Pt();
740     Float_t lPhiVis = lVVec1.Phi();
741     Float_t lSumEtVis = iPt1 + iPt2;
742 pharris 1.7 Met lPFRec = fRecoilTools->pfRecoil (lPtVis,lPhiVis,lSumEtVis,iCands);
743 pharris 1.2 Met lTKRec = fRecoilTools->trackRecoil(lPtVis,lPhiVis,lSumEtVis,iCands,iVertex);
744 pharris 1.3 Met lNPRec = fRecoilTools->NoPURecoil (lPtVis,lPhiVis,lSumEtVis,iJets,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2);
745     Met lPCRec = fRecoilTools->PUCRecoil (lPtVis,lPhiVis,lSumEtVis,iJets,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2);
746     Met lPUMet = fRecoilTools->PUMet ( iJets,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2);
747 pharris 1.2
748     Double_t lPt0 = 0; const PFJet *lLead = 0;
749     Double_t lPt1 = 0; const PFJet *l2nd = 0;
750     int lNAllJet = 0;
751     int lNJet = 0;
752     for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
753     const PFJet *pJet = iJets->At(i0);
754     Double_t pPt = pJet->Pt();
755     if( pJet->TrackCountingHighEffBJetTagsDisc() == -100 && pPt < 10.) continue;
756     double pDEta1 = pJet->Eta() - iEta1;
757     double pDPhi1 = fabs(pJet->Phi() - iPhi1); if(pDPhi1 > 2.*TMath::Pi()-pDPhi1) pDPhi1 = 2.*TMath::Pi()-pDPhi1;
758     double pDR1 = sqrt(pDEta1*pDEta1 + pDPhi1*pDPhi1);
759     if(pDR1 < 0.5) continue;
760     double pDEta2 = pJet->Eta() - iEta2;
761     double pDPhi2 = fabs(pJet->Phi() - iPhi2); if(pDPhi2 > 2.*TMath::Pi()-pDPhi2) pDPhi2 = 2.*TMath::Pi()-pDPhi2;
762     double pDR2 = sqrt(pDEta2*pDEta2 + pDPhi2*pDPhi2);
763     if(pDR2 < 0.5) continue;
764     if(!JetTools::passPFLooseId(pJet)) continue;
765     lNAllJet++;
766     if(pPt > 30.) lNJet++;
767     if(lPt0 < pPt) {lPt0 = pPt; lLead = pJet; continue;}
768     if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
769     }
770 pharris 1.5 fSumEt = lPFRec.SumEt();
771 pharris 1.2 fU = lPFRec.Pt();
772     fUPhi = lPFRec.Phi();
773     fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
774     fTKU = lTKRec.Pt();
775     fTKUPhi = lTKRec.Phi();
776     fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
777     fNPU = lNPRec.Pt();
778     fNPUPhi = lNPRec.Phi();
779     fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
780     fPUMet = lPUMet.Pt();
781     fPUMetPhi= lPUMet.Phi();
782     fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
783     fPCU = lPCRec.Pt() ;
784     fPCUPhi = lPCRec.Phi() ;
785     fJSPt1 = lPt0;
786     fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
787     fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
788     fJSPt2 = lPt1;
789     fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
790     fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
791     fNJet = lNJet ;
792     fNAllJet = lNAllJet;
793     fNPV = iNPV ;
794    
795     Float_t lMVA = evaluatePhi();//fPhiReader ->EvaluateMVA( fPhiMethodName );
796 pharris 1.1
797     if(!iPhi) fUPhiMVA = fUPhi + lMVA;
798 pharris 1.4 //fTKSumEt /= lPFRec.SumEt();
799     //fNPSumEt /= lPFRec.SumEt();
800     //fPUSumEt /= lPFRec.SumEt();
801     //fPCSumEt /= lPFRec.SumEt();
802 pharris 1.7 if(!iPhi) lMVA = evaluateU1();//fU1Reader ->EvaluateMVA( fU1MethodName );
803     fUMVA = lMVA*fU;
804    
805     TLorentzVector lUVec (0,0,0,0); lUVec .SetPtEtaPhiM(fUMVA ,0,fUPhiMVA,0);
806 pharris 1.1 TLorentzVector lVVec (0,0,0,0); lVVec .SetPtEtaPhiM(lPtVis ,0,lPhiVis ,0);
807 pharris 1.7 if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
808 pharris 1.1 lUVec -= lVVec;
809     Met lMet(lUVec.Px(),lUVec.Py());
810 pharris 1.4 //Cov matrix
811     double lCovU1 = evaluateCovU1();
812     double lCovU2 = evaluateCovU2();
813    
814     double lCos2 = cos(fUPhiMVA)*cos(fUPhiMVA);
815     double lSin2 = sin(fUPhiMVA)*sin(fUPhiMVA);
816     //Now Compute teh covariance matrix in X and Y
817 pharris 1.5 (*fCov)(0,0) = lCovU1*lCos2+lCovU2*lSin2;
818     (*fCov)(1,0) = -lCovU1*sin(fUPhiMVA)*cos(fUPhiMVA)+lCovU2*sin(fUPhiMVA)*cos(fUPhiMVA);
819     (*fCov)(0,1) = (*fCov)(1,0);
820     (*fCov)(1,1) = lCovU1*lSin2+lCovU2*lCos2;
821 pharris 1.10 TMatrixD lInv(2,2);
822     lInv(0,0) = (*fCov)(0,0); lInv(1,1) = (*fCov)(1,1); lInv(1,0) = (*fCov)(1,0); lInv(0,1) = (*fCov)(0,1);
823     if(lInv.Determinant() != 0) lInv.Invert();
824     fSignificance = TMath::Sqrt(lUVec.Px()*lUVec.Px()*(lInv)(0,0) + 2.*lUVec.Px()*lUVec.Py()*(lInv)(1,0) + lUVec.Py()*lUVec.Py()*(lInv)(1,1));
825     fUncertainty = sqrt(lCovU1+lCovU2);
826 pharris 1.7
827 pharris 1.1 if (printDebug == kTRUE) {
828 ceballos 1.11 std::cout << "Debug Met MVA: "
829 pharris 1.1 << fU << " : "
830     << fUPhi << " : "
831     << fTKSumEt << " : "
832     << fTKU << " : "
833     << fTKUPhi << " : "
834     << fNPSumEt << " : "
835     << fNPU << " : "
836     << fNPUPhi << " : "
837     << fPUSumEt << " : "
838     << fPUMet << " : "
839     << fPUMetPhi << " : "
840     << fPCSumEt << " : "
841     << fPCU << " : "
842     << fPCUPhi << " : "
843     << fJSPt1 << " : "
844     << fJSEta1 << " : "
845     << fJSPhi1 << " : "
846     << fJSPt2 << " : "
847     << fJSEta2 << " : "
848     << fJSPhi2 << " : "
849     << fNJet << " : "
850     << fNAllJet << " : "
851     << fNPV << " : "
852     << " === : === "
853     << std::endl;
854     }
855    
856     return lMet;
857     }
858