ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MVAMet.cc
Revision: 1.1
Committed: Wed Mar 21 18:56:26 2012 UTC (13 years, 1 month ago) by pharris
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_025e, Mit_025d
Log Message:
Adding MET regression

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     #include "TMVA/Tools.h"
8     #include "TMVA/Reader.h"
9    
10    
11     ClassImp(mithep::MVAMet)
12    
13     using namespace mithep;
14    
15     //--------------------------------------------------------------------------------------------------
16     MVAMet::MVAMet() :
17     fRecoilTools(0),
18     fPhiMethodName("PhiMVAMet"),
19     fU1MethodName ("U1MVAMet"),
20     fIsInitialized(kFALSE),
21     fU (0),
22     fUPhi (0),
23     fTKSumEt(0),
24     fTKU (0),
25     fTKUPhi (0),
26     fNPSumEt(0),
27     fNPU (0),
28     fNPUPhi (0),
29     fPUSumEt(0),
30     fPUMet (0),
31     fPUMetPhi(0),
32     fPCSumEt(0),
33     fPCU (0),
34     fPCUPhi (0),
35     fJSPt1 (0),
36     fJSEta1 (0),
37     fJSPhi1 (0),
38     fJSPt2 (0),
39     fJSEta2 (0),
40     fJSPhi2 (0),
41     fNJet (0),
42     fNAllJet(0),
43     fNPV (0),
44     fPhiReader(0),
45     fU1Reader(0)
46     { }
47     //--------------------------------------------------------------------------------------------------
48     MVAMet::~MVAMet()
49     {
50     fPhiReader = 0;
51     fU1Reader = 0;
52     }
53     //--------------------------------------------------------------------------------------------------
54     void MVAMet::setVariables(TMVA::Reader *iReader,bool iScale) {
55     iReader->AddVariable( "npv" , &fNPV );
56     iReader->AddVariable( "u" , &fU );
57     iReader->AddVariable( "uphi" , &fUPhi );
58     iReader->AddVariable( "chsumet/sumet" , &fTKSumEt );
59     iReader->AddVariable( "tku" , &fTKU );
60     iReader->AddVariable( "tkuphi" , &fTKUPhi );
61     iReader->AddVariable( "nopusumet/sumet" , &fNPSumEt );
62     iReader->AddVariable( "nopuu" , &fNPU );
63     iReader->AddVariable( "nopuuphi" , &fNPUPhi );
64     iReader->AddVariable( "pusumet/sumet" , &fPUSumEt );
65     iReader->AddVariable( "pumet" , &fPUMet );
66     iReader->AddVariable( "pumetphi" , &fPUMetPhi);
67     iReader->AddVariable( "pucsumet/sumet" , &fPCSumEt );
68     iReader->AddVariable( "pucu" , &fPCU );
69     iReader->AddVariable( "pucuphi" , &fPCUPhi );
70     iReader->AddVariable( "jspt_1" , &fJSPt1 );
71     iReader->AddVariable( "jseta_1" , &fJSEta1 );
72     iReader->AddVariable( "jsphi_1" , &fJSPhi1 );
73     iReader->AddVariable( "jspt_2" , &fJSPt2 );
74     iReader->AddVariable( "jseta_2" , &fJSEta2 );
75     iReader->AddVariable( "jsphi_2" , &fJSPhi2 );
76     iReader->AddVariable( "nalljet" , &fNAllJet );
77     iReader->AddVariable( "njet" , &fNJet );
78     if(iScale) iReader->AddVariable( "uphi_mva" , &fUPhiMVA );
79     }
80     //--------------------------------------------------------------------------------------------------
81     void MVAMet::Initialize( TString iU1MethodName,
82     TString iPhiMethodName,
83     TString iJetMVAFile,
84     TString iU1Weights,
85     TString iPhiWeights,
86     MVAMet::MVAType iType) {
87    
88     fIsInitialized = kTRUE;
89     fRecoilTools = new RecoilTools(iJetMVAFile);
90    
91     fU1MethodName = iU1MethodName;
92     fPhiMethodName = iPhiMethodName;
93     fType = iType;
94     fPhiReader = new TMVA::Reader( "!Color:!Silent:Error" );
95     fU1Reader = new TMVA::Reader( "!Color:!Silent:Error" );
96     if (fType == kBaseline) {
97     setVariables(fPhiReader,false);
98     setVariables(fU1Reader ,true);
99     }
100     fPhiReader->BookMVA(fPhiMethodName , iPhiWeights );
101     fU1Reader ->BookMVA(fU1MethodName , iU1Weights );
102     std::cout << "Jet ID MVA Initialization\n";
103     std::cout << "Phi Method Name : " << fPhiMethodName << " , type == " << iType << std::endl;
104     std::cout << "U1 Method Name : " << fU1MethodName << " , type == " << iType << std::endl;
105     }
106    
107     //--------------------------------------------------------------------------------------------------
108     Double_t MVAMet::MVAValue( Bool_t iPhi,
109     Float_t iPFSumEt,
110     Float_t iU ,
111     Float_t iUPhi ,
112     Float_t iTKSumEt,
113     Float_t iTKU ,
114     Float_t iTKUPhi ,
115     Float_t iNPSumEt,
116     Float_t iNPU ,
117     Float_t iNPUPhi ,
118     Float_t iPUSumEt,
119     Float_t iPUMet ,
120     Float_t iPUMetPhi,
121     Float_t iPCSumEt,
122     Float_t iPCU ,
123     Float_t iPCUPhi ,
124     Float_t iJSPt1 ,
125     Float_t iJSEta1 ,
126     Float_t iJSPhi1 ,
127     Float_t iJSPt2 ,
128     Float_t iJSEta2 ,
129     Float_t iJSPhi2 ,
130     Float_t iNAllJet,
131     Float_t iNJet ,
132     Float_t iNPV
133     ){
134    
135     if (!fIsInitialized) {
136     std::cout << "Error: MVA Met not properly initialized.\n";
137     return -9999;
138     }
139    
140     fU = iU ;
141     fUPhi = iUPhi ;
142     fTKSumEt = iTKSumEt/iPFSumEt;
143     fTKU = iTKU ;
144     fTKUPhi = iTKUPhi ;
145     fNPSumEt = iNPSumEt/iPFSumEt;
146     fNPU = iNPU ;
147     fNPUPhi = iNPUPhi ;
148     fPUSumEt = iPUSumEt/iPFSumEt;
149     fPUMet = iPUMet ;
150     fPUMetPhi = iPUMetPhi;
151     fPCSumEt = iPCSumEt/iPFSumEt;
152     fPCU = iPCU ;
153     fPCUPhi = iPCUPhi ;
154     fJSPt1 = iJSPt1 ;
155     fJSEta1 = iJSEta1 ;
156     fJSPhi1 = iJSPhi1 ;
157     fJSPt2 = iJSPt2 ;
158     fJSEta2 = iJSEta2 ;
159     fJSPhi2 = iJSPhi2 ;
160     fNAllJet = iNAllJet;
161     fNJet = iNJet ;
162     fNPV = iNPV ;
163    
164     Double_t lMVA = -9999;
165     lMVA = fPhiReader->EvaluateMVA( fPhiMethodName );
166     if(!iPhi) fUPhiMVA = iUPhi+lMVA;
167     //Not no nice feature of teh training
168     fTKSumEt /= iPFSumEt;
169     fNPSumEt /= iPFSumEt;
170     fPUSumEt /= iPFSumEt;
171     fPCSumEt /= iPFSumEt;
172     if(!iPhi) lMVA = fU1Reader ->EvaluateMVA( fU1MethodName );
173    
174     return lMVA;
175     }
176     //--------------------------------------------------------------------------------------------------
177     //====> Please not that the jet collection must be cleaned => all jets near leptons must be removed
178     Met MVAMet::GetMet( Bool_t iPhi,Float_t iPtVis,Float_t iPhiVis,Float_t iSumEtVis,
179     const PFMet *iMet ,
180     const PFCandidateCol *iCands,const Vertex *iVertex,
181     const PFJetCol *iJets ,
182     FactorizedJetCorrector *iJetCorrector,
183     const PileupEnergyDensityCol *iPUEnergyDensity,
184     int iNPV,
185     Bool_t printDebug) {
186    
187     Met lPFRec = fRecoilTools->pfRecoil (iPtVis,iPhiVis,iSumEtVis,iMet);
188     Met lTKRec = fRecoilTools->trackRecoil(iPtVis,iPhiVis,iSumEtVis,iCands,iVertex);
189     Met lNPRec = fRecoilTools->NoPURecoil (iPtVis,iPhiVis,iSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex);
190     Met lPCRec = fRecoilTools->PUCRecoil (iPtVis,iPhiVis,iSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex);
191     Met lPUMet = fRecoilTools->PUMet (iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex);
192    
193     Double_t lPt0 = 0; const PFJet *lLead = 0;
194     Double_t lPt1 = 0; const PFJet *l2nd = 0;
195     int lNAllJet = 0;
196     int lNJet = 0;
197     for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
198     const PFJet *pJet = iJets->At(i0);
199     Double_t pPt = fRecoilTools->fJetIDMVA->correctedPt(pJet,iJetCorrector,iPUEnergyDensity);
200     if( pJet->TrackCountingHighEffBJetTagsDisc() == -100 && pPt < 10.) continue;
201     if(!JetTools::passPFLooseId(pJet)) continue;
202     lNAllJet++;
203     if(pPt > 30) lNJet++;
204     if(lPt0 < pPt) {lPt0 = pPt; lLead = pJet; continue;}
205     if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
206     }
207    
208     fU = lPFRec.Pt();
209     fUPhi = lPFRec.Phi();
210     fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
211     fTKU = lTKRec.Pt();
212     fTKUPhi = lTKRec.Phi();
213     fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
214     fNPU = lNPRec.Pt();
215     fNPUPhi = lNPRec.Phi();
216     fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
217     fPUMet = lPUMet.Pt();
218     fPUMetPhi = lPUMet.Phi();
219     fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
220     fPCU = lPCRec.Pt() ;
221     fPCUPhi = lPCRec.Phi() ;
222     fJSPt1 = lPt0;
223     fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
224     fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
225     fJSPt2 = lPt1;
226     fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
227     fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
228     fNJet = lNJet ;
229     fNAllJet = lNAllJet;
230     fNPV = iNPV ;
231    
232     Float_t lMVA = fPhiReader ->EvaluateMVA( fPhiMethodName );
233    
234     if(!iPhi) fUPhiMVA = fUPhi + lMVA;
235     //Not no nice feature of teh training
236     fTKSumEt /= lPFRec.SumEt();
237     fNPSumEt /= lPFRec.SumEt();
238     fPUSumEt /= lPFRec.SumEt();
239     fPCSumEt /= lPFRec.SumEt();
240     if(!iPhi) lMVA = fU1Reader ->EvaluateMVA( fU1MethodName );
241    
242     TLorentzVector lUVec(0,0,0,0); lUVec .SetPtEtaPhiM(fU*lMVA,0,fUPhiMVA,0);
243     TLorentzVector lVVec(0,0,0,0); lVVec .SetPtEtaPhiM(iPtVis ,0,iPhiVis ,0);
244     if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
245     lUVec -= lVVec;
246     Met lMet(lUVec.Px(),lUVec.Py());
247    
248     if (printDebug == kTRUE) {
249     std::cout << "Debug Jet MVA: "
250     << fU << " : "
251     << fUPhi << " : "
252     << fTKSumEt << " : "
253     << fTKU << " : "
254     << fTKUPhi << " : "
255     << fNPSumEt << " : "
256     << fNPU << " : "
257     << fNPUPhi << " : "
258     << fPUSumEt << " : "
259     << fPUMet << " : "
260     << fPUMetPhi << " : "
261     << fPCUPhi << " : "
262     << fPCSumEt << " : "
263     << fPCU << " : "
264     << fPCUPhi << " : "
265     << fJSPt1 << " : "
266     << fJSEta1 << " : "
267     << fJSPhi1 << " : "
268     << fJSPt2 << " : "
269     << fJSEta2 << " : "
270     << fJSPhi2 << " : "
271     << fNJet << " : "
272     << fNAllJet << " : "
273     << fNPV << " : "
274     << " === : === "
275     << std::endl;
276     }
277    
278     return lMet;
279     }
280     //Interms of the two candidates
281     Met MVAMet::GetMet( Bool_t iPhi,
282     Float_t iPt1,Float_t iPhi1,Float_t iEta1,
283     Float_t iPt2,Float_t iPhi2,Float_t iEta2,
284     const PFMet *iMet ,
285     const PFCandidateCol *iCands,const Vertex *iVertex,
286     const PFJetCol *iJets ,
287     FactorizedJetCorrector *iJetCorrector,
288     const PileupEnergyDensityCol *iPUEnergyDensity,
289     int iNPV,
290     Bool_t printDebug) {
291    
292     TLorentzVector lVVec1(0,0,0,0); lVVec1.SetPtEtaPhiM(iPt1,0,iPhi1 ,0);
293     TLorentzVector lVVec2(0,0,0,0); lVVec2.SetPtEtaPhiM(iPt2,0,iPhi2 ,0);
294     lVVec1+=lVVec2;
295     Float_t lPtVis = lVVec1.Pt();
296     Float_t lPhiVis = lVVec1.Phi();
297     Float_t lSumEtVis = iPt1 + iPt2;
298     Met lPFRec = fRecoilTools->pfRecoil (lPtVis,lPhiVis,lSumEtVis,iMet);
299     Met lTKRec = fRecoilTools->trackRecoil(lPtVis,lPhiVis,lSumEtVis,iCands,iVertex);
300     Met lNPRec = fRecoilTools->NoPURecoil (lPtVis,lPhiVis,lSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iPhi1,iEta1,iPhi2,iEta2);
301     Met lPCRec = fRecoilTools->PUCRecoil (lPtVis,lPhiVis,lSumEtVis,iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iPhi1,iEta1,iPhi2,iEta2);
302     Met lPUMet = fRecoilTools->PUMet ( iJets,iJetCorrector,iPUEnergyDensity,iCands,iVertex,iPhi1,iEta1,iPhi2,iEta2);
303    
304     Double_t lPt0 = 0; const PFJet *lLead = 0;
305     Double_t lPt1 = 0; const PFJet *l2nd = 0;
306     int lNAllJet = 0;
307     int lNJet = 0;
308     for(unsigned int i0 = 0; i0 < iJets->GetEntries(); i0++) {
309     const PFJet *pJet = iJets->At(i0);
310     Double_t pPt = fRecoilTools->fJetIDMVA->correctedPt(pJet,iJetCorrector,iPUEnergyDensity);
311     if( pJet->TrackCountingHighEffBJetTagsDisc() == -100 && pPt < 10.) continue;
312     double pDEta1 = pJet->Eta() - iEta1;
313     double pDPhi1 = fabs(pJet->Phi() - iPhi1); if(pDPhi1 > 2.*TMath::Pi()-pDPhi1) pDPhi1 = 2.*TMath::Pi()-pDPhi1;
314     double pDR1 = sqrt(pDEta1*pDEta1 + pDPhi1*pDPhi1);
315     if(pDR1 < 0.5) continue;
316     double pDEta2 = pJet->Eta() - iEta2;
317     double pDPhi2 = fabs(pJet->Phi() - iPhi2); if(pDPhi2 > 2.*TMath::Pi()-pDPhi2) pDPhi2 = 2.*TMath::Pi()-pDPhi2;
318     double pDR2 = sqrt(pDEta2*pDEta2 + pDPhi2*pDPhi2);
319     if(pDR2 < 0.5) continue;
320     if(!JetTools::passPFLooseId(pJet)) continue;
321     lNAllJet++;
322     if(pPt > 30.) lNJet++;
323     if(lPt0 < pPt) {lPt0 = pPt; lLead = pJet; continue;}
324     if(lPt1 < pPt) {lPt1 = pPt; l2nd = pJet; continue;}
325     }
326    
327     fU = lPFRec.Pt();
328     fUPhi = lPFRec.Phi();
329     fTKSumEt = lTKRec.SumEt()/lPFRec.SumEt();
330     fTKU = lTKRec.Pt();
331     fTKUPhi = lTKRec.Phi();
332     fNPSumEt = lNPRec.SumEt()/lPFRec.SumEt();
333     fNPU = lNPRec.Pt();
334     fNPUPhi = lNPRec.Phi();
335     fPUSumEt = lPUMet.SumEt()/lPFRec.SumEt();
336     fPUMet = lPUMet.Pt();
337     fPUMetPhi= lPUMet.Phi();
338     fPCSumEt = lPCRec.SumEt()/lPFRec.SumEt();
339     fPCU = lPCRec.Pt() ;
340     fPCUPhi = lPCRec.Phi() ;
341     fJSPt1 = lPt0;
342     fJSEta1 = 0; if(lLead != 0) fJSEta1 = lLead->Eta();
343     fJSPhi1 = 0; if(lLead != 0) fJSPhi1 = lLead->Phi();
344     fJSPt2 = lPt1;
345     fJSEta2 = 0; if(l2nd != 0) fJSEta2 = l2nd ->Eta();
346     fJSPhi2 = 0; if(l2nd != 0) fJSPhi2 = l2nd ->Phi();
347     fNJet = lNJet ;
348     fNAllJet = lNAllJet;
349     fNPV = iNPV ;
350    
351     Float_t lMVA = fPhiReader ->EvaluateMVA( fPhiMethodName );
352    
353     if(!iPhi) fUPhiMVA = fUPhi + lMVA;
354     fTKSumEt /= lPFRec.SumEt();
355     fNPSumEt /= lPFRec.SumEt();
356     fPUSumEt /= lPFRec.SumEt();
357     fPCSumEt /= lPFRec.SumEt();
358     if(!iPhi) lMVA = fU1Reader ->EvaluateMVA( fU1MethodName );
359    
360     TLorentzVector lUVec (0,0,0,0); lUVec .SetPtEtaPhiM(fU*lMVA,0,fUPhiMVA,0);
361     TLorentzVector lVVec (0,0,0,0); lVVec .SetPtEtaPhiM(lPtVis ,0,lPhiVis ,0);
362     if(lMVA < 0) lUVec .RotateZ(TMath::Pi());
363     lUVec -= lVVec;
364     Met lMet(lUVec.Px(),lUVec.Py());
365    
366     if (printDebug == kTRUE) {
367     std::cout << "Debug Jet MVA: "
368     << fU << " : "
369     << fUPhi << " : "
370     << fTKSumEt << " : "
371     << fTKU << " : "
372     << fTKUPhi << " : "
373     << fNPSumEt << " : "
374     << fNPU << " : "
375     << fNPUPhi << " : "
376     << fPUSumEt << " : "
377     << fPUMet << " : "
378     << fPUMetPhi << " : "
379     << fPCSumEt << " : "
380     << fPCU << " : "
381     << fPCUPhi << " : "
382     << fJSPt1 << " : "
383     << fJSEta1 << " : "
384     << fJSPhi1 << " : "
385     << fJSPt2 << " : "
386     << fJSEta2 << " : "
387     << fJSPhi2 << " : "
388     << fNJet << " : "
389     << fNAllJet << " : "
390     << fNPV << " : "
391     << " === : === "
392     << std::endl;
393     }
394    
395     return lMet;
396     }
397