ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/JetIDMVA.cc
Revision: 1.4
Committed: Sat Apr 7 09:36:33 2012 UTC (13 years ago) by pharris
Content type: text/plain
Branch: MAIN
Changes since 1.3: +211 -266 lines
Log Message:
Update of Met regresion

File Contents

# User Rev Content
1 pharris 1.4 #include "FWCore/PythonParameterSet/interface/MakeParameterSets.h"
2     #include "FWCore/ParameterSet/interface/ParameterSet.h"
3 pharris 1.1 #include "MitPhysics/Utils/interface/JetIDMVA.h"
4     #include "MitPhysics/Utils/interface/JetTools.h"
5     #include "MitAna/DataTree/interface/StableData.h"
6     #include <TFile.h>
7     #include <TRandom3.h>
8     #include "TMVA/Tools.h"
9     #include "TMVA/Reader.h"
10    
11    
12     ClassImp(mithep::JetIDMVA)
13    
14     using namespace mithep;
15    
16     //--------------------------------------------------------------------------------------------------
17     JetIDMVA::JetIDMVA() :
18 pharris 1.4 fJetPtMin(0) , //We need to lower this
19     fDZCut (0.2),
20     fLowPtMethodName ("JetIDMVALowPt" ),
21     fHighPtMethodName("JetIDMVAHighPt"),
22 pharris 1.1 fIsInitialized(kFALSE),
23 pharris 1.4 fNVtx (0),
24     fJPt1 (0),
25     fJEta1 (0),
26     fJPhi1 (0),
27     fJD01 (0),
28     fJDZ1 (0),
29     fBeta (0),
30     fBetaStar (0),
31     fNCharged (0),
32     fNNeutrals(0),
33     fDRMean (0),
34     fFrac01 (0),
35     fFrac02 (0),
36     fFrac03 (0),
37     fFrac04 (0),
38     fFrac05 (0)
39 pharris 1.1 {
40     fReader = 0;
41     }
42     //--------------------------------------------------------------------------------------------------
43 pharris 1.4 JetIDMVA::~JetIDMVA() {
44    
45 pharris 1.1 fReader = 0;
46     }
47    
48     //--------------------------------------------------------------------------------------------------
49 pharris 1.4 void JetIDMVA::Initialize( JetIDMVA::CutType iCutType,
50     TString iLowPtWeights,
51     TString iHighPtWeights,
52     JetIDMVA::MVAType iType,
53     TString iCutFileName)
54 pharris 1.1 {
55    
56     fIsInitialized = kTRUE;
57     fType = iType;
58 pharris 1.4 fCutType = iCutType;
59 pharris 1.1 fReader = 0;
60     fReader = new TMVA::Reader( "!Color:!Silent:Error" );
61     if (fType == kBaseline) {
62 pharris 1.4 fReader->AddVariable( "nvtx" , &fNVtx );
63     fReader->AddVariable( "jetPt" , &fJPt1 );
64     fReader->AddVariable( "jetEta" , &fJEta1 );
65     fReader->AddVariable( "jetPhi" , &fJPhi1 );
66     fReader->AddVariable( "dZ" , &fJDZ1 );
67     fReader->AddVariable( "d0" , &fJD01 );
68     fReader->AddVariable( "beta" , &fBeta );
69     fReader->AddVariable( "betaStar" , &fBetaStar );
70     fReader->AddVariable( "nCharged" , &fNCharged );
71     fReader->AddVariable( "nNeutrals", &fNNeutrals );
72     fReader->AddVariable( "dRMean" , &fDRMean );
73     fReader->AddVariable( "frac01" , &fFrac01 );
74     fReader->AddVariable( "frac02" , &fFrac02 );
75     fReader->AddVariable( "frac03" , &fFrac03 );
76     fReader->AddVariable( "frac04" , &fFrac04 );
77     fReader->AddVariable( "frac05" , &fFrac05 );
78 pharris 1.1 }
79 pharris 1.4 fReader->BookMVA(fLowPtMethodName , iLowPtWeights );
80     fReader->BookMVA(fHighPtMethodName , iHighPtWeights );
81 pharris 1.1 std::cout << "Jet ID MVA Initialization\n";
82 pharris 1.4 std::cout << "MethodName : " << fLowPtMethodName << " , type == " << fType << std::endl;
83    
84     //Load Cut Matrix
85     edm::ParameterSet lConfig = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>("JetIdParams");
86     std::string lCutType = "Tight";
87     if(fCutType == kMedium) lCutType = "Medium";
88     if(fCutType == kLoose ) lCutType = "Loose";
89     std::vector<double> lPt010 = lConfig.getParameter<std::vector<double> >(("Pt010_" +lCutType).c_str());
90     std::vector<double> lPt1020 = lConfig.getParameter<std::vector<double> >(("Pt1020_"+lCutType).c_str());
91     std::vector<double> lPt2030 = lConfig.getParameter<std::vector<double> >(("Pt2030_"+lCutType).c_str());
92     std::vector<double> lPt3050 = lConfig.getParameter<std::vector<double> >(("Pt3050_"+lCutType).c_str());
93     for(int i0 = 0; i0 < 4; i0++) fMVACut[0][i0] = lPt010 [i0];
94     for(int i0 = 0; i0 < 4; i0++) fMVACut[1][i0] = lPt1020[i0];
95     for(int i0 = 0; i0 < 4; i0++) fMVACut[2][i0] = lPt2030[i0];
96     for(int i0 = 0; i0 < 4; i0++) fMVACut[3][i0] = lPt3050[i0];
97 pharris 1.1 }
98    
99     //--------------------------------------------------------------------------------------------------
100     Double_t JetIDMVA::MVAValue(
101     Float_t iNPV ,
102     Float_t iJPt1 ,
103     Float_t iJEta1 ,
104     Float_t iJPhi1 ,
105     Float_t iJD01 ,
106     Float_t iJDZ1 ,
107 pharris 1.4 Float_t iBeta ,
108     Float_t iBetaStar,
109     Float_t iNCharged,
110     Float_t iNNeutrals,
111     Float_t iDRMean ,
112     Float_t iFrac01 ,
113     Float_t iFrac02 ,
114     Float_t iFrac03 ,
115     Float_t iFrac04 ,
116     Float_t iFrac05
117 pharris 1.1 ){
118    
119 pharris 1.4 if(!fIsInitialized) {
120 pharris 1.1 std::cout << "Error: JetIDMVA not properly initialized.\n";
121     return -9999;
122     }
123    
124 pharris 1.4 fNVtx = iNPV;
125     fJPt1 = iJPt1;
126     fJEta1 = iJEta1;
127     fJPhi1 = fJPhi1;
128     fJD01 = iJD01;
129     fJDZ1 = iJDZ1;
130     fBeta = iBeta;
131     fBetaStar = iBetaStar;
132     fNCharged = iNCharged;
133     fNNeutrals = iNNeutrals;
134     fDRMean = iDRMean;
135     fFrac01 = iFrac01;
136     fFrac02 = iFrac02;
137     fFrac03 = iFrac03;
138     fFrac04 = iFrac04;
139     fFrac05 = iFrac05;
140 pharris 1.1
141     Double_t lMVA = -9999;
142 pharris 1.4 if(iJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName );
143     if(iJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
144 pharris 1.1
145     return lMVA;
146     }
147     //--------------------------------------------------------------------------------------------------
148 pharris 1.4 Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices,
149 pharris 1.1 FactorizedJetCorrector *iJetCorrector,
150     const PileupEnergyDensityCol *iPileupEnergyDensity) {
151 pharris 1.4
152 pharris 1.1 if(!JetTools::passPFLooseId(iJet)) return false;
153 pharris 1.4 if(iJet->Pt() < fJetPtMin) return false;
154     if(iJet->Pt() > 50) return true; //==> we can raise this
155     if(fabs(iJet->Eta()) > 4.99) return false;//Castor
156    
157     double lMVA = MVAValue (iJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity);
158     double lPt = correctedPt(iJet, iJetCorrector,iPileupEnergyDensity);
159    
160     int lPtId = 0;
161     if(lPt > 10 && iJet->Pt() < 20) lPtId = 1;
162     if(lPt > 20 && iJet->Pt() < 30) lPtId = 2;
163     if(lPt > 30 ) lPtId = 3;
164    
165     int lEtaId = 0;
166     if(fabs(iJet->Eta()) > 2.5 && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
167     if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
168     if(fabs(iJet->Eta()) > 3.0 && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
169    
170     double lMVACut = fMVACut[lPtId][lEtaId];
171     if(lMVA < lMVACut) return false;
172 pharris 1.1 return true;
173 pharris 1.4 //if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
174     //if(correctedPt(iJet,iJetCorrector,iPileupEnergyDensity) < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false;
175     //This line is a bug in the Met training
176     //if(lMVA < -0.8) return false;
177     //if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
178 pharris 1.1 }
179     //--------------------------------------------------------------------------------------------------
180 pharris 1.4 Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices) {
181 pharris 1.2 if(!JetTools::passPFLooseId(iJet)) return false;
182 pharris 1.4 if(iJet->Pt() < fJetPtMin) return false;
183     if(iJet->Pt() > 50) return true; //==> we can raise this
184     if(fabs(iJet->Eta()) > 4.99) return false;//Castor
185    
186     double lMVA = MVAValue(iJet,iVertex,iVertices);
187    
188     int lPtId = 0;
189     if(iJet->Pt() > 10 && iJet->Pt() < 20) lPtId = 1;
190     if(iJet->Pt() > 20 && iJet->Pt() < 30) lPtId = 2;
191     if(iJet->Pt() > 30 ) lPtId = 3;
192    
193     int lEtaId = 0;
194     if(fabs(iJet->Eta()) > 2.5 && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
195     if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
196     if(fabs(iJet->Eta()) > 3.0 && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
197    
198     double lMVACut = fMVACut[lPtId][lEtaId];
199     if(lMVA < lMVACut) return false;
200 pharris 1.2 return true;
201 pharris 1.4 //if(lMVA < -0.8) return false;
202     //if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
203     //if(iJet->Pt() < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false; //This line is a bug in the Met training
204     //if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
205 pharris 1.2 }
206     //--------------------------------------------------------------------------------------------------
207 pharris 1.4 Double_t JetIDMVA::MVAValue(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices, //Vertex here is the PV
208 pharris 1.1 FactorizedJetCorrector *iJetCorrector,
209     const PileupEnergyDensityCol *iPileupEnergyDensity,
210     Bool_t printDebug) {
211    
212     if (!fIsInitialized) {
213     std::cout << "Error: JetIDMVA not properly initialized.\n";
214     return -9999;
215     }
216     if(!JetTools::passPFLooseId(iJet)) return -2.;
217    
218     //set all input variables
219 pharris 1.4 fNVtx = iVertices->GetEntries();
220 pharris 1.1 fJPt1 = correctedPt(iJet,iJetCorrector,iPileupEnergyDensity);
221     fJEta1 = iJet->RawMom().Eta();
222     fJPhi1 = iJet->RawMom().Phi();
223 pharris 1.4 fJD01 = JetTools::impactParameter(iJet,iVertex);
224     fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
225     fBeta = JetTools::Beta(iJet,iVertex,fDZCut);
226     fBetaStar = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
227     fNCharged = iJet->ChargedMultiplicity();
228     fNNeutrals = iJet->NeutralMultiplicity();
229    
230     fDRMean = JetTools::dRMean(iJet,-1);
231     fFrac01 = JetTools::frac (iJet,0.1,0. ,-1);
232     fFrac02 = JetTools::frac (iJet,0.2,0.1,-1);
233     fFrac03 = JetTools::frac (iJet,0.3,0.2,-1);
234     fFrac04 = JetTools::frac (iJet,0.4,0.3,-1);
235     fFrac05 = JetTools::frac (iJet,0.5,0.4,-1);
236    
237     double lMVA = 0;
238     if(fJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName );
239     if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
240 pharris 1.1
241     if (printDebug == kTRUE) {
242     std::cout << "Debug Jet MVA: "
243 pharris 1.4 << fNVtx << " "
244     << fJPt1 << " "
245     << fJEta1 << " "
246     << fJPhi1 << " "
247     << fJD01 << " "
248     << fJDZ1 << " "
249     << fBeta << " "
250     << fBetaStar << " "
251     << fNCharged << " "
252     << fNNeutrals << " "
253     << fDRMean << " "
254     << fFrac01 << " "
255     << fFrac02 << " "
256     << fFrac03 << " "
257     << fFrac04 << " "
258     << fFrac05
259 pharris 1.1 << " === : === "
260     << lMVA << " "
261     << std::endl;
262     }
263    
264     return lMVA;
265     }
266 pharris 1.4 Double_t JetIDMVA::MVAValue(const PFJet *iJet,const Vertex *iVertex, const VertexCol *iVertices,//Vertex here is the PV
267 pharris 1.2 Bool_t printDebug) {
268    
269     if (!fIsInitialized) {
270     std::cout << "Error: JetIDMVA not properly initialized.\n";
271     return -9999;
272     }
273     if(!JetTools::passPFLooseId(iJet)) return -2.;
274    
275     //set all input variables
276 pharris 1.4 fNVtx = iVertices->GetEntries();
277 pharris 1.2 fJPt1 = iJet->Pt();
278     fJEta1 = iJet->RawMom().Eta();
279     fJPhi1 = iJet->RawMom().Phi();
280 pharris 1.4 fJD01 = JetTools::impactParameter(iJet,iVertex);
281     fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
282     fBeta = JetTools::Beta(iJet,iVertex,fDZCut);
283     fBetaStar = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
284     fNCharged = iJet->ChargedMultiplicity();
285     fNNeutrals = iJet->NeutralMultiplicity();
286    
287     fDRMean = JetTools::dRMean(iJet,-1);
288     fFrac01 = JetTools::frac (iJet,0.1,0. ,-1);
289     fFrac02 = JetTools::frac (iJet,0.2,0.1,-1);
290     fFrac03 = JetTools::frac (iJet,0.3,0.2,-1);
291     fFrac04 = JetTools::frac (iJet,0.4,0.3,-1);
292     fFrac05 = JetTools::frac (iJet,0.5,0.4,-1);
293    
294     double lMVA = 0;
295     if(fJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName );
296     if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
297 pharris 1.2
298     if (printDebug == kTRUE) {
299     std::cout << "Debug Jet MVA: "
300 pharris 1.4 << fNVtx << " "
301     << fJPt1 << " "
302     << fJEta1 << " "
303     << fJPhi1 << " "
304     << fJD01 << " "
305     << fJDZ1 << " "
306     << fBeta << " "
307     << fBetaStar << " "
308     << fNCharged << " "
309     << fNNeutrals << " "
310     << fDRMean << " "
311     << fFrac01 << " "
312     << fFrac02 << " "
313     << fFrac03 << " "
314     << fFrac04 << " "
315     << fFrac05
316 pharris 1.2 << " === : === "
317     << lMVA << " "
318     << std::endl;
319     }
320    
321     return lMVA;
322     }
323 pharris 1.1 Double_t JetIDMVA::correctedPt(const PFJet *iJet, FactorizedJetCorrector *iJetCorrector,const PileupEnergyDensityCol *iPUEnergyDensity) {
324     const FourVectorM rawMom = iJet->RawMom();
325     iJetCorrector->setJetEta(rawMom.Eta());
326     iJetCorrector->setJetPt (rawMom.Pt());
327     iJetCorrector->setJetPhi(rawMom.Phi());
328     iJetCorrector->setJetE (rawMom.E());
329     iJetCorrector->setRho (iPUEnergyDensity->At(0)->RhoHighEta());
330     iJetCorrector->setJetA (iJet->JetArea());
331     iJetCorrector->setJetEMF(-99.0);
332     Double_t correction = iJetCorrector->getCorrection();
333     return rawMom.Pt()*correction;
334     }