[ViewVC] Diff of: cvsroot/UserCode/MitPhysics/Utils/src/JetIDMVA.cc

Comparing UserCode/MitPhysics/Utils/src/JetIDMVA.cc (file contents):
Revision 1.6 by pharris, Mon Apr 23 15:38:00 2012 UTC vs.
Revision 1.21 by pharris, Tue Jul 16 22:29:52 2013 UTC

#	Line 15 \| Line 15 \| using namespace mithep;
15
16		//--------------------------------------------------------------------------------------------------
17		JetIDMVA::JetIDMVA() :
18	<	fJetPtMin(0) , //We need to lower this
18	>	fJetPtMin(0.) , //We need to lower this
19		fDZCut (0.2),
20		fLowPtMethodName ("JetIDMVALowPt" ),
21		fHighPtMethodName("JetIDMVAHighPt"),
#	Line 31 \| Line 31 \| JetIDMVA::JetIDMVA() :
31		fNCharged (0),
32		fNNeutrals(0),
33		fDRMean (0),
34	+	fPtD (0),
35		fFrac01 (0),
36		fFrac02 (0),
37		fFrac03 (0),
38		fFrac04 (0),
39	<	fFrac05 (0)
39	>	fFrac05 (0),
40	>	fDR2Mean (0)
41		{
42		fReader = 0;
43		}
44		//--------------------------------------------------------------------------------------------------
45		JetIDMVA::~JetIDMVA() {
46
47	<	fReader = 0;
47	>	delete fReader;
48	>	delete fLowPtReader;
49		}
50
51		//--------------------------------------------------------------------------------------------------
#	Line 50 \| Line 53 \| void JetIDMVA::Initialize( JetIDMVA::Cut
53		TString iLowPtWeights,
54		TString iHighPtWeights,
55		JetIDMVA::MVAType iType,
56	<	TString iCutFileName)
56	>	TString iCutFileName,bool i42)
57		{
58
59		fIsInitialized = kTRUE;
60		fType = iType;
61		fCutType = iCutType;
62	+	f42 = i42;
63	+
64	+	std::string lCutId = "JetIdParams";
65	+	if(fType == k42) lCutId = "PuJetIdOptMVA_wp";
66	+	if(fType == k52) lCutId = "full_5x_wp";
67	+	if(fType == kCut) lCutId = "PuJetIdCutBased_wp";
68	+	if(fType == k53) lCutId = "full_53x_wp";
69	+	if(fType == k53MET) lCutId = "met_53x_wp";
70	+
71	+	//Load Cut Matrix
72	+	edm::ParameterSet lDefConfig = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>("JetIdParams");
73	+	edm::ParameterSet lConfig = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>(lCutId);
74	+	std::string lCutType = "Tight";
75	+	if(fCutType == kMedium) lCutType = "Medium";
76	+	if(fCutType == kLoose ) lCutType = "Loose";
77	+	if(fCutType == kMET ) lCutType = "MET";
78	+	if(fType != kCut) {
79	+	std::string lLowPtCut = "MET";
80	+	std::vector<double> lPt010 = lDefConfig.getParameter<std::vector<double> >(("Pt010_" +lLowPtCut).c_str());
81	+	std::vector<double> lPt1020 = lConfig.getParameter<std::vector<double> >(("Pt1020_"+lCutType).c_str());
82	+	std::vector<double> lPt2030 = lConfig.getParameter<std::vector<double> >(("Pt2030_"+lCutType).c_str());
83	+	std::vector<double> lPt3050 = lConfig.getParameter<std::vector<double> >(("Pt3050_"+lCutType).c_str());
84	+	for(int i0 = 0; i0 < 4; i0++) fMVACut[0][i0] = lPt010 [i0];
85	+	for(int i0 = 0; i0 < 4; i0++) fMVACut[1][i0] = lPt1020[i0];
86	+	for(int i0 = 0; i0 < 4; i0++) fMVACut[2][i0] = lPt2030[i0];
87	+	for(int i0 = 0; i0 < 4; i0++) fMVACut[3][i0] = lPt3050[i0];
88	+	}
89	+	if(fType == kCut) {
90	+	for(int i0 = 0; i0 < 2; i0++) {
91	+	std::string lFullCutType = lCutType;
92	+	if(i0 == 0) lFullCutType = "BetaStar"+ lCutType;
93	+	if(i0 == 1) lFullCutType = "RMS" + lCutType;
94	+	std::vector<double> pt010 = lConfig.getParameter<std::vector<double> >(("Pt010_" +lFullCutType).c_str());
95	+	std::vector<double> pt1020 = lConfig.getParameter<std::vector<double> >(("Pt1020_"+lFullCutType).c_str());
96	+	std::vector<double> pt2030 = lConfig.getParameter<std::vector<double> >(("Pt2030_"+lFullCutType).c_str());
97	+	std::vector<double> pt3050 = lConfig.getParameter<std::vector<double> >(("Pt3050_"+lFullCutType).c_str());
98	+	if(i0 == 0) {
99	+	for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[0][i2] = pt010 [i2];
100	+	for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[1][i2] = pt1020[i2];
101	+	for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[2][i2] = pt2030[i2];
102	+	for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[3][i2] = pt3050[i2];
103	+	}
104	+	if(i0 == 1) {
105	+	for(int i2 = 0; i2 < 4; i2++) fRMSCut[0][i2] = pt010 [i2];
106	+	for(int i2 = 0; i2 < 4; i2++) fRMSCut[1][i2] = pt1020[i2];
107	+	for(int i2 = 0; i2 < 4; i2++) fRMSCut[2][i2] = pt2030[i2];
108	+	for(int i2 = 0; i2 < 4; i2++) fRMSCut[3][i2] = pt3050[i2];
109	+	}
110	+	}
111	+	return;
112	+	}
113	+
114	+
115	+	fLowPtReader = 0;
116	+	fLowPtReader = new TMVA::Reader( "!Color:!Silent:Error" );
117	+	fLowPtReader->AddVariable( "nvtx" , &fNVtx );
118	+	if(fType != k53) fLowPtReader->AddVariable( "jetPt" , &fJPt1 );
119	+	if(fType != k53) fLowPtReader->AddVariable( "jetEta" , &fJEta1 );
120	+	if(fType != k53) fLowPtReader->AddVariable( "jetPhi" , &fJPhi1 );
121	+	fLowPtReader->AddVariable( "dZ" , &fJDZ1 );
122	+	if(fType != k53 && fType != k53MET) fLowPtReader->AddVariable( "d0" , &fJD01 );
123	+	fLowPtReader->AddVariable( "beta" , &fBeta );
124	+	fLowPtReader->AddVariable( "betaStar" , &fBetaStar );
125	+	fLowPtReader->AddVariable( "nCharged" , &fNCharged );
126	+	fLowPtReader->AddVariable( "nNeutrals", &fNNeutrals );
127	+	if(fType != k53 && fType != k53MET) fLowPtReader->AddVariable( "dRMean" , &fDRMean );
128	+	if(fType == k53 \|\| fType == k53MET) fLowPtReader->AddVariable( "dR2Mean" , &fDR2Mean );
129	+	if(fType == k53 \|\| fType == k53MET) fLowPtReader->AddVariable( "ptD" , &fPtD );
130	+	fLowPtReader->AddVariable( "frac01" , &fFrac01 );
131	+	fLowPtReader->AddVariable( "frac02" , &fFrac02 );
132	+	fLowPtReader->AddVariable( "frac03" , &fFrac03 );
133	+	fLowPtReader->AddVariable( "frac04" , &fFrac04 );
134	+	fLowPtReader->AddVariable( "frac05" , &fFrac05 );
135	+	if(fType == k53) fLowPtReader->AddSpectator( "jetPt" , &fJPt1 );
136	+	if(fType == k53) fLowPtReader->AddSpectator( "jetEta" , &fJEta1 );
137	+	if(fType == k53) fLowPtReader->AddSpectator( "jetPhi" , &fJPhi1 );
138	+
139		fReader = 0;
140		fReader = new TMVA::Reader( "!Color:!Silent:Error" );
141		if (fType == kBaseline) {
#	Line 76 \| Line 156 \| void JetIDMVA::Initialize( JetIDMVA::Cut
156		fReader->AddVariable( "frac04" , &fFrac04 );
157		fReader->AddVariable( "frac05" , &fFrac05 );
158		}
159	<	fReader->BookMVA(fLowPtMethodName , iLowPtWeights );
159	>	if (fType == k42) {
160	>	fReader->AddVariable( "frac01" , &fFrac01 );
161	>	fReader->AddVariable( "frac02" , &fFrac02 );
162	>	fReader->AddVariable( "frac03" , &fFrac03 );
163	>	fReader->AddVariable( "frac04" , &fFrac04 );
164	>	fReader->AddVariable( "frac05" , &fFrac05 );
165	>	fReader->AddVariable( "nvtx" , &fNVtx );
166	>	fReader->AddVariable( "nNeutrals", &fNNeutrals );
167	>	fReader->AddVariable( "beta" , &fBeta );
168	>	fReader->AddVariable( "betaStar" , &fBetaStar );
169	>	fReader->AddVariable( "dZ" , &fJDZ1 );
170	>	fReader->AddVariable( "nCharged" , &fNCharged );
171	>	fReader->AddSpectator( "jetPt" , &fJPt1 );
172	>	fReader->AddSpectator( "jetEta" , &fJEta1 );
173	>	}
174	>	if (fType == k52) {
175	>	fReader->AddVariable( "frac01" , &fFrac01 );
176	>	fReader->AddVariable( "frac02" , &fFrac02 );
177	>	fReader->AddVariable( "frac03" , &fFrac03 );
178	>	fReader->AddVariable( "frac04" , &fFrac04 );
179	>	fReader->AddVariable( "frac05" , &fFrac05 );
180	>	fReader->AddVariable( "dR2Mean" , &fDR2Mean );
181	>	fReader->AddVariable( "nvtx" , &fNVtx );
182	>	fReader->AddVariable( "nNeutrals", &fNNeutrals );
183	>	fReader->AddVariable( "beta" , &fBeta );
184	>	fReader->AddVariable( "betaStar" , &fBetaStar );
185	>	fReader->AddVariable( "dZ" , &fJDZ1 );
186	>	fReader->AddVariable( "nCharged" , &fNCharged );
187	>	fReader->AddSpectator( "jetPt" , &fJPt1 );
188	>	fReader->AddSpectator( "jetEta" , &fJEta1 );
189	>	}
190	>	if (fType == k53) {
191	>	fReader->AddVariable( "nvtx" , &fNVtx );
192	>	fReader->AddVariable( "dZ" , &fJDZ1 );
193	>	fReader->AddVariable( "beta" , &fBeta );
194	>	fReader->AddVariable( "betaStar" , &fBetaStar );
195	>	fReader->AddVariable( "nCharged" , &fNCharged );
196	>	fReader->AddVariable( "nNeutrals", &fNNeutrals );
197	>	fReader->AddVariable( "dR2Mean" , &fDR2Mean );
198	>	fReader->AddVariable( "ptD" , &fPtD );
199	>	fReader->AddVariable( "frac01" , &fFrac01 );
200	>	fReader->AddVariable( "frac02" , &fFrac02 );
201	>	fReader->AddVariable( "frac03" , &fFrac03 );
202	>	fReader->AddVariable( "frac04" , &fFrac04 );
203	>	fReader->AddVariable( "frac05" , &fFrac05 );
204	>	fReader->AddSpectator( "jetPt" , &fJPt1 );
205	>	fReader->AddSpectator( "jetEta" , &fJEta1 );
206	>	fReader->AddSpectator( "jetPhi" , &fJPhi1 );
207	>	}
208	>	if (fType == k53MET) {
209	>	fReader->AddVariable( "nvtx" , &fNVtx );
210	>	fReader->AddVariable( "jetPt" , &fJPt1 );
211	>	fReader->AddVariable( "jetEta" , &fJEta1 );
212	>	fReader->AddVariable( "jetPhi" , &fJPhi1 );
213	>	fReader->AddVariable( "dZ" , &fJDZ1 );
214	>	fReader->AddVariable( "beta" , &fBeta );
215	>	fReader->AddVariable( "betaStar" , &fBetaStar );
216	>	fReader->AddVariable( "nCharged" , &fNCharged );
217	>	fReader->AddVariable( "nNeutrals", &fNNeutrals );
218	>	fReader->AddVariable( "dR2Mean" , &fDR2Mean );
219	>	fReader->AddVariable( "ptD" , &fPtD );
220	>	fReader->AddVariable( "frac01" , &fFrac01 );
221	>	fReader->AddVariable( "frac02" , &fFrac02 );
222	>	fReader->AddVariable( "frac03" , &fFrac03 );
223	>	fReader->AddVariable( "frac04" , &fFrac04 );
224	>	fReader->AddVariable( "frac05" , &fFrac05 );
225	>	}
226	>	if (fType == kQGP) {
227	>	fReader->AddVariable( "nvtx" , &fNVtx );
228	>	fReader->AddVariable( "jetPt" , &fJPt1 );
229	>	fReader->AddVariable( "jetEta" , &fJEta1 );
230	>	fReader->AddVariable( "jetPhi" , &fJPhi1 );
231	>	fReader->AddVariable( "beta" , &fBeta );
232	>	fReader->AddVariable( "betaStar" , &fBetaStar );
233	>	fReader->AddVariable( "nParticles", &fNNeutrals );
234	>	fReader->AddVariable( "nCharged" , &fNCharged );
235	>	fReader->AddVariable( "dRMean" , &fDRMean );
236	>	fReader->AddVariable( "ptD" , &fPtD );
237	>	fReader->AddVariable( "frac01" , &fFrac01 );
238	>	fReader->AddVariable( "frac02" , &fFrac02 );
239	>	fReader->AddVariable( "frac03" , &fFrac03 );
240	>	fReader->AddVariable( "frac04" , &fFrac04 );
241	>	fReader->AddVariable( "frac05" , &fFrac05 );
242	>	}
243	>
244	>	fLowPtReader->BookMVA(fLowPtMethodName , iLowPtWeights );
245		fReader->BookMVA(fHighPtMethodName , iHighPtWeights );
246		std::cout << "Jet ID MVA Initialization\n";
247		std::cout << "MethodName : " << fLowPtMethodName << " , type == " << fType << std::endl;
248
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	–	if(fCutType == kMET ) lCutType = "MET";
90	–	std::vector<double> lPt010 = lConfig.getParameter<std::vector<double> >(("Pt010_" +lCutType).c_str());
91	–	std::vector<double> lPt1020 = lConfig.getParameter<std::vector<double> >(("Pt1020_"+lCutType).c_str());
92	–	std::vector<double> lPt2030 = lConfig.getParameter<std::vector<double> >(("Pt2030_"+lCutType).c_str());
93	–	std::vector<double> lPt3050 = lConfig.getParameter<std::vector<double> >(("Pt3050_"+lCutType).c_str());
94	–	for(int i0 = 0; i0 < 4; i0++) fMVACut[0][i0] = lPt010 [i0];
95	–	for(int i0 = 0; i0 < 4; i0++) fMVACut[1][i0] = lPt1020[i0];
96	–	for(int i0 = 0; i0 < 4; i0++) fMVACut[2][i0] = lPt2030[i0];
97	–	for(int i0 = 0; i0 < 4; i0++) fMVACut[3][i0] = lPt3050[i0];
98	–	//std::cout << " Working Points : << " << std::endl;
99	–	//for(int i0 = 0; i0 < 4; i0++) for(int i1 = 0; i1 < 4; i1++)
100	–	// std::cout << " ==> " << i0 << " -- " << i1 << " -- " << fMVACut[i0][i1] << std::endl;
249		}
250
251		//--------------------------------------------------------------------------------------------------
#	Line 117 \| Line 265 \| Double_t JetIDMVA::MVAValue(
265		Float_t iFrac02 ,
266		Float_t iFrac03 ,
267		Float_t iFrac04 ,
268	<	Float_t iFrac05
268	>	Float_t iFrac05 ,
269	>	Float_t iDR2Mean ,
270	>	Float_t iPtD
271		){
272
273		if(!fIsInitialized) {
#	Line 136 \| Line 286 \| Double_t JetIDMVA::MVAValue(
286		fNCharged = iNCharged;
287		fNNeutrals = iNNeutrals;
288		fDRMean = iDRMean;
289	+	fPtD = iPtD;
290		fFrac01 = iFrac01;
291		fFrac02 = iFrac02;
292		fFrac03 = iFrac03;
293		fFrac04 = iFrac04;
294		fFrac05 = iFrac05;
295	+	fDR2Mean = iDR2Mean;
296
297		Double_t lMVA = -9999;
298	<	if(iJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName );
298	>	if(iJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName );
299		if(iJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
148	–
300		return lMVA;
301		}
302		//--------------------------------------------------------------------------------------------------
303	+	Bool_t JetIDMVA::passPt(const PFJet iJet, FactorizedJetCorrector iJetCorrector,
304	+	const PileupEnergyDensityCol *iPileupEnergyDensity,
305	+	RhoUtilities::RhoType type) {
306	+	if(iJetCorrector == 0) return (iJet->Pt() > fJetPtMin);
307	+	double lPt = correctedPt(iJet, iJetCorrector,iPileupEnergyDensity,type);
308	+	if(lPt < fJetPtMin) return false;
309	+	return true;
310	+	}
311	+	//--------------------------------------------------------------------------------------------------
312		Bool_t JetIDMVA::pass(const PFJet iJet,const Vertex iVertex,const VertexCol *iVertices,
313		FactorizedJetCorrector *iJetCorrector,
314	<	const PileupEnergyDensityCol *iPileupEnergyDensity) {
314	>	const PileupEnergyDensityCol *iPileupEnergyDensity,
315	>	RhoUtilities::RhoType type) {
316
317		if(!JetTools::passPFLooseId(iJet)) return false;
318	<	if(iJet->Pt() < fJetPtMin) return false;
158	<	if(fabs(iJet->Eta()) > 4.99) return true; //==> Castor
159	<	//if(iJet->Pt() > 50) return true; //==> we can raise this
318	>	if(fabs(iJet->Eta()) > 4.99) return false;
319
320		double lMVA = MVAValue (iJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity);
321	<	double lPt = correctedPt(iJet, iJetCorrector,iPileupEnergyDensity);
322	<
321	>	double lPt = correctedPt(iJet, iJetCorrector,iPileupEnergyDensity,type);
322	>	if(lPt < fJetPtMin) return false;
323	>
324		int lPtId = 0;
325	<	if(lPt > 10 && iJet->Pt() < 20) lPtId = 1;
326	<	if(lPt > 20 && iJet->Pt() < 30) lPtId = 2;
325	>	if(lPt > 10 && lPt < 20) lPtId = 1;
326	>	if(lPt > 20 && lPt < 30) lPtId = 2;
327		if(lPt > 30 ) lPtId = 3;
328
329		int lEtaId = 0;
#	Line 174 \| Line 334 \| *Bool_t JetIDMVA::pass(const PFJet iJet,**
334		double lMVACut = fMVACut[lPtId][lEtaId];
335		if(lMVA < lMVACut) return false;
336		return true;
337	<	//if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
338	<	//if(correctedPt(iJet,iJetCorrector,iPileupEnergyDensity) < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false;
339	<	//This line is a bug in the Met training
340	<	//if(lMVA < -0.8) return false;
341	<	//if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
337	>	}
338	>	//--------------------------------------------------------------------------------------------------
339	>	Bool_t JetIDMVA::passCut(const PFJet iJet,const Vertex iVertex,const VertexCol *iVertices) {
340	>	if(!JetTools::passPFLooseId(iJet)) return false;
341	>	if(iJet->Pt() < fJetPtMin) return false;
342	>	if(fabs(iJet->Eta()) > 4.99) return false;
343	>	//if(fType == kCut) passCut(iJet,iVertex,iVertices);
344	>
345	>	double lPt = iJet->Pt();
346	>	int lPtId = 0;
347	>	if(lPt > 10 && lPt < 20) lPtId = 1;
348	>	if(lPt > 20 && lPt < 30) lPtId = 2;
349	>	if(lPt > 30 ) lPtId = 3;
350	>
351	>	int lEtaId = 0;
352	>	if(fabs(iJet->Eta()) > 2.5 && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
353	>	if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
354	>	if(fabs(iJet->Eta()) > 3.0 && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
355	>	float betaStarModified = JetTools::betaStarClassic(iJet,iVertex,iVertices)/log(iVertices ->GetEntries()-0.64);
356	>	float dR2Mean = JetTools::dR2Mean(iJet,-1);
357	>
358	>	if(betaStarModified < fBetaStarCut[lPtId][lEtaId] &&
359	>	dR2Mean < fRMSCut [lPtId][lEtaId]
360	>	) return true;
361	>
362	>	return false;
363		}
364		//--------------------------------------------------------------------------------------------------
365		Bool_t JetIDMVA::pass(const PFJet iJet,const Vertex iVertex,const VertexCol *iVertices) {
366		if(!JetTools::passPFLooseId(iJet)) return false;
367		if(iJet->Pt() < fJetPtMin) return false;
368	<	//if(iJet->Pt() > 50) return true; //==> we can raise this
369	<	if(fabs(iJet->Eta()) > 4.99) return true; //==> Castor
368	>	if(fabs(iJet->Eta()) > 4.99) return false;
369	>	if(fType == kCut) return passCut(iJet,iVertex,iVertices);
370		double lMVA = MVAValue(iJet,iVertex,iVertices);
371
372		int lPtId = 0;
#	Line 231 \| Line 412 \| Double_t JetIDMVA::MVAValue(const PFJet
412		fNNeutrals = iJet->NeutralMultiplicity();
413
414		fDRMean = JetTools::dRMean(iJet,-1);
415	+	fDR2Mean = JetTools::dR2Mean(iJet,-1);
416	+	fPtD = JetTools::W(iJet,-1,0);
417		fFrac01 = JetTools::frac (iJet,0.1,0. ,-1);
418		fFrac02 = JetTools::frac (iJet,0.2,0.1,-1);
419		fFrac03 = JetTools::frac (iJet,0.3,0.2,-1);
#	Line 238 \| Line 421 \| Double_t JetIDMVA::MVAValue(const PFJet
421		fFrac05 = JetTools::frac (iJet,0.5,0.4,-1);
422
423		double lMVA = 0;
424	<	if(fJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName );
424	>	if(fJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName );
425		if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
426		if (printDebug == kTRUE) {
427		std::cout << "Debug Jet MVA: "
#	Line 253 \| Line 436 \| Double_t JetIDMVA::MVAValue(const PFJet
436		<< fNCharged << " "
437		<< fNNeutrals << " "
438		<< fDRMean << " "
439	+	<< fPtD << " "
440		<< fFrac01 << " "
441		<< fFrac02 << " "
442		<< fFrac03 << " "
443		<< fFrac04 << " "
444	<	<< fFrac05
444	>	<< fFrac05 << " "
445	>	<< fDR2Mean
446		<< " === : === "
447		<< lMVA << " "
448	<	<< std::endl;
448	>	<< " -----> raw pt " << iJet->Pt() << std::endl;
449		}
450	<
450	>	//std::cout << " === " << iJet->Pt() << " -- " << iJet->Eta() << " -- " << fPtD << " -- " << lMVA << std::endl;
451		return lMVA;
452		}
453		Double_t JetIDMVA::MVAValue(const PFJet iJet,const Vertex iVertex, const VertexCol *iVertices,//Vertex here is the PV
#	Line 285 \| Line 470 \| Double_t JetIDMVA::MVAValue(const PFJet
470		fBetaStar = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
471		fNCharged = iJet->ChargedMultiplicity();
472		fNNeutrals = iJet->NeutralMultiplicity();
473	+	fPtD = JetTools::W(iJet,-1,0);
474	+	fDRMean = JetTools::dRMean (iJet,-1);
475	+	fDR2Mean = JetTools::dR2Mean(iJet,-1);
476	+	fFrac01 = JetTools::frac (iJet,0.1,0. ,-1);
477	+	fFrac02 = JetTools::frac (iJet,0.2,0.1,-1);
478	+	fFrac03 = JetTools::frac (iJet,0.3,0.2,-1);
479	+	fFrac04 = JetTools::frac (iJet,0.4,0.3,-1);
480	+	fFrac05 = JetTools::frac (iJet,0.5,0.4,-1);
481	+
482	+	double lMVA = 0;
483	+	if(fJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName );
484	+	if(fJPt1 > 10) lMVA = fReader ->EvaluateMVA( fHighPtMethodName );
485	+
486	+	if (printDebug == kTRUE) {
487	+	std::cout << "Debug Jet MVA: "
488	+	<< fNVtx << " "
489	+	<< fJPt1 << " "
490	+	<< fJEta1 << " "
491	+	<< fJPhi1 << " "
492	+	<< fJD01 << " "
493	+	<< fJDZ1 << " "
494	+	<< fBeta << " "
495	+	<< fBetaStar << " "
496	+	<< fNCharged << " "
497	+	<< fNNeutrals << " "
498	+	<< fDRMean << " "
499	+	<< fPtD << " "
500	+	<< fFrac01 << " "
501	+	<< fFrac02 << " "
502	+	<< fFrac03 << " "
503	+	<< fFrac04 << " "
504	+	<< fFrac05 << " "
505	+	<< fDR2Mean
506	+	<< " === : === "
507	+	<< lMVA << " "
508	+	<< std::endl;
509	+	}
510	+
511	+	return lMVA;
512	+	}
513	+	//--------------------------------------------------------------------------------------------------
514	+	Double_t* JetIDMVA::QGValue(const PFJet iJet,const Vertex iVertex,const VertexCol *iVertices, //Vertex here is the PV
515	+	FactorizedJetCorrector *iJetCorrector,
516	+	const PileupEnergyDensityCol *iPileupEnergyDensity,
517	+	Bool_t printDebug) {
518	+
519	+	Double_t *lId = new double[3];
520	+	lId[0] = -2;
521	+	lId[1] = -2;
522	+	lId[2] = -2;
523	+	if (!fIsInitialized) {
524	+	std::cout << "Error: JetIDMVA not properly initialized.\n";
525	+	return lId;
526	+	}
527	+	if(!JetTools::passPFLooseId(iJet)) return lId;
528	+
529	+	fJPt1 = correctedPt(iJet,iJetCorrector,iPileupEnergyDensity);
530	+	if(fJPt1 < 20) return lId;
531	+
532	+	//set all input variables
533	+	fNVtx = iVertices->GetEntries();
534	+	fJEta1 = iJet->RawMom().Eta();
535	+	fJPhi1 = iJet->RawMom().Phi();
536	+	fJD01 = JetTools::impactParameter(iJet,iVertex);
537	+	fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
538	+	fBeta = JetTools::Beta(iJet,iVertex,fDZCut);
539	+	fBetaStar = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
540	+	fNCharged = iJet->ChargedMultiplicity();
541	+	fNNeutrals = iJet->NeutralMultiplicity();
542	+	fNParticles = fNCharged+fNNeutrals;
543	+	fPtD = JetTools::W(iJet,-1,0);
544
545		fDRMean = JetTools::dRMean(iJet,-1);
546	+	fDR2Mean = JetTools::dR2Mean(iJet,-1);
547		fFrac01 = JetTools::frac (iJet,0.1,0. ,-1);
548		fFrac02 = JetTools::frac (iJet,0.2,0.1,-1);
549		fFrac03 = JetTools::frac (iJet,0.3,0.2,-1);
#	Line 294 \| Line 551 \| Double_t JetIDMVA::MVAValue(const PFJet
551		fFrac05 = JetTools::frac (iJet,0.5,0.4,-1);
552
553		double lMVA = 0;
554	<	if(fJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName );
555	<	if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
556	<
554	>	lId[0] = fReader->EvaluateMulticlass( fHighPtMethodName )[0];
555	>	lId[1] = fReader->EvaluateMulticlass( fHighPtMethodName )[1];
556	>	lId[2] = fReader->EvaluateMulticlass( fHighPtMethodName )[2];
557		if (printDebug == kTRUE) {
558		std::cout << "Debug Jet MVA: "
559		<< fNVtx << " "
#	Line 314 \| Line 571 \| Double_t JetIDMVA::MVAValue(const PFJet
571		<< fFrac02 << " "
572		<< fFrac03 << " "
573		<< fFrac04 << " "
574	<	<< fFrac05
574	>	<< fFrac05 << " "
575	>	<< fDRMean
576		<< " === : === "
577		<< lMVA << " "
578		<< std::endl;
579		}
580	<
323	<	return lMVA;
580	>	return lId;
581		}
582	<	Double_t JetIDMVA::correctedPt(const PFJet iJet, FactorizedJetCorrector iJetCorrector,const PileupEnergyDensityCol *iPUEnergyDensity) {
582	>	Double_t JetIDMVA::correctedPt(const PFJet iJet, FactorizedJetCorrector iJetCorrector,
583	>	const PileupEnergyDensityCol *iPUEnergyDensity,
584	>	RhoUtilities::RhoType type,int iId) {
585	>	Double_t Rho = 0.0;
586	>	switch(type) {
587	>	case RhoUtilities::MIT_RHO_VORONOI_HIGH_ETA:
588	>	Rho = iPUEnergyDensity->At(0)->Rho();
589	>	break;
590	>	case RhoUtilities::MIT_RHO_VORONOI_LOW_ETA:
591	>	Rho = iPUEnergyDensity->At(0)->RhoLowEta();
592	>	break;
593	>	case RhoUtilities::MIT_RHO_RANDOM_HIGH_ETA:
594	>	Rho = iPUEnergyDensity->At(0)->RhoRandom();
595	>	break;
596	>	case RhoUtilities::MIT_RHO_RANDOM_LOW_ETA:
597	>	Rho = iPUEnergyDensity->At(0)->RhoRandomLowEta();
598	>	break;
599	>	case RhoUtilities::CMS_RHO_RHOKT6PFJETS:
600	>	Rho = iPUEnergyDensity->At(0)->RhoKt6PFJets();
601	>	break;
602	>	default:
603	>	// use the old default
604	>	f42 ? Rho = iPUEnergyDensity->At(0)->RhoRandom() : Rho = iPUEnergyDensity->At(0)->Rho();
605	>	break;
606	>	}
607	>
608		const FourVectorM rawMom = iJet->RawMom();
609		iJetCorrector->setJetEta(rawMom.Eta());
610		iJetCorrector->setJetPt (rawMom.Pt());
611		iJetCorrector->setJetPhi(rawMom.Phi());
612		iJetCorrector->setJetE (rawMom.E());
613	<	iJetCorrector->setRho (iPUEnergyDensity->At(0)->RhoHighEta());
613	>	iJetCorrector->setRho (Rho);
614		iJetCorrector->setJetA (iJet->JetArea());
615		iJetCorrector->setJetEMF(-99.0);
616	<	Double_t correction = iJetCorrector->getCorrection();
616	>	Double_t correction = 1.;
617	>	if(iId < 0 \|\| iId == 100) correction = iJetCorrector->getCorrection();
618	>	std::vector<float> lCorrections; if(iId != -1 && iId != 100) lCorrections = iJetCorrector->getSubCorrections();
619	>	if(iId > -1 && iId < int(lCorrections.size())) correction = lCorrections[iId];
620	>	if(iId == 100) {
621	>	iJetCorrector->setJetEta(rawMom.Eta());
622	>	iJetCorrector->setJetPt (rawMom.Pt());
623	>	iJetCorrector->setJetPhi(rawMom.Phi());
624	>	iJetCorrector->setJetE (rawMom.E());
625	>	iJetCorrector->setRho (Rho);
626	>	iJetCorrector->setJetA (iJet->JetArea());
627	>	iJetCorrector->setJetEMF(-99.0);
628	>	lCorrections = iJetCorrector->getSubCorrections();
629	>	double correction2 = 1;
630	>	correction2 *= lCorrections[0];
631	>	correction = correction-correction2;
632	>	}
633	>
634		return rawMom.Pt()*correction;
635		}

Diff Legend

-–
+Removed lines
-+
+Added lines
-<
+Changed lines
->
+Changed lines

Comparing UserCode/MitPhysics/Utils/src/JetIDMVA.cc (file contents): Revision 1.6 by pharris, Mon Apr 23 15:38:00 2012 UTC vs. Revision 1.21 by pharris, Tue Jul 16 22:29:52 2013 UTC

Diff Legend

Comparing UserCode/MitPhysics/Utils/src/JetIDMVA.cc (file contents):
Revision 1.6 by pharris, Mon Apr 23 15:38:00 2012 UTC vs.
Revision 1.21 by pharris, Tue Jul 16 22:29:52 2013 UTC