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.24 by dnisson, Tue Dec 15 03:33:35 2009 UTC

#	Line 1 | Line 1
1	<	// TODO find why energies suddenly bigger
2	<	/* A simple jet-finding analyzer */
1	>	/* A JetFitAnalyzer that makes histograms with smearing */
2
3		#include "UserCode/JetFitAnalyzer/interface/JetFitAnalyzer.h"
4
#	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>
25		#include <cmath>
26		#include <cstdlib>
27	+	#include <iostream>
28		#include <fstream>
29		#include <sstream>
30
#	Line 30 | Line 36
36		using namespace std;
37		using namespace fastjet;
38
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	–
39		class JetFinderAnalyzer : public JetFitAnalyzer {
40		public:
90	–	struct jet {
91	–	double energy;
92	–	double eta;
93	–	double phi;
94	–	};
95	–
41		explicit JetFinderAnalyzer( const edm::ParameterSet&);
42		~JetFinderAnalyzer() {}
43
44		private:
45	<	static map<TH2 *, vector< vector<jet> > > unique_jets;
46	<
47	<	static double phi_cutoff_;
103	<
104	<	static double g2int(double xlo, double xhi, double ylo, double yhi,
105	<	double *pval) {
106	<	double sum1 = 0.0;
107	<	double sum2 = 0.0;
108	<	double xmid = 0.5 * (xlo + xhi);
109	<	double ymid = 0.5 * (ylo + yhi);
110	<	double xstep = (xhi - xlo) / 50.0;
111	<	double ystep = (yhi - ylo) / 50.0;
112	<	for (int i = 0; i < 50; i++) {
113	<	double x = (static_cast<double>(i) + 0.5) * xstep + xlo;
114	<	sum1 += xstep * jetfit::fit_fcn(x, ymid, pval);
115	<	}
116	<	for (int i = 0; i < 50; i++) {
117	<	double y = (static_cast<double>(i) + 0.5) * ystep + ylo;
118	<	sum2 += ystep * jetfit::fit_fcn(xmid, y, pval);
119	<	}
120	<	return sum1 * sum2;
121	<	}
122	<
123	<	static void jetfinder(TMinuit *gMinuit, TH2 *hist, int ngauss) {
124	<	double dist_sq = numeric_limits<double>::infinity();
125	<	unique_jets[hist].resize(ngauss);
126	<	int nbinsX = hist->GetXaxis()->GetNbins();
127	<	int nbinsY = hist->GetYaxis()->GetNbins();
128	<	double XbinSize = (hist->GetXaxis()->GetXmax()
129	<	- hist->GetXaxis()->GetXmin())
130	<	/ static_cast<double>(nbinsX);
131	<	double YbinSize = (hist->GetYaxis()->GetXmax()
132	<	- hist->GetYaxis()->GetXmin())
133	<	/ static_cast<double>(nbinsY);
134	<	for (int i = 0; i < ngauss; i++) {
135	<	double N, mu_x, mu_y, sig, err, lo, hi;
136	<	int iuint;
137	<	TString name;
138	<	gMinuit->mnpout(4*i, name, N, err, lo, hi, iuint);
139	<	gMinuit->mnpout(4*i + 1, name, mu_x, err, lo, hi, iuint);
140	<	gMinuit->mnpout(4*i + 2, name, mu_y, err, lo, hi, iuint);
141	<	gMinuit->mnpout(4*i + 3, name, sig, err, lo, hi, iuint);
142	<	for (int j = 0; j < i; j++) {
143	<	double N2, mu_x2, mu_y2, sig2;
144	<	gMinuit->mnpout(4*j, name, N2, err, lo, hi, iuint);
145	<	gMinuit->mnpout(4*j + 1, name, mu_x2, err, lo, hi, iuint);
146	<	gMinuit->mnpout(4*j + 2, name, mu_y2, err, lo, hi, iuint);
147	<	gMinuit->mnpout(4*j + 3, name, sig2, err, lo, hi, iuint);
148	<	double _dist_sq = (mu_x2 - mu_x)*(mu_x2 - mu_x)
149	<	+ (mu_y2 - mu_y)*(mu_y2 - mu_y);
150	<	if (_dist_sq < dist_sq)
151	<	dist_sq = _dist_sq;
152	<	}
153	<
154	<	jet j;
155	<	j.energy = N;
156	<	j.eta = mu_x; j.phi = mu_y;
157	<	unique_jets[hist][ngauss-1].push_back(j);
158	<	}
159	<	}
160	<
161	<	virtual void beginJob(const edm::EventSetup&);
162	<	virtual TH2D* make_histo(const edm::Event&, const edm::EventSetup&);
163	<	virtual jetfit::model_def& make_model_def(const edm::Event&,
45	>	virtual void beginJob(const edm::EventSetup &es);
46	>	virtual TH2D * make_histo(const edm::Event &, const edm::EventSetup &);
47	>	virtual HistoFitter::ModelDefinition& make_model_def(const edm::Event&,
48		const edm::EventSetup&,
49		TH2 *);
50	<	virtual void analyze_results(jetfit::results r,
51	<	std::vector<jetfit::trouble> t, TH2 *);
52	<	vector<GenericParticle *> get_particles(const edm::Event&);
169	<	void fetchCandidateCollection(edm::Handle<reco::PFCandidateCollection>&,
170	<	const edm::InputTag&, const edm::Event&) const;
50	>	virtual void analyze_results(HistoFitter::FitResults r,
51	>	std::vector<HistoFitter::Trouble> t, TH2 *);
52	>	vector<reco::Candidate *> get_particles(const edm::Event&);
53
54		fstream ofs;
55		edm::InputTag inputTagPFCandidates_;
56	+	edm::InputTag inputTagGenParticles_;
57		int info_type_;
58		double smear_;
59		int smear_coord_;
60	+	string jet_algo_;
61	+	string reclustered_jet_algo_;
62		};
63
179	–	map<TH2 *, vector< vector< JetFinderAnalyzer::jet > > >
180	–	JetFinderAnalyzer::unique_jets;
181	–
64		JetFinderAnalyzer::JetFinderAnalyzer(const edm::ParameterSet &pSet)
65		: JetFitAnalyzer(pSet) // this is important!
66		{
67		info_type_ = pSet.getUntrackedParameter("info_type", 0);
68
69	+	if (info_type_ == 0) {
70	+	inputTagGenParticles_ = pSet.getParameter<edm::InputTag>("GenParticles");
71	+	}
72		if (info_type_ == 1) {
73		inputTagPFCandidates_ = pSet.getParameter<edm::InputTag>("PFCandidates");
74		}
#	Line 192 | Line 77 | JetFinderAnalyzer::JetFinderAnalyzer(con
77		smear_coord_ = pSet.getUntrackedParameter("smear_coord", 0);
78		// 0 = eta-phi smear
79		// 1 = proper angle smear
80	<	set_user_minuit(jetfinder);
80	>	jet_algo_ = pSet.getParameter<string>("jet_algo");
81		}
82
83	<	void
84	<	JetFinderAnalyzer::fetchCandidateCollection(edm::Handle<reco::PFCandidateCollection>& c,
200	<	const edm::InputTag& tag,
201	<	const edm::Event& iEvent) const {
202	<
203	<	bool found = iEvent.getByLabel(tag, c);
204	<
205	<	if(!found ) {
206	<	ostringstream  err;
207	<	err<<" cannot get PFCandidates: "
208	<	<<tag<<endl;
209	<	edm::LogError("PFCandidates")<<err.str();
210	<	throw cms::Exception( "MissingProduct", err.str());
211	<	}
212	<
213	<	}
83	>	TH2D * JetFinderAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&) {
84	>	const reco::Jet * highest_e_jet = 0;
85
215	–	vector<GenericParticle *> JetFinderAnalyzer::get_particles(const edm::Event &evt) {
216	–	// fill unreduced histo
217	–	edm::Handle<edm::HepMCProduct> hRaw;
218	–	edm::Handle<reco::PFCandidateCollection> hPFlow;
86		if (info_type_ == 0) {
87	<	evt.getByLabel("source", hRaw);
88	<	}
89	<	if (info_type_ == 1) {
90	<	fetchCandidateCollection(hPFlow,
91	<	inputTagPFCandidates_,
92	<	evt);
87	>	edm::Handle< vector<reco::GenJet> > evtJets;
88	>	evt.getByLabel(jet_algo_, evtJets);
89	>	for (unsigned i = 0; i < evtJets->size(); i++) {
90	>	if (highest_e_jet == 0
91	>	|| (*evtJets)[i].energy() > highest_e_jet->energy()) {
92	>	highest_e_jet = &((*evtJets)[i]);
93	>	}
94	>	}
95		}
96	<
97	<	vector<GenericParticle *> particles;
98	<
99	<	switch (info_type_) {
100	<	case 0:
101	<	const HepMC::GenEvent *hmcEvt = hRaw->GetEvent();
102	<	for (HepMC::GenEvent::particle_const_iterator
234	<	pit = hmcEvt->particles_begin(); pit != hmcEvt->particles_end();
235	<	pit++) {
236	<	if ((*pit)->status() == 1) {
237	<	particles.push_back(new GenericParticle(**pit));
96	>	else if (info_type_ == 1) {
97	>	edm::Handle< vector<reco::PFJet> > evtJets;
98	>	evt.getByLabel(jet_algo_, evtJets);
99	>	for (unsigned i = 0; i < evtJets->size(); i++) {
100	>	if (highest_e_jet == 0
101	>	|| (*evtJets)[i].energy() > highest_e_jet->energy()) {
102	>	highest_e_jet = &((*evtJets)[i]);
103		}
104		}
105	+	}
106
107	<	break;
108	<	case 1:
109	<	for (unsigned i = 0; i < hPFlow->size(); i++) {
244	<	particles.push_back(new GenericParticle((*hPFlow)[i]));
245	<	}
246	<	break;
247	<	default:
248	<	cerr << "Unknown event type" << endl; // TODO use MessageLogger
107	>	if (highest_e_jet == 0) {
108	>	cerr << "No fat jets found!" << endl;
109	>	return 0;
110		}
111
112	<	return particles;
113	<	}
114	<
254	<	TH2D * JetFinderAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&) {
112	>	vector<const reco::Candidate *> particles =
113	>	highest_e_jet->getJetConstituentsQuick();
114	>
115		ostringstream oss;
116	<	oss << "eta_phi_energy_unred"<<evt.id().event() << flush;
257	<	TH2D *unred_histo = new TH2D(oss.str().c_str(), oss.str().c_str(),
258	<	600, -2.5, 2.5, 600, -PI, PI);
259	<
260	<	vector<GenericParticle *> particles = get_particles(evt);
261	<	for (int i = 0; i < particles.size(); i++) {
262	<	unred_histo->Fill(particles[i]->eta(),
263	<	particles[i]->phi(),
264	<	particles[i]->e());
265	<	}
266	<
267	<	// reduce histo
268	<	ostringstream oss2;
269	<	oss2 << "eta_phi_energy_red"<<evt.id().event() << flush;
116	>	oss << "eta_phi_energy"<<evt.id().event() << flush;
117		edm::Service<TFileService> fs;
118	<	// draw cone of radius 0.5 around highest energy bin, reduce
119	<	double maxE = 0.0;
120	<	int max_i = 29, max_j = 29;
121	<	for (int i = 0; i < unred_histo->GetNbinsX(); i++) {
122	<	for (int j = 0; j < unred_histo->GetNbinsY(); j++) {
123	<	double E = unred_histo->GetBinContent(i+1, j+1);
124	<	if (E > maxE) {
125	<	maxE = E;
126	<	max_i = i;
127	<	max_j = j;
128	<	}
129	<	}
130	<	}
131	<
132	<	double rcone = 0.5;
133	<	double Xlo = unred_histo->GetXaxis()->GetXmin();
134	<	double Xhi = unred_histo->GetXaxis()->GetXmax();
135	<	double Ylo = unred_histo->GetYaxis()->GetXmin();
136	<	double Yhi = unred_histo->GetYaxis()->GetXmax();
137	<	double XbinSize = (Xhi - Xlo) /
138	<	static_cast<double>(unred_histo->GetXaxis()->GetNbins());
139	<	double YbinSize = (Yhi - Ylo) /
293	<	static_cast<double>(unred_histo->GetYaxis()->GetNbins());
294	<	double max_x = (static_cast<double>(max_i) + 0.5) * XbinSize + Xlo;
295	<	double max_y = (static_cast<double>(max_j) + 0.5) * YbinSize + Ylo;
296	<	TH2D *histo = fs->make<TH2D>(oss2.str().c_str(), oss2.str().c_str(),
297	<	60, max_x-rcone, max_x+rcone,
298	<	60, max_y-rcone, max_y+rcone);
299	<
300	<	// create an unsmeared reduced histo
301	<	TH2D *histo_unsmeared = fs->make<TH2D>((oss2.str()+"_unsmeared").c_str(),
302	<	(oss2.str()+"_unsmeared").c_str(),
303	<	60, max_x-rcone, max_x+rcone,
304	<	60, max_y-rcone, max_y+rcone);
305	<	for (int i = 0; i < particles.size(); i++) {
306	<	double N = particles[i]->e();
307	<	double x = particles[i]->eta();
308	<	double y = particles[i]->phi();
309	<	histo_unsmeared->Fill(x, y, N);
310	<	}
311	<
312	<	// create a smeared reduced histo
313	<	// create a temporary 2D vector for smeared energies
314	<	XbinSize = (histo->GetXaxis()->GetXmax()
315	<	- histo->GetXaxis()->GetXmin()) /
316	<	static_cast<double>(histo->GetXaxis()->GetNbins());
317	<	YbinSize = (histo->GetYaxis()->GetXmax()
318	<	- histo->GetYaxis()->GetXmin()) /
319	<	static_cast<double>(histo->GetYaxis()->GetNbins());
320	<	vector< vector<double> > smeared(60, vector<double>(60, 0.0) );
321	<	switch (smear_coord_) {
322	<	case 1:
118	>	TH2D *histo = fs->make<TH2D>(oss.str().c_str(), oss.str().c_str(),
119	>	30,
120	>	highest_e_jet->eta()-0.5,
121	>	highest_e_jet->eta()+0.5,
122	>	30,
123	>	highest_e_jet->phi()-0.5,
124	>	highest_e_jet->phi()+0.5);
125	>
126	>	makeSmearedHisto
127	>	if (smear_ > 0.0) {
128	>	// create a temporary 2D vector for smeared energies
129	>	double XbinSize = (histo->GetXaxis()->GetXmax()
130	>	- histo->GetXaxis()->GetXmin()) /
131	>	static_cast<double>(histo->GetXaxis()->GetNbins());
132	>	double YbinSize = (histo->GetYaxis()->GetXmax()
133	>	- histo->GetYaxis()->GetXmin()) /
134	>	static_cast<double>(histo->GetYaxis()->GetNbins());
135	>	double Xlo = histo->GetXaxis()->GetXmin();
136	>	double Xhi = histo->GetXaxis()->GetXmax();
137	>	double Ylo = histo->GetYaxis()->GetXmin();
138	>	double Yhi = histo->GetYaxis()->GetXmax();
139	>	vector< vector<double> > smeared(30, vector<double>(30, 0.0) );
140		for (int i = 0; i < particles.size(); i++) {
141	<	double N = particles[i]->e();
141	>	double N = particles[i]->energy();
142		double x = particles[i]->eta();
143		double y = particles[i]->phi();
144	+	if (y < Ylo) {
145	+	y += 2.0*PI;
146	+	}
147	+	if (y > Yhi) {
148	+	y -= 2.0*PI;
149	+	}
150		// loop over bins and add Gaussian in proper angle to smeared
151	<	for (vector< vector<double> >::size_type i2 = 0; i2 < 60; i2++) {
152	<	for (vector< double >::size_type j2 = 0; j2 < 60; j2++) {
151	>	for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) {
152	>	for (vector< double >::size_type j2 = 0; j2 < 30; j2++) {
153		double eta = static_cast<double>((signed int)i2) * XbinSize +
154	<	max_x - rcone - x;
155	<	double phi = acos(cos(static_cast<double>((signed int)j2) * YbinSize +
156	<	max_y - rcone - y));
334	<	phi = sin(phi) > 0 ? phi : -phi;
335	<
336	<	// transform eta, phi to proper angle
337	<	double theta = 2.0*atan(exp(-eta));
338	<	double iota = asin(sin(theta)*sin(phi));
154	>	Xlo - x;
155	>	double phi = static_cast<double>((signed int)j2) * YbinSize +
156	>	Ylo - y;
157
158	<	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
159	<	* exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_));
158	>	double theta, iota;
159	>	switch (smear_coord_) {
160	>	case 1:
161	>	// transform eta, phi to proper angle
162	>	theta = 2.0*atan(exp(-eta));
163	>	iota = asin(sin(theta)*sin(phi));
164	>
165	>	// add smeared Gaussian
166	>	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
167	>	* exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_));
168	>	break;
169	>	case 0:
170	>	default:
171	>	// eta-phi smear--just add Gaussian centered around particle'
172	>	// in eta and phi
173	>	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
174	>	* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_));
175	>	}
176		}
177		}
178	<	}
179	<	break;
180	<	case 0:
181	<	default:
348	<	for (int i = 0; i < particles.size(); i++) {
349	<	double N = particles[i]->e();
350	<	double x = particles[i]->eta();
351	<	double y = particles[i]->phi();
352	<	// loop over bins and add Gaussian to smeared
353	<	for (vector< vector<double> >::size_type i2 = 0; i2 < 60; i2++) {
354	<	for (vector< double >::size_type j2 = 0; j2 < 60; j2++) {
355	<	double eta = static_cast<double>((signed int)i2) * XbinSize
356	<	+ max_x - rcone - x;
357	<	double phi = acos(cos(static_cast<double>((signed int)j2) * YbinSize
358	<	+ max_y - rcone - y));
359	<	phi = sin(phi) > 0 ? phi : -phi;
360	<	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
361	<	* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_));
178	>	// set histogram to match smear vector
179	>	for (int i = 1; i <= 30; i++) {
180	>	for (int j = 1; j <= 30; j++) {
181	>	histo->SetBinContent(i, j, smeared[i-1][j-1]);
182		}
183		}
184	<	}
184	>	}
185		}
186	<	// set histogram to match smear vector
187	<	for (int i = 1; i <= 60; i++) {
188	<	for (int j = 1; j <= 60; j++) {
189	<	histo->SetBinContent(i, j, smeared[i-1][j-1]);
186	>	else {
187	>	// don't smear--just fill with particles
188	>	for (int i = 0; i < particles.size(); i++) {
189	>	histo->Fill(particles[i]->eta(), particles[i]->phi(),
190	>	particles[i]->energy());
191		}
192		}
193	<
193	>
194		return histo;
195		}
196	<
197	<	void seed_with_CA(vector<GenericParticle *> gParticles, TH2 *histo,
198	<	jetfit::model_def &_mdef) {
199	<	// create a PseudoJet vector
200	<	vector<PseudoJet> particles;
201	<	for (unsigned i = 0; i < gParticles.size(); i++) {
202	<	double x_max = (histo->GetXaxis()->GetXmax()
203	<	+ histo->GetXaxis()->GetXmin()) / 2.0;
204	<	double y_max = (histo->GetYaxis()->GetXmax()
205	<	+ histo->GetYaxis()->GetXmin()) / 2.0;
206	<	valarray<double> pmom(4);
207	<	pmom[0] = gParticles[i]->px();
208	<	pmom[1] = gParticles[i]->py();
209	<	pmom[2] = gParticles[i]->pz();
210	<	pmom[3] = gParticles[i]->e();
211	<	double eta = gParticles[i]->eta();
212	<	double phi = gParticles[i]->phi();
213	<	if ((eta - x_max)*(eta - x_max) + (phi - y_max)*(phi - y_max) < 0.25) {
214	<	PseudoJet j(pmom);
215	<	particles.push_back(j);
216	<	}
217	<	}
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
414	<	int ijset = 0;
415	<	for (unsigned ij = 0; ij < jets.size(); ij++) {
416	<	double N = jets[ij].e();
417	<	if (N > 50.0) {
418	<	_mdef.set_special_par(ijset, 0, N, _mdef.chisquare_error(N)*0.1,
419	<	0.0, 1.0e6);
420	<	_mdef.set_special_par(ijset, 1, jets[ij].eta(), 0.01,
421	<	0.0, 0.0);
422	<	double mdef_phi = jets[ij].phi() > PI ? jets[ij].phi() - 2*PI
423	<	: jets[ij].phi();
424	<	_mdef.set_special_par(ijset, 2, mdef_phi, 0.01,
425	<	0.0, 0.0);
426	<	_mdef.set_special_par(ijset, 3, 0.1, 0.001, 0.0, 0.0);
427	<	ijset++;
428	<	}
429	<	}
196	>
197	>	template <class Jet>
198	>	void seed_with_jetcoll(vector<Jet> jets,
199	>	HistoFitter::ModelDefinition &_mdef, double phi_lo = -PI,
200	>	double phi_hi = PI) {
201	>	// seed with jet collection
202	>	double N = jets[0].energy();
203	>	double eta = jets[0].eta();
204	>	double phi = jets[0].phi();
205	>	_mdef.setIndivParameter(0, string("N0"), N, _mdef.chisquareSigma(N)*0.1,
206	>	0.0, 1.0e6);
207	>	_mdef.setIndivParameter(1, string("mu_x0"), eta, 0.01,
208	>	0.0, 0.0);
209	>	if (phi < phi_lo) {
210	>	phi += 2.0*PI;
211	>	}
212	>	if (phi > phi_hi) {
213	>	phi -= 2.0*PI;
214	>	}
215	>	_mdef.setIndivParameter(2, string("mu_y0"), phi, 0.01,
216	>	0.0, 0.0);
217	>	_mdef.setIndivParameter(3, string("sig0"), 0.1, 0.001, 0.0, 0.0);
218		}
219
220	<	jetfit::model_def& JetFinderAnalyzer::make_model_def(const edm::Event& evt,
220	>	HistoFitter::ModelDefinition& JetFinderAnalyzer::make_model_def(const edm::Event& evt,
221		const edm::EventSetup&,
222		TH2 *histo) {
223	<	class jf_model_def : public jetfit::model_def {
223	>	class jf_model_def : public HistoFitter::ModelDefinition {
224		public:
225	<	virtual double chisquare_error(double E) {
226	<	return 0.97*E + 14.0;
439	<	// study from 08-27-09
225	>	virtual double chisquareSigma(double E) {
226	>	return 0.298 - 0.9557*E; // study 11-09-09
227		}
228		};
229
230		jf_model_def *_mdef = new jf_model_def();
231		TFormula *formula = new TFormula("gaus2d",
232		"[0]*exp(-0.5*((x-[1])**2 + (y-[2])**2)/([3]**2))/(2*pi*[3]**2)");
233	<	_mdef->set_formula(formula);
234	<	_mdef->set_indiv_max_E(0);
235	<	_mdef->set_indiv_max_x(1);
236	<	_mdef->set_indiv_max_y(2);
237	<	_mdef->set_indiv_par(0, string("N"), 0.0, 0.0, 0.0, 1.0e6);
238	<	_mdef->set_indiv_par(1, string("mu_x"), 0.0, 0.0, 0.0, 0.0);
239	<	_mdef->set_indiv_par(2, string("mu_y"), 0.0, 0.0, 0.0, 0.0);
240	<	_mdef->set_indiv_par(3, string("sig"), 0.1, 0.001, 0.0, 0.0);
233	>	_mdef->setFormula(formula);
234	>	_mdef->setIndivEnergy(0);
235	>	_mdef->setIndivMeanX(1);
236	>	_mdef->setIndivMeanY(2);
237	>	_mdef->setIndivParameter(0, string("N"), 0.0, 0.0, 0.0, 1.0e6);
238	>	_mdef->setIndivParameter(1, string("mu_x"), 0.0, 0.0, 0.0, 0.0);
239	>	_mdef->setIndivParameter(2, string("mu_y"), 0.0, 0.0, 0.0, 0.0);
240	>	_mdef->setIndivParameter(3, string("sig"), 0.1, 0.001, 0.0, 0.0);
241
242	<	seed_with_CA(get_particles(evt), histo, *_mdef);
242	>	_mdef->setMaxNumberOfGaussians(3);
243
244	<	jetfit::set_model_def(_mdef);
245	<
246	<	// generate initial fit histogram
247	<	edm::Service<TFileService> fs;
248	<	TH2D *init_fit_histo = fs->make<TH2D>(("init_fit_"+string(histo->GetName()))
249	<	.c_str(),
250	<	("Initial fit for "
251	<	+string(histo->GetName())).c_str(),
252	<	histo->GetXaxis()->GetNbins(),
253	<	histo->GetXaxis()->GetXmin(),
254	<	histo->GetXaxis()->GetXmax(),
255	<	histo->GetXaxis()->GetNbins(),
256	<	histo->GetXaxis()->GetXmin(),
257	<	histo->GetXaxis()->GetXmax());
258	<	double XbinSize = (histo->GetXaxis()->GetXmax()
259	<	- histo->GetXaxis()->GetXmin()) /
260	<	static_cast<double>(histo->GetXaxis()->GetNbins());
261	<	double YbinSize = (histo->GetYaxis()->GetXmax()
262	<	- histo->GetYaxis()->GetXmin()) /
263	<	static_cast<double>(histo->GetYaxis()->GetNbins());
264	<	double Xlo = histo->GetXaxis()->GetXmin();
265	<	double Xhi = histo->GetXaxis()->GetXmax();
266	<	double Ylo = histo->GetYaxis()->GetXmin();
267	<	double Yhi = histo->GetYaxis()->GetXmax();
268	<
269	<	for (int i = 0; i < 60; i++) {
270	<	for (int j = 0; j < 60; j++) {
271	<	double x = (static_cast<double>(i) + 0.5)*XbinSize + Xlo;
272	<	double y = (static_cast<double>(j) + 0.5)*YbinSize + Ylo;
486	<	double pval[256];
487	<	if (_mdef->get_n_special_par_sets() > 64) {
488	<	cerr << "Parameter overload" << endl;
489	<	return *_mdef;
490	<	}
491	<	else {
492	<	for (int is = 0; is < _mdef->get_n_special_par_sets(); is++) {
493	<	for (int ii = 0; ii < 4; ii++) {
494	<	double spval, sperr, splo, sphi;
495	<	_mdef->get_special_par(is, ii, spval, sperr, splo, sphi);
496	<	pval[4*is + ii] = spval;
497	<	}
244	>	// get jetcoll from event file and select highest e jet
245	>	if (info_type_ == 0) {
246	>	edm::Handle< vector<reco::GenJet> > jet_collection;
247	>	evt.getByLabel(jet_algo_, jet_collection);
248	>	reco::GenJet highest_e_jet;
249	>	bool found_jet = false;
250	>	for (unsigned i = 0; i < jet_collection->size(); i++) {
251	>	if (!found_jet || (*jet_collection)[i].energy() > highest_e_jet.energy())
252	>	{
253	>	highest_e_jet = (*jet_collection)[i];
254	>	found_jet = true;
255	>	}
256	>	}
257	>	vector<reco::GenJet> highest_e_jet_coll;
258	>	highest_e_jet_coll.push_back(highest_e_jet);
259	>	seed_with_jetcoll(highest_e_jet_coll, *_mdef,
260	>	histo->GetYaxis()->GetXmin(),
261	>	histo->GetYaxis()->GetXmax());
262	>	}
263	>	if (info_type_ == 1) {
264	>	edm::Handle< vector<reco::PFJet> > jet_collection;
265	>	evt.getByLabel(jet_algo_, jet_collection);
266	>	reco::PFJet highest_e_jet;
267	>	bool found_jet = false;
268	>	for (unsigned i = 0; i < jet_collection->size(); i++) {
269	>	if (!found_jet || (*jet_collection)[i].energy() > highest_e_jet.energy())
270	>	{
271	>	highest_e_jet = (*jet_collection)[i];
272	>	found_jet = true;
273		}
499	–	}
500	–	jetfit::set_ngauss(_mdef->get_n_special_par_sets());
501	–	init_fit_histo->SetBinContent(i+1, j+1,
502	–	jetfit::fit_fcn(x, y, pval));
274		}
275	+	vector<reco::PFJet> highest_e_jet_coll;
276	+	highest_e_jet_coll.push_back(highest_e_jet);
277	+	seed_with_jetcoll(highest_e_jet_coll, *_mdef,
278	+	histo->GetYaxis()->GetXmin(),
279	+	histo->GetYaxis()->GetXmax());
280		}
281
282		return *_mdef;
283		}
284
285		void JetFinderAnalyzer::beginJob(const edm::EventSetup &es) {
510	–	ofs.open("jetfindlog.txt", ios::out);
511	–	if (ofs.fail()) {
512	–	cerr << "Opening jetfindlog.txt FAILED" << endl;
513	–	}
514	–	ofs << "Jetfinder log" << endl
515	–	<< "=============" << endl << endl;
516	–	}
517	–
518	–	ostream& operator<<(ostream &out, jetfit::trouble t) {
519	–	string action, error_string;
520	–
521	–	if (t.istat != 3) {
522	–	switch(t.occ) {
523	–	case jetfit::T_NULL:
524	–	action = "Program"; break;
525	–	case jetfit::T_SIMPLEX:
526	–	action = "SIMPLEX"; break;
527	–	case jetfit::T_MIGRAD:
528	–	action = "MIGRAD"; break;
529	–	case jetfit::T_MINOS:
530	–	action = "MINOS"; break;
531	–	default:
532	–	action = "Program"; break;
533	–	}
534	–
535	–	switch (t.istat) {
536	–	case 0:
537	–	error_string = "Unable to calculate error matrix"; break;
538	–	case 1:
539	–	error_string = "Error matrix a diagonal approximation"; break;
540	–	case 2:
541	–	error_string = "Error matrix not positive definite"; break;
542	–	case 3:
543	–	error_string = "Converged successfully"; break;
544	–	default:
545	–	ostringstream oss;
546	–	oss<<"Unknown status code "<<t.istat << flush;
547	–	error_string = oss.str(); break;
548	–	}
286
550	–	if (t.occ != jetfit::T_NULL)
551	–	out << action<<" trouble: "<<error_string;
552	–	else
553	–	out << "Not calculated" << endl;
554	–	}
555	–
556	–	return out;
287		}
288
289	<	void JetFinderAnalyzer::analyze_results(jetfit::results r,
290	<	std::vector<jetfit::trouble> t,
289	>	void JetFinderAnalyzer::analyze_results(HistoFitter::FitResults r,
290	>	std::vector<HistoFitter::Trouble> t,
291		TH2 *hist_orig) {
292	<	ofs << "Histogram "<<hist_orig->GetName() << endl;
563	<	for (int i = 0; i < unique_jets[hist_orig].size(); i++) {
564	<	ofs << "For "<<i+1<<" gaussians: " << endl
565	<	<< t.at(i) << endl
566	<	<< unique_jets[hist_orig][i].size()<<" unique jets found" << endl;
567	<	for (int j = 0; j < unique_jets[hist_orig][i].size(); j++) {
568	<	jet _jet = unique_jets[hist_orig][i][j];
569	<	ofs << "Jet "<<j<<": Energy = "<<_jet.energy<<", eta = "<<_jet.eta
570	<	<< ", phi = "<<_jet.phi << endl;
571	<	}
572	<	ofs << endl;
573	<	}
574	<	ofs << endl;
575	<
576	<	// save fit function histograms to root file
577	<	edm::Service<TFileService> fs;
578	<	for (vector< vector<double> >::size_type i = 0;
579	<	i < r.pval.size(); i++) {
580	<	jetfit::set_ngauss(r.pval[i].size() / 4);
581	<	TF2 *tf2 = new TF2("fit_func", jetfit::fit_fcn_TF2,
582	<	hist_orig->GetXaxis()->GetXmin(),
583	<	hist_orig->GetXaxis()->GetXmax(),
584	<	hist_orig->GetYaxis()->GetXmin(),
585	<	hist_orig->GetYaxis()->GetXmax(),
586	<	r.pval[i].size());
587	<	for (vector<double>::size_type j = 0; j < r.pval[i].size(); j++) {
588	<	tf2->SetParameter(j, r.pval[i][j]);
589	<	}
590	<	ostringstream fit_histo_oss;
591	<	fit_histo_oss << hist_orig->GetName()<<"_fit_"<<i << flush;
592	<	tf2->SetNpx(hist_orig->GetXaxis()->GetNbins());
593	<	tf2->SetNpy(hist_orig->GetYaxis()->GetNbins());
594	<	TH2D *fit_histo = fs->make<TH2D>(fit_histo_oss.str().c_str(),
595	<	fit_histo_oss.str().c_str(),
596	<	hist_orig->GetXaxis()->GetNbins(),
597	<	hist_orig->GetXaxis()->GetXmin(),
598	<	hist_orig->GetXaxis()->GetXmax(),
599	<	hist_orig->GetYaxis()->GetNbins(),
600	<	hist_orig->GetYaxis()->GetXmin(),
601	<	hist_orig->GetYaxis()->GetXmax());
602	<	TH1 *tf2_histo = tf2->CreateHistogram();
603	<	double XbinSize = (fit_histo->GetXaxis()->GetXmax()
604	<	- fit_histo->GetXaxis()->GetXmin())
605	<	/ static_cast<double>(fit_histo->GetXaxis()->GetNbins());
606	<	double YbinSize = (fit_histo->GetYaxis()->GetXmax()
607	<	- fit_histo->GetYaxis()->GetXmin())
608	<	/ static_cast<double>(fit_histo->GetYaxis()->GetNbins());
609	<	for (int ih = 0; ih < tf2->GetNpx(); ih++) {
610	<	for (int jh = 0; jh < tf2->GetNpy(); jh++) {
611	<	fit_histo->SetBinContent(ih+1, jh+1,
612	<	tf2_histo->GetBinContent(ih+1, jh+1)
613	<	* XbinSize * YbinSize);
614	<	}
615	<	}
616	<	}
617	<
618	<	// save results to file
619	<	ostringstream res_tree_oss, rt_title_oss;
620	<	res_tree_oss << hist_orig->GetName()<<"_results" << flush;
621	<	rt_title_oss << "Fit results for "<<hist_orig->GetName() << flush;
292	>	// TODO perform analysis of fit results
293		}
294
295		DEFINE_FWK_MODULE(JetFinderAnalyzer);

Diff Legend

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