ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/JetIDMVA.cc
Revision: 1.1
Committed: Wed Mar 21 18:56:25 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/JetIDMVA.h"
2     #include "MitPhysics/Utils/interface/JetTools.h"
3     #include "MitAna/DataTree/interface/StableData.h"
4     #include <TFile.h>
5     #include <TRandom3.h>
6     #include "TMVA/Tools.h"
7     #include "TMVA/Reader.h"
8    
9    
10     ClassImp(mithep::JetIDMVA)
11    
12     using namespace mithep;
13    
14     //--------------------------------------------------------------------------------------------------
15     JetIDMVA::JetIDMVA() :
16     fMethodName("JetIDMA"),
17     fIsInitialized(kFALSE),
18     fJetPtMin(10), //We need to lower this
19     fNPV (0),
20     fJPt1 (0),
21     fJEta1 (0),
22     fJPhi1 (0),
23     fJD01 (0),
24     fJDZ1 (0),
25     fJM1 (0),
26     fNPart1 (0),
27     fLPt1 (0),
28     fLEta1 (0),
29     fLPhi1 (0),
30     fSPt1 (0),
31     fSEta1 (0),
32     fSPhi1 (0),
33     fNEPt1 (0),
34     fNEEta1 (0),
35     fNEPhi1 (0),
36     fEMPt1 (0),
37     fEMEta1 (0),
38     fEMPhi1 (0),
39     fChPt1 (0),
40     fChPhi1 (0),
41     fLFr1 (0),
42     fDRLC1 (0),
43     fDRLS1 (0),
44     fDRM1 (0),
45     fDRNE1 (0),
46     fDREM1 (0),
47     fDRCH1 (0)
48     {
49     fReader = 0;
50     }
51     //--------------------------------------------------------------------------------------------------
52     JetIDMVA::~JetIDMVA()
53     {
54     fReader = 0;
55     }
56    
57     //--------------------------------------------------------------------------------------------------
58     void JetIDMVA::Initialize( TString iMethodName,
59     TString iWeights,
60     JetIDMVA::MVAType iType)
61     {
62    
63     fIsInitialized = kTRUE;
64     fMethodName = iMethodName;
65     fType = iType;
66     fReader = 0;
67     fReader = new TMVA::Reader( "!Color:!Silent:Error" );
68     if (fType == kBaseline) {
69     fReader->AddSpectator( "npv" , &fNPV );
70     fReader->AddVariable( "jspt_1" , &fJPt1 );
71     fReader->AddVariable( "jseta_1" , &fJEta1 );
72     fReader->AddVariable( "jsphi_1" , &fJPhi1 );
73     fReader->AddVariable( "jd0_1" , &fJD01 );
74     fReader->AddVariable( "jdZ_1" , &fJDZ1 );
75     fReader->AddVariable( "jm_1" , &fJM1 );
76     fReader->AddVariable( "npart_1" , &fNPart1 );
77     fReader->AddVariable( "lpt_1" , &fLPt1 );
78     fReader->AddVariable( "leta_1" , &fLEta1 );
79     fReader->AddVariable( "lphi_1" , &fLPhi1 );
80     fReader->AddVariable( "spt_1" , &fSPt1 );
81     fReader->AddVariable( "seta_1" , &fSEta1 );
82     fReader->AddVariable( "sphi_1" , &fSPhi1 );
83     fReader->AddVariable( "lnept_1" , &fNEPt1 );
84     fReader->AddVariable( "lneeta_1" , &fNEEta1 );
85     fReader->AddVariable( "lnephi_1" , &fNEPhi1 );
86     fReader->AddVariable( "lempt_1" , &fEMPt1 );
87     fReader->AddVariable( "lemeta_1" , &fEMEta1 );
88     fReader->AddVariable( "lemphi_1" , &fEMPhi1 );
89     fReader->AddVariable( "lchpt_1" , &fChPt1 );
90     fReader->AddVariable( "lchphi_1" , &fChPhi1 );
91     fReader->AddVariable( "lLfr_1" , &fLFr1 );
92     fReader->AddVariable( "drlc_1" , &fDRLC1 );
93     fReader->AddVariable( "drls_1" , &fDRLS1 );
94     fReader->AddVariable( "drm_1" , &fDRM1 );
95     fReader->AddVariable( "drmne_1" , &fDRNE1 );
96     fReader->AddVariable( "drem_1" , &fDREM1 );
97     fReader->AddVariable( "drch_1" , &fDRCH1 );
98     }
99     fReader->BookMVA(fMethodName , iWeights );
100     std::cout << "Jet ID MVA Initialization\n";
101     std::cout << "MethodName : " << fMethodName << " , type == " << fType << std::endl;
102     }
103    
104     //--------------------------------------------------------------------------------------------------
105     Double_t JetIDMVA::MVAValue(
106     Float_t iNPV ,
107     Float_t iJPt1 ,
108     Float_t iJEta1 ,
109     Float_t iJPhi1 ,
110     Float_t iJD01 ,
111     Float_t iJDZ1 ,
112     Float_t iJM1 ,
113     Float_t iNPart1 ,
114     Float_t iLPt1 ,
115     Float_t iLEta1 ,
116     Float_t iLPhi1 ,
117     Float_t iSPt1 ,
118     Float_t iSEta1 ,
119     Float_t iSPhi1 ,
120     Float_t iNEPt1 ,
121     Float_t iNEEta1 ,
122     Float_t iNEPhi1 ,
123     Float_t iEMPt1 ,
124     Float_t iEMEta1 ,
125     Float_t iEMPhi1 ,
126     Float_t iChPt1 ,
127     Float_t iChPhi1 ,
128     Float_t iLFr1 ,
129     Float_t iDRLC1 ,
130     Float_t iDRLS1 ,
131     Float_t iDRM1 ,
132     Float_t iDRNE1 ,
133     Float_t iDREM1 ,
134     Float_t iDRCH1
135     ){
136    
137     if (!fIsInitialized) {
138     std::cout << "Error: JetIDMVA not properly initialized.\n";
139     return -9999;
140     }
141    
142     fNPV = iNPV;
143     fJPt1 = iJPt1;
144     fJEta1 = iJEta1;
145     fJPhi1 = fJPhi1;
146     fJD01 = iJD01;
147     fJDZ1 = iJDZ1;
148     fJM1 = iJM1 ;
149     fNPart1 = iNPart1;
150     fLPt1 = iLPt1;
151     fLEta1 = iLEta1;
152     fLPhi1 = iLPhi1;
153     fSPt1 = iSPt1;
154     fSEta1 = iSEta1;
155     fSPhi1 = iSPhi1;
156     fNEPt1 = iNEPt1;
157     fNEEta1 = iNEEta1;
158     fNEPhi1 = iNEPhi1;
159     fEMPt1 = iEMPt1;
160     fEMEta1 = iEMEta1;
161     fEMPhi1 = iEMPhi1;
162     fChPt1 = iChPt1;
163     fChPhi1 = iChPhi1;
164     fLFr1 = iLFr1;
165     fDRLC1 = iDRLC1;
166     fDRLS1 = iDRLS1;
167     fDRM1 = iDRM1;
168     fDRNE1 = iDRNE1;
169     fDREM1 = iDREM1;
170     fDRCH1 = iDRCH1;
171    
172     Double_t lMVA = -9999;
173     lMVA = fReader->EvaluateMVA( fMethodName );
174    
175     return lMVA;
176     }
177     //--------------------------------------------------------------------------------------------------
178     Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,
179     FactorizedJetCorrector *iJetCorrector,
180     const PileupEnergyDensityCol *iPileupEnergyDensity) {
181     if(!JetTools::passPFLooseId(iJet)) return false;
182     if(correctedPt(iJet,iJetCorrector,iPileupEnergyDensity) < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false; //This line is a bug in the Met training
183     if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
184     double lMVA = MVAValue(iJet,iVertex,iJetCorrector,iPileupEnergyDensity);
185     if(lMVA < -0.8) return false;
186     if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
187     return true;
188     }
189     //--------------------------------------------------------------------------------------------------
190     Double_t JetIDMVA::MVAValue(const PFJet *iJet,const Vertex *iVertex, //Vertex here is the PV
191     FactorizedJetCorrector *iJetCorrector,
192     const PileupEnergyDensityCol *iPileupEnergyDensity,
193     Bool_t printDebug) {
194    
195     if (!fIsInitialized) {
196     std::cout << "Error: JetIDMVA not properly initialized.\n";
197     return -9999;
198     }
199     if(!JetTools::passPFLooseId(iJet)) return -2.;
200    
201     //set all input variables
202     fJPt1 = correctedPt(iJet,iJetCorrector,iPileupEnergyDensity);
203     fJEta1 = iJet->RawMom().Eta();
204     fJPhi1 = iJet->RawMom().Phi();
205     fJM1 = iJet->Mass();
206    
207     const mithep::PFCandidate *lLead = JetTools::leadCand(iJet,-1);
208     const mithep::PFCandidate *lSecond = JetTools::leadCand(iJet,-1,true);
209     const mithep::PFCandidate *lLeadNeut = JetTools::leadCand(iJet ,5);
210     const mithep::PFCandidate *lLeadEm = JetTools::leadCand(iJet ,4);
211     const mithep::PFCandidate *lLeadCh = JetTools::leadCand(iJet ,1);
212    
213     fJD01 = JetTools::impactParameter(iJet,iVertex);
214     fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
215     fNPart1 = iJet->NPFCands();
216     fLPt1 = lLead ->Pt();
217     fLEta1 = lLead ->Eta();
218     fLPhi1 = lLead ->Phi();
219     fSPt1 = lSecond ->Pt();
220     fSEta1 = lSecond ->Eta();
221     fSPhi1 = lSecond ->Phi();
222     fNEPt1 = lLeadNeut->Pt();
223     fNEEta1 = lLeadNeut->Eta();
224     fNEPhi1 = lLeadNeut->Phi();
225     fEMPt1 = lLeadEm ->Pt();
226     fEMEta1 = lLeadEm ->Eta();
227     fEMPhi1 = lLeadEm ->Phi();
228     fChPt1 = lLeadCh ->Pt();
229     //fChEta1 = lLeadCh ->Eta();
230     fChPhi1 = lLeadCh ->Phi();
231     fLFr1 = lLead->Pt()/iJet->Pt();
232    
233     fDRLC1 = MathUtils::DeltaR(iJet->Mom(),lLead ->Mom());
234     fDRLS1 = MathUtils::DeltaR(iJet->Mom(),lSecond->Mom());
235     fDRM1 = JetTools::dRMean (iJet,-1);
236     fDRNE1 = JetTools::dRMean (iJet, 5);
237     fDREM1 = JetTools::dRMean (iJet, 4);
238     fDRCH1 = JetTools::dRMean (iJet, 1);
239    
240     double lMVA = fReader->EvaluateMVA( fMethodName );
241    
242     if (printDebug == kTRUE) {
243     std::cout << "Debug Jet MVA: "
244     << fNPV << " "
245     << fJPt1 << " "
246     << fJEta1 << " "
247     << fJPhi1 << " "
248     << fJD01 << " "
249     << fJDZ1 << " "
250     << fJM1 << " "
251     << fNPart1 << " "
252     << fLPt1 << " "
253     << fLEta1 << " "
254     << fLPhi1 << " "
255     << fSPt1 << " "
256     << fSEta1 << " "
257     << fSPhi1 << " "
258     << fNEPt1 << " "
259     << fNEEta1 << " "
260     << fNEPhi1 << " "
261     << fEMPt1 << " "
262     << fEMEta1 << " "
263     << fEMPhi1 << " "
264     << fChPt1 << " "
265     << fChPhi1 << " "
266     << fLFr1 << " "
267     << fDRLC1 << " "
268     << fDRLS1 << " "
269     << fDRM1 << " "
270     << fDRNE1 << " "
271     << fDREM1 << " "
272     << fDRCH1 << " "
273     << " === : === "
274     << lMVA << " "
275     << std::endl;
276     }
277    
278     return lMVA;
279     }
280     Double_t JetIDMVA::correctedPt(const PFJet *iJet, FactorizedJetCorrector *iJetCorrector,const PileupEnergyDensityCol *iPUEnergyDensity) {
281     const FourVectorM rawMom = iJet->RawMom();
282     iJetCorrector->setJetEta(rawMom.Eta());
283     iJetCorrector->setJetPt (rawMom.Pt());
284     iJetCorrector->setJetPhi(rawMom.Phi());
285     iJetCorrector->setJetE (rawMom.E());
286     iJetCorrector->setRho (iPUEnergyDensity->At(0)->RhoHighEta());
287     iJetCorrector->setJetA (iJet->JetArea());
288     iJetCorrector->setJetEMF(-99.0);
289     Double_t correction = iJetCorrector->getCorrection();
290     return rawMom.Pt()*correction;
291     }