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.11 by pharris, Mon May 14 09:42:27 2012 UTC vs.
Revision 1.16 by pharris, Tue Sep 25 15:39:15 2012 UTC

#	Line 58 | Line 58 | void JetIDMVA::Initialize( JetIDMVA::Cut
58		fIsInitialized = kTRUE;
59		fType          = iType;
60		fCutType       = iCutType;
61	+
62	+	std::string lCutId = "JetIdParams";
63	+	if(fType == k42)  lCutId = "PuJetIdOptMVA_wp";
64	+	if(fType == k52)  lCutId = "full_5x_wp";
65	+	if(fType == kCut) lCutId = "PuJetIdCutBased_wp";
66	+	//Load Cut Matrix
67	+	edm::ParameterSet lDefConfig = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>("JetIdParams");
68	+	edm::ParameterSet lConfig    = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>(lCutId);
69	+	std::string lCutType = "Tight";
70	+	if(fCutType == kMedium) lCutType = "Medium";
71	+	if(fCutType == kLoose ) lCutType = "Loose";
72	+	if(fCutType == kMET   ) lCutType = "MET";
73	+	if(fType != kCut) {
74	+	std::string lLowPtCut = "MET";
75	+	std::vector<double> lPt010  = lDefConfig.getParameter<std::vector<double> >(("Pt010_" +lLowPtCut).c_str());
76	+	std::vector<double> lPt1020 = lConfig.getParameter<std::vector<double> >(("Pt1020_"+lCutType).c_str());
77	+	std::vector<double> lPt2030 = lConfig.getParameter<std::vector<double> >(("Pt2030_"+lCutType).c_str());
78	+	std::vector<double> lPt3050 = lConfig.getParameter<std::vector<double> >(("Pt3050_"+lCutType).c_str());
79	+	for(int i0 = 0; i0 < 4; i0++) fMVACut[0][i0] = lPt010 [i0];
80	+	for(int i0 = 0; i0 < 4; i0++) fMVACut[1][i0] = lPt1020[i0];
81	+	for(int i0 = 0; i0 < 4; i0++) fMVACut[2][i0] = lPt2030[i0];
82	+	for(int i0 = 0; i0 < 4; i0++) fMVACut[3][i0] = lPt3050[i0];
83	+	}
84	+	if(fType == kCut) {
85	+	for(int i0 = 0; i0 < 2; i0++) {
86	+	std::string lFullCutType = lCutType;
87	+	if(i0 == 0) lFullCutType = "BetaStar"+ lCutType;
88	+	if(i0 == 1) lFullCutType = "RMS"     + lCutType;
89	+	std::vector<double> pt010  = lConfig.getParameter<std::vector<double> >(("Pt010_" +lFullCutType).c_str());
90	+	std::vector<double> pt1020 = lConfig.getParameter<std::vector<double> >(("Pt1020_"+lFullCutType).c_str());
91	+	std::vector<double> pt2030 = lConfig.getParameter<std::vector<double> >(("Pt2030_"+lFullCutType).c_str());
92	+	std::vector<double> pt3050 = lConfig.getParameter<std::vector<double> >(("Pt3050_"+lFullCutType).c_str());
93	+	if(i0 == 0) {
94	+	for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[0][i2] = pt010 [i2];
95	+	for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[1][i2] = pt1020[i2];
96	+	for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[2][i2] = pt2030[i2];
97	+	for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[3][i2] = pt3050[i2];
98	+	}
99	+	if(i0 == 1) {
100	+	for(int i2 = 0; i2 < 4; i2++) fRMSCut[0][i2] = pt010 [i2];
101	+	for(int i2 = 0; i2 < 4; i2++) fRMSCut[1][i2] = pt1020[i2];
102	+	for(int i2 = 0; i2 < 4; i2++) fRMSCut[2][i2] = pt2030[i2];
103	+	for(int i2 = 0; i2 < 4; i2++) fRMSCut[3][i2] = pt3050[i2];
104	+	}
105	+	}
106	+	return;
107	+	}
108	+
109	+
110		fLowPtReader   = 0;
111		fLowPtReader   = new TMVA::Reader( "!Color:!Silent:Error" );
112		fLowPtReader->AddVariable( "nvtx"     , &fNVtx      );
#	Line 128 | Line 177 | void JetIDMVA::Initialize( JetIDMVA::Cut
177		fReader->AddSpectator( "jetPt"    , &fJPt1      );
178		fReader->AddSpectator( "jetEta"   , &fJEta1     );
179		}
180	+	if (fType == kQGP) {
181	+	fReader->AddVariable( "nvtx"      , &fNVtx       );
182	+	fReader->AddVariable( "jetPt"     , &fJPt1       );
183	+	fReader->AddVariable( "jetEta"    , &fJEta1      );
184	+	fReader->AddVariable( "jetPhi"    , &fJPhi1      );
185	+	fReader->AddVariable( "beta"      , &fBeta      );
186	+	fReader->AddVariable( "betaStar"  , &fBetaStar  );
187	+	fReader->AddVariable( "nParticles", &fNNeutrals );
188	+	fReader->AddVariable( "nCharged"  , &fNCharged  );
189	+	fReader->AddVariable( "dRMean"    , &fDRMean    );
190	+	fReader->AddVariable( "ptD"       , &fPtD       );
191	+	fReader->AddVariable( "frac01"    , &fFrac01    );
192	+	fReader->AddVariable( "frac02"    , &fFrac02    );
193	+	fReader->AddVariable( "frac03"    , &fFrac03    );
194	+	fReader->AddVariable( "frac04"    , &fFrac04    );
195	+	fReader->AddVariable( "frac05"    , &fFrac05    );
196	+	}
197
198		fLowPtReader->BookMVA(fLowPtMethodName  , iLowPtWeights );
199		fReader->BookMVA(fHighPtMethodName , iHighPtWeights );
200		std::cout << "Jet ID MVA Initialization\n";
201		std::cout << "MethodName : " << fLowPtMethodName << " , type == " << fType << std::endl;
202
137	–	std::string lCutId = "JetIdParams";
138	–	if(fType == k42) lCutId = "PuJetIdOptMVA_wp";
139	–	if(fType == k52) lCutId = "full_5x_wp";
140	–	//Load Cut Matrix
141	–	edm::ParameterSet lDefConfig = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>("JetIdParams");
142	–	edm::ParameterSet lConfig    = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>(lCutId);
143	–	std::string lCutType = "Tight";
144	–	if(fCutType == kMedium) lCutType = "Medium";
145	–	if(fCutType == kLoose ) lCutType = "Loose";
146	–	if(fCutType == kMET   ) lCutType = "MET";
147	–	std::string lLowPtCut = "MET";
148	–	std::vector<double> lPt010  = lDefConfig.getParameter<std::vector<double> >(("Pt010_" +lLowPtCut).c_str());
149	–	std::vector<double> lPt1020 = lConfig.getParameter<std::vector<double> >(("Pt1020_"+lCutType).c_str());
150	–	std::vector<double> lPt2030 = lConfig.getParameter<std::vector<double> >(("Pt2030_"+lCutType).c_str());
151	–	std::vector<double> lPt3050 = lConfig.getParameter<std::vector<double> >(("Pt3050_"+lCutType).c_str());
152	–	for(int i0 = 0; i0 < 4; i0++) fMVACut[0][i0] = lPt010 [i0];
153	–	for(int i0 = 0; i0 < 4; i0++) fMVACut[1][i0] = lPt1020[i0];
154	–	for(int i0 = 0; i0 < 4; i0++) fMVACut[2][i0] = lPt2030[i0];
155	–	for(int i0 = 0; i0 < 4; i0++) fMVACut[3][i0] = lPt3050[i0];
156	–	//std::cout << " Working Points : << " << std::endl;
157	–	//for(int i0 = 0; i0 < 4; i0++) for(int i1 = 0; i1 < 4; i1++)
158	–	//  std::cout << " ==> " << i0 << " -- " << i1 << " -- " << fMVACut[i0][i1] << std::endl;
203		}
204
205		//--------------------------------------------------------------------------------------------------
#	Line 209 | Line 253 | Double_t JetIDMVA::MVAValue(
253		return lMVA;
254		}
255		//--------------------------------------------------------------------------------------------------
256	+	Bool_t JetIDMVA::passPt(const PFJet *iJet, FactorizedJetCorrector *iJetCorrector,
257	+	const PileupEnergyDensityCol *iPileupEnergyDensity,
258	+	RhoUtilities::RhoType type) {
259	+	if(iJetCorrector == 0) return (iJet->Pt() > fJetPtMin);
260	+	double lPt  = correctedPt(iJet,                  iJetCorrector,iPileupEnergyDensity,type);
261	+	if(lPt < fJetPtMin)                         return false;
262	+	return true;
263	+	}
264	+	//--------------------------------------------------------------------------------------------------
265		Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices,
266		FactorizedJetCorrector *iJetCorrector,
267		const PileupEnergyDensityCol *iPileupEnergyDensity,
#	Line 236 | Line 289 | Bool_t JetIDMVA::pass(const PFJet *iJet,
289		return true;
290		}
291		//--------------------------------------------------------------------------------------------------
292	+	Bool_t JetIDMVA::passCut(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices) {
293	+	if(!JetTools::passPFLooseId(iJet))                 return false;
294	+	if(iJet->Pt()        < fJetPtMin) return false;
295	+	if(fabs(iJet->Eta()) > 4.99)      return false;
296	+	//if(fType == kCut) passCut(iJet,iVertex,iVertices);
297	+
298	+	double lPt = iJet->Pt();
299	+	int lPtId = 0;
300	+	if(lPt > 10 && lPt < 20) lPtId = 1;
301	+	if(lPt > 20 && lPt < 30) lPtId = 2;
302	+	if(lPt > 30                   ) lPtId = 3;
303	+
304	+	int lEtaId = 0;
305	+	if(fabs(iJet->Eta()) > 2.5  && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
306	+	if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
307	+	if(fabs(iJet->Eta()) > 3.0  && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
308	+	float betaStarModified = JetTools::betaStarClassic(iJet,iVertex,iVertices)/log(iVertices ->GetEntries()-0.64);
309	+	float dR2Mean          = JetTools::dR2Mean(iJet,-1);
310	+
311	+	if(betaStarModified < fBetaStarCut[lPtId][lEtaId] &&
312	+	dR2Mean          < fRMSCut     [lPtId][lEtaId]
313	+	) return true;
314	+
315	+	return false;
316	+	}
317	+	//--------------------------------------------------------------------------------------------------
318		Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices) {
319		if(!JetTools::passPFLooseId(iJet))                 return false;
320		if(iJet->Pt()        < fJetPtMin) return false;
321		if(fabs(iJet->Eta()) > 4.99)      return false;
322	+	if(fType == kCut) return passCut(iJet,iVertex,iVertices);
323		double lMVA = MVAValue(iJet,iVertex,iVertices);
324
325		int lPtId = 0;
#	Line 380 | Line 460 | Double_t JetIDMVA::MVAValue(const PFJet
460
461		return lMVA;
462		}
463	+	//--------------------------------------------------------------------------------------------------
464	+	Double_t* JetIDMVA::QGValue(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices, //Vertex here is the PV
465	+	FactorizedJetCorrector *iJetCorrector,
466	+	const PileupEnergyDensityCol *iPileupEnergyDensity,
467	+	Bool_t printDebug) {
468	+
469	+	Double_t *lId = new double[3];
470	+	lId[0] = -2;
471	+	lId[1] = -2;
472	+	lId[2] = -2;
473	+	if (!fIsInitialized) {
474	+	std::cout << "Error: JetIDMVA not properly initialized.\n";
475	+	return lId;
476	+	}
477	+	if(!JetTools::passPFLooseId(iJet)) return lId;
478	+
479	+	fJPt1       = correctedPt(iJet,iJetCorrector,iPileupEnergyDensity);
480	+	if(fJPt1 < 20) return lId;
481	+
482	+	//set all input variables
483	+	fNVtx       = iVertices->GetEntries();
484	+	fJEta1      = iJet->RawMom().Eta();
485	+	fJPhi1      = iJet->RawMom().Phi();
486	+	fJD01       = JetTools::impactParameter(iJet,iVertex);
487	+	fJDZ1       = JetTools::impactParameter(iJet,iVertex,true);
488	+	fBeta       = JetTools::Beta(iJet,iVertex,fDZCut);
489	+	fBetaStar   = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
490	+	fNCharged   = iJet->ChargedMultiplicity();
491	+	fNNeutrals  = iJet->NeutralMultiplicity();
492	+	fNParticles = fNCharged+fNNeutrals;
493	+	fPtD        = JetTools::W(iJet,-1,0);
494	+
495	+	fDRMean    = JetTools::dRMean(iJet,-1);
496	+	fDR2Mean   = JetTools::dR2Mean(iJet,-1);
497	+	fFrac01    = JetTools::frac  (iJet,0.1,0. ,-1);
498	+	fFrac02    = JetTools::frac  (iJet,0.2,0.1,-1);
499	+	fFrac03    = JetTools::frac  (iJet,0.3,0.2,-1);
500	+	fFrac04    = JetTools::frac  (iJet,0.4,0.3,-1);
501	+	fFrac05    = JetTools::frac  (iJet,0.5,0.4,-1);
502	+
503	+	double lMVA = 0;
504	+	lId[0] = fReader->EvaluateMulticlass( fHighPtMethodName )[0];
505	+	lId[1] = fReader->EvaluateMulticlass( fHighPtMethodName )[1];
506	+	lId[2] = fReader->EvaluateMulticlass( fHighPtMethodName )[2];
507	+	if (printDebug == kTRUE) {
508	+	std::cout << "Debug Jet MVA: "
509	+	<< fNVtx      << " "
510	+	<< fJPt1      << " "
511	+	<< fJEta1     << " "
512	+	<< fJPhi1     << " "
513	+	<< fJD01      << " "
514	+	<< fJDZ1      << " "
515	+	<< fBeta      << " "
516	+	<< fBetaStar  << " "
517	+	<< fNCharged  << " "
518	+	<< fNNeutrals << " "
519	+	<< fDRMean    << " "
520	+	<< fFrac01    << " "
521	+	<< fFrac02    << " "
522	+	<< fFrac03    << " "
523	+	<< fFrac04    << " "
524	+	<< fFrac05    << " "
525	+	<< fDRMean
526	+	<< " === : === "
527	+	<< lMVA << " "
528	+	<< std::endl;
529	+	}
530	+	return lId;
531	+	}
532		Double_t JetIDMVA::correctedPt(const PFJet *iJet, FactorizedJetCorrector *iJetCorrector,
533		const PileupEnergyDensityCol *iPUEnergyDensity,
534	<	RhoUtilities::RhoType type) {
534	>	RhoUtilities::RhoType type,int iId) {
535		Double_t Rho = 0.0;
536		switch(type) {
537		case RhoUtilities::MIT_RHO_VORONOI_HIGH_ETA:
#	Line 414 | Line 563 | Double_t JetIDMVA::correctedPt(const PFJ
563		iJetCorrector->setRho   (Rho);
564		iJetCorrector->setJetA  (iJet->JetArea());
565		iJetCorrector->setJetEMF(-99.0);
566	<	Double_t correction = iJetCorrector->getCorrection();
566	>	Double_t correction = 1.;
567	>	if(iId < 0 || iId == 100) correction = iJetCorrector->getCorrection();
568	>	std::vector<float> lCorrections; if(iId != -1 && iId != 100) lCorrections = iJetCorrector->getSubCorrections();
569	>	if(iId > -1 && iId < int(lCorrections.size())) correction = lCorrections[iId];
570	>	if(iId == 100) {
571	>	iJetCorrector->setJetEta(rawMom.Eta());
572	>	iJetCorrector->setJetPt (rawMom.Pt());
573	>	iJetCorrector->setJetPhi(rawMom.Phi());
574	>	iJetCorrector->setJetE  (rawMom.E());
575	>	iJetCorrector->setRho   (Rho);
576	>	iJetCorrector->setJetA  (iJet->JetArea());
577	>	iJetCorrector->setJetEMF(-99.0);
578	>	lCorrections = iJetCorrector->getSubCorrections();
579	>	double correction2 = 1;
580	>	correction2 *= lCorrections[0];
581	>	correction = correction-correction2;
582	>	}
583	>
584		return rawMom.Pt()*correction;
585		}

Diff Legend

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