ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/JetFitAnalyzer/test/JetFinderAnalyzer.cc

Comparing UserCode/JetFitAnalyzer/test/JetFinderAnalyzer.cc (file contents):
Revision 1.14 by dnisson, Wed Sep 23 03:31:30 2009 UTC vs.
Revision 1.19 by dnisson, Tue Nov 10 01:25:49 2009 UTC

#	Line 24 | Line 24
24		#include <limits>
25		#include <cmath>
26		#include <cstdlib>
27	+	#include <iostream>
28		#include <fstream>
29		#include <sstream>
30
#	Line 70 | Line 71 | private:
71		return sum1 * sum2;
72		}
73
74	<	static void jetfinder(TMinuit *gMinuit, TH2 *hist, int ngauss) {
75	<	double dist_sq = numeric_limits<double>::infinity();
75	<	unique_jets[hist].resize(ngauss);
76	<	int nbinsX = hist->GetXaxis()->GetNbins();
77	<	int nbinsY = hist->GetYaxis()->GetNbins();
78	<	double XbinSize = (hist->GetXaxis()->GetXmax()
79	<	- hist->GetXaxis()->GetXmin())
80	<	/ static_cast<double>(nbinsX);
81	<	double YbinSize = (hist->GetYaxis()->GetXmax()
82	<	- hist->GetYaxis()->GetXmin())
83	<	/ static_cast<double>(nbinsY);
84	<	for (int i = 0; i < ngauss; i++) {
85	<	double N, mu_x, mu_y, sig, err, lo, hi;
86	<	int iuint;
87	<	TString name;
88	<	gMinuit->mnpout(4*i, name, N, err, lo, hi, iuint);
89	<	gMinuit->mnpout(4*i + 1, name, mu_x, err, lo, hi, iuint);
90	<	gMinuit->mnpout(4*i + 2, name, mu_y, err, lo, hi, iuint);
91	<	gMinuit->mnpout(4*i + 3, name, sig, err, lo, hi, iuint);
92	<	for (int j = 0; j < i; j++) {
93	<	double N2, mu_x2, mu_y2, sig2;
94	<	gMinuit->mnpout(4*j, name, N2, err, lo, hi, iuint);
95	<	gMinuit->mnpout(4*j + 1, name, mu_x2, err, lo, hi, iuint);
96	<	gMinuit->mnpout(4*j + 2, name, mu_y2, err, lo, hi, iuint);
97	<	gMinuit->mnpout(4*j + 3, name, sig2, err, lo, hi, iuint);
98	<	double _dist_sq = (mu_x2 - mu_x)*(mu_x2 - mu_x)
99	<	+ (mu_y2 - mu_y)*(mu_y2 - mu_y);
100	<	if (_dist_sq < dist_sq)
101	<	dist_sq = _dist_sq;
102	<	}
103	<
104	<	jet j;
105	<	j.energy = N;
106	<	j.eta = mu_x; j.phi = mu_y;
107	<	unique_jets[hist][ngauss-1].push_back(j);
108	<	}
109	<	}
110	<
111	<	virtual void beginJob(const edm::EventSetup&);
112	<	virtual TH2D* make_histo(const edm::Event&, const edm::EventSetup&);
74	>	virtual void beginJob(const edm::EventSetup &es);
75	>	virtual TH2D * make_histo(const edm::Event &, const edm::EventSetup &);
76		virtual jetfit::model_def& make_model_def(const edm::Event&,
77		const edm::EventSetup&,
78		TH2 *);
79		virtual void analyze_results(jetfit::results r,
80		std::vector<jetfit::trouble> t, TH2 *);
81		vector<reco::Candidate *> get_particles(const edm::Event&);
119	–	void fetchCandidateCollection(edm::Handle<reco::PFCandidateCollection>&,
120	–	const edm::InputTag&, const edm::Event&) const;
121	–	void fetchCandidateCollection(edm::Handle<reco::GenParticleCollection>&,
122	–	const edm::InputTag&, const edm::Event&) const;
82
83		fstream ofs;
84		edm::InputTag inputTagPFCandidates_;
#	Line 128 | Line 87 | private:
87		double smear_;
88		int smear_coord_;
89		string jet_algo_;
90	+	string reclustered_jet_algo_;
91		};
92
93		map<TH2 *, vector< vector< JetFinderAnalyzer::jet > > >
#	Line 150 | Line 110 | JetFinderAnalyzer::JetFinderAnalyzer(con
110		// 0 = eta-phi smear
111		// 1 = proper angle smear
112		jet_algo_ = pSet.getParameter<string>("jet_algo");
113	<	set_user_minuit(jetfinder);
154	<	}
155	<
156	<	void
157	<	JetFinderAnalyzer::fetchCandidateCollection(edm::Handle<reco::PFCandidateCollection>& c,
158	<	const edm::InputTag& tag,
159	<	const edm::Event& iEvent) const {
160	<
161	<	bool found = iEvent.getByLabel(tag, c);
162	<
163	<	if(!found ) {
164	<	ostringstream  err;
165	<	err<<" cannot get PFCandidates: "
166	<	<<tag<<endl;
167	<	edm::LogError("PFCandidates")<<err.str();
168	<	throw cms::Exception( "MissingProduct", err.str());
169	<	}
170	<
171	<	}
172	<
173	<	void
174	<	JetFinderAnalyzer::fetchCandidateCollection(edm::Handle<reco::GenParticleCollection>& c,
175	<	const edm::InputTag& tag,
176	<	const edm::Event& iEvent) const {
177	<
178	<	bool found = iEvent.getByLabel(tag, c);
179	<
180	<	if(!found ) {
181	<	ostringstream  err;
182	<	err<<" cannot get GenParticles: "
183	<	<<tag<<endl;
184	<	edm::LogError("GenParticles")<<err.str();
185	<	throw cms::Exception( "MissingProduct", err.str());
186	<	}
187	<
188	<	}
189	<
190	<	vector<reco::Candidate *> JetFinderAnalyzer::get_particles(const edm::Event &evt) {
191	<	// fill unreduced histo
192	<	edm::Handle<reco::GenParticleCollection> hRaw;
193	<	edm::Handle<reco::PFCandidateCollection> hPFlow;
194	<	if (info_type_ == 0) {
195	<	fetchCandidateCollection(hRaw,
196	<	inputTagGenParticles_,
197	<	evt);
198	<	}
199	<	if (info_type_ == 1) {
200	<	fetchCandidateCollection(hPFlow,
201	<	inputTagPFCandidates_,
202	<	evt);
203	<	}
204	<
205	<	vector<reco::Candidate *> particles;
206	<
207	<	switch (info_type_) {
208	<	case 0:
209	<	for (unsigned i = 0; i < hRaw->size(); i++) {
210	<	if ((*hRaw)[i].status() == 1) {
211	<	particles.push_back((reco::Candidate *)&((*hRaw)[i]));
212	<	}
213	<	}
214	<	break;
215	<	case 1:
216	<	for (unsigned i = 0; i < hPFlow->size(); i++) {
217	<	particles.push_back((reco::Candidate *)&((*hPFlow)[i]));
218	<	}
219	<	break;
220	<	default:
221	<	cerr << "Unknown event type" << endl; // TODO use MessageLogger
222	<	}
223	<
224	<	return particles;
113	>	set_user_minuit(0); // we don't need this feature--it's deprecated
114		}
115
116		TH2D * JetFinderAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&) {
117	<	const reco::Jet * highest_e_jet;
117	>	const reco::Jet * highest_e_jet = 0;
118
119		if (info_type_ == 0) {
120		edm::Handle< vector<reco::GenJet> > evtJets;
121		evt.getByLabel(jet_algo_, evtJets);
122		for (unsigned i = 0; i < evtJets->size(); i++) {
123	<	if (i == 0 || (*evtJets)[i].energy() > highest_e_jet->energy()) {
123	>	if (highest_e_jet == 0
124	>	|| (*evtJets)[i].energy() > highest_e_jet->energy()) {
125		highest_e_jet = &((*evtJets)[i]);
126		}
127		}
#	Line 240 | Line 130 | TH2D * JetFinderAnalyzer::make_histo(con
130		edm::Handle< vector<reco::PFJet> > evtJets;
131		evt.getByLabel(jet_algo_, evtJets);
132		for (unsigned i = 0; i < evtJets->size(); i++) {
133	<	if (i == 0 || (*evtJets)[i].energy() > highest_e_jet->energy()) {
133	>	if (highest_e_jet == 0
134	>	|| (*evtJets)[i].energy() > highest_e_jet->energy()) {
135		highest_e_jet = &((*evtJets)[i]);
136		}
137		}
138		}
248	–	cout << "found highest e jet" << endl;
139
140		if (highest_e_jet == 0) {
141	<	cerr << "No fatjets found!" << endl; exit(1);
141	>	cerr << "No fat jets found!" << endl;
142	>	return 0;
143		}
144
145		vector<const reco::Candidate *> particles =
146		highest_e_jet->getJetConstituentsQuick();
147	<	cout << "found highest e jet" << endl;
257	<
147	>
148		ostringstream oss;
149		oss << "eta_phi_energy"<<evt.id().event() << flush;
150	<	TH2D *histo = new TH2D(oss.str().c_str(), oss.str().c_str(),
150	>	edm::Service<TFileService> fs;
151	>	TH2D *histo = fs->make<TH2D>(oss.str().c_str(), oss.str().c_str(),
152		30,
153		highest_e_jet->eta()-0.5,
154		highest_e_jet->eta()+0.5,
#	Line 265 | Line 156 | TH2D * JetFinderAnalyzer::make_histo(con
156		highest_e_jet->phi()-0.5,
157		highest_e_jet->phi()+0.5);
158
159	<	for (int i = 0; i < particles.size(); i++) {
269	<	histo->Fill(particles[i]->eta(),
270	<	particles[i]->phi(),
271	<	particles[i]->energy());
272	<	}
273	<
159	>	if (smear_ > 0.0) {
160		// create a smeared histo
161		// create a temporary 2D vector for smeared energies
162		double XbinSize = (histo->GetXaxis()->GetXmax()
#	Line 280 | Line 166 | TH2D * JetFinderAnalyzer::make_histo(con
166		- histo->GetYaxis()->GetXmin()) /
167		static_cast<double>(histo->GetYaxis()->GetNbins());
168		double Xlo = histo->GetXaxis()->GetXmin();
169	+	double Xhi = histo->GetXaxis()->GetXmax();
170		double Ylo = histo->GetYaxis()->GetXmin();
171	+	double Yhi = histo->GetYaxis()->GetXmax();
172		vector< vector<double> > smeared(30, vector<double>(30, 0.0) );
173		switch (smear_coord_) {
174		case 1:
#	Line 288 | Line 176 | TH2D * JetFinderAnalyzer::make_histo(con
176		double N = particles[i]->energy();
177		double x = particles[i]->eta();
178		double y = particles[i]->phi();
179	+	if (y < Ylo) {
180	+	y += 2.0*PI;
181	+	}
182	+	if (y > Yhi) {
183	+	y -= 2.0*PI;
184	+	}
185		// loop over bins and add Gaussian in proper angle to smeared
186		for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) {
187		for (vector< double >::size_type j2 = 0; j2 < 30; j2++) {
188		double eta = static_cast<double>((signed int)i2) * XbinSize +
189	<	Xlo;
190	<	double phi0 = static_cast<double>((signed int)j2) * YbinSize +
191	<	Ylo;
192	<	double phi = acos(cos(phi0));
299	<	phi = sin(phi0) > 0 ? phi : -phi;
300	<
189	>	Xlo - x;
190	>	double phi = static_cast<double>((signed int)j2) * YbinSize +
191	>	Ylo - y;
192	>
193		// transform eta, phi to proper angle
194		double theta = 2.0*atan(exp(-eta));
195		double iota = asin(sin(theta)*sin(phi));
196
197		smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
198		* exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_));
307	–
308	–	// correct for out-of-range phi
309	–	double transform = 0.0;
310	–	if (histo->GetYaxis()->GetXmin() < -PI) {
311	–	transform = -2.0*PI;
312	–	phi += transform;
313	–	iota = asin(sin(theta)*sin(phi));
314	–	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
315	–	* exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_));
316	–	}
317	–	if (histo->GetYaxis()->GetXmax() > PI) {
318	–	transform = 2.0*PI;
319	–	phi += transform;
320	–	iota = asin(sin(theta)*sin(phi));
321	–	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
322	–	* exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_));
323	–	}
199		}
200		}
201		}
#	Line 331 | Line 206 | TH2D * JetFinderAnalyzer::make_histo(con
206		double N = particles[i]->energy();
207		double x = particles[i]->eta();
208		double y = particles[i]->phi();
209	+	if (y < Ylo) {
210	+	y += 2.0*PI;
211	+	}
212	+	if (y > Yhi) {
213	+	y -= 2.0*PI;
214	+	}
215		// loop over bins and add Gaussian to smeared
216		for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) {
217		for (vector< double >::size_type j2 = 0; j2 < 30; j2++) {
218		double eta = static_cast<double>((signed int)i2) * XbinSize
219	<	+ Xlo;
220	<	double phi0 = static_cast<double>((signed int)j2) * YbinSize +
221	<	Ylo;
341	<	double phi = acos(cos(phi0));
342	<	phi = sin(phi0) > 0 ? phi : -phi;
343	<	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
344	<	* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_));
219	>	+ Xlo - x;
220	>	double phi = static_cast<double>((signed int)j2) * YbinSize +
221	>	Ylo - y;
222
223	<	// correct for out-of-range phi
347	<	double transform = 0.0;
348	<	if (histo->GetYaxis()->GetXmin() < -PI) {
349	<	transform = -2.0*PI;
350	<	phi += transform;
351	<	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
352	<	* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_));
353	<	}
354	<	if (histo->GetYaxis()->GetXmax() > PI) {
355	<	transform = 2.0*PI;
356	<	phi += transform;
357	<	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
223	>	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
224		* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_));
359	–	}
225		}
226		}
227		}
#	Line 367 | Line 232 | TH2D * JetFinderAnalyzer::make_histo(con
232		histo->SetBinContent(i, j, smeared[i-1][j-1]);
233		}
234		}
235	+	}
236	+	else {
237	+	// don't smear--just fill with particles
238	+	for (int i = 0; i < particles.size(); i++) {
239	+	histo->Fill(particles[i]->eta(), particles[i]->phi(),
240	+	particles[i]->energy());
241	+	}
242	+	}
243
244		return histo;
245		}
246
247		template <class Jet>
248		void seed_with_jetcoll(vector<Jet> jets,
249	<	jetfit::model_def &_mdef) {
249	>	jetfit::model_def &_mdef, double phi_lo = -PI,
250	>	double phi_hi = PI) {
251		// seed with jet collection
252		int ijset = 0;
253		for (unsigned ij = 0; ij < jets.size(); ij++) {
#	Line 385 | Line 259 | void seed_with_jetcoll(vector<Jet> jets,
259		0.0, 1.0e6);
260		_mdef.set_special_par(ijset, 1, eta, 0.01,
261		0.0, 0.0);
262	<	double mdef_phi = phi > PI ? phi - 2*PI
263	<	: phi;
264	<	_mdef.set_special_par(ijset, 2, mdef_phi, 0.01,
262	>	if (phi < phi_lo) {
263	>	phi += 2.0*PI;
264	>	}
265	>	if (phi > phi_hi) {
266	>	phi -= 2.0*PI;
267	>	}
268	>	_mdef.set_special_par(ijset, 2, phi, 0.01,
269		0.0, 0.0);
270		_mdef.set_special_par(ijset, 3, 0.1, 0.001, 0.0, 0.0);
271		ijset++;
#	Line 401 | Line 279 | jetfit::model_def& JetFinderAnalyzer::ma
279		class jf_model_def : public jetfit::model_def {
280		public:
281		virtual double chisquare_error(double E) {
282	<	return exp(-(4.2 + 0.11*E));
405	<	// study from 09-04-09
282	>	return 0.298 - 0.9557*E; // study 11-09-09
283		}
284		};
285
#	Line 428 | Line 305 | jetfit::model_def& JetFinderAnalyzer::ma
305		if (!found_jet || (*jet_collection)[i].energy() > highest_e_jet.energy())
306		{
307		highest_e_jet = (*jet_collection)[i];
308	+	found_jet = true;
309		}
310		}
311		vector<reco::GenJet> highest_e_jet_coll;
312		highest_e_jet_coll.push_back(highest_e_jet);
313	<	seed_with_jetcoll(highest_e_jet_coll, *_mdef);
313	>	seed_with_jetcoll(highest_e_jet_coll, *_mdef,
314	>	histo->GetYaxis()->GetXmin(),
315	>	histo->GetYaxis()->GetXmax());
316		}
317		if (info_type_ == 1) {
318		edm::Handle< vector<reco::PFJet> > jet_collection;
#	Line 443 | Line 323 | jetfit::model_def& JetFinderAnalyzer::ma
323		if (!found_jet || (*jet_collection)[i].energy() > highest_e_jet.energy())
324		{
325		highest_e_jet = (*jet_collection)[i];
326	+	found_jet = true;
327		}
328		}
329		vector<reco::PFJet> highest_e_jet_coll;
330		highest_e_jet_coll.push_back(highest_e_jet);
331	<	seed_with_jetcoll(highest_e_jet_coll, *_mdef);
331	>	seed_with_jetcoll(highest_e_jet_coll, *_mdef,
332	>	histo->GetYaxis()->GetXmin(),
333	>	histo->GetYaxis()->GetXmax());
334		}
335
336		jetfit::set_model_def(_mdef);
#	Line 475 | Line 358 | jetfit::model_def& JetFinderAnalyzer::ma
358		double Ylo = histo->GetYaxis()->GetXmin();
359		double Yhi = histo->GetYaxis()->GetXmax();
360
361	<	for (int i = 0; i < 60; i++) {
362	<	for (int j = 0; j < 60; j++) {
361	>	for (int i = 0; i < 30; i++) {
362	>	for (int j = 0; j < 30; j++) {
363		double x = (static_cast<double>(i) + 0.5)*XbinSize + Xlo;
364		double y = (static_cast<double>(j) + 0.5)*YbinSize + Ylo;
365		double pval[256];
#	Line 603 | Line 486 | void JetFinderAnalyzer::analyze_results(
486		hist_orig->GetYaxis()->GetNbins(),
487		hist_orig->GetYaxis()->GetXmin(),
488		hist_orig->GetYaxis()->GetXmax());
489	<	TH1 *tf2_histo = tf2->CreateHistogram();
490	<	double XbinSize = (fit_histo->GetXaxis()->GetXmax()
491	<	- fit_histo->GetXaxis()->GetXmin())
492	<	/ static_cast<double>(fit_histo->GetXaxis()->GetNbins());
493	<	double YbinSize = (fit_histo->GetYaxis()->GetXmax()
494	<	- fit_histo->GetYaxis()->GetXmin())
495	<	/ static_cast<double>(fit_histo->GetYaxis()->GetNbins());
496	<	for (int ih = 0; ih < tf2->GetNpx(); ih++) {
497	<	for (int jh = 0; jh < tf2->GetNpy(); jh++) {
498	<	fit_histo->SetBinContent(ih+1, jh+1,
499	<	tf2_histo->GetBinContent(ih+1, jh+1)
500	<	* XbinSize * YbinSize);
489	>	double XbinSize = (hist_orig->GetXaxis()->GetXmax()
490	>	- hist_orig->GetXaxis()->GetXmin())
491	>	/ hist_orig->GetXaxis()->GetNbins();
492	>	double YbinSize = (hist_orig->GetYaxis()->GetXmax()
493	>	- hist_orig->GetYaxis()->GetXmin())
494	>	/ hist_orig->GetYaxis()->GetNbins();
495	>	for (int i = 1; i <= fit_histo->GetXaxis()->GetNbins(); i++) {
496	>	for (int j = 1; j <= fit_histo->GetYaxis()->GetNbins(); j++) {
497	>	double x = static_cast<double>(i - 1)
498	>	* (fit_histo->GetXaxis()->GetXmax()
499	>	- fit_histo->GetXaxis()->GetXmin())
500	>	/ static_cast<double>(fit_histo->GetXaxis()->GetNbins())
501	>	+ fit_histo->GetXaxis()->GetXmin();;
502	>	double y = static_cast<double>(j - 1)
503	>	* (fit_histo->GetYaxis()->GetXmax()
504	>	- fit_histo->GetYaxis()->GetXmin())
505	>	/ static_cast<double>(fit_histo->GetYaxis()->GetNbins())
506	>	+ fit_histo->GetYaxis()->GetXmin();
507	>	fit_histo->SetBinContent(i, j, tf2->Eval(x, y) * XbinSize * YbinSize);
508		}
509		}
510		}

Diff Legend

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