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.5 by dnisson, Fri Sep 4 15:30:38 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"
#	Line 30 | Line 29
29		using namespace std;
30		using namespace fastjet;
31
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	–
32		class JetFinderAnalyzer : public JetFitAnalyzer {
33		public:
34		struct jet {
#	Line 165 | Line 109 | private:
109		TH2 *);
110		virtual void analyze_results(jetfit::results r,
111		std::vector<jetfit::trouble> t, TH2 *);
112	<	vector<GenericParticle *> get_particles(const edm::Event&);
112	>	vector<reco::Candidate *> get_particles(const edm::Event&);
113		void fetchCandidateCollection(edm::Handle<reco::PFCandidateCollection>&,
114		const edm::InputTag&, const edm::Event&) const;
115	+	void fetchCandidateCollection(edm::Handle<reco::GenParticleCollection>&,
116	+	const edm::InputTag&, const edm::Event&) const;
117
118		fstream ofs;
119		edm::InputTag inputTagPFCandidates_;
120	+	edm::InputTag inputTagGenParticles_;
121		int info_type_;
122		double smear_;
123		int smear_coord_;
#	Line 184 | Line 131 | JetFinderAnalyzer::JetFinderAnalyzer(con
131		{
132		info_type_ = pSet.getUntrackedParameter("info_type", 0);
133
134	+	if (info_type_ == 0) {
135	+	inputTagGenParticles_ = pSet.getParameter<edm::InputTag>("GenParticles");
136	+	}
137		if (info_type_ == 1) {
138		inputTagPFCandidates_ = pSet.getParameter<edm::InputTag>("PFCandidates");
139		}
#	Line 212 | Line 162 | JetFinderAnalyzer::fetchCandidateCollect
162
163		}
164
165	<	vector<GenericParticle *> JetFinderAnalyzer::get_particles(const edm::Event &evt) {
165	>	void
166	>	JetFinderAnalyzer::fetchCandidateCollection(edm::Handle<reco::GenParticleCollection>& c,
167	>	const edm::InputTag& tag,
168	>	const edm::Event& iEvent) const {
169	>
170	>	bool found = iEvent.getByLabel(tag, c);
171	>
172	>	if(!found ) {
173	>	ostringstream  err;
174	>	err<<" cannot get GenParticles: "
175	>	<<tag<<endl;
176	>	edm::LogError("GenParticles")<<err.str();
177	>	throw cms::Exception( "MissingProduct", err.str());
178	>	}
179	>
180	>	}
181	>
182	>	vector<reco::Candidate *> JetFinderAnalyzer::get_particles(const edm::Event &evt) {
183		// fill unreduced histo
184	<	edm::Handle<edm::HepMCProduct> hRaw;
184	>	edm::Handle<reco::GenParticleCollection> hRaw;
185		edm::Handle<reco::PFCandidateCollection> hPFlow;
186		if (info_type_ == 0) {
187	<	evt.getByLabel("source", hRaw);
187	>	fetchCandidateCollection(hRaw,
188	>	inputTagGenParticles_,
189	>	evt);
190		}
191		if (info_type_ == 1) {
192		fetchCandidateCollection(hPFlow,
#	Line 225 | Line 194 | vector<GenericParticle *> JetFinderAnaly
194		evt);
195		}
196
197	<	vector<GenericParticle *> particles;
197	>	vector<reco::Candidate *> particles;
198
199		switch (info_type_) {
200		case 0:
201	<	const HepMC::GenEvent *hmcEvt = hRaw->GetEvent();
202	<	for (HepMC::GenEvent::particle_const_iterator
203	<	pit = hmcEvt->particles_begin(); pit != hmcEvt->particles_end();
235	<	pit++) {
236	<	if ((*pit)->status() == 1) {
237	<	particles.push_back(new GenericParticle(**pit));
201	>	for (unsigned i = 0; i < hRaw->size(); i++) {
202	>	if ((*hRaw)[i].status() == 1) {
203	>	particles.push_back((reco::Candidate *)&((*hRaw)[i]));
204		}
205		}
240	–
206		break;
207		case 1:
208		for (unsigned i = 0; i < hPFlow->size(); i++) {
209	<	particles.push_back(new GenericParticle((*hPFlow)[i]));
209	>	particles.push_back((reco::Candidate *)&((*hPFlow)[i]));
210		}
211		break;
212		default:
#	Line 257 | Line 222 | TH2D * JetFinderAnalyzer::make_histo(con
222		TH2D *unred_histo = new TH2D(oss.str().c_str(), oss.str().c_str(),
223		600, -2.5, 2.5, 600, -PI, PI);
224
225	<	vector<GenericParticle *> particles = get_particles(evt);
225	>	vector<reco::Candidate *> particles = get_particles(evt);
226		for (int i = 0; i < particles.size(); i++) {
227		unred_histo->Fill(particles[i]->eta(),
228		particles[i]->phi(),
229	<	particles[i]->e());
229	>	particles[i]->energy());
230		}
231
232		// reduce histo
#	Line 303 | Line 268 | TH2D * JetFinderAnalyzer::make_histo(con
268		60, max_x-rcone, max_x+rcone,
269		60, max_y-rcone, max_y+rcone);
270		for (int i = 0; i < particles.size(); i++) {
271	<	double N = particles[i]->e();
271	>	double N = particles[i]->energy();
272		double x = particles[i]->eta();
273		double y = particles[i]->phi();
274		histo_unsmeared->Fill(x, y, N);
#	Line 321 | Line 286 | TH2D * JetFinderAnalyzer::make_histo(con
286		switch (smear_coord_) {
287		case 1:
288		for (int i = 0; i < particles.size(); i++) {
289	<	double N = particles[i]->e();
289	>	double N = particles[i]->energy();
290		double x = particles[i]->eta();
291		double y = particles[i]->phi();
292		// loop over bins and add Gaussian in proper angle to smeared
#	Line 346 | Line 311 | TH2D * JetFinderAnalyzer::make_histo(con
311		case 0:
312		default:
313		for (int i = 0; i < particles.size(); i++) {
314	<	double N = particles[i]->e();
314	>	double N = particles[i]->energy();
315		double x = particles[i]->eta();
316		double y = particles[i]->phi();
317		// loop over bins and add Gaussian to smeared
#	Line 373 | Line 338 | TH2D * JetFinderAnalyzer::make_histo(con
338		return histo;
339		}
340
341	<	void seed_with_CA(vector<GenericParticle *> gParticles, TH2 *histo,
341	>	void seed_with_CA(vector<reco::Candidate *> gParticles, TH2 *histo,
342		jetfit::model_def &_mdef) {
343		// create a PseudoJet vector
344		vector<PseudoJet> particles;
#	Line 386 | Line 351 | void seed_with_CA(vector<GenericParticle
351		pmom[0] = gParticles[i]->px();
352		pmom[1] = gParticles[i]->py();
353		pmom[2] = gParticles[i]->pz();
354	<	pmom[3] = gParticles[i]->e();
354	>	pmom[3] = gParticles[i]->energy();
355		double eta = gParticles[i]->eta();
356		double phi = gParticles[i]->phi();
357		if ((eta - x_max)*(eta - x_max) + (phi - y_max)*(phi - y_max) < 0.25) {
#	Line 415 | Line 380 | void seed_with_CA(vector<GenericParticle
380		for (unsigned ij = 0; ij < jets.size(); ij++) {
381		double N = jets[ij].e();
382		if (N > 50.0) {
383	+	cout << N << endl;
384		_mdef.set_special_par(ijset, 0, N, _mdef.chisquare_error(N)*0.1,
385		0.0, 1.0e6);
386		_mdef.set_special_par(ijset, 1, jets[ij].eta(), 0.01,
#	Line 435 | Line 401 | jetfit::model_def& JetFinderAnalyzer::ma
401		class jf_model_def : public jetfit::model_def {
402		public:
403		virtual double chisquare_error(double E) {
404	<	return 0.97*E + 14.0;
405	<	// study from 08-27-09
404	>	return E > 2.7 ? 0.55*E - 1.0 : 0.5;
405	>	// study from 09-03-09
406		}
407		};
408

Diff Legend

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