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.13 by dnisson, Wed Sep 23 03:13:07 2009 UTC vs.
Revision 1.23 by dnisson, Thu Nov 26 21:06:02 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 37 | Line 38 | using namespace fastjet;
38
39		class JetFinderAnalyzer : public JetFitAnalyzer {
40		public:
40	–	struct jet {
41	–	double energy;
42	–	double eta;
43	–	double phi;
44	–	};
45	–
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_;
53	<
54	<	static double g2int(double xlo, double xhi, double ylo, double yhi,
55	<	double *pval) {
56	<	double sum1 = 0.0;
57	<	double sum2 = 0.0;
58	<	double xmid = 0.5 * (xlo + xhi);
59	<	double ymid = 0.5 * (ylo + yhi);
60	<	double xstep = (xhi - xlo) / 50.0;
61	<	double ystep = (yhi - ylo) / 50.0;
62	<	for (int i = 0; i < 50; i++) {
63	<	double x = (static_cast<double>(i) + 0.5) * xstep + xlo;
64	<	sum1 += xstep * jetfit::fit_fcn(x, ymid, pval);
65	<	}
66	<	for (int i = 0; i < 50; i++) {
67	<	double y = (static_cast<double>(i) + 0.5) * ystep + ylo;
68	<	sum2 += ystep * jetfit::fit_fcn(xmid, y, pval);
69	<	}
70	<	return sum1 * sum2;
71	<	}
72	<
73	<	static void jetfinder(TMinuit *gMinuit, TH2 *hist, int ngauss) {
74	<	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&);
113	<	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 *);
50	>	virtual void analyze_results(HistoFitter::FitResults r,
51	>	std::vector<HistoFitter::Trouble> t, TH2 *);
52		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;
53
54		fstream ofs;
55		edm::InputTag inputTagPFCandidates_;
#	Line 128 | Line 58 | private:
58		double smear_;
59		int smear_coord_;
60		string jet_algo_;
61	+	string reclustered_jet_algo_;
62		};
63
133	–	map<TH2 *, vector< vector< JetFinderAnalyzer::jet > > >
134	–	JetFinderAnalyzer::unique_jets;
135	–
64		JetFinderAnalyzer::JetFinderAnalyzer(const edm::ParameterSet &pSet)
65		: JetFitAnalyzer(pSet) // this is important!
66		{
#	Line 150 | Line 78 | JetFinderAnalyzer::JetFinderAnalyzer(con
78		// 0 = eta-phi smear
79		// 1 = proper angle smear
80		jet_algo_ = pSet.getParameter<string>("jet_algo");
153	–	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;
81		}
82
83		TH2D * JetFinderAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&) {
84	<	const reco::Jet * highest_e_jet;
84	>	const reco::Jet * highest_e_jet = 0;
85
86		if (info_type_ == 0) {
87		edm::Handle< vector<reco::GenJet> > evtJets;
88		evt.getByLabel(jet_algo_, evtJets);
89		for (unsigned i = 0; i < evtJets->size(); i++) {
90	<	if (i == 0 || (*evtJets)[i].energy() > highest_e_jet->energy()) {
90	>	if (highest_e_jet == 0
91	>	|| (*evtJets)[i].energy() > highest_e_jet->energy()) {
92		highest_e_jet = &((*evtJets)[i]);
93		}
94		}
#	Line 240 | Line 97 | TH2D * JetFinderAnalyzer::make_histo(con
97		edm::Handle< vector<reco::PFJet> > evtJets;
98		evt.getByLabel(jet_algo_, evtJets);
99		for (unsigned i = 0; i < evtJets->size(); i++) {
100	<	if (i == 0 || (*evtJets)[i].energy() > highest_e_jet->energy()) {
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		if (highest_e_jet == 0) {
108	<	cerr << "No fatjets found!" << endl; exit(1);
108	>	cerr << "No fat jets found!" << endl;
109	>	return 0;
110		}
111
112		vector<const reco::Candidate *> particles =
113		highest_e_jet->getJetConstituentsQuick();
114	<	//  cout << "found highest e jet" << endl;
256	<
114	>
115		ostringstream oss;
116		oss << "eta_phi_energy"<<evt.id().event() << flush;
117	<	TH2D *histo = new TH2D(oss.str().c_str(), oss.str().c_str(),
118	<	600,
117	>	edm::Service<TFileService> fs;
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	<	600,
122	>	30,
123		highest_e_jet->phi()-0.5,
124		highest_e_jet->phi()+0.5);
125
126	<	for (int i = 0; i < particles.size(); i++) {
268	<	histo->Fill(particles[i]->eta(),
269	<	particles[i]->phi(),
270	<	particles[i]->energy());
271	<	}
272	<
126	>	if (smear_ > 0.0) {
127		// create a smeared histo
128		// create a temporary 2D vector for smeared energies
129		double XbinSize = (histo->GetXaxis()->GetXmax()
#	Line 279 | Line 133 | TH2D * JetFinderAnalyzer::make_histo(con
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	<	vector< vector<double> > smeared(60, vector<double>(60, 0.0) );
138	>	double Yhi = histo->GetYaxis()->GetXmax();
139	>	vector< vector<double> > smeared(30, vector<double>(30, 0.0) );
140		switch (smear_coord_) {
141		case 1:
142		for (int i = 0; i < particles.size(); i++) {
143		double N = particles[i]->energy();
144		double x = particles[i]->eta();
145		double y = particles[i]->phi();
146	+	if (y < Ylo) {
147	+	y += 2.0*PI;
148	+	}
149	+	if (y > Yhi) {
150	+	y -= 2.0*PI;
151	+	}
152		// loop over bins and add Gaussian in proper angle to smeared
153	<	for (vector< vector<double> >::size_type i2 = 0; i2 < 60; i2++) {
154	<	for (vector< double >::size_type j2 = 0; j2 < 60; j2++) {
153	>	for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) {
154	>	for (vector< double >::size_type j2 = 0; j2 < 30; j2++) {
155		double eta = static_cast<double>((signed int)i2) * XbinSize +
156	<	Xlo;
157	<	double phi0 = static_cast<double>((signed int)j2) * YbinSize +
158	<	Ylo;
159	<	double phi = acos(cos(phi0));
298	<	phi = sin(phi0) > 0 ? phi : -phi;
299	<
156	>	Xlo - x;
157	>	double phi = static_cast<double>((signed int)j2) * YbinSize +
158	>	Ylo - y;
159	>
160		// transform eta, phi to proper angle
161		double theta = 2.0*atan(exp(-eta));
162		double iota = asin(sin(theta)*sin(phi));
163
164		smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
165		* exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_));
306	–
307	–	// correct for out-of-range phi
308	–	double transform = 0.0;
309	–	if (histo->GetYaxis()->GetXmin() < -PI) {
310	–	transform = -2.0*PI;
311	–	phi += transform;
312	–	iota = asin(sin(theta)*sin(phi));
313	–	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
314	–	* exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_));
315	–	}
316	–	if (histo->GetYaxis()->GetXmax() > PI) {
317	–	transform = 2.0*PI;
318	–	phi += transform;
319	–	iota = asin(sin(theta)*sin(phi));
320	–	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
321	–	* exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_));
322	–	}
166		}
167		}
168		}
#	Line 330 | Line 173 | TH2D * JetFinderAnalyzer::make_histo(con
173		double N = particles[i]->energy();
174		double x = particles[i]->eta();
175		double y = particles[i]->phi();
176	+	if (y < Ylo) {
177	+	y += 2.0*PI;
178	+	}
179	+	if (y > Yhi) {
180	+	y -= 2.0*PI;
181	+	}
182		// loop over bins and add Gaussian to smeared
183	<	for (vector< vector<double> >::size_type i2 = 0; i2 < 60; i2++) {
184	<	for (vector< double >::size_type j2 = 0; j2 < 60; j2++) {
183	>	for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) {
184	>	for (vector< double >::size_type j2 = 0; j2 < 30; j2++) {
185		double eta = static_cast<double>((signed int)i2) * XbinSize
186	<	+ Xlo;
187	<	double phi0 = static_cast<double>((signed int)j2) * YbinSize +
188	<	Ylo;
340	<	double phi = acos(cos(phi0));
341	<	phi = sin(phi0) > 0 ? phi : -phi;
342	<	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
343	<	* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_));
186	>	+ Xlo - x;
187	>	double phi = static_cast<double>((signed int)j2) * YbinSize +
188	>	Ylo - y;
189
190	<	// correct for out-of-range phi
346	<	double transform = 0.0;
347	<	if (histo->GetYaxis()->GetXmin() < -PI) {
348	<	transform = -2.0*PI;
349	<	phi += transform;
350	<	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
351	<	* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_));
352	<	}
353	<	if (histo->GetYaxis()->GetXmax() > PI) {
354	<	transform = 2.0*PI;
355	<	phi += transform;
356	<	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
190	>	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
191		* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_));
358	–	}
192		}
193		}
194		}
195		}
196		// set histogram to match smear vector
197	<	for (int i = 1; i <= 60; i++) {
198	<	for (int j = 1; j <= 60; j++) {
197	>	for (int i = 1; i <= 30; i++) {
198	>	for (int j = 1; j <= 30; j++) {
199		histo->SetBinContent(i, j, smeared[i-1][j-1]);
200		}
201		}
202	+	}
203	+	else {
204	+	// don't smear--just fill with particles
205	+	for (int i = 0; i < particles.size(); i++) {
206	+	histo->Fill(particles[i]->eta(), particles[i]->phi(),
207	+	particles[i]->energy());
208	+	}
209	+	}
210
211		return histo;
212		}
213
214		template <class Jet>
215		void seed_with_jetcoll(vector<Jet> jets,
216	<	jetfit::model_def &_mdef) {
216	>	HistoFitter::ModelDefinition &_mdef, double phi_lo = -PI,
217	>	double phi_hi = PI) {
218		// seed with jet collection
219	<	int ijset = 0;
220	<	for (unsigned ij = 0; ij < jets.size(); ij++) {
221	<	double N = jets[ij].energy();
222	<	double eta = jets[ij].eta();
223	<	double phi = jets[ij].phi();
224	<	if (N > 10.0) {
225	<	_mdef.set_special_par(ijset, 0, N, _mdef.chisquare_error(N)*0.1,
226	<	0.0, 1.0e6);
227	<	_mdef.set_special_par(ijset, 1, eta, 0.01,
228	<	0.0, 0.0);
229	<	double mdef_phi = phi > PI ? phi - 2*PI
230	<	: phi;
231	<	_mdef.set_special_par(ijset, 2, mdef_phi, 0.01,
232	<	0.0, 0.0);
233	<	_mdef.set_special_par(ijset, 3, 0.1, 0.001, 0.0, 0.0);
234	<	ijset++;
393	<	}
394	<	}
219	>	double N = jets[0].energy();
220	>	double eta = jets[0].eta();
221	>	double phi = jets[0].phi();
222	>	_mdef.setIndivParameter(0, string("N0"), N, _mdef.chisquareSigma(N)*0.1,
223	>	0.0, 1.0e6);
224	>	_mdef.setIndivParameter(1, string("mu_x0"), eta, 0.01,
225	>	0.0, 0.0);
226	>	if (phi < phi_lo) {
227	>	phi += 2.0*PI;
228	>	}
229	>	if (phi > phi_hi) {
230	>	phi -= 2.0*PI;
231	>	}
232	>	_mdef.setIndivParameter(2, string("mu_y0"), phi, 0.01,
233	>	0.0, 0.0);
234	>	_mdef.setIndivParameter(3, string("sig0"), 0.1, 0.001, 0.0, 0.0);
235		}
236
237	<	jetfit::model_def& JetFinderAnalyzer::make_model_def(const edm::Event& evt,
237	>	HistoFitter::ModelDefinition& JetFinderAnalyzer::make_model_def(const edm::Event& evt,
238		const edm::EventSetup&,
239		TH2 *histo) {
240	<	class jf_model_def : public jetfit::model_def {
240	>	class jf_model_def : public HistoFitter::ModelDefinition {
241		public:
242	<	virtual double chisquare_error(double E) {
243	<	return exp(-(4.2 + 0.11*E));
404	<	// study from 09-04-09
242	>	virtual double chisquareSigma(double E) {
243	>	return 0.298 - 0.9557*E; // study 11-09-09
244		}
245		};
246
247		jf_model_def *_mdef = new jf_model_def();
248		TFormula *formula = new TFormula("gaus2d",
249		"[0]*exp(-0.5*((x-[1])**2 + (y-[2])**2)/([3]**2))/(2*pi*[3]**2)");
250	<	_mdef->set_formula(formula);
251	<	_mdef->set_indiv_max_E(0);
252	<	_mdef->set_indiv_max_x(1);
253	<	_mdef->set_indiv_max_y(2);
254	<	_mdef->set_indiv_par(0, string("N"), 0.0, 0.0, 0.0, 1.0e6);
255	<	_mdef->set_indiv_par(1, string("mu_x"), 0.0, 0.0, 0.0, 0.0);
256	<	_mdef->set_indiv_par(2, string("mu_y"), 0.0, 0.0, 0.0, 0.0);
257	<	_mdef->set_indiv_par(3, string("sig"), 0.1, 0.001, 0.0, 0.0);
250	>	_mdef->setFormula(formula);
251	>	_mdef->setIndivEnergy(0);
252	>	_mdef->setIndivMeanX(1);
253	>	_mdef->setIndivMeanY(2);
254	>	_mdef->setIndivParameter(0, string("N"), 0.0, 0.0, 0.0, 1.0e6);
255	>	_mdef->setIndivParameter(1, string("mu_x"), 0.0, 0.0, 0.0, 0.0);
256	>	_mdef->setIndivParameter(2, string("mu_y"), 0.0, 0.0, 0.0, 0.0);
257	>	_mdef->setIndivParameter(3, string("sig"), 0.1, 0.001, 0.0, 0.0);
258	>
259	>	_mdef->setMaxNumberOfGaussians(3);
260
261		// get jetcoll from event file and select highest e jet
262		if (info_type_ == 0) {
#	Line 427 | Line 268 | jetfit::model_def& JetFinderAnalyzer::ma
268		if (!found_jet || (*jet_collection)[i].energy() > highest_e_jet.energy())
269		{
270		highest_e_jet = (*jet_collection)[i];
271	+	found_jet = true;
272		}
273		}
274		vector<reco::GenJet> highest_e_jet_coll;
275		highest_e_jet_coll.push_back(highest_e_jet);
276	<	seed_with_jetcoll(highest_e_jet_coll, *_mdef);
276	>	seed_with_jetcoll(highest_e_jet_coll, *_mdef,
277	>	histo->GetYaxis()->GetXmin(),
278	>	histo->GetYaxis()->GetXmax());
279		}
280		if (info_type_ == 1) {
281		edm::Handle< vector<reco::PFJet> > jet_collection;
#	Line 442 | Line 286 | jetfit::model_def& JetFinderAnalyzer::ma
286		if (!found_jet || (*jet_collection)[i].energy() > highest_e_jet.energy())
287		{
288		highest_e_jet = (*jet_collection)[i];
289	+	found_jet = true;
290		}
291		}
292		vector<reco::PFJet> highest_e_jet_coll;
293		highest_e_jet_coll.push_back(highest_e_jet);
294	<	seed_with_jetcoll(highest_e_jet_coll, *_mdef);
295	<	}
296	<
452	<	jetfit::set_model_def(_mdef);
453	<
454	<	// generate initial fit histogram
455	<	edm::Service<TFileService> fs;
456	<	TH2D *init_fit_histo = fs->make<TH2D>(("init_fit_"+string(histo->GetName()))
457	<	.c_str(),
458	<	("Initial fit for "
459	<	+string(histo->GetName())).c_str(),
460	<	histo->GetXaxis()->GetNbins(),
461	<	histo->GetXaxis()->GetXmin(),
462	<	histo->GetXaxis()->GetXmax(),
463	<	histo->GetXaxis()->GetNbins(),
464	<	histo->GetXaxis()->GetXmin(),
465	<	histo->GetXaxis()->GetXmax());
466	<	double XbinSize = (histo->GetXaxis()->GetXmax()
467	<	- histo->GetXaxis()->GetXmin()) /
468	<	static_cast<double>(histo->GetXaxis()->GetNbins());
469	<	double YbinSize = (histo->GetYaxis()->GetXmax()
470	<	- histo->GetYaxis()->GetXmin()) /
471	<	static_cast<double>(histo->GetYaxis()->GetNbins());
472	<	double Xlo = histo->GetXaxis()->GetXmin();
473	<	double Xhi = histo->GetXaxis()->GetXmax();
474	<	double Ylo = histo->GetYaxis()->GetXmin();
475	<	double Yhi = histo->GetYaxis()->GetXmax();
476	<
477	<	for (int i = 0; i < 60; i++) {
478	<	for (int j = 0; j < 60; j++) {
479	<	double x = (static_cast<double>(i) + 0.5)*XbinSize + Xlo;
480	<	double y = (static_cast<double>(j) + 0.5)*YbinSize + Ylo;
481	<	double pval[256];
482	<	if (_mdef->get_n_special_par_sets() > 64) {
483	<	cerr << "Parameter overload" << endl;
484	<	return *_mdef;
485	<	}
486	<	else {
487	<	for (int is = 0; is < _mdef->get_n_special_par_sets(); is++) {
488	<	for (int ii = 0; ii < 4; ii++) {
489	<	double spval, sperr, splo, sphi;
490	<	_mdef->get_special_par(is, ii, spval, sperr, splo, sphi);
491	<	pval[4*is + ii] = spval;
492	<	}
493	<	}
494	<	}
495	<	jetfit::set_ngauss(_mdef->get_n_special_par_sets());
496	<	init_fit_histo->SetBinContent(i+1, j+1,
497	<	jetfit::fit_fcn(x, y, pval)
498	<	* XbinSize * YbinSize);
499	<	}
294	>	seed_with_jetcoll(highest_e_jet_coll, *_mdef,
295	>	histo->GetYaxis()->GetXmin(),
296	>	histo->GetYaxis()->GetXmax());
297		}
298
299		return *_mdef;
300		}
301
302		void JetFinderAnalyzer::beginJob(const edm::EventSetup &es) {
506	–	ofs.open("jetfindlog.txt", ios::out);
507	–	if (ofs.fail()) {
508	–	cerr << "Opening jetfindlog.txt FAILED" << endl;
509	–	}
510	–	ofs << "Jetfinder log" << endl
511	–	<< "=============" << endl << endl;
512	–	}
513	–
514	–	ostream& operator<<(ostream &out, jetfit::trouble t) {
515	–	string action, error_string;
516	–
517	–	if (t.istat != 3) {
518	–	switch(t.occ) {
519	–	case jetfit::T_NULL:
520	–	action = "Program"; break;
521	–	case jetfit::T_SIMPLEX:
522	–	action = "SIMPLEX"; break;
523	–	case jetfit::T_MIGRAD:
524	–	action = "MIGRAD"; break;
525	–	case jetfit::T_MINOS:
526	–	action = "MINOS"; break;
527	–	default:
528	–	action = "Program"; break;
529	–	}
303
531	–	switch (t.istat) {
532	–	case 0:
533	–	error_string = "Unable to calculate error matrix"; break;
534	–	case 1:
535	–	error_string = "Error matrix a diagonal approximation"; break;
536	–	case 2:
537	–	error_string = "Error matrix not positive definite"; break;
538	–	case 3:
539	–	error_string = "Converged successfully"; break;
540	–	default:
541	–	ostringstream oss;
542	–	oss<<"Unknown status code "<<t.istat << flush;
543	–	error_string = oss.str(); break;
544	–	}
545	–
546	–	if (t.occ != jetfit::T_NULL)
547	–	out << action<<" trouble: "<<error_string;
548	–	else
549	–	out << "Not calculated" << endl;
550	–	}
551	–	else {
552	–	out << "Error matrix accurate" << endl;
553	–	}
554	–
555	–	return out;
304		}
305
306	<	void JetFinderAnalyzer::analyze_results(jetfit::results r,
307	<	std::vector<jetfit::trouble> t,
306	>	void JetFinderAnalyzer::analyze_results(HistoFitter::FitResults r,
307	>	std::vector<HistoFitter::Trouble> t,
308		TH2 *hist_orig) {
309	<	ofs << "Histogram "<<hist_orig->GetName() << endl;
562	<	for (int i = unique_jets[hist_orig].size() - r.chisquare.size();
563	<	i < unique_jets[hist_orig].size(); i++) {
564	<	int ir = i - unique_jets[hist_orig].size() + r.chisquare.size();
565	<	ofs << "For "<<i+1<<" gaussians: " << endl
566	<	<< t.at(i) << endl;
567	<	ofs << "chisquare="<<r.chisquare.at(ir)
568	<	<< endl;
569	<	ofs << unique_jets[hist_orig][i].size()<<" unique jets found" << endl;
570	<	for (int j = 0; j < unique_jets[hist_orig][i].size(); j++) {
571	<	jet _jet = unique_jets[hist_orig][i][j];
572	<	ofs << "Jet "<<j<<": Energy = "<<_jet.energy<<", eta = "<<_jet.eta
573	<	<< ", phi = "<<_jet.phi << endl;
574	<	}
575	<	ofs << endl;
576	<	}
577	<	ofs << endl;
578	<
579	<	// save fit function histograms to root file
580	<	edm::Service<TFileService> fs;
581	<	for (vector< vector<double> >::size_type i = 0;
582	<	i < r.pval.size(); i++) {
583	<	jetfit::set_ngauss(r.pval[i].size() / 4);
584	<	TF2 *tf2 = new TF2("fit_func", jetfit::fit_fcn_TF2,
585	<	hist_orig->GetXaxis()->GetXmin(),
586	<	hist_orig->GetXaxis()->GetXmax(),
587	<	hist_orig->GetYaxis()->GetXmin(),
588	<	hist_orig->GetYaxis()->GetXmax(),
589	<	r.pval[i].size());
590	<	for (vector<double>::size_type j = 0; j < r.pval[i].size(); j++) {
591	<	tf2->SetParameter(j, r.pval[i][j]);
592	<	}
593	<	ostringstream fit_histo_oss;
594	<	fit_histo_oss << hist_orig->GetName()<<"_fit_"<<i << flush;
595	<	tf2->SetNpx(hist_orig->GetXaxis()->GetNbins());
596	<	tf2->SetNpy(hist_orig->GetYaxis()->GetNbins());
597	<	TH2D *fit_histo = fs->make<TH2D>(fit_histo_oss.str().c_str(),
598	<	fit_histo_oss.str().c_str(),
599	<	hist_orig->GetXaxis()->GetNbins(),
600	<	hist_orig->GetXaxis()->GetXmin(),
601	<	hist_orig->GetXaxis()->GetXmax(),
602	<	hist_orig->GetYaxis()->GetNbins(),
603	<	hist_orig->GetYaxis()->GetXmin(),
604	<	hist_orig->GetYaxis()->GetXmax());
605	<	TH1 *tf2_histo = tf2->CreateHistogram();
606	<	double XbinSize = (fit_histo->GetXaxis()->GetXmax()
607	<	- fit_histo->GetXaxis()->GetXmin())
608	<	/ static_cast<double>(fit_histo->GetXaxis()->GetNbins());
609	<	double YbinSize = (fit_histo->GetYaxis()->GetXmax()
610	<	- fit_histo->GetYaxis()->GetXmin())
611	<	/ static_cast<double>(fit_histo->GetYaxis()->GetNbins());
612	<	for (int ih = 0; ih < tf2->GetNpx(); ih++) {
613	<	for (int jh = 0; jh < tf2->GetNpy(); jh++) {
614	<	fit_histo->SetBinContent(ih+1, jh+1,
615	<	tf2_histo->GetBinContent(ih+1, jh+1)
616	<	* XbinSize * YbinSize);
617	<	}
618	<	}
619	<	}
620	<
621	<	// save results to file
622	<	ostringstream res_tree_oss, rt_title_oss;
623	<	res_tree_oss << hist_orig->GetName()<<"_results" << flush;
624	<	rt_title_oss << "Fit results for "<<hist_orig->GetName() << flush;
309	>	// perform analysis of fit results
310		}
311
312		DEFINE_FWK_MODULE(JetFinderAnalyzer);

Diff Legend

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