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.1 by dnisson, Wed Sep 2 21:48:13 2009 UTC vs.
Revision 1.13 by dnisson, Wed Sep 23 03:13:07 2009 UTC

#	Line 2 | Line 2
2
3		#include "UserCode/JetFitAnalyzer/interface/JetFitAnalyzer.h"
4
5	–	#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
5		#include "fastjet/ClusterSequence.hh"
6		#include "FWCore/ServiceRegistry/interface/Service.h"
7	+	#include "FWCore/MessageLogger/interface/MessageLogger.h"
8		#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
9
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 103 | Line 115 | private:
115		TH2 *);
116		virtual void analyze_results(jetfit::results r,
117		std::vector<jetfit::trouble> t, TH2 *);
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 115 | Line 136 | JetFinderAnalyzer::unique_jets;
136		JetFinderAnalyzer::JetFinderAnalyzer(const edm::ParameterSet &pSet)
137		: JetFitAnalyzer(pSet) // this is important!
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	+	}
147	+
148		smear_ = pSet.getUntrackedParameter("smear", 0.02);
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
156	<	TH2D * JetFinderAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&) {
157	<	ostringstream oss;
158	<	oss << "eta_phi_energy_unred"<<evt.id().event() << flush;
159	<	TH2D *unred_histo = new TH2D(oss.str().c_str(), oss.str().c_str(),
160	<	600, -2.5, 2.5, 600, -PI, PI);
161	<
162	<	// fill unreduced histo
163	<	edm::Handle<edm::HepMCProduct> hRaw;
164	<	evt.getByLabel("source",
165	<	hRaw);
166	<	vector<HepMC::GenParticle *> particles;
167	<	const HepMC::GenEvent *hmcEvt = hRaw->GetEvent();
168	<	for (HepMC::GenEvent::particle_const_iterator
138	<	pit = hmcEvt->particles_begin(); pit != hmcEvt->particles_end();
139	<	pit++) {
140	<	if ((*pit)->status() == 1) {
141	<	particles.push_back(*pit);
142	<	}
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	<	for (int i = 0; i < particles.size(); i++) {
174	<	unred_histo->Fill(particles[i]->momentum().eta(),
175	<	particles[i]->momentum().phi(),
176	<	particles[i]->momentum().e());
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	<	// reduce histo
191	<	ostringstream oss2;
192	<	oss2 << "eta_phi_energy_red"<<evt.id().event() << flush;
193	<	edm::Service<TFileService> fs;
194	<	// draw cone of radius 0.5 around highest energy bin, reduce
195	<	double maxE = 0.0;
196	<	int max_i = 29, max_j = 29;
197	<	for (int i = 0; i < unred_histo->GetNbinsX(); i++) {
198	<	for (int j = 0; j < unred_histo->GetNbinsY(); j++) {
199	<	double E = unred_histo->GetBinContent(i+1, j+1);
200	<	if (E > maxE) {
201	<	maxE = E;
202	<	max_i = i;
203	<	max_j = j;
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	<	double rcone = 0.5;
225	<	double Xlo = unred_histo->GetXaxis()->GetXmin();
226	<	double Xhi = unred_histo->GetXaxis()->GetXmax();
227	<	double Ylo = unred_histo->GetYaxis()->GetXmin();
228	<	double Yhi = unred_histo->GetYaxis()->GetXmax();
229	<	double XbinSize = (Xhi - Xlo) /
230	<	static_cast<double>(unred_histo->GetXaxis()->GetNbins());
231	<	double YbinSize = (Yhi - Ylo) /
232	<	static_cast<double>(unred_histo->GetYaxis()->GetNbins());
233	<	double max_x = (static_cast<double>(max_i) + 0.5) * XbinSize + Xlo;
234	<	double max_y = (static_cast<double>(max_j) + 0.5) * YbinSize + Ylo;
235	<	TH2D *histo = fs->make<TH2D>(oss2.str().c_str(), oss2.str().c_str(),
236	<	60, max_x-rcone, max_x+rcone,
237	<	60, max_y-rcone, max_y+rcone);
238	<
239	<	// create an unsmeared reduced histo
240	<	TH2D *histo_unsmeared = fs->make<TH2D>((oss2.str()+"_unsmeared").c_str(),
241	<	(oss2.str()+"_unsmeared").c_str(),
242	<	60, max_x-rcone, max_x+rcone,
243	<	60, max_y-rcone, max_y+rcone);
223	>
224	>	return particles;
225	>	}
226	>
227	>	TH2D * JetFinderAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&) {
228	>	const reco::Jet * highest_e_jet;
229	>
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	>	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	>	if (highest_e_jet == 0) {
250	>	cerr << "No fatjets found!" << endl; exit(1);
251	>	}
252	>
253	>	vector<const reco::Candidate *> particles =
254	>	highest_e_jet->getJetConstituentsQuick();
255	>	//  cout << "found highest e jet" << endl;
256	>
257	>	ostringstream oss;
258	>	oss << "eta_phi_energy"<<evt.id().event() << flush;
259	>	TH2D *histo = new TH2D(oss.str().c_str(), oss.str().c_str(),
260	>	600,
261	>	highest_e_jet->eta()-0.5,
262	>	highest_e_jet->eta()+0.5,
263	>	600,
264	>	highest_e_jet->phi()-0.5,
265	>	highest_e_jet->phi()+0.5);
266	>
267		for (int i = 0; i < particles.size(); i++) {
268	<	double N = particles[i]->momentum().e();
269	<	double x = particles[i]->momentum().eta();
270	<	double y = particles[i]->momentum().phi();
193	<	histo_unsmeared->Fill(x, y, N);
268	>	histo->Fill(particles[i]->eta(),
269	>	particles[i]->phi(),
270	>	particles[i]->energy());
271		}
272
273	<	// create a smeared reduced histo
273	>	// create a smeared histo
274		// create a temporary 2D vector for smeared energies
275	<	XbinSize = (histo->GetXaxis()->GetXmax()
275	>	double XbinSize = (histo->GetXaxis()->GetXmax()
276		- histo->GetXaxis()->GetXmin()) /
277		static_cast<double>(histo->GetXaxis()->GetNbins());
278	<	YbinSize = (histo->GetYaxis()->GetXmax()
278	>	double YbinSize = (histo->GetYaxis()->GetXmax()
279		- histo->GetYaxis()->GetXmin()) /
280		static_cast<double>(histo->GetYaxis()->GetNbins());
281	+	double Xlo = histo->GetXaxis()->GetXmin();
282	+	double Ylo = histo->GetYaxis()->GetXmin();
283		vector< vector<double> > smeared(60, vector<double>(60, 0.0) );
284		switch (smear_coord_) {
285		case 1:
286		for (int i = 0; i < particles.size(); i++) {
287	<	double N = particles[i]->momentum().e();
288	<	double x = particles[i]->momentum().eta();
289	<	double y = particles[i]->momentum().phi();
287	>	double N = particles[i]->energy();
288	>	double x = particles[i]->eta();
289	>	double y = particles[i]->phi();
290		// loop over bins and add Gaussian in proper angle to smeared
291		for (vector< vector<double> >::size_type i2 = 0; i2 < 60; i2++) {
292		for (vector< double >::size_type j2 = 0; j2 < 60; j2++) {
293		double eta = static_cast<double>((signed int)i2) * XbinSize +
294	<	max_x - rcone - x;
295	<	double phi = acos(cos(static_cast<double>((signed int)j2) * YbinSize +
296	<	max_y - rcone - y));
297	<	phi = sin(phi) > 0 ? phi : -phi;
294	>	Xlo;
295	>	double phi0 = static_cast<double>((signed int)j2) * YbinSize +
296	>	Ylo;
297	>	double phi = acos(cos(phi0));
298	>	phi = sin(phi0) > 0 ? phi : -phi;
299
300		// transform eta, phi to proper angle
301		double theta = 2.0*atan(exp(-eta));
#	Line 223 | Line 303 | TH2D * JetFinderAnalyzer::make_histo(con
303
304		smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
305		* 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	+	}
323		}
324		}
325		}
#	Line 230 | Line 327 | TH2D * JetFinderAnalyzer::make_histo(con
327		case 0:
328		default:
329		for (int i = 0; i < particles.size(); i++) {
330	<	double N = particles[i]->momentum().e();
331	<	double x = particles[i]->momentum().eta();
332	<	double y = particles[i]->momentum().phi();
330	>	double N = particles[i]->energy();
331	>	double x = particles[i]->eta();
332	>	double y = particles[i]->phi();
333		// loop over bins and add Gaussian to smeared
334		for (vector< vector<double> >::size_type i2 = 0; i2 < 60; i2++) {
335		for (vector< double >::size_type j2 = 0; j2 < 60; j2++) {
336		double eta = static_cast<double>((signed int)i2) * XbinSize
337	<	+ max_x - rcone - x;
338	<	double phi = acos(cos(static_cast<double>((signed int)j2) * YbinSize
339	<	+ max_y - rcone - y));
340	<	phi = sin(phi) > 0 ? phi : -phi;
337	>	+ Xlo;
338	>	double phi0 = static_cast<double>((signed int)j2) * YbinSize +
339	>	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_));
344	+
345	+	// 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_))
357	+	* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_));
358	+	}
359		}
360		}
361		}
#	Line 257 | Line 370 | TH2D * JetFinderAnalyzer::make_histo(con
370		return histo;
371		}
372
373	<	void seed_with_CA(const edm::Event &evt, TH2 *histo,
373	>	template <class Jet>
374	>	void seed_with_jetcoll(vector<Jet> jets,
375		jetfit::model_def &_mdef) {
376	<	edm::Handle<edm::HepMCProduct> hMC;
263	<	evt.getByLabel("source", hMC);
264	<
265	<	// create a PseudoJet vector
266	<	vector<PseudoJet> particles;
267	<	const HepMC::GenEvent *hmcEvt = hMC->GetEvent();
268	<	for (HepMC::GenEvent::particle_const_iterator
269	<	pit = hmcEvt->particles_begin(); pit != hmcEvt->particles_end();
270	<	pit++) {
271	<	if ((*pit)->status() == 1) {
272	<	double x_max = (histo->GetXaxis()->GetXmax()
273	<	+ histo->GetXaxis()->GetXmin()) / 2.0;
274	<	double y_max = (histo->GetYaxis()->GetXmax()
275	<	+ histo->GetYaxis()->GetXmin()) / 2.0;
276	<	valarray<double> pmom(4);
277	<	pmom[0] = (*pit)->momentum().px();
278	<	pmom[1] = (*pit)->momentum().py();
279	<	pmom[2] = (*pit)->momentum().pz();
280	<	pmom[3] = (*pit)->momentum().e();
281	<	double eta = (*pit)->momentum().eta();
282	<	double phi = (*pit)->momentum().phi();
283	<	if ((eta - x_max)*(eta - x_max) + (phi - y_max)*(phi - y_max) < 0.25) {
284	<	PseudoJet j(pmom);
285	<	particles.push_back(j);
286	<	}
287	<	}
288	<	}
289	<
290	<	// choose a jet definition
291	<	double R = 0.2;
292	<	JetDefinition jet_def(cambridge_algorithm, R);
293	<
294	<	// run clustering and extract the jets
295	<	ClusterSequence cs(particles, jet_def);
296	<	vector<PseudoJet> jets = cs.inclusive_jets();
297	<
298	<	double XbinSize = (histo->GetXaxis()->GetXmax()
299	<	- histo->GetXaxis()->GetXmin()) /
300	<	static_cast<double>(histo->GetXaxis()->GetNbins());
301	<	double YbinSize = (histo->GetYaxis()->GetXmax()
302	<	- histo->GetYaxis()->GetXmin()) /
303	<	static_cast<double>(histo->GetYaxis()->GetNbins());
304	<
305	<	// seed with C-A jets
376	>	// seed with jet collection
377		int ijset = 0;
378		for (unsigned ij = 0; ij < jets.size(); ij++) {
379	<	double N = jets[ij].e();
379	>	double N = jets[ij].energy();
380	>	double eta = jets[ij].eta();
381	>	double phi = jets[ij].phi();
382		if (N > 10.0) {
383		_mdef.set_special_par(ijset, 0, N, _mdef.chisquare_error(N)*0.1,
384		0.0, 1.0e6);
385	<	_mdef.set_special_par(ijset, 1, jets[ij].eta(), 0.01,
385	>	_mdef.set_special_par(ijset, 1, eta, 0.01,
386		0.0, 0.0);
387	<	double mdef_phi = jets[ij].phi() > PI ? jets[ij].phi() - 2*PI
388	<	: jets[ij].phi();
387	>	double mdef_phi = phi > PI ? phi - 2*PI
388	>	: phi;
389		_mdef.set_special_par(ijset, 2, mdef_phi, 0.01,
390		0.0, 0.0);
391		_mdef.set_special_par(ijset, 3, 0.1, 0.001, 0.0, 0.0);
#	Line 327 | Line 400 | jetfit::model_def& JetFinderAnalyzer::ma
400		class jf_model_def : public jetfit::model_def {
401		public:
402		virtual double chisquare_error(double E) {
403	<	return 0.97*E + 14.0;
404	<	// study from 08-27-09
403	>	return exp(-(4.2 + 0.11*E));
404	>	// study from 09-04-09
405		}
406		};
407
#	Line 344 | Line 417 | jetfit::model_def& JetFinderAnalyzer::ma
417		_mdef->set_indiv_par(2, string("mu_y"), 0.0, 0.0, 0.0, 0.0);
418		_mdef->set_indiv_par(3, string("sig"), 0.1, 0.001, 0.0, 0.0);
419
420	<	seed_with_CA(evt, histo, *_mdef);
420	>	// get jetcoll from event file and select highest e jet
421	>	if (info_type_ == 0) {
422	>	edm::Handle< vector<reco::GenJet> > jet_collection;
423	>	evt.getByLabel(jet_algo_, jet_collection);
424	>	reco::GenJet highest_e_jet;
425	>	bool found_jet = false;
426	>	for (unsigned i = 0; i < jet_collection->size(); i++) {
427	>	if (!found_jet || (*jet_collection)[i].energy() > highest_e_jet.energy())
428	>	{
429	>	highest_e_jet = (*jet_collection)[i];
430	>	}
431	>	}
432	>	vector<reco::GenJet> highest_e_jet_coll;
433	>	highest_e_jet_coll.push_back(highest_e_jet);
434	>	seed_with_jetcoll(highest_e_jet_coll, *_mdef);
435	>	}
436	>	if (info_type_ == 1) {
437	>	edm::Handle< vector<reco::PFJet> > jet_collection;
438	>	evt.getByLabel(jet_algo_, jet_collection);
439	>	reco::PFJet highest_e_jet;
440	>	bool found_jet = false;
441	>	for (unsigned i = 0; i < jet_collection->size(); i++) {
442	>	if (!found_jet || (*jet_collection)[i].energy() > highest_e_jet.energy())
443	>	{
444	>	highest_e_jet = (*jet_collection)[i];
445	>	}
446	>	}
447	>	vector<reco::PFJet> highest_e_jet_coll;
448	>	highest_e_jet_coll.push_back(highest_e_jet);
449	>	seed_with_jetcoll(highest_e_jet_coll, *_mdef);
450	>	}
451
452		jetfit::set_model_def(_mdef);
453
#	Line 391 | Line 494 | jetfit::model_def& JetFinderAnalyzer::ma
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));
497	>	jetfit::fit_fcn(x, y, pval)
498	>	* XbinSize * YbinSize);
499		}
500		}
501
#	Line 444 | Line 548 | ostream& operator<<(ostream &out, jetfit
548		else
549		out << "Not calculated" << endl;
550		}
551	+	else {
552	+	out << "Error matrix accurate" << endl;
553	+	}
554
555		return out;
556		}
#	Line 452 | Line 559 | void JetFinderAnalyzer::analyze_results(
559		std::vector<jetfit::trouble> t,
560		TH2 *hist_orig) {
561		ofs << "Histogram "<<hist_orig->GetName() << endl;
562	<	for (int i = 0; i < unique_jets[hist_orig].size(); i++) {
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	<	<< unique_jets[hist_orig][i].size()<<" unique jets found" << 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

Diff Legend

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