ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/JetFitAnalyzer/src/JetFitAnalyzer.cc

Comparing UserCode/JetFitAnalyzer/src/JetFitAnalyzer.cc (file contents):
Revision 1.2 by dnisson, Wed Oct 14 02:46:26 2009 UTC vs.
Revision 1.23 by dnisson, Thu Apr 8 05:01:28 2010 UTC

#	Line 1 | Line 1
1	<	// -*- C++ -*-
1	>	//*- C++ -*-
2		//
3		// Package:    JetFitAnalyzer
4		// Class:      JetFitAnalyzer
5		//
6		/**\class JetFitAnalyzer JetFitAnalyzer.cc UserCode/JetFitAnalyzer/src/JetFitAnalyzer.cc
7
8	<	Description: Base class for using jetfit with CMSSW
8	>	Defaription: Base class for using jetfit with CMSSW
9
10		Implementation:
11		o Method make_histo is used to prepare a histogram for analysis;
#	Line 25 | Line 25
25		//
26		//
27
28	+	#include <HepMC/GenEvent.h>
29	+	#include <HepMC/GenVertex.h>
30	+	#include <HepMC/GenParticle.h>
31	+
32		#include "UserCode/JetFitAnalyzer/interface/JetFitAnalyzer.h"
33
34	<	#include <vector>
34	>	#include "DataFormats/JetReco/interface/PFJetCollection.h"
35	>	#include "DataFormats/JetReco/interface/GenJetCollection.h"
36	>	#include "DataFormats/JetReco/interface/PFJet.h"
37	>	#include "DataFormats/JetReco/interface/GenJet.h"
38	>	#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
39	>
40	>	#include "FWCore/ServiceRegistry/interface/Service.h"
41	>	#include "FWCore/MessageLogger/interface/MessageLogger.h"
42	>	#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
43	>
44	>	#include "TFormula.h"
45	>	#define PI 3.141592
46	>
47	>	#include <iostream>
48	>	#include <sstream>
49	>	#include <string>
50	>	#include <fstream>
51	>	// to save the W decay product energy, mass, eta, phi
52	>
53	>	using namespace std;
54
55		JetFitAnalyzer::JetFitAnalyzer(const edm::ParameterSet& iConfig)
33	–	: user_minuit(0)
56		{
57		// get parameters from iConfig
36	–	ignorezero_ = iConfig.getUntrackedParameter("ignorezero", false);
37	–	rebinX_ = iConfig.getUntrackedParameter("rebinX", 1);
38	–	rebinY_ = iConfig.getUntrackedParameter("rebinY", 1);
58		P_cutoff_val_ = iConfig.getUntrackedParameter("P_cutoff_val", 0.5);
59	<
60	<	jetfit::set_ignorezero(ignorezero_);
59	>	infoType_ = iConfig.getUntrackedParameter("infoType", 1);
60	>	smear_ = iConfig.getUntrackedParameter("smear", 0.05);
61	>	jet_algo_ = iConfig.getParameter<string>("jet_algo");
62	>	maxGaussians_ = iConfig.getUntrackedParameter("maxGaussians", 4);
63	>	gausSpread_ = iConfig.getUntrackedParameter("gausSpread", 0.1);
64		}
65
66
#	Line 61 | Line 83 | JetFitAnalyzer::analyze(const edm::Event
83		{
84		using namespace edm;
85
86	<	jetfit::results r; std::vector<jetfit::trouble> t;
86	>	HistoFitter histoFitter;
87	>	std::vector<HistoFitter::Trouble> t; // record troubles fitting
88
89	<	TH2D *histo = make_histo(iEvent, iSetup);
90	<	if (histo != 0) {
91	<	jetfit::model_def &_mdef = make_model_def(iEvent, iSetup, histo);
92	<	jetfit::set_model_def(&_mdef);
93	<	r = jetfit::fit_histo(histo, t, user_minuit, _mdef.get_n_special_par_sets(),
94	<	rebinX_, rebinY_, P_cutoff_val_);
95	<	analyze_results(r, t, histo);
96	<	}
97	<	else {
98	<	cerr << "Fitting not performed" << endl;
89	>	std::vector<unsigned> numberOfParticles;
90	>	std::vector<TH2D *> histos = make_histo(iEvent, iSetup, numberOfParticles);
91	>	// create the histogrms
92	>	for (unsigned i= 0; i < histos.size(); i++) {
93	>	TH2D *histo = histos[i]; // obtain the histogram
94	>	if (histo != 0) {
95	>	HistoFitter::ModelDefinition &_mdef = make_model_def(iEvent, iSetup, histo);
96	>	histoFitter.set_ModelDefinition(&_mdef);
97	>	HistoFitter::FitResults r(histoFitter.fit_histo(histo, t, 1,
98	>	P_cutoff_val_,
99	>	numberOfParticles[i]));
100	>	// fit the histogram
101	>	analyze_results(r, t, histo, iEvent); // perform user analysis
102	>	}
103	>	else {
104	>	std::cerr << "Fitting not performed" << std::endl;
105	>	}
106		}
107		}
108
109	+	/* create histograms of the particles in all jets found */
110	+	std::vector<TH2D *> JetFitAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&, std::vector<unsigned> &numberOfParticles) {
111	+	std::vector<TH2D*> histos;
112	+
113	+	std::vector<const reco::Jet *> jets;
114	+	if (infoType_ == 0) {
115	+	edm::Handle< vector<reco::GenJet> > evtJets;
116	+	evt.getByLabel(jet_algo_, evtJets);
117	+	jets.resize(evtJets->size());
118	+	for (unsigned i = 0; i < evtJets->size(); i++) {
119	+	jets[i] = &((*evtJets)[i]);
120	+	}
121	+	}
122	+	else if (infoType_ == 1) {
123	+	edm::Handle< vector<reco::PFJet> > evtJets;
124	+	evt.getByLabel(jet_algo_, evtJets);
125	+	jets.resize(evtJets->size());
126	+	for (unsigned i = 0; i < evtJets->size(); i++) {
127	+	jets[i] = &((*evtJets)[i]);
128	+	}
129	+	}
130	+
131	+	histos.resize(jets.size());
132	+	numberOfParticles.resize(jets.size());
133	+
134	+	for (unsigned i = 0; i < jets.size(); i++) {
135	+	/* Get the particles in the fat jet */
136	+	vector<const reco::Candidate *> particles =
137	+	jets[i]->getJetConstituentsQuick();
138	+
139	+	/* Create the histogram name */
140	+	ostringstream oss;
141	+	oss << "eta_phi_energy"<<evt.id().event() <<"jet"<<i << flush;
142	+
143	+	edm::Service<TFileService> fs;
144	+
145	+	/* Initialize the histogram */
146	+	TH2D *histo = fs->make<TH2D>(oss.str().c_str(), oss.str().c_str(),
147	+	30,
148	+	jets[i]->eta()-0.5,
149	+	jets[i]->eta()+0.5,
150	+	30,
151	+	jets[i]->phi()-0.5,
152	+	jets[i]->phi()+0.5);
153	+	numberOfParticles[i] = particles.size();
154	+	if (smear_ > 0.0) {
155	+	// get the limits of the histogram
156	+	double Xlo = histo->GetXaxis()->GetXmin();
157	+	double Xhi = histo->GetXaxis()->GetXmax();
158	+	double Ylo = histo->GetYaxis()->GetXmin();
159	+	double Yhi = histo->GetYaxis()->GetXmax();
160	+
161	+	// get the bin sizes
162	+	double XbinSize = (Xhi - Xlo) / static_cast<double>(histo->GetNbinsX());
163	+	double YbinSize = (Yhi - Ylo) / static_cast<double>(histo->GetNbinsY());
164	+
165	+	for (int i = 0; i < particles.size(); i++) {
166	+	// compute the parameters of the Gaussian to be added
167	+	double N = particles[i]->energy();
168	+	double x = particles[i]->eta();
169	+	double y = particles[i]->phi();
170	+
171	+	// Compute correct phi for this histogram if it is out of range
172	+	if (y < Ylo) {
173	+	y += 2.0*PI;
174	+	}
175	+	if (y > Yhi) {
176	+	y -= 2.0*PI;
177	+	}
178	+
179	+	// loop over bins and add Gaussian in eta and phi to histogram
180	+	for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) {
181	+	for (vector< double >::size_type j2 = 0; j2 < 30; j2++) {
182	+	// compute the distance of the bin from the particle eta, phi
183	+	double eta = static_cast<double>((signed int)i2) * XbinSize +
184	+	Xlo - x;
185	+	double phi = static_cast<double>((signed int)j2) * YbinSize +
186	+	Ylo - y;
187	+
188	+	/* Fill the histogram with a Gaussian centered around particle */
189	+	histo->SetBinContent(i2+1, j2+1, histo->GetBinContent(i2+1, j2+1)
190	+	+ (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
191	+	* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_)));
192	+	}
193	+	}
194	+	}
195	+	}
196	+	else {
197	+	// don't smear--just fill with particles
198	+	for (int i = 0; i < particles.size(); i++) {
199	+	histo->Fill(particles[i]->eta(), particles[i]->phi(),
200	+	particles[i]->energy());
201	+	}
202	+	}
203	+	histos[i] = histo;
204	+	}
205	+
206	+	return histos;
207	+	}
208	+
209	+	/* Create a model definition based on the histogram and event information */
210	+	HistoFitter::ModelDefinition& JetFitAnalyzer::make_model_def(const edm::Event& evt,
211	+	const edm::EventSetup&,
212	+	TH2 *histo) {
213	+	class DefaultModelDef : public HistoFitter::ModelDefinition {
214	+	public:
215	+	virtual double chisquareSigma(double E, double Emax) {
216	+	double sigma = 0.14*Emax*(1.0 - exp(-2.3*E/Emax));
217	+
218	+	return sigma;
219	+	// based on Nisson Elog entry 33
220	+	}
221	+	};
222	+
223	+	DefaultModelDef *_mdef = new DefaultModelDef();
224	+
225	+	// 2D Gaussian formula
226	+	TFormula *formula = new TFormula("gaus2d",
227	+	"[0]*exp(-0.5*((x-[1])**2 + (y-[2])**2)/([3]**2))/(2*pi*[3]**2)");
228	+	_mdef->setFormula(formula);
229	+
230	+	_mdef->setIndivEnergy(0);
231	+	_mdef->setIndivMeanX(1);
232	+	_mdef->setIndivMeanY(2);
233	+	_mdef->setIndivParameter(0, string("N"), 0.0, 0.0, 0.0, 1.0e6);
234	+	_mdef->setIndivParameter(1, string("mu_x"), 0.0, 0.0, 0.0, 0.0);
235	+	_mdef->setIndivParameter(2, string("mu_y"), 0.0, 0.0, 0.0, 0.0);
236	+
237	+	// set fixed Gaussian sigma equal to gausSpread
238	+	_mdef->setIndivParameter(3, string("sig"), gausSpread_, 0.001, 0.0, 0.0);
239	+
240	+	_mdef->setMaxNumberOfGaussians(maxGaussians_);
241	+
242	+	return *_mdef;
243	+	}
244
245		// ------------ method called once each job just before starting event loop  ------------
246		void

Diff Legend

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