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.2 by dnisson, Thu Sep 3 18:00:45 2009 UTC vs.
Revision 1.11 by dnisson, Fri Sep 11 20:03:24 2009 UTC

#	Line 1 | Line 1
1	–	// TODO find why energies suddenly bigger
1		/* A simple jet-finding analyzer */
2
3		#include "UserCode/JetFitAnalyzer/interface/JetFitAnalyzer.h"
#	Line 8 | Line 7
7		#include "FWCore/MessageLogger/interface/MessageLogger.h"
8		#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
9
10	<	#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
11	<	#include "DataFormats/Candidate/interface/Particle.h"
10	>	#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
11	>	#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
12		#include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
13		#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
14		#include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
15
16	+	/* Jet reco stuff */
17	+	#include "DataFormats/JetReco/interface/PFJetCollection.h"
18	+	#include "DataFormats/JetReco/interface/GenJetCollection.h"
19	+	#include "DataFormats/JetReco/interface/PFJet.h"
20	+	#include "DataFormats/JetReco/interface/GenJet.h"
21	+
22		#include <map>
23		#include <vector>
24		#include <limits>
#	Line 30 | Line 35
35		using namespace std;
36		using namespace fastjet;
37
33	–	// Class to represent a "generic" particle, whether raw or reconstructed
34	–	class GenericParticle {
35	–	public:
36	–	GenericParticle(double __px, double __py, double __pz, double __e,
37	–	double __charge = 0.0)
38	–	: _px(__px), _py(__py), _pz(__pz), _e(__e), _charge(__charge) {
39	–	}
40	–	GenericParticle(const HepMC::GenParticle &genParticle)
41	–	: _charge(0.0) {
42	–	_px = genParticle.momentum().px();
43	–	_py = genParticle.momentum().py();
44	–	_pz = genParticle.momentum().pz();
45	–	_e = genParticle.momentum().e();
46	–	}
47	–	GenericParticle(const reco::PFCandidate &pfCandidate)
48	–	: _charge(0.0) {
49	–	_px = pfCandidate.px();
50	–	_py = pfCandidate.py();
51	–	_pz = pfCandidate.pz();
52	–	_e = pfCandidate.energy();
53	–	}
54	–	double px() {
55	–	return _px;
56	–	}
57	–	double py() {
58	–	return _py;
59	–	}
60	–	double pz() {
61	–	return _pz;
62	–	}
63	–	double e() {
64	–	return _e;
65	–	}
66	–	double charge() {
67	–	return _charge;
68	–	}
69	–	double eta() {
70	–	double theta = acos(_pz/sqrt(_px*_px + _py*_py + _pz*_pz));
71	–	return -log(tan(theta*0.5));
72	–	}
73	–	double phi() {
74	–	double phi0 = acos(_px/sqrt(_px*_px + _py*_py));
75	–	return _py > 0.0 ? phi0 : -phi0;
76	–	}
77	–	double m() {
78	–	return sqrt(_e*_e - _px*_px - _py*_py - _pz*_pz);
79	–	}
80	–	private:
81	–	double _px;
82	–	double _py;
83	–	double _pz;
84	–	double _e;
85	–	double _charge;
86	–	};
87	–
38		class JetFinderAnalyzer : public JetFitAnalyzer {
39		public:
40		struct jet {
#	Line 165 | Line 115 | private:
115		TH2 *);
116		virtual void analyze_results(jetfit::results r,
117		std::vector<jetfit::trouble> t, TH2 *);
118	<	vector<GenericParticle *> get_particles(const edm::Event&);
118	>	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;
123
124		fstream ofs;
125		edm::InputTag inputTagPFCandidates_;
126	+	edm::InputTag inputTagGenParticles_;
127		int info_type_;
128		double smear_;
129		int smear_coord_;
130	+	string jet_algo_;
131		};
132
133		map<TH2 *, vector< vector< JetFinderAnalyzer::jet > > >
#	Line 184 | Line 138 | JetFinderAnalyzer::JetFinderAnalyzer(con
138		{
139		info_type_ = pSet.getUntrackedParameter("info_type", 0);
140
141	+	if (info_type_ == 0) {
142	+	inputTagGenParticles_ = pSet.getParameter<edm::InputTag>("GenParticles");
143	+	}
144		if (info_type_ == 1) {
145		inputTagPFCandidates_ = pSet.getParameter<edm::InputTag>("PFCandidates");
146		}
#	Line 192 | Line 149 | JetFinderAnalyzer::JetFinderAnalyzer(con
149		smear_coord_ = pSet.getUntrackedParameter("smear_coord", 0);
150		// 0 = eta-phi smear
151		// 1 = proper angle smear
152	+	jet_algo_ = pSet.getParameter<string>("jet_algo");
153		set_user_minuit(jetfinder);
154		}
155
#	Line 212 | Line 170 | JetFinderAnalyzer::fetchCandidateCollect
170
171		}
172
173	<	vector<GenericParticle *> JetFinderAnalyzer::get_particles(const edm::Event &evt) {
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<edm::HepMCProduct> hRaw;
192	>	edm::Handle<reco::GenParticleCollection> hRaw;
193		edm::Handle<reco::PFCandidateCollection> hPFlow;
194		if (info_type_ == 0) {
195	<	evt.getByLabel("source", hRaw);
195	>	fetchCandidateCollection(hRaw,
196	>	inputTagGenParticles_,
197	>	evt);
198		}
199		if (info_type_ == 1) {
200		fetchCandidateCollection(hPFlow,
#	Line 225 | Line 202 | vector<GenericParticle *> JetFinderAnaly
202		evt);
203		}
204
205	<	vector<GenericParticle *> particles;
205	>	vector<reco::Candidate *> particles;
206
207		switch (info_type_) {
208		case 0:
209	<	const HepMC::GenEvent *hmcEvt = hRaw->GetEvent();
210	<	for (HepMC::GenEvent::particle_const_iterator
211	<	pit = hmcEvt->particles_begin(); pit != hmcEvt->particles_end();
235	<	pit++) {
236	<	if ((*pit)->status() == 1) {
237	<	particles.push_back(new GenericParticle(**pit));
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		}
240	–
214		break;
215		case 1:
216		for (unsigned i = 0; i < hPFlow->size(); i++) {
217	<	particles.push_back(new GenericParticle((*hPFlow)[i]));
217	>	particles.push_back((reco::Candidate *)&((*hPFlow)[i]));
218		}
219		break;
220		default:
#	Line 257 | Line 230 | TH2D * JetFinderAnalyzer::make_histo(con
230		TH2D *unred_histo = new TH2D(oss.str().c_str(), oss.str().c_str(),
231		600, -2.5, 2.5, 600, -PI, PI);
232
233	<	vector<GenericParticle *> particles = get_particles(evt);
233	>	vector<reco::Candidate *> particles = get_particles(evt);
234		for (int i = 0; i < particles.size(); i++) {
235		unred_histo->Fill(particles[i]->eta(),
236		particles[i]->phi(),
237	<	particles[i]->e());
237	>	particles[i]->energy());
238		}
239
240		// reduce histo
#	Line 303 | Line 276 | TH2D * JetFinderAnalyzer::make_histo(con
276		60, max_x-rcone, max_x+rcone,
277		60, max_y-rcone, max_y+rcone);
278		for (int i = 0; i < particles.size(); i++) {
279	<	double N = particles[i]->e();
279	>	double N = particles[i]->energy();
280		double x = particles[i]->eta();
281		double y = particles[i]->phi();
282		histo_unsmeared->Fill(x, y, N);
#	Line 321 | Line 294 | TH2D * JetFinderAnalyzer::make_histo(con
294		switch (smear_coord_) {
295		case 1:
296		for (int i = 0; i < particles.size(); i++) {
297	<	double N = particles[i]->e();
297	>	double N = particles[i]->energy();
298		double x = particles[i]->eta();
299		double y = particles[i]->phi();
300		// loop over bins and add Gaussian in proper angle to smeared
#	Line 329 | Line 302 | TH2D * JetFinderAnalyzer::make_histo(con
302		for (vector< double >::size_type j2 = 0; j2 < 60; j2++) {
303		double eta = static_cast<double>((signed int)i2) * XbinSize +
304		max_x - rcone - x;
305	<	double phi = acos(cos(static_cast<double>((signed int)j2) * YbinSize +
306	<	max_y - rcone - y));
307	<	phi = sin(phi) > 0 ? phi : -phi;
305	>	double phi0 = static_cast<double>((signed int)j2) * YbinSize +
306	>	max_y - rcone - y;
307	>	double phi = acos(cos(phi0));
308	>	phi = sin(phi0) > 0 ? phi : -phi;
309
310		// transform eta, phi to proper angle
311		double theta = 2.0*atan(exp(-eta));
#	Line 339 | Line 313 | TH2D * JetFinderAnalyzer::make_histo(con
313
314		smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
315		* exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_));
316	+
317	+	// correct for out-of-range phi
318	+	double transform = 0.0;
319	+	if (histo->GetYaxis()->GetXmin() < -PI) {
320	+	transform = -2.0*PI;
321	+	phi += transform;
322	+	iota = asin(sin(theta)*sin(phi));
323	+	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
324	+	* exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_));
325	+	}
326	+	if (histo->GetYaxis()->GetXmax() > PI) {
327	+	transform = 2.0*PI;
328	+	phi += transform;
329	+	iota = asin(sin(theta)*sin(phi));
330	+	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
331	+	* exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_));
332	+	}
333		}
334		}
335		}
#	Line 346 | Line 337 | TH2D * JetFinderAnalyzer::make_histo(con
337		case 0:
338		default:
339		for (int i = 0; i < particles.size(); i++) {
340	<	double N = particles[i]->e();
340	>	double N = particles[i]->energy();
341		double x = particles[i]->eta();
342		double y = particles[i]->phi();
343		// loop over bins and add Gaussian to smeared
#	Line 354 | Line 345 | TH2D * JetFinderAnalyzer::make_histo(con
345		for (vector< double >::size_type j2 = 0; j2 < 60; j2++) {
346		double eta = static_cast<double>((signed int)i2) * XbinSize
347		+ max_x - rcone - x;
348	<	double phi = acos(cos(static_cast<double>((signed int)j2) * YbinSize
349	<	+ max_y - rcone - y));
350	<	phi = sin(phi) > 0 ? phi : -phi;
348	>	double phi0 = static_cast<double>((signed int)j2) * YbinSize +
349	>	max_y - rcone - y;
350	>	double phi = acos(cos(phi0));
351	>	phi = sin(phi0) > 0 ? phi : -phi;
352		smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
353		* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_));
354	+
355	+	// correct for out-of-range phi
356	+	double transform = 0.0;
357	+	if (histo->GetYaxis()->GetXmin() < -PI) {
358	+	transform = -2.0*PI;
359	+	phi += transform;
360	+	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
361	+	* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_));
362	+	}
363	+	if (histo->GetYaxis()->GetXmax() > PI) {
364	+	transform = 2.0*PI;
365	+	phi += transform;
366	+	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
367	+	* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_));
368	+	}
369		}
370		}
371		}
#	Line 373 | Line 380 | TH2D * JetFinderAnalyzer::make_histo(con
380		return histo;
381		}
382
383	<	void seed_with_CA(vector<GenericParticle *> gParticles, TH2 *histo,
383	>	void seed_with_jetcoll(vector<reco::Jet> jets,
384		jetfit::model_def &_mdef) {
385	<	// create a PseudoJet vector
379	<	vector<PseudoJet> particles;
380	<	for (unsigned i = 0; i < gParticles.size(); i++) {
381	<	double x_max = (histo->GetXaxis()->GetXmax()
382	<	+ histo->GetXaxis()->GetXmin()) / 2.0;
383	<	double y_max = (histo->GetYaxis()->GetXmax()
384	<	+ histo->GetYaxis()->GetXmin()) / 2.0;
385	<	valarray<double> pmom(4);
386	<	pmom[0] = gParticles[i]->px();
387	<	pmom[1] = gParticles[i]->py();
388	<	pmom[2] = gParticles[i]->pz();
389	<	pmom[3] = gParticles[i]->e();
390	<	double eta = gParticles[i]->eta();
391	<	double phi = gParticles[i]->phi();
392	<	if ((eta - x_max)*(eta - x_max) + (phi - y_max)*(phi - y_max) < 0.25) {
393	<	PseudoJet j(pmom);
394	<	particles.push_back(j);
395	<	}
396	<	}
397	<
398	<	// choose a jet definition
399	<	double R = 0.2;
400	<	JetDefinition jet_def(cambridge_algorithm, R);
401	<
402	<	// run clustering and extract the jets
403	<	ClusterSequence cs(particles, jet_def);
404	<	vector<PseudoJet> jets = cs.inclusive_jets();
405	<
406	<	double XbinSize = (histo->GetXaxis()->GetXmax()
407	<	- histo->GetXaxis()->GetXmin()) /
408	<	static_cast<double>(histo->GetXaxis()->GetNbins());
409	<	double YbinSize = (histo->GetYaxis()->GetXmax()
410	<	- histo->GetYaxis()->GetXmin()) /
411	<	static_cast<double>(histo->GetYaxis()->GetNbins());
412	<
413	<	// seed with C-A jets
385	>	// seed with jet collection
386		int ijset = 0;
387		for (unsigned ij = 0; ij < jets.size(); ij++) {
388	<	double N = jets[ij].e();
389	<	if (N > 50.0) {
388	>	double N = jets[ij].energy();
389	>	double eta = jets[ij].eta();
390	>	double phi = jets[ij].phi();
391	>	if (N > 10.0) {
392		_mdef.set_special_par(ijset, 0, N, _mdef.chisquare_error(N)*0.1,
393		0.0, 1.0e6);
394	<	_mdef.set_special_par(ijset, 1, jets[ij].eta(), 0.01,
394	>	_mdef.set_special_par(ijset, 1, eta, 0.01,
395		0.0, 0.0);
396	<	double mdef_phi = jets[ij].phi() > PI ? jets[ij].phi() - 2*PI
397	<	: jets[ij].phi();
396	>	double mdef_phi = phi > PI ? phi - 2*PI
397	>	: phi;
398		_mdef.set_special_par(ijset, 2, mdef_phi, 0.01,
399		0.0, 0.0);
400		_mdef.set_special_par(ijset, 3, 0.1, 0.001, 0.0, 0.0);
#	Line 435 | Line 409 | jetfit::model_def& JetFinderAnalyzer::ma
409		class jf_model_def : public jetfit::model_def {
410		public:
411		virtual double chisquare_error(double E) {
412	<	return 0.97*E + 14.0;
413	<	// study from 08-27-09
412	>	return exp(-(4.2 + 0.11*E));
413	>	// study from 09-04-09
414		}
415		};
416
#	Line 452 | Line 426 | jetfit::model_def& JetFinderAnalyzer::ma
426		_mdef->set_indiv_par(2, string("mu_y"), 0.0, 0.0, 0.0, 0.0);
427		_mdef->set_indiv_par(3, string("sig"), 0.1, 0.001, 0.0, 0.0);
428
429	<	seed_with_CA(get_particles(evt), histo, *_mdef);
429	>	// get jetcoll from event file
430	>	edm::Handle< vector<reco::Jet> > jet_collection;
431	>	evt.getByLabel(jet_algo_, jet_collection);
432	>	seed_with_jetcoll(*jet_collection, *_mdef);
433
434		jetfit::set_model_def(_mdef);
435
#	Line 499 | Line 476 | jetfit::model_def& JetFinderAnalyzer::ma
476		}
477		jetfit::set_ngauss(_mdef->get_n_special_par_sets());
478		init_fit_histo->SetBinContent(i+1, j+1,
479	<	jetfit::fit_fcn(x, y, pval));
479	>	jetfit::fit_fcn(x, y, pval)
480	>	* XbinSize * YbinSize);
481		}
482		}
483
#	Line 552 | Line 530 | ostream& operator<<(ostream &out, jetfit
530		else
531		out << "Not calculated" << endl;
532		}
533	+	else {
534	+	out << "Error matrix accurate" << endl;
535	+	}
536
537		return out;
538		}
#	Line 560 | Line 541 | void JetFinderAnalyzer::analyze_results(
541		std::vector<jetfit::trouble> t,
542		TH2 *hist_orig) {
543		ofs << "Histogram "<<hist_orig->GetName() << endl;
544	<	for (int i = 0; i < unique_jets[hist_orig].size(); i++) {
544	>	for (int i = unique_jets[hist_orig].size() - r.chisquare.size();
545	>	i < unique_jets[hist_orig].size(); i++) {
546	>	int ir = i - unique_jets[hist_orig].size() + r.chisquare.size();
547		ofs << "For "<<i+1<<" gaussians: " << endl
548	<	<< t.at(i) << endl
549	<	<< unique_jets[hist_orig][i].size()<<" unique jets found" << endl;
548	>	<< t.at(i) << endl;
549	>	ofs << "chisquare="<<r.chisquare.at(ir)
550	>	<< endl;
551	>	ofs << unique_jets[hist_orig][i].size()<<" unique jets found" << endl;
552		for (int j = 0; j < unique_jets[hist_orig][i].size(); j++) {
553		jet _jet = unique_jets[hist_orig][i][j];
554		ofs << "Jet "<<j<<": Energy = "<<_jet.energy<<", eta = "<<_jet.eta

Diff Legend

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