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.5 by dnisson, Fri Sep 4 15:30:38 2009 UTC vs.
Revision 1.14 by dnisson, Wed Sep 23 03:31:30 2009 UTC

#	Line 13 | Line 13
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 121 | Line 127 | private:
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 142 | 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 217 | Line 225 | vector<reco::Candidate *> JetFinderAnaly
225		}
226
227		TH2D * JetFinderAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&) {
228	<	ostringstream oss;
221	<	oss << "eta_phi_energy_unred"<<evt.id().event() << flush;
222	<	TH2D *unred_histo = new TH2D(oss.str().c_str(), oss.str().c_str(),
223	<	600, -2.5, 2.5, 600, -PI, PI);
228	>	const reco::Jet * highest_e_jet;
229
230	<	vector<reco::Candidate *> particles = get_particles(evt);
231	<	for (int i = 0; i < particles.size(); i++) {
232	<	unred_histo->Fill(particles[i]->eta(),
233	<	particles[i]->phi(),
234	<	particles[i]->energy());
230	>	if (info_type_ == 0) {
231	>	edm::Handle< vector<reco::GenJet> > evtJets;
232	>	evt.getByLabel(jet_algo_, evtJets);
233	>	for (unsigned i = 0; i < evtJets->size(); i++) {
234	>	if (i == 0 || (*evtJets)[i].energy() > highest_e_jet->energy()) {
235	>	highest_e_jet = &((*evtJets)[i]);
236	>	}
237	>	}
238		}
239	<
240	<	// reduce histo
241	<	ostringstream oss2;
242	<	oss2 << "eta_phi_energy_red"<<evt.id().event() << flush;
243	<	edm::Service<TFileService> fs;
244	<	// draw cone of radius 0.5 around highest energy bin, reduce
237	<	double maxE = 0.0;
238	<	int max_i = 29, max_j = 29;
239	<	for (int i = 0; i < unred_histo->GetNbinsX(); i++) {
240	<	for (int j = 0; j < unred_histo->GetNbinsY(); j++) {
241	<	double E = unred_histo->GetBinContent(i+1, j+1);
242	<	if (E > maxE) {
243	<	maxE = E;
244	<	max_i = i;
245	<	max_j = j;
239	>	else if (info_type_ == 1) {
240	>	edm::Handle< vector<reco::PFJet> > evtJets;
241	>	evt.getByLabel(jet_algo_, evtJets);
242	>	for (unsigned i = 0; i < evtJets->size(); i++) {
243	>	if (i == 0 || (*evtJets)[i].energy() > highest_e_jet->energy()) {
244	>	highest_e_jet = &((*evtJets)[i]);
245		}
246		}
247		}
248	<
249	<	double rcone = 0.5;
250	<	double Xlo = unred_histo->GetXaxis()->GetXmin();
251	<	double Xhi = unred_histo->GetXaxis()->GetXmax();
252	<	double Ylo = unred_histo->GetYaxis()->GetXmin();
253	<	double Yhi = unred_histo->GetYaxis()->GetXmax();
254	<	double XbinSize = (Xhi - Xlo) /
255	<	static_cast<double>(unred_histo->GetXaxis()->GetNbins());
256	<	double YbinSize = (Yhi - Ylo) /
257	<	static_cast<double>(unred_histo->GetYaxis()->GetNbins());
258	<	double max_x = (static_cast<double>(max_i) + 0.5) * XbinSize + Xlo;
259	<	double max_y = (static_cast<double>(max_j) + 0.5) * YbinSize + Ylo;
260	<	TH2D *histo = fs->make<TH2D>(oss2.str().c_str(), oss2.str().c_str(),
261	<	60, max_x-rcone, max_x+rcone,
262	<	60, max_y-rcone, max_y+rcone);
263	<
264	<	// create an unsmeared reduced histo
265	<	TH2D *histo_unsmeared = fs->make<TH2D>((oss2.str()+"_unsmeared").c_str(),
266	<	(oss2.str()+"_unsmeared").c_str(),
267	<	60, max_x-rcone, max_x+rcone,
269	<	60, max_y-rcone, max_y+rcone);
248	>	cout << "found highest e jet" << endl;
249	>
250	>	if (highest_e_jet == 0) {
251	>	cerr << "No fatjets found!" << endl; exit(1);
252	>	}
253	>
254	>	vector<const reco::Candidate *> particles =
255	>	highest_e_jet->getJetConstituentsQuick();
256	>	cout << "found highest e jet" << endl;
257	>
258	>	ostringstream oss;
259	>	oss << "eta_phi_energy"<<evt.id().event() << flush;
260	>	TH2D *histo = new TH2D(oss.str().c_str(), oss.str().c_str(),
261	>	30,
262	>	highest_e_jet->eta()-0.5,
263	>	highest_e_jet->eta()+0.5,
264	>	30,
265	>	highest_e_jet->phi()-0.5,
266	>	highest_e_jet->phi()+0.5);
267	>
268		for (int i = 0; i < particles.size(); i++) {
269	<	double N = particles[i]->energy();
270	<	double x = particles[i]->eta();
271	<	double y = particles[i]->phi();
274	<	histo_unsmeared->Fill(x, y, N);
269	>	histo->Fill(particles[i]->eta(),
270	>	particles[i]->phi(),
271	>	particles[i]->energy());
272		}
273
274	<	// create a smeared reduced histo
274	>	// create a smeared histo
275		// create a temporary 2D vector for smeared energies
276	<	XbinSize = (histo->GetXaxis()->GetXmax()
276	>	double XbinSize = (histo->GetXaxis()->GetXmax()
277		- histo->GetXaxis()->GetXmin()) /
278		static_cast<double>(histo->GetXaxis()->GetNbins());
279	<	YbinSize = (histo->GetYaxis()->GetXmax()
279	>	double YbinSize = (histo->GetYaxis()->GetXmax()
280		- histo->GetYaxis()->GetXmin()) /
281		static_cast<double>(histo->GetYaxis()->GetNbins());
282	<	vector< vector<double> > smeared(60, vector<double>(60, 0.0) );
282	>	double Xlo = histo->GetXaxis()->GetXmin();
283	>	double Ylo = histo->GetYaxis()->GetXmin();
284	>	vector< vector<double> > smeared(30, vector<double>(30, 0.0) );
285		switch (smear_coord_) {
286		case 1:
287		for (int i = 0; i < particles.size(); i++) {
#	Line 290 | Line 289 | TH2D * JetFinderAnalyzer::make_histo(con
289		double x = particles[i]->eta();
290		double y = particles[i]->phi();
291		// loop over bins and add Gaussian in proper angle to smeared
292	<	for (vector< vector<double> >::size_type i2 = 0; i2 < 60; i2++) {
293	<	for (vector< double >::size_type j2 = 0; j2 < 60; j2++) {
292	>	for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) {
293	>	for (vector< double >::size_type j2 = 0; j2 < 30; j2++) {
294		double eta = static_cast<double>((signed int)i2) * XbinSize +
295	<	max_x - rcone - x;
296	<	double phi = acos(cos(static_cast<double>((signed int)j2) * YbinSize +
297	<	max_y - rcone - y));
298	<	phi = sin(phi) > 0 ? phi : -phi;
295	>	Xlo;
296	>	double phi0 = static_cast<double>((signed int)j2) * YbinSize +
297	>	Ylo;
298	>	double phi = acos(cos(phi0));
299	>	phi = sin(phi0) > 0 ? phi : -phi;
300
301		// transform eta, phi to proper angle
302		double theta = 2.0*atan(exp(-eta));
#	Line 304 | Line 304 | TH2D * JetFinderAnalyzer::make_histo(con
304
305		smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
306		* 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	+	}
324		}
325		}
326		}
#	Line 315 | Line 332 | TH2D * JetFinderAnalyzer::make_histo(con
332		double x = particles[i]->eta();
333		double y = particles[i]->phi();
334		// loop over bins and add Gaussian to smeared
335	<	for (vector< vector<double> >::size_type i2 = 0; i2 < 60; i2++) {
336	<	for (vector< double >::size_type j2 = 0; j2 < 60; j2++) {
335	>	for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) {
336	>	for (vector< double >::size_type j2 = 0; j2 < 30; j2++) {
337		double eta = static_cast<double>((signed int)i2) * XbinSize
338	<	+ max_x - rcone - x;
339	<	double phi = acos(cos(static_cast<double>((signed int)j2) * YbinSize
340	<	+ max_y - rcone - y));
341	<	phi = sin(phi) > 0 ? phi : -phi;
338	>	+ Xlo;
339	>	double phi0 = static_cast<double>((signed int)j2) * YbinSize +
340	>	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_));
345	+
346	+	// 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_))
358	+	* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_));
359	+	}
360		}
361		}
362		}
363		}
364		// set histogram to match smear vector
365	<	for (int i = 1; i <= 60; i++) {
366	<	for (int j = 1; j <= 60; j++) {
365	>	for (int i = 1; i <= 30; i++) {
366	>	for (int j = 1; j <= 30; j++) {
367		histo->SetBinContent(i, j, smeared[i-1][j-1]);
368		}
369		}
#	Line 338 | Line 371 | TH2D * JetFinderAnalyzer::make_histo(con
371		return histo;
372		}
373
374	<	void seed_with_CA(vector<reco::Candidate *> gParticles, TH2 *histo,
374	>	template <class Jet>
375	>	void seed_with_jetcoll(vector<Jet> jets,
376		jetfit::model_def &_mdef) {
377	<	// create a PseudoJet vector
344	<	vector<PseudoJet> particles;
345	<	for (unsigned i = 0; i < gParticles.size(); i++) {
346	<	double x_max = (histo->GetXaxis()->GetXmax()
347	<	+ histo->GetXaxis()->GetXmin()) / 2.0;
348	<	double y_max = (histo->GetYaxis()->GetXmax()
349	<	+ histo->GetYaxis()->GetXmin()) / 2.0;
350	<	valarray<double> pmom(4);
351	<	pmom[0] = gParticles[i]->px();
352	<	pmom[1] = gParticles[i]->py();
353	<	pmom[2] = gParticles[i]->pz();
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) {
358	<	PseudoJet j(pmom);
359	<	particles.push_back(j);
360	<	}
361	<	}
362	<
363	<	// choose a jet definition
364	<	double R = 0.2;
365	<	JetDefinition jet_def(cambridge_algorithm, R);
366	<
367	<	// run clustering and extract the jets
368	<	ClusterSequence cs(particles, jet_def);
369	<	vector<PseudoJet> jets = cs.inclusive_jets();
370	<
371	<	double XbinSize = (histo->GetXaxis()->GetXmax()
372	<	- histo->GetXaxis()->GetXmin()) /
373	<	static_cast<double>(histo->GetXaxis()->GetNbins());
374	<	double YbinSize = (histo->GetYaxis()->GetXmax()
375	<	- histo->GetYaxis()->GetXmin()) /
376	<	static_cast<double>(histo->GetYaxis()->GetNbins());
377	<
378	<	// seed with C-A jets
377	>	// seed with jet collection
378		int ijset = 0;
379		for (unsigned ij = 0; ij < jets.size(); ij++) {
380	<	double N = jets[ij].e();
381	<	if (N > 50.0) {
382	<	cout << N << endl;
380	>	double N = jets[ij].energy();
381	>	double eta = jets[ij].eta();
382	>	double phi = jets[ij].phi();
383	>	if (N > 10.0) {
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,
386	>	_mdef.set_special_par(ijset, 1, eta, 0.01,
387		0.0, 0.0);
388	<	double mdef_phi = jets[ij].phi() > PI ? jets[ij].phi() - 2*PI
389	<	: jets[ij].phi();
388	>	double mdef_phi = phi > PI ? phi - 2*PI
389	>	: phi;
390		_mdef.set_special_par(ijset, 2, mdef_phi, 0.01,
391		0.0, 0.0);
392		_mdef.set_special_par(ijset, 3, 0.1, 0.001, 0.0, 0.0);
#	Line 401 | 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 E > 2.7 ? 0.55*E - 1.0 : 0.5;
405	<	// study from 09-03-09
404	>	return exp(-(4.2 + 0.11*E));
405	>	// study from 09-04-09
406		}
407		};
408
#	Line 418 | Line 418 | jetfit::model_def& JetFinderAnalyzer::ma
418		_mdef->set_indiv_par(2, string("mu_y"), 0.0, 0.0, 0.0, 0.0);
419		_mdef->set_indiv_par(3, string("sig"), 0.1, 0.001, 0.0, 0.0);
420
421	<	seed_with_CA(get_particles(evt), histo, *_mdef);
421	>	// get jetcoll from event file and select highest e jet
422	>	if (info_type_ == 0) {
423	>	edm::Handle< vector<reco::GenJet> > jet_collection;
424	>	evt.getByLabel(jet_algo_, jet_collection);
425	>	reco::GenJet highest_e_jet;
426	>	bool found_jet = false;
427	>	for (unsigned i = 0; i < jet_collection->size(); i++) {
428	>	if (!found_jet || (*jet_collection)[i].energy() > highest_e_jet.energy())
429	>	{
430	>	highest_e_jet = (*jet_collection)[i];
431	>	}
432	>	}
433	>	vector<reco::GenJet> highest_e_jet_coll;
434	>	highest_e_jet_coll.push_back(highest_e_jet);
435	>	seed_with_jetcoll(highest_e_jet_coll, *_mdef);
436	>	}
437	>	if (info_type_ == 1) {
438	>	edm::Handle< vector<reco::PFJet> > jet_collection;
439	>	evt.getByLabel(jet_algo_, jet_collection);
440	>	reco::PFJet highest_e_jet;
441	>	bool found_jet = false;
442	>	for (unsigned i = 0; i < jet_collection->size(); i++) {
443	>	if (!found_jet || (*jet_collection)[i].energy() > highest_e_jet.energy())
444	>	{
445	>	highest_e_jet = (*jet_collection)[i];
446	>	}
447	>	}
448	>	vector<reco::PFJet> highest_e_jet_coll;
449	>	highest_e_jet_coll.push_back(highest_e_jet);
450	>	seed_with_jetcoll(highest_e_jet_coll, *_mdef);
451	>	}
452
453		jetfit::set_model_def(_mdef);
454
#	Line 465 | Line 495 | jetfit::model_def& JetFinderAnalyzer::ma
495		}
496		jetfit::set_ngauss(_mdef->get_n_special_par_sets());
497		init_fit_histo->SetBinContent(i+1, j+1,
498	<	jetfit::fit_fcn(x, y, pval));
498	>	jetfit::fit_fcn(x, y, pval)
499	>	* XbinSize * YbinSize);
500		}
501		}
502
#	Line 518 | Line 549 | ostream& operator<<(ostream &out, jetfit
549		else
550		out << "Not calculated" << endl;
551		}
552	+	else {
553	+	out << "Error matrix accurate" << endl;
554	+	}
555
556		return out;
557		}
#	Line 526 | Line 560 | void JetFinderAnalyzer::analyze_results(
560		std::vector<jetfit::trouble> t,
561		TH2 *hist_orig) {
562		ofs << "Histogram "<<hist_orig->GetName() << endl;
563	<	for (int i = 0; i < unique_jets[hist_orig].size(); i++) {
563	>	for (int i = unique_jets[hist_orig].size() - r.chisquare.size();
564	>	i < unique_jets[hist_orig].size(); i++) {
565	>	int ir = i - unique_jets[hist_orig].size() + r.chisquare.size();
566		ofs << "For "<<i+1<<" gaussians: " << endl
567	<	<< t.at(i) << endl
568	<	<< unique_jets[hist_orig][i].size()<<" unique jets found" << endl;
567	>	<< t.at(i) << endl;
568	>	ofs << "chisquare="<<r.chisquare.at(ir)
569	>	<< endl;
570	>	ofs << unique_jets[hist_orig][i].size()<<" unique jets found" << endl;
571		for (int j = 0; j < unique_jets[hist_orig][i].size(); j++) {
572		jet _jet = unique_jets[hist_orig][i][j];
573		ofs << "Jet "<<j<<": Energy = "<<_jet.energy<<", eta = "<<_jet.eta

Diff Legend

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