ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MVAMet.cc
Revision: 1.8
Committed: Tue Apr 24 23:26:17 2012 UTC (13 years ago) by pharris
Content type: text/plain
Branch: MAIN
Changes since 1.7: +5 -5 lines
Log Message:
Fixed significance calculationo

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    
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.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));
545 pharris 1.7 fSignificance = lMet.Pt()/fUncertainty;
546 pharris 1.2
547     if (printDebug == kTRUE) {
548     std::cout << "Debug Jet MVA: "
549     << fU << " : "
550     << fUPhi << " : "
551     << fTKSumEt << " : "
552     << fTKU << " : "
553     << fTKUPhi << " : "
554     << fNPSumEt << " : "
555     << fNPU << " : "
556     << fNPUPhi << " : "
557     << fPUSumEt << " : "
558     << fPUMet << " : "
559     << fPUMetPhi << " : "
560     << fPCUPhi << " : "
561     << fPCSumEt << " : "
562     << fPCU << " : "
563     << fPCUPhi << " : "
564     << fJSPt1 << " : "
565     << fJSEta1 << " : "
566     << fJSPhi1 << " : "
567     << fJSPt2 << " : "
568     << fJSEta2 << " : "
569     << fJSPhi2 << " : "
570     << fNJet << " : "
571     << fNAllJet << " : "
572     << fNPV << " : "
573     << " === : === "
574     << std::endl;
575     }
576    
577     return lMet;
578     }
579     //--------------------------------------------------------------------------------------------------
580 pharris 1.1 //Interms of the two candidates
581     Met MVAMet::GetMet( Bool_t iPhi,
582     Float_t iPt1,Float_t iPhi1,Float_t iEta1,
583     Float_t iPt2,Float_t iPhi2,Float_t iEta2,
584     const PFMet *iMet ,
585 pharris 1.3 const PFCandidateCol *iCands,
586     const Vertex *iVertex,const VertexCol *iVertices,
587 pharris 1.1 const PFJetCol *iJets ,
588     FactorizedJetCorrector *iJetCorrector,
589     const PileupEnergyDensityCol *iPUEnergyDensity,
590     int iNPV,
591     Bool_t printDebug) {
592    
593     TLorentzVector lVVec1(0,0,0,0); lVVec1.SetPtEtaPhiM(iPt1,0,iPhi1 ,0);
594     TLorentzVector lVVec2(0,0,0,0); lVVec2.SetPtEtaPhiM(iPt2,0,iPhi2 ,0);
595     lVVec1+=lVVec2;
596     Float_t lPtVis = lVVec1.Pt();
597     Float_t lPhiVis = lVVec1.Phi();
598     Float_t lSumEtVis = iPt1 + iPt2;
599 pharris 1.7 Met lPFRec = fRecoilTools->pfRecoil (lPtVis,lPhiVis,lSumEtVis,iCands);
600 pharris 1.1 Met lTKRec = fRecoilTools->trackRecoil(lPtVis,lPhiVis,lSumEtVis,iCands,iVertex);
601 pharris 1.3 Met lNPRec = fRecoilTools->NoPURecoil (lPtVis,lPhiVis,lSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices,
602     iPhi1,iEta1,iPhi2,iEta2);
603     Met lPCRec = fRecoilTools->PUCRecoil (lPtVis,lPhiVis,lSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices,
604     iPhi1,iEta1,iPhi2,iEta2);
605     Met lPUMet = fRecoilTools->PUMet ( iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iVertices,
606     iPhi1,iEta1,iPhi2,iEta2);
607 pharris 1.1
608     Double_t lPt0 = 0; const PFJet *lLead = 0;
609     Double_t lPt1 = 0; const PFJet *l2nd = 0;
610     int lNAllJet = 0;
611     int lNJet = 0;
612     for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
613     const PFJet *pJet = iJets->At(i0);
614     Double_t pPt = fRecoilTools->fJetIDMVA->correctedPt(pJet,iJetCorrector,iPUEnergyDensity);
615     double pDEta1 = pJet->Eta() - iEta1;
616     double pDPhi1 = fabs(pJet->Phi() - iPhi1); if(pDPhi1 > 2.*TMath::Pi()-pDPhi1) pDPhi1 = 2.*TMath::Pi()-pDPhi1;
617     double pDR1 = sqrt(pDEta1*pDEta1 + pDPhi1*pDPhi1);
618     if(pDR1 < 0.5) continue;
619     double pDEta2 = pJet->Eta() - iEta2;
620     double pDPhi2 = fabs(pJet->Phi() - iPhi2); if(pDPhi2 > 2.*TMath::Pi()-pDPhi2) pDPhi2 = 2.*TMath::Pi()-pDPhi2;
621     double pDR2 = sqrt(pDEta2*pDEta2 + pDPhi2*pDPhi2);
622     if(pDR2 < 0.5) continue;
623     if(!JetTools::passPFLooseId(pJet)) continue;
624     lNAllJet++;
625     if(pPt > 30.) lNJet++;
626     if(lPt0 < pPt) {lPt0 = pPt; lLead = pJet; continue;}
627     if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
628     }
629    
630 pharris 1.5 fSumEt = lPFRec.SumEt();
631 pharris 1.1 fU = lPFRec.Pt();
632     fUPhi = lPFRec.Phi();
633     fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
634     fTKU = lTKRec.Pt();
635     fTKUPhi = lTKRec.Phi();
636     fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
637     fNPU = lNPRec.Pt();
638     fNPUPhi = lNPRec.Phi();
639     fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
640     fPUMet = lPUMet.Pt();
641     fPUMetPhi= lPUMet.Phi();
642     fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
643     fPCU = lPCRec.Pt() ;
644     fPCUPhi = lPCRec.Phi() ;
645     fJSPt1 = lPt0;
646     fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
647     fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
648     fJSPt2 = lPt1;
649     fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
650     fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
651     fNJet = lNJet ;
652     fNAllJet = lNAllJet;
653     fNPV = iNPV ;
654 pharris 1.5
655 pharris 1.2 Float_t lMVA = evaluatePhi();//fPhiReader ->EvaluateMVA( fPhiMethodName );
656     if(!iPhi) fUPhiMVA = fUPhi + lMVA;
657 pharris 1.4 //fTKSumEt /= lPFRec.SumEt();
658     //fNPSumEt /= lPFRec.SumEt();
659     //fPUSumEt /= lPFRec.SumEt();
660     //fPCSumEt /= lPFRec.SumEt();
661 pharris 1.7 if(!iPhi) lMVA = evaluateU1();//fU1Reader ->EvaluateMVA( fU1MethodName );
662     fUMVA = fU*lMVA;
663    
664     TLorentzVector lUVec (0,0,0,0); lUVec .SetPtEtaPhiM(fUMVA ,0,fUPhiMVA,0);
665 pharris 1.2 TLorentzVector lVVec (0,0,0,0); lVVec .SetPtEtaPhiM(lPtVis ,0,lPhiVis ,0);
666     if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
667     lUVec -= lVVec;
668     Met lMet(lUVec.Px(),lUVec.Py());
669 pharris 1.4 //Cov matrix
670     double lCovU1 = evaluateCovU1();
671     double lCovU2 = evaluateCovU2();
672    
673     double lCos2 = cos(fUPhiMVA)*cos(fUPhiMVA);
674     double lSin2 = sin(fUPhiMVA)*sin(fUPhiMVA);
675    
676     //Now Compute teh covariance matrix in X and Y
677 pharris 1.5 (*fCov)(0,0) = lCovU1*lCos2+lCovU2*lSin2;
678     (*fCov)(1,0) = -lCovU1*sin(fUPhiMVA)*cos(fUPhiMVA)+lCovU2*sin(fUPhiMVA)*cos(fUPhiMVA);
679     (*fCov)(0,1) = (*fCov)(1,0);
680     (*fCov)(1,1) = lCovU1*lSin2+lCovU2*lCos2;
681 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));
682 pharris 1.7 fSignificance = lMet.Pt()/fUncertainty;
683 pharris 1.2
684     if (printDebug == kTRUE) {
685     std::cout << "Debug Jet MVA: "
686     << fU << " : "
687     << fUPhi << " : "
688     << fTKSumEt << " : "
689     << fTKU << " : "
690     << fTKUPhi << " : "
691     << fNPSumEt << " : "
692     << fNPU << " : "
693     << fNPUPhi << " : "
694     << fPUSumEt << " : "
695     << fPUMet << " : "
696     << fPUMetPhi << " : "
697     << fPCSumEt << " : "
698     << fPCU << " : "
699     << fPCUPhi << " : "
700     << fJSPt1 << " : "
701     << fJSEta1 << " : "
702     << fJSPhi1 << " : "
703     << fJSPt2 << " : "
704     << fJSEta2 << " : "
705     << fJSPhi2 << " : "
706     << fNJet << " : "
707     << fNAllJet << " : "
708     << fNPV << " : "
709     << " === : === "
710     << std::endl;
711     }
712    
713     return lMet;
714     }
715     //--------------------------------------------------------------------------------------------------
716     //Interms of the two candidates => corrected Jets
717     Met MVAMet::GetMet( Bool_t iPhi,
718     Float_t iPt1,Float_t iPhi1,Float_t iEta1,
719     Float_t iPt2,Float_t iPhi2,Float_t iEta2,
720     const PFMet *iMet ,
721 pharris 1.3 const PFCandidateCol *iCands,
722 pharris 1.7 const Vertex *iVertex,const VertexCol *iVertices,Double_t iRho,
723 pharris 1.2 const PFJetCol *iJets ,
724     int iNPV,
725     Bool_t printDebug) {
726    
727     TLorentzVector lVVec1(0,0,0,0); lVVec1.SetPtEtaPhiM(iPt1,0,iPhi1 ,0);
728     TLorentzVector lVVec2(0,0,0,0); lVVec2.SetPtEtaPhiM(iPt2,0,iPhi2 ,0);
729     lVVec1+=lVVec2;
730     Float_t lPtVis = lVVec1.Pt();
731     Float_t lPhiVis = lVVec1.Phi();
732     Float_t lSumEtVis = iPt1 + iPt2;
733 pharris 1.7 Met lPFRec = fRecoilTools->pfRecoil (lPtVis,lPhiVis,lSumEtVis,iCands);
734 pharris 1.2 Met lTKRec = fRecoilTools->trackRecoil(lPtVis,lPhiVis,lSumEtVis,iCands,iVertex);
735 pharris 1.3 Met lNPRec = fRecoilTools->NoPURecoil (lPtVis,lPhiVis,lSumEtVis,iJets,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2);
736     Met lPCRec = fRecoilTools->PUCRecoil (lPtVis,lPhiVis,lSumEtVis,iJets,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2);
737     Met lPUMet = fRecoilTools->PUMet ( iJets,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2);
738 pharris 1.2
739     Double_t lPt0 = 0; const PFJet *lLead = 0;
740     Double_t lPt1 = 0; const PFJet *l2nd = 0;
741     int lNAllJet = 0;
742     int lNJet = 0;
743     for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
744     const PFJet *pJet = iJets->At(i0);
745     Double_t pPt = pJet->Pt();
746     if( pJet->TrackCountingHighEffBJetTagsDisc() == -100 && pPt < 10.) continue;
747     double pDEta1 = pJet->Eta() - iEta1;
748     double pDPhi1 = fabs(pJet->Phi() - iPhi1); if(pDPhi1 > 2.*TMath::Pi()-pDPhi1) pDPhi1 = 2.*TMath::Pi()-pDPhi1;
749     double pDR1 = sqrt(pDEta1*pDEta1 + pDPhi1*pDPhi1);
750     if(pDR1 < 0.5) continue;
751     double pDEta2 = pJet->Eta() - iEta2;
752     double pDPhi2 = fabs(pJet->Phi() - iPhi2); if(pDPhi2 > 2.*TMath::Pi()-pDPhi2) pDPhi2 = 2.*TMath::Pi()-pDPhi2;
753     double pDR2 = sqrt(pDEta2*pDEta2 + pDPhi2*pDPhi2);
754     if(pDR2 < 0.5) continue;
755     if(!JetTools::passPFLooseId(pJet)) continue;
756     lNAllJet++;
757     if(pPt > 30.) lNJet++;
758     if(lPt0 < pPt) {lPt0 = pPt; lLead = pJet; continue;}
759     if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
760     }
761 pharris 1.5 fSumEt = lPFRec.SumEt();
762 pharris 1.2 fU = lPFRec.Pt();
763     fUPhi = lPFRec.Phi();
764     fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
765     fTKU = lTKRec.Pt();
766     fTKUPhi = lTKRec.Phi();
767     fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
768     fNPU = lNPRec.Pt();
769     fNPUPhi = lNPRec.Phi();
770     fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
771     fPUMet = lPUMet.Pt();
772     fPUMetPhi= lPUMet.Phi();
773     fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
774     fPCU = lPCRec.Pt() ;
775     fPCUPhi = lPCRec.Phi() ;
776     fJSPt1 = lPt0;
777     fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
778     fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
779     fJSPt2 = lPt1;
780     fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
781     fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
782     fNJet = lNJet ;
783     fNAllJet = lNAllJet;
784     fNPV = iNPV ;
785    
786     Float_t lMVA = evaluatePhi();//fPhiReader ->EvaluateMVA( fPhiMethodName );
787 pharris 1.1
788     if(!iPhi) fUPhiMVA = fUPhi + lMVA;
789 pharris 1.4 //fTKSumEt /= lPFRec.SumEt();
790     //fNPSumEt /= lPFRec.SumEt();
791     //fPUSumEt /= lPFRec.SumEt();
792     //fPCSumEt /= lPFRec.SumEt();
793 pharris 1.7 if(!iPhi) lMVA = evaluateU1();//fU1Reader ->EvaluateMVA( fU1MethodName );
794     fUMVA = lMVA*fU;
795    
796     TLorentzVector lUVec (0,0,0,0); lUVec .SetPtEtaPhiM(fUMVA ,0,fUPhiMVA,0);
797 pharris 1.1 TLorentzVector lVVec (0,0,0,0); lVVec .SetPtEtaPhiM(lPtVis ,0,lPhiVis ,0);
798 pharris 1.7 if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
799 pharris 1.1 lUVec -= lVVec;
800     Met lMet(lUVec.Px(),lUVec.Py());
801 pharris 1.4 //Cov matrix
802     double lCovU1 = evaluateCovU1();
803     double lCovU2 = evaluateCovU2();
804    
805     double lCos2 = cos(fUPhiMVA)*cos(fUPhiMVA);
806     double lSin2 = sin(fUPhiMVA)*sin(fUPhiMVA);
807     //Now Compute teh covariance matrix in X and Y
808 pharris 1.5 (*fCov)(0,0) = lCovU1*lCos2+lCovU2*lSin2;
809     (*fCov)(1,0) = -lCovU1*sin(fUPhiMVA)*cos(fUPhiMVA)+lCovU2*sin(fUPhiMVA)*cos(fUPhiMVA);
810     (*fCov)(0,1) = (*fCov)(1,0);
811     (*fCov)(1,1) = lCovU1*lSin2+lCovU2*lCos2;
812 pharris 1.8
813     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));
814 pharris 1.7 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