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.14 by dnisson, Wed Sep 23 03:31:30 2009 UTC vs.
Revision 1.29 by dnisson, Wed Apr 7 23:43:32 2010 UTC

#	Line 1 | Line 1
1	<	/* A simple jet-finding analyzer */
1	>	/* A JetFitAnalyze_hat makes histograms with smearing */
2
3		#include "UserCode/JetFitAnalyzer/interface/JetFitAnalyzer.h"
4
5	–	#include "fastjet/ClusterSequence.hh"
5		#include "FWCore/ServiceRegistry/interface/Service.h"
6		#include "FWCore/MessageLogger/interface/MessageLogger.h"
7		#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
8	+	#include "DataFormats/JetReco/interface/Jet.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>
25	<	#include <cmath>
26	<	#include <cstdlib>
27	<	#include <fstream>
10	>	#include <iostream>
11		#include <sstream>
12
30	–	#include "TFormula.h"
13		#include "TF2.h"
14	+	#include "TNtuple.h"
15	+	#include "TH1.h"
16
17	<	#define PI 3.141593
17	>	#define PI 3.141593
18
19		using namespace std;
36	–	using namespace fastjet;
20
21		class JetFinderAnalyzer : public JetFitAnalyzer {
22		public:
40	–	struct jet {
41	–	double energy;
42	–	double eta;
43	–	double phi;
44	–	};
45	–
23		explicit JetFinderAnalyzer( const edm::ParameterSet&);
24		~JetFinderAnalyzer() {}
25
26		private:
27	<	static map<TH2 *, vector< vector<jet> > > unique_jets;
28	<
29	<	static double phi_cutoff_;
27	>	virtual void beginJob(const edm::EventSetup &es);
28	>	virtual void analyze_results(HistoFitter::FitResults, std::vector<HistoFitter::Trouble>,
29	>	TH2 *);
30
31	<	static double g2int(double xlo, double xhi, double ylo, double yhi,
32	<	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&,
114	<	const edm::EventSetup&,
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_;
31	>	TH1D *pairMassHist;
32	>	TH1D *topMassHist;
33		};
34
133	–	map<TH2 *, vector< vector< JetFinderAnalyzer::jet > > >
134	–	JetFinderAnalyzer::unique_jets;
135	–
35		JetFinderAnalyzer::JetFinderAnalyzer(const edm::ParameterSet &pSet)
36		: JetFitAnalyzer(pSet) // this is important!
37		{
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	–	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	–
38		}
39
40	<	void
41	<	JetFinderAnalyzer::fetchCandidateCollection(edm::Handle<reco::GenParticleCollection>& c,
42	<	const edm::InputTag& tag,
43	<	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	<	}
40	>	void JetFinderAnalyzer::beginJob(const edm::EventSetup &es) {
41	>	edm::Service<TFileService> fs;
42	>	pairMassHist = fs->make<TH1D>("pairMass", "Mass of Jet Pairs",
43	>	100, 0.0, 250.0);
44
45		}
46
47	<	vector<reco::Candidate *> JetFinderAnalyzer::get_particles(const edm::Event &evt) {
48	<	// fill unreduced histo
49	<	edm::Handle<reco::GenParticleCollection> hRaw;
50	<	edm::Handle<reco::PFCandidateCollection> hPFlow;
51	<	if (info_type_ == 0) {
52	<	fetchCandidateCollection(hRaw,
53	<	inputTagGenParticles_,
54	<	evt);
55	<	}
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;
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	<	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	<	histo->Fill(particles[i]->eta(),
270	<	particles[i]->phi(),
271	<	particles[i]->energy());
272	<	}
273	<
274	<	// create a smeared histo
275	<	// create a temporary 2D vector for smeared energies
276	<	double XbinSize = (histo->GetXaxis()->GetXmax()
277	<	- histo->GetXaxis()->GetXmin()) /
278	<	static_cast<double>(histo->GetXaxis()->GetNbins());
279	<	double YbinSize = (histo->GetYaxis()->GetXmax()
280	<	- histo->GetYaxis()->GetXmin()) /
281	<	static_cast<double>(histo->GetYaxis()->GetNbins());
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++) {
288	<	double N = particles[i]->energy();
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 < 30; i2++) {
293	<	for (vector< double >::size_type j2 = 0; j2 < 30; j2++) {
294	<	double eta = static_cast<double>((signed int)i2) * XbinSize +
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));
303	<	double iota = asin(sin(theta)*sin(phi));
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	<	}
327	<	break;
328	<	case 0:
329	<	default:
330	<	for (int i = 0; i < particles.size(); i++) {
331	<	double N = particles[i]->energy();
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 < 30; i2++) {
336	<	for (vector< double >::size_type j2 = 0; j2 < 30; j2++) {
337	<	double eta = static_cast<double>((signed int)i2) * XbinSize
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 <= 30; i++) {
366	<	for (int j = 1; j <= 30; j++) {
367	<	histo->SetBinContent(i, j, smeared[i-1][j-1]);
368	<	}
369	<	}
47	>	double evalFitFunction(HistoFitter::FitResults r, double x, double y,
48	>	int i) {
49	>	unsigned nGauss = r.pars[i].size() / 4;
50	>	double fitVal = 0.0;
51	>	for (unsigned j = 0; j < nGauss; j++) {
52	>	double N = r.pval[i][4*j];
53	>	double mu_x = r.pval[i][4*j + 1];
54	>	double mu_y = r.pval[i][4*j + 2];
55	>	double sig = r.pval[i][4*j + 3];
56
57	<	return histo;
58	<	}
59	<
374	<	template <class Jet>
375	<	void seed_with_jetcoll(vector<Jet> jets,
376	<	jetfit::model_def &_mdef) {
377	<	// seed with jet collection
378	<	int ijset = 0;
379	<	for (unsigned ij = 0; ij < jets.size(); ij++) {
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, eta, 0.01,
387	<	0.0, 0.0);
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);
393	<	ijset++;
394	<	}
57	>	double rel_x = x - mu_x; double rel_y = y - mu_y;
58	>	fitVal += (N / 2.0 / M_PI / sig / sig)
59	>	* exp(-(rel_x * rel_x + rel_y * rel_y)/2.0/sig/sig);
60		}
61	+	return fitVal;
62		}
63
64	<	jetfit::model_def& JetFinderAnalyzer::make_model_def(const edm::Event& evt,
65	<	const edm::EventSetup&,
66	<	TH2 *histo) {
67	<	class jf_model_def : public jetfit::model_def {
402	<	public:
403	<	virtual double chisquare_error(double E) {
404	<	return exp(-(4.2 + 0.11*E));
405	<	// study from 09-04-09
406	<	}
407	<	};
408	<
409	<	jf_model_def *_mdef = new jf_model_def();
410	<	TFormula *formula = new TFormula("gaus2d",
411	<	"[0]*exp(-0.5*((x-[1])**2 + (y-[2])**2)/([3]**2))/(2*pi*[3]**2)");
412	<	_mdef->set_formula(formula);
413	<	_mdef->set_indiv_max_E(0);
414	<	_mdef->set_indiv_max_x(1);
415	<	_mdef->set_indiv_max_y(2);
416	<	_mdef->set_indiv_par(0, string("N"), 0.0, 0.0, 0.0, 1.0e6);
417	<	_mdef->set_indiv_par(1, string("mu_x"), 0.0, 0.0, 0.0, 0.0);
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	<	// 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	<
455	<	// generate initial fit histogram
64	>	void JetFinderAnalyzer::analyze_results(HistoFitter::FitResults r,
65	>	std::vector<HistoFitter::Trouble> t,
66	>	TH2 *hist_orig) {
67	>	// perform analysis of fit results
68		edm::Service<TFileService> fs;
69	<	TH2D *init_fit_histo = fs->make<TH2D>(("init_fit_"+string(histo->GetName()))
70	<	.c_str(),
71	<	("Initial fit for "
72	<	+string(histo->GetName())).c_str(),
73	<	histo->GetXaxis()->GetNbins(),
74	<	histo->GetXaxis()->GetXmin(),
75	<	histo->GetXaxis()->GetXmax(),
76	<	histo->GetXaxis()->GetNbins(),
77	<	histo->GetXaxis()->GetXmin(),
78	<	histo->GetXaxis()->GetXmax());
79	<	double XbinSize = (histo->GetXaxis()->GetXmax()
80	<	- histo->GetXaxis()->GetXmin()) /
81	<	static_cast<double>(histo->GetXaxis()->GetNbins());
82	<	double YbinSize = (histo->GetYaxis()->GetXmax()
83	<	- histo->GetYaxis()->GetXmin()) /
84	<	static_cast<double>(histo->GetYaxis()->GetNbins());
85	<	double Xlo = histo->GetXaxis()->GetXmin();
86	<	double Xhi = histo->GetXaxis()->GetXmax();
87	<	double Ylo = histo->GetYaxis()->GetXmin();
88	<	double Yhi = histo->GetYaxis()->GetXmax();
89	<
90	<	for (int i = 0; i < 60; i++) {
91	<	for (int j = 0; j < 60; j++) {
92	<	double x = (static_cast<double>(i) + 0.5)*XbinSize + Xlo;
93	<	double y = (static_cast<double>(j) + 0.5)*YbinSize + Ylo;
482	<	double pval[256];
483	<	if (_mdef->get_n_special_par_sets() > 64) {
484	<	cerr << "Parameter overload" << endl;
485	<	return *_mdef;
486	<	}
487	<	else {
488	<	for (int is = 0; is < _mdef->get_n_special_par_sets(); is++) {
489	<	for (int ii = 0; ii < 4; ii++) {
490	<	double spval, sperr, splo, sphi;
491	<	_mdef->get_special_par(is, ii, spval, sperr, splo, sphi);
492	<	pval[4*is + ii] = spval;
493	<	}
494	<	}
69	>	for (unsigned i = 0; i < r.pars.size(); i++) {
70	>	ostringstream fitHistoOss;
71	>	fitHistoOss << hist_orig->GetName() << "_fit" << i << flush;
72	>	TH2D *fitHisto = fs->make<TH2D>(fitHistoOss.str().c_str(),
73	>	("Fitted distribution to "
74	>	+std::string(hist_orig->GetName())).c_str(),
75	>	hist_orig->GetNbinsX(),
76	>	hist_orig->GetXaxis()->GetXmin(),
77	>	hist_orig->GetXaxis()->GetXmax(),
78	>	hist_orig->GetNbinsY(),
79	>	hist_orig->GetYaxis()->GetXmin(),
80	>	hist_orig->GetYaxis()->GetXmax());
81	>
82	>	double Xlo = fitHisto->GetXaxis()->GetXmin();
83	>	double Xhi = fitHisto->GetXaxis()->GetXmax();
84	>	double Ylo = fitHisto->GetYaxis()->GetXmin();
85	>	double Yhi = fitHisto->GetYaxis()->GetXmax();
86	>	double XbinSize = (Xhi - Xlo) / static_cast<double>(fitHisto->GetNbinsX());
87	>	double YbinSize = (Yhi - Ylo) / static_cast<double>(fitHisto->GetNbinsY());
88	>
89	>	for (int ib = 1; ib <= fitHisto->GetNbinsX(); ib++) {
90	>	for (int jb = 1; jb <= fitHisto->GetNbinsY(); jb++) {
91	>	double x = (static_cast<double>(ib) - 0.5) * XbinSize + Xlo;
92	>	double y = (static_cast<double>(jb) - 0.5) * YbinSize + Ylo;
93	>	fitHisto->SetBinContent(ib, jb, evalFitFunction(r, x, y, i) * XbinSize * YbinSize);
94		}
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)
499	–	* XbinSize * YbinSize);
95		}
96	<	}
97	<
98	<	return *_mdef;
99	<	}
100	<
101	<	void JetFinderAnalyzer::beginJob(const edm::EventSetup &es) {
102	<	ofs.open("jetfindlog.txt", ios::out);
103	<	if (ofs.fail()) {
104	<	cerr << "Opening jetfindlog.txt FAILED" << endl;
105	<	}
106	<	ofs << "Jetfinder log" << endl
512	<	<< "=============" << endl << endl;
513	<	}
514	<
515	<	ostream& operator<<(ostream &out, jetfit::trouble t) {
516	<	string action, error_string;
517	<
518	<	if (t.istat != 3) {
519	<	switch(t.occ) {
520	<	case jetfit::T_NULL:
521	<	action = "Program"; break;
522	<	case jetfit::T_SIMPLEX:
523	<	action = "SIMPLEX"; break;
524	<	case jetfit::T_MIGRAD:
525	<	action = "MIGRAD"; break;
526	<	case jetfit::T_MINOS:
527	<	action = "MINOS"; break;
528	<	default:
529	<	action = "Program"; break;
530	<	}
531	<
532	<	switch (t.istat) {
533	<	case 0:
534	<	error_string = "Unable to calculate error matrix"; break;
535	<	case 1:
536	<	error_string = "Error matrix a diagonal approximation"; break;
537	<	case 2:
538	<	error_string = "Error matrix not positive definite"; break;
539	<	case 3:
540	<	error_string = "Converged successfully"; break;
541	<	default:
542	<	ostringstream oss;
543	<	oss<<"Unknown status code "<<t.istat << flush;
544	<	error_string = oss.str(); break;
96	>
97	>	// save fit results to an ntuple
98	>	ostringstream fitResultsOss;
99	>	fitResultsOss << hist_orig->GetName() << "_results" << i << flush;
100	>	TNtuple *rNtuple = fs->make<TNtuple>(fitResultsOss.str().c_str(),
101	>	("Fit results for "+std::string(hist_orig->GetName())).c_str(),
102	>	"N:mu_x:mu_y:sigma");
103	>	unsigned nGauss = r.pval[i].size() / 4;
104	>	for (unsigned j = 0; j < nGauss; j++) {
105	>	rNtuple->Fill(r.pval[i][4*j], r.pval[i][4*j+1], r.pval[i][4*j+2],
106	>	r.pval[i][4*j+3]);
107		}
108	<
109	<	if (t.occ != jetfit::T_NULL)
110	<	out << action<<" trouble: "<<error_string;
111	<	else
112	<	out << "Not calculated" << endl;
113	<	}
114	<	else {
115	<	out << "Error matrix accurate" << endl;
116	<	}
117	<
118	<	return out;
119	<	}
120	<
121	<	void JetFinderAnalyzer::analyze_results(jetfit::results r,
122	<	std::vector<jetfit::trouble> t,
123	<	TH2 *hist_orig) {
124	<	ofs << "Histogram "<<hist_orig->GetName() << endl;
125	<	for (int i = unique_jets[hist_orig].size() - r.chisquare.size();
126	<	i < unique_jets[hist_orig].size(); i++) {
127	<	int ir = i - unique_jets[hist_orig].size() + r.chisquare.size();
128	<	ofs << "For "<<i+1<<" gaussians: " << endl
129	<	<< t.at(i) << endl;
130	<	ofs << "chisquare="<<r.chisquare.at(ir)
131	<	<< endl;
132	<	ofs << unique_jets[hist_orig][i].size()<<" unique jets found" << endl;
133	<	for (int j = 0; j < unique_jets[hist_orig][i].size(); j++) {
134	<	jet _jet = unique_jets[hist_orig][i][j];
135	<	ofs << "Jet "<<j<<": Energy = "<<_jet.energy<<", eta = "<<_jet.eta
136	<	<< ", phi = "<<_jet.phi << endl;
108	>
109	>	// save jets found
110	>	vector<reco::Jet::LorentzVector> jets;
111	>	for (unsigned j = 0; j < nGauss; j++) {
112	>	double energy = r.pval[i][4*j];
113	>	double eta = r.pval[i][4*j + 1];
114	>	double phi = r.pval[i][4*j + 2];
115	>	double width = r.pval[i][4*j + 3];
116	>	//true jet width
117	>	width = sqrt(width*width - 0.0025);
118	>
119	>	reco::LeafCandidate::LorentzVector P4;
120	>	reco::LeafCandidate::Point vertex(0.0, 0.0, 0.0);
121	>	P4.SetE(energy);
122	>	// SetEta and SetPhi don't seem to be working
123	>	//P4.SetEta(eta);
124	>	//P4.SetPhi(phi);
125	>	double theta = 2.0*atan(exp(-eta));
126	>	P4.SetPx(energy*sin(theta)*cos(phi));
127	>	P4.SetPy(energy*sin(theta)*sin(phi));
128	>	P4.SetPz(energy*cos(theta));
129	>	jets.push_back(P4);
130	>	}
131	>	// compute pair masses
132	>	for (unsigned s1 = 0; s1 < jets.size() - 1; s1++) {
133	>	for (unsigned s2 = s1 + 1; s2 < jets.size(); s2++) {
134	>	// create TOTAL MOMENTUM VECTOR
135	>	reco::Jet::LorentzVector P4 = jets[s1] + jets[s2];
136	>	pairMassHist->Fill(static_cast<double>(P4.M()));
137	>	}
138	>	}
139	>
140	>	// save chisquares to ntuple
141	>	for (unsigned j = 0; j < r.chisquare.size(); j++) {
142	>	ostringstream csNtupleName, csNtupleTitle;
143	>	csNtupleName << hist_orig->GetName() << "_chi2_" << j << flush;
144	>	csNtupleTitle << "Chisquare "<<j<<" for histo "<<hist_orig->GetName()
145	>	<< flush;
146	>	TNtuple *csNtuple = fs->make<TNtuple>(csNtupleName.str().c_str(),
147	>	csNtupleTitle.str().c_str(),
148	>	"chisq");
149	>	csNtuple->Fill(r.chisquare[j]);
150		}
576	–	ofs << endl;
577	–	}
578	–	ofs << endl;
151
152	<	// save fit function histograms to root file
153	<	edm::Service<TFileService> fs;
582	<	for (vector< vector<double> >::size_type i = 0;
583	<	i < r.pval.size(); i++) {
584	<	jetfit::set_ngauss(r.pval[i].size() / 4);
585	<	TF2 *tf2 = new TF2("fit_func", jetfit::fit_fcn_TF2,
586	<	hist_orig->GetXaxis()->GetXmin(),
587	<	hist_orig->GetXaxis()->GetXmax(),
588	<	hist_orig->GetYaxis()->GetXmin(),
589	<	hist_orig->GetYaxis()->GetXmax(),
590	<	r.pval[i].size());
591	<	for (vector<double>::size_type j = 0; j < r.pval[i].size(); j++) {
592	<	tf2->SetParameter(j, r.pval[i][j]);
593	<	}
594	<	ostringstream fit_histo_oss;
595	<	fit_histo_oss << hist_orig->GetName()<<"_fit_"<<i << flush;
596	<	tf2->SetNpx(hist_orig->GetXaxis()->GetNbins());
597	<	tf2->SetNpy(hist_orig->GetYaxis()->GetNbins());
598	<	TH2D *fit_histo = fs->make<TH2D>(fit_histo_oss.str().c_str(),
599	<	fit_histo_oss.str().c_str(),
600	<	hist_orig->GetXaxis()->GetNbins(),
601	<	hist_orig->GetXaxis()->GetXmin(),
602	<	hist_orig->GetXaxis()->GetXmax(),
603	<	hist_orig->GetYaxis()->GetNbins(),
604	<	hist_orig->GetYaxis()->GetXmin(),
605	<	hist_orig->GetYaxis()->GetXmax());
606	<	TH1 *tf2_histo = tf2->CreateHistogram();
607	<	double XbinSize = (fit_histo->GetXaxis()->GetXmax()
608	<	- fit_histo->GetXaxis()->GetXmin())
609	<	/ static_cast<double>(fit_histo->GetXaxis()->GetNbins());
610	<	double YbinSize = (fit_histo->GetYaxis()->GetXmax()
611	<	- fit_histo->GetYaxis()->GetXmin())
612	<	/ static_cast<double>(fit_histo->GetYaxis()->GetNbins());
613	<	for (int ih = 0; ih < tf2->GetNpx(); ih++) {
614	<	for (int jh = 0; jh < tf2->GetNpy(); jh++) {
615	<	fit_histo->SetBinContent(ih+1, jh+1,
616	<	tf2_histo->GetBinContent(ih+1, jh+1)
617	<	* XbinSize * YbinSize);
618	<	}
619	<	}
152	>	// save jets found
153	>
154		}
621	–
622	–	// save results to file
623	–	ostringstream res_tree_oss, rt_title_oss;
624	–	res_tree_oss << hist_orig->GetName()<<"_results" << flush;
625	–	rt_title_oss << "Fit results for "<<hist_orig->GetName() << flush;
155		}
156
157		DEFINE_FWK_MODULE(JetFinderAnalyzer);

Diff Legend

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