ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/MitPhysics/Utils/src/JetIDMVA.cc

Comparing UserCode/MitPhysics/Utils/src/JetIDMVA.cc (file contents):
Revision 1.5 by pharris, Fri Apr 13 14:24:00 2012 UTC vs.
Revision 1.12 by pharris, Fri May 25 20:01:40 2012 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 35 | Line 35 | JetIDMVA::JetIDMVA() :
35		fFrac02   (0),
36		fFrac03   (0),
37		fFrac04   (0),
38	<	fFrac05   (0)
38	>	fFrac05   (0),
39	>	fDR2Mean  (0)
40		{
41		fReader = 0;
42		}
43		//--------------------------------------------------------------------------------------------------
44		JetIDMVA::~JetIDMVA() {
45
46	<	fReader = 0;
46	>	delete fReader;
47	>	delete fLowPtReader;
48		}
49
50		//--------------------------------------------------------------------------------------------------
#	Line 56 | Line 58 | void JetIDMVA::Initialize( JetIDMVA::Cut
58		fIsInitialized = kTRUE;
59		fType          = iType;
60		fCutType       = iCutType;
61	+	fLowPtReader   = 0;
62	+	fLowPtReader   = new TMVA::Reader( "!Color:!Silent:Error" );
63	+	fLowPtReader->AddVariable( "nvtx"     , &fNVtx      );
64	+	fLowPtReader->AddVariable( "jetPt"    , &fJPt1      );
65	+	fLowPtReader->AddVariable( "jetEta"   , &fJEta1     );
66	+	fLowPtReader->AddVariable( "jetPhi"   , &fJPhi1     );
67	+	fLowPtReader->AddVariable( "dZ"       , &fJDZ1      );
68	+	fLowPtReader->AddVariable( "d0"       , &fJD01      );
69	+	fLowPtReader->AddVariable( "beta"     , &fBeta      );
70	+	fLowPtReader->AddVariable( "betaStar" , &fBetaStar  );
71	+	fLowPtReader->AddVariable( "nCharged" , &fNCharged  );
72	+	fLowPtReader->AddVariable( "nNeutrals", &fNNeutrals );
73	+	fLowPtReader->AddVariable( "dRMean"   , &fDRMean    );
74	+	fLowPtReader->AddVariable( "frac01"   , &fFrac01    );
75	+	fLowPtReader->AddVariable( "frac02"   , &fFrac02    );
76	+	fLowPtReader->AddVariable( "frac03"   , &fFrac03    );
77	+	fLowPtReader->AddVariable( "frac04"   , &fFrac04    );
78	+	fLowPtReader->AddVariable( "frac05"   , &fFrac05    );
79	+
80		fReader        = 0;
81		fReader        = new TMVA::Reader( "!Color:!Silent:Error" );
82		if (fType == kBaseline) {
#	Line 76 | Line 97 | void JetIDMVA::Initialize( JetIDMVA::Cut
97		fReader->AddVariable( "frac04"   , &fFrac04    );
98		fReader->AddVariable( "frac05"   , &fFrac05    );
99		}
100	<	fReader->BookMVA(fLowPtMethodName  , iLowPtWeights );
100	>	if (fType == k42) {
101	>	fReader->AddVariable( "frac01"   , &fFrac01    );
102	>	fReader->AddVariable( "frac02"   , &fFrac02    );
103	>	fReader->AddVariable( "frac03"   , &fFrac03    );
104	>	fReader->AddVariable( "frac04"   , &fFrac04    );
105	>	fReader->AddVariable( "frac05"   , &fFrac05    );
106	>	fReader->AddVariable( "nvtx"     , &fNVtx      );
107	>	fReader->AddVariable( "nNeutrals", &fNNeutrals );
108	>	fReader->AddVariable( "beta"     , &fBeta      );
109	>	fReader->AddVariable( "betaStar" , &fBetaStar  );
110	>	fReader->AddVariable( "dZ"       , &fJDZ1      );
111	>	fReader->AddVariable( "nCharged" , &fNCharged  );
112	>	fReader->AddSpectator( "jetPt"    , &fJPt1      );
113	>	fReader->AddSpectator( "jetEta"   , &fJEta1     );
114	>	}
115	>	if (fType == k52) {
116	>	fReader->AddVariable( "frac01"   , &fFrac01    );
117	>	fReader->AddVariable( "frac02"   , &fFrac02    );
118	>	fReader->AddVariable( "frac03"   , &fFrac03    );
119	>	fReader->AddVariable( "frac04"   , &fFrac04    );
120	>	fReader->AddVariable( "frac05"   , &fFrac05    );
121	>	fReader->AddVariable( "dR2Mean"  , &fDR2Mean   );
122	>	fReader->AddVariable( "nvtx"     , &fNVtx      );
123	>	fReader->AddVariable( "nNeutrals", &fNNeutrals );
124	>	fReader->AddVariable( "beta"     , &fBeta      );
125	>	fReader->AddVariable( "betaStar" , &fBetaStar  );
126	>	fReader->AddVariable( "dZ"       , &fJDZ1      );
127	>	fReader->AddVariable( "nCharged" , &fNCharged  );
128	>	fReader->AddSpectator( "jetPt"    , &fJPt1      );
129	>	fReader->AddSpectator( "jetEta"   , &fJEta1     );
130	>	}
131	>
132	>	fLowPtReader->BookMVA(fLowPtMethodName  , iLowPtWeights );
133		fReader->BookMVA(fHighPtMethodName , iHighPtWeights );
134		std::cout << "Jet ID MVA Initialization\n";
135		std::cout << "MethodName : " << fLowPtMethodName << " , type == " << fType << std::endl;
136
137	+	std::string lCutId = "JetIdParams";
138	+	if(fType == k42)  lCutId = "PuJetIdOptMVA_wp";
139	+	if(fType == k52)  lCutId = "full_5x_wp";
140	+	if(fType == kCut) lCutId = "PuJetIdCutBased_wp";
141		//Load Cut Matrix
142	<	edm::ParameterSet lConfig = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>("JetIdParams");
142	>	edm::ParameterSet lDefConfig = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>("JetIdParams");
143	>	edm::ParameterSet lConfig    = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>(lCutId);
144		std::string lCutType = "Tight";
145		if(fCutType == kMedium) lCutType = "Medium";
146		if(fCutType == kLoose ) lCutType = "Loose";
147	<	std::vector<double> lPt010  = lConfig.getParameter<std::vector<double> >(("Pt010_" +lCutType).c_str());
148	<	std::vector<double> lPt1020 = lConfig.getParameter<std::vector<double> >(("Pt1020_"+lCutType).c_str());
149	<	std::vector<double> lPt2030 = lConfig.getParameter<std::vector<double> >(("Pt2030_"+lCutType).c_str());
150	<	std::vector<double> lPt3050 = lConfig.getParameter<std::vector<double> >(("Pt3050_"+lCutType).c_str());
151	<	for(int i0 = 0; i0 < 4; i0++) fMVACut[0][i0] = lPt010 [i0];
152	<	for(int i0 = 0; i0 < 4; i0++) fMVACut[1][i0] = lPt1020[i0];
153	<	for(int i0 = 0; i0 < 4; i0++) fMVACut[2][i0] = lPt2030[i0];
154	<	for(int i0 = 0; i0 < 4; i0++) fMVACut[3][i0] = lPt3050[i0];
155	<	//std::cout << " Working Points : << " << std::endl;
156	<	//for(int i0 = 0; i0 < 4; i0++) for(int i1 = 0; i1 < 4; i1++)
157	<	//  std::cout << " ==> " << i0 << " -- " << i1 << " -- " << fMVACut[i0][i1] << std::endl;
147	>	if(fCutType == kMET   ) lCutType = "MET";
148	>	if(fType != kCut) {
149	>	std::string lLowPtCut = "MET";
150	>	std::vector<double> lPt010  = lDefConfig.getParameter<std::vector<double> >(("Pt010_" +lLowPtCut).c_str());
151	>	std::vector<double> lPt1020 = lConfig.getParameter<std::vector<double> >(("Pt1020_"+lCutType).c_str());
152	>	std::vector<double> lPt2030 = lConfig.getParameter<std::vector<double> >(("Pt2030_"+lCutType).c_str());
153	>	std::vector<double> lPt3050 = lConfig.getParameter<std::vector<double> >(("Pt3050_"+lCutType).c_str());
154	>	for(int i0 = 0; i0 < 4; i0++) fMVACut[0][i0] = lPt010 [i0];
155	>	for(int i0 = 0; i0 < 4; i0++) fMVACut[1][i0] = lPt1020[i0];
156	>	for(int i0 = 0; i0 < 4; i0++) fMVACut[2][i0] = lPt2030[i0];
157	>	for(int i0 = 0; i0 < 4; i0++) fMVACut[3][i0] = lPt3050[i0];
158	>	}
159	>	if(fType == kCut) {
160	>	for(int i0 = 0; i0 < 2; i0++) {
161	>	std::string lFullCutType = lCutType;
162	>	if(i0 == 0) lFullCutType = "BetaStar"+ lCutType;
163	>	if(i0 == 1) lFullCutType = "RMS"     + lCutType;
164	>	std::vector<double> pt010  = lConfig.getParameter<std::vector<double> >(("Pt010_" +lFullCutType).c_str());
165	>	std::vector<double> pt1020 = lConfig.getParameter<std::vector<double> >(("Pt1020_"+lFullCutType).c_str());
166	>	std::vector<double> pt2030 = lConfig.getParameter<std::vector<double> >(("Pt2030_"+lFullCutType).c_str());
167	>	std::vector<double> pt3050 = lConfig.getParameter<std::vector<double> >(("Pt3050_"+lFullCutType).c_str());
168	>	if(i0 == 0) {
169	>	for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[0][i2] = pt010 [i2];
170	>	for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[1][i2] = pt1020[i2];
171	>	for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[2][i2] = pt2030[i2];
172	>	for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[3][i2] = pt3050[i2];
173	>	}
174	>	if(i0 == 1) {
175	>	for(int i2 = 0; i2 < 4; i2++) fRMSCut[0][i2] = pt010 [i2];
176	>	for(int i2 = 0; i2 < 4; i2++) fRMSCut[1][i2] = pt1020[i2];
177	>	for(int i2 = 0; i2 < 4; i2++) fRMSCut[2][i2] = pt2030[i2];
178	>	for(int i2 = 0; i2 < 4; i2++) fRMSCut[3][i2] = pt3050[i2];
179	>	}
180	>	}
181	>	}
182		}
183
184		//--------------------------------------------------------------------------------------------------
#	Line 116 | Line 198 | Double_t JetIDMVA::MVAValue(
198		Float_t iFrac02  ,
199		Float_t iFrac03  ,
200		Float_t iFrac04  ,
201	<	Float_t iFrac05
201	>	Float_t iFrac05  ,
202	>	Float_t iDR2Mean
203		){
204
205		if(!fIsInitialized) {
#	Line 140 | Line 223 | Double_t JetIDMVA::MVAValue(
223		fFrac03    = iFrac03;
224		fFrac04    = iFrac04;
225		fFrac05    = iFrac05;
226	+	fDR2Mean   = iDR2Mean;
227
228		Double_t lMVA = -9999;
229	<	if(iJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName  );
229	>	if(iJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName  );
230		if(iJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
231
232		return lMVA;
#	Line 150 | Line 234 | Double_t JetIDMVA::MVAValue(
234		//--------------------------------------------------------------------------------------------------
235		Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices,
236		FactorizedJetCorrector *iJetCorrector,
237	<	const PileupEnergyDensityCol *iPileupEnergyDensity) {
237	>	const PileupEnergyDensityCol *iPileupEnergyDensity,
238	>	RhoUtilities::RhoType type) {
239
240		if(!JetTools::passPFLooseId(iJet))                 return false;
241	<	if(iJet->Pt() < fJetPtMin)                         return false;
157	<	if(fabs(iJet->Eta()) > 4.99)                       return true; //==> Castor
158	<	//if(iJet->Pt() > 50)                                return true; //==> we can raise this
241	>	if(fabs(iJet->Eta()) > 4.99)                       return false;
242
243		double lMVA = MVAValue   (iJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity);
244	<	double lPt  = correctedPt(iJet,                  iJetCorrector,iPileupEnergyDensity);
245	<
244	>	double lPt  = correctedPt(iJet,                  iJetCorrector,iPileupEnergyDensity,type);
245	>	if(lPt < fJetPtMin)                         return false;
246	>
247		int lPtId = 0;
248	<	if(lPt > 10 && iJet->Pt() < 20) lPtId = 1;
249	<	if(lPt > 20 && iJet->Pt() < 30) lPtId = 2;
248	>	if(lPt > 10 && lPt < 20) lPtId = 1;
249	>	if(lPt > 20 && lPt < 30) lPtId = 2;
250		if(lPt > 30                   ) lPtId = 3;
251
252		int lEtaId = 0;
#	Line 173 | Line 257 | Bool_t JetIDMVA::pass(const PFJet *iJet,
257		double lMVACut = fMVACut[lPtId][lEtaId];
258		if(lMVA < lMVACut) return false;
259		return true;
260	<	//if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
261	<	//if(correctedPt(iJet,iJetCorrector,iPileupEnergyDensity) < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false;
262	<	//This line is a bug in the Met training
263	<	//if(lMVA < -0.8)                            return false;
264	<	//if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
260	>	}
261	>	//--------------------------------------------------------------------------------------------------
262	>	Bool_t JetIDMVA::passCut(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices) {
263	>	if(!JetTools::passPFLooseId(iJet))                 return false;
264	>	if(iJet->Pt()        < fJetPtMin) return false;
265	>	if(fabs(iJet->Eta()) > 4.99)      return false;
266	>	if(fType == kCut) passCut(iJet,iVertex,iVertices);
267	>
268	>	double lPt = iJet->Pt();
269	>	int lPtId = 0;
270	>	if(lPt > 10 && lPt < 20) lPtId = 1;
271	>	if(lPt > 20 && lPt < 30) lPtId = 2;
272	>	if(lPt > 30                   ) lPtId = 3;
273	>
274	>	int lEtaId = 0;
275	>	if(fabs(iJet->Eta()) > 2.5  && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
276	>	if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
277	>	if(fabs(iJet->Eta()) > 3.0  && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
278	>	float betaStarModified = JetTools::betaStarClassic(iJet,iVertex,iVertices)/log(iVertices ->GetEntries()-0.64);
279	>	float dR2Mean          = JetTools::dR2Mean(iJet,-1);
280	>
281	>	if(betaStarModified < fBetaStarCut[lPtId][lEtaId] &&
282	>	dR2Mean          < fRMSCut     [lPtId][lEtaId]
283	>	) return true;
284	>
285	>	return false;
286		}
287		//--------------------------------------------------------------------------------------------------
288		Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices) {
289		if(!JetTools::passPFLooseId(iJet))                 return false;
290		if(iJet->Pt()        < fJetPtMin) return false;
291	<	//if(iJet->Pt()        > 50)        return true; //==> we can raise this
292	<	if(fabs(iJet->Eta()) > 4.99)     return true; //==> Castor
291	>	if(fabs(iJet->Eta()) > 4.99)      return false;
292	>	if(fType == kCut) passCut(iJet,iVertex,iVertices);
293		double lMVA = MVAValue(iJet,iVertex,iVertices);
294
295		int lPtId = 0;
#	Line 230 | Line 335 | Double_t JetIDMVA::MVAValue(const PFJet
335		fNNeutrals = iJet->NeutralMultiplicity();
336
337		fDRMean    = JetTools::dRMean(iJet,-1);
338	+	fDR2Mean   = JetTools::dR2Mean(iJet,-1);
339		fFrac01    = JetTools::frac  (iJet,0.1,0. ,-1);
340		fFrac02    = JetTools::frac  (iJet,0.2,0.1,-1);
341		fFrac03    = JetTools::frac  (iJet,0.3,0.2,-1);
#	Line 237 | Line 343 | Double_t JetIDMVA::MVAValue(const PFJet
343		fFrac05    = JetTools::frac  (iJet,0.5,0.4,-1);
344
345		double lMVA = 0;
346	<	if(fJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName  );
346	>	if(fJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName  );
347		if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
348		if (printDebug == kTRUE) {
349		std::cout << "Debug Jet MVA: "
#	Line 256 | Line 362 | Double_t JetIDMVA::MVAValue(const PFJet
362		<< fFrac02    << " "
363		<< fFrac03    << " "
364		<< fFrac04    << " "
365	<	<< fFrac05
365	>	<< fFrac05    << " "
366	>	<< fDRMean
367		<< " === : === "
368		<< lMVA << " "
369		<< std::endl;
#	Line 285 | Line 392 | Double_t JetIDMVA::MVAValue(const PFJet
392		fNCharged  = iJet->ChargedMultiplicity();
393		fNNeutrals = iJet->NeutralMultiplicity();
394
395	<	fDRMean    = JetTools::dRMean(iJet,-1);
396	<	fFrac01    = JetTools::frac  (iJet,0.1,0. ,-1);
397	<	fFrac02    = JetTools::frac  (iJet,0.2,0.1,-1);
398	<	fFrac03    = JetTools::frac  (iJet,0.3,0.2,-1);
399	<	fFrac04    = JetTools::frac  (iJet,0.4,0.3,-1);
400	<	fFrac05    = JetTools::frac  (iJet,0.5,0.4,-1);
395	>	fDRMean    = JetTools::dRMean (iJet,-1);
396	>	fDR2Mean   = JetTools::dR2Mean(iJet,-1);
397	>	fFrac01    = JetTools::frac   (iJet,0.1,0. ,-1);
398	>	fFrac02    = JetTools::frac   (iJet,0.2,0.1,-1);
399	>	fFrac03    = JetTools::frac   (iJet,0.3,0.2,-1);
400	>	fFrac04    = JetTools::frac   (iJet,0.4,0.3,-1);
401	>	fFrac05    = JetTools::frac   (iJet,0.5,0.4,-1);
402
403		double lMVA = 0;
404	<	if(fJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName  );
405	<	if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
404	>	if(fJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName  );
405	>	if(fJPt1 > 10) lMVA = fReader     ->EvaluateMVA( fHighPtMethodName );
406
407		if (printDebug == kTRUE) {
408		std::cout << "Debug Jet MVA: "
#	Line 313 | Line 421 | Double_t JetIDMVA::MVAValue(const PFJet
421		<< fFrac02    << " "
422		<< fFrac03    << " "
423		<< fFrac04    << " "
424	<	<< fFrac05
424	>	<< fFrac05    << " "
425	>	<< fDRMean
426		<< " === : === "
427		<< lMVA << " "
428		<< std::endl;
#	Line 321 | Line 430 | Double_t JetIDMVA::MVAValue(const PFJet
430
431		return lMVA;
432		}
433	<	Double_t JetIDMVA::correctedPt(const PFJet *iJet, FactorizedJetCorrector *iJetCorrector,const PileupEnergyDensityCol *iPUEnergyDensity) {
433	>	Double_t JetIDMVA::correctedPt(const PFJet *iJet, FactorizedJetCorrector *iJetCorrector,
434	>	const PileupEnergyDensityCol *iPUEnergyDensity,
435	>	RhoUtilities::RhoType type) {
436	>	Double_t Rho = 0.0;
437	>	switch(type) {
438	>	case RhoUtilities::MIT_RHO_VORONOI_HIGH_ETA:
439	>	Rho = iPUEnergyDensity->At(0)->Rho();
440	>	break;
441	>	case RhoUtilities::MIT_RHO_VORONOI_LOW_ETA:
442	>	Rho = iPUEnergyDensity->At(0)->RhoLowEta();
443	>	break;
444	>	case RhoUtilities::MIT_RHO_RANDOM_HIGH_ETA:
445	>	Rho = iPUEnergyDensity->At(0)->RhoRandom();
446	>	break;
447	>	case RhoUtilities::MIT_RHO_RANDOM_LOW_ETA:
448	>	Rho = iPUEnergyDensity->At(0)->RhoRandomLowEta();
449	>	break;
450	>	case RhoUtilities::CMS_RHO_RHOKT6PFJETS:
451	>	Rho = iPUEnergyDensity->At(0)->RhoKt6PFJets();
452	>	break;
453	>	default:
454	>	// use the old default
455	>	Rho = iPUEnergyDensity->At(0)->Rho();
456	>	break;
457	>	}
458	>
459		const FourVectorM rawMom = iJet->RawMom();
460		iJetCorrector->setJetEta(rawMom.Eta());
461		iJetCorrector->setJetPt (rawMom.Pt());
462		iJetCorrector->setJetPhi(rawMom.Phi());
463		iJetCorrector->setJetE  (rawMom.E());
464	<	iJetCorrector->setRho   (iPUEnergyDensity->At(0)->RhoHighEta());
464	>	iJetCorrector->setRho   (Rho);
465		iJetCorrector->setJetA  (iJet->JetArea());
466		iJetCorrector->setJetEMF(-99.0);
467		Double_t correction = iJetCorrector->getCorrection();

Diff Legend

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