ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/dasu/UltraFastSim/UltraFastSim.cc

Comparing UserCode/dasu/UltraFastSim/UltraFastSim.cc (file contents):
Revision 1.1.1.1 by dasu, Wed Feb 9 05:39:29 2011 UTC vs.
Revision 1.22 by dasu, Thu Jun 13 15:46:28 2013 UTC

#	Line 1 | Line 1
1		#include "UltraFastSim.h"
2
3		#include "Event.h"
4	–
4		#include "Analysis.h"
5
6	+	#include "TParticle.h"
7	+	#include "TLorentzVector.h"
8	+
9		#include "fastjet/PseudoJet.hh"
10		#include "fastjet/JetDefinition.hh"
11		#include "fastjet/ClusterSequence.hh"
12
13	+	#include "UFSDataStore.h"
14	+	#include "UFSDataStore.cc"
15	+
16		using namespace Pythia8;
17		using namespace fastjet;
18		using namespace std;
19
20	<	UltraFastSim::UltraFastSim(Rndm* r) : rndmPtr(r),
21	<	jetDefPtr(0),
22	<	cs(0),
18	<	trackerResolution(0.001),
19	<	ecalResolution(0.01),
20	<	ecalConstantTerm(0.01),
21	<	hcalResolution(1.),
22	<	hcalConstantTerm(0.1)
23	<	{
24	<	jetDefPtr = new JetDefinition(antikt_algorithm, 0.5);
25	<	return;
26	<	}
20	>	Rndm *rndmPtr;
21	>
22	>	bool UltraFastSim::run(Event &event, Rndm *r) {
23
28	–	bool UltraFastSim::run(Event &event) {
24
25		// Clear the previous event
26
27		clear();
28	+	rndmPtr = r;
29
30		// Select particles of interest for later analysis
31
32	+	primaryVertex.SetXYZT(-1000000., -1000000., -1000000., -1000000.);
33		for (int i = 0; i < event.size(); ++i) {
34		Particle& particle = event[i];
35
36	<	// Select generated b quarks in detector acceptance
36	>	// Select generated b and c quarks, and taus with some basic cuts
37	>	// Ignore photon and guon bremmstrahlung
38
39	<	if(abs(particle.id()) == 5 && particle.pT() > 10. && abs(particle.eta()) < 2.5) {
40	<	selectedBQuarks.push_back(PseudoJet(particle.px(), particle.py(), particle.pz(), particle.e()));
39	>	if(abs(particle.id()) == 5 &&
40	>	abs(event[particle.daughter1()].id()) != 5 &&
41	>	particle.pT() > 5. &&
42	>	abs(particle.eta()) < 10.) {
43	>	bQuarks.push_back(TParticle(particle.id(),
44	>	particle.status(),
45	>	particle.mother1(),
46	>	i,   // Mother 2 is mostly useless; Store particle location for later navigation of partial pythia list
47	>	particle.daughter1(),
48	>	particle.daughter2(),
49	>	particle.px(),
50	>	particle.py(),
51	>	particle.pz(),
52	>	particle.e(),
53	>	particle.xProd(),
54	>	particle.yProd(),
55	>	particle.zProd(),
56	>	particle.tProd())
57	>	);
58		}
59
60	<	// Select generated taus in detector acceptance
61	<
62	<	if(abs(particle.id()) == 15 && particle.pT() > 10. && abs(particle.eta()) < 2.5) {
63	<	selectedTaus.push_back(PseudoJet(particle.px(), particle.py(), particle.pz(), particle.e()));
60	>	if(abs(particle.id()) == 4 &&
61	>	abs(event[particle.daughter1()].id()) != 4 &&
62	>	particle.pT() > 5. &&
63	>	abs(particle.eta()) < 10.) {
64	>	cQuarks.push_back(TParticle(particle.id(),
65	>	particle.status(),
66	>	particle.mother1(),
67	>	i,   // Mother 2 is mostly useless; Store particle location for later navigation of partial pythia list
68	>	particle.daughter1(),
69	>	particle.daughter2(),
70	>	particle.px(),
71	>	particle.py(),
72	>	particle.pz(),
73	>	particle.e(),
74	>	particle.xProd(),
75	>	particle.yProd(),
76	>	particle.zProd(),
77	>	particle.tProd())
78	>	);
79	>	}
80	>
81	>	if(abs(particle.id()) == 15 &&
82	>	abs(event[particle.daughter1()].id()) != 15 &&
83	>	particle.pT() > 5. &&
84	>	abs(particle.eta()) < 10.) {
85	>	genTaus.push_back(TParticle(particle.id(),
86	>	particle.status(),
87	>	particle.mother1(),
88	>	i,   // Mother 2 is mostly useless; Store particle location for later navigation of partial pythia list
89	>	particle.daughter1(),
90	>	particle.daughter2(),
91	>	particle.px(),
92	>	particle.py(),
93	>	particle.pz(),
94	>	particle.e(),
95	>	particle.xProd(),
96	>	particle.yProd(),
97	>	particle.zProd(),
98	>	particle.tProd())
99	>	);
100	>	TLorentzVector visP;
101	>	for(int k = particle.daughter1(); k < particle.daughter2(); k++) {
102	>	if(abs(event[k].isVisible())) {
103	>	TLorentzVector kP(event[k].px(), event[k].py(), event[k].pz(), event[k].e());
104	>	visP += kP;
105	>	}
106	>	}
107	>	visTaus.push_back(TParticle(particle.id(),
108	>	particle.status(),
109	>	particle.mother1(),
110	>	i,   // Mother 2 is mostly useless; Store particle location for later navigation of partial pythia list
111	>	particle.daughter1(),
112	>	particle.daughter2(),
113	>	visP.Px(),
114	>	visP.Py(),
115	>	visP.Pz(),
116	>	visP.E(),
117	>	particle.xProd(),
118	>	particle.yProd(),
119	>	particle.zProd(),
120	>	particle.tProd())
121	>	);
122		}
123
124		// Consider only particles with good status
#	Line 55 | Line 128 | bool UltraFastSim::run(Event &event) {
128		// Ignore soft particles and those outside the detector
129		if(particle.pT() > 1.0 && abs(particle.eta()) < 5.0) {
130
131	+	if(primaryVertex.X() < -9999.) {
132	+	primaryVertex.SetXYZT(particle.xProd(), particle.yProd(), particle.zProd(), particle.tProd());
133	+	}
134	+
135	+	TParticle smearedParticle;
136	+	setCommon(particle, smearedParticle);
137	+
138		// Select electrons within detector acceptance
139	<	if(abs(particle.id()) == 11 && particle.pT() > 10. && abs(particle.eta()) < 2.5) {
140	<	selectedElectrons.push_back(PseudoJet(particle.px(), particle.py(), particle.pz(), particle.e()));
139	>	if(abs(particle.id()) == 11 && particle.pT() > 5. && abs(particle.eta()) < 2.5) {
140	>	emSmear(particle, smearedParticle);
141	>	electrons.push_back(smearedParticle);
142		}
143
144		// Select muons within detector acceptance
145	<	if(abs(particle.id()) == 13 && particle.pT() > 10. && abs(particle.eta()) < 2.5) {
146	<	selectedMuons.push_back(PseudoJet(particle.px(), particle.py(), particle.pz(), particle.e()));
145	>	else if(abs(particle.id()) == 13 && particle.pT() > 5. && abs(particle.eta()) < 2.5) {
146	>	float randomNumber = rndmPtr->flat();
147	>	tkSmear(particle, smearedParticle);
148	>	muons.push_back(smearedParticle);
149		}
150	<
150	>
151		// Select other charged tracks and smear them
152	+	else if(abs(particle.charge()) != 0) {
153	+	tkSmear(particle, smearedParticle);
154	+	chargedHadrons.push_back(smearedParticle);
155	+	}
156	+
157		// Select photons and smear them
158	+	else if(particle.id() == 22) {
159	+	emSmear(particle, smearedParticle);
160	+	photons.push_back(smearedParticle);
161	+	}
162	+
163		// Select other neutral particles and smear them
164
165	<	// Select for making jets using fastjet
166	<	selectedParticles.push_back(PseudoJet(particle.px(), particle.py(), particle.pz(), particle.e()));
165	>	else {
166	>	hdSmear(particle, smearedParticle);
167	>	neutralHadrons.push_back(smearedParticle);
168	>	}
169
170		}
171		}
#	Line 79 | Line 174 | bool UltraFastSim::run(Event &event) {
174
175		makeTaus();
176		makeJets();
177	<	makeBJets();
83	<
84	<	cout << "Number of Particles = " << selectedParticles.size() << endl;
85	<	cout << "Number of Gen Elecs = " << selectedElectrons.size() << endl;
86	<	cout << "Number of Gen Muons = " << selectedMuons.size() << endl;
87	<	cout << "Number of Gen Taus =  " << selectedTaus.size() << endl;
88	<	cout << "Number of b Quarks =  " << selectedBQuarks.size() << endl;
89	<	cout << "Number of Jets =      " << sortedJets.size() << endl;
90	<	cout << "Number of bJets =     " << bJets.size() << endl;
91	<	cout << "Number of taus =      " << taus.size() << endl;
92	<	cout << endl;
177	>	makeETSums();
178
179		return true;
180
181		}
182
183		void UltraFastSim::clear() {
184	<	delete cs;
185	<	selectedParticles.clear();
186	<	selectedElectrons.clear();
187	<	selectedMuons.clear();
188	<	selectedTaus.clear();
189	<	selectedBQuarks.clear();
190	<	jets.clear();
106	<	sortedJets.clear();
107	<	bJets.clear();
184	>	genTaus.clear();
185	>	bQuarks.clear();
186	>	cQuarks.clear();
187	>	photons.clear();
188	>	electrons.clear();
189	>	muons.clear();
190	>	visTaus.clear();
191		taus.clear();
192	+	chargedHadrons.clear();
193	+	neutralHadrons.clear();
194	+	jets.clear();
195	+	ET = 0;
196	+	HT = 0;
197	+	MET = TLorentzVector();
198	+	MHT = TLorentzVector();
199		}
200
201		void UltraFastSim::makeJets() {
202	<	cs = new ClusterSequence(selectedParticles, *jetDefPtr);
202	>	vector<PseudoJet> particles;
203	>	for(unsigned int i = 0; i < photons.size(); i++) {
204	>	// Take all neutral particles
205	>	particles.push_back(PseudoJet(photons[i].Px(), photons[i].Py(), photons[i].Pz(), photons[i].Energy()));
206	>	}
207	>	for(unsigned int i = 0; i < electrons.size(); i++) {
208	>	// Take in all electrons within 1 mm of the zVertex
209	>	if(abs(electrons[i].Vz() - primaryVertex.Z()) < 1) {
210	>	particles.push_back(PseudoJet(electrons[i].Px(), electrons[i].Py(), electrons[i].Pz(), electrons[i].Energy()));
211	>	}
212	>	}
213	>	for(unsigned int i = 0; i < chargedHadrons.size(); i++) {
214	>	// Take in all chargedHadrons within 1 mm of the primary vertex, (x,y,z)=0
215	>	if(abs(chargedHadrons[i].Vz() - primaryVertex.Z()) < 1) {
216	>	particles.push_back(PseudoJet(chargedHadrons[i].Px(), chargedHadrons[i].Py(), chargedHadrons[i].Pz(), chargedHadrons[i].Energy()));
217	>	}
218	>	}
219	>	for(unsigned int i = 0; i < neutralHadrons.size(); i++) {
220	>	// Take all neutral particles
221	>	particles.push_back(PseudoJet(neutralHadrons[i].Px(), neutralHadrons[i].Py(), neutralHadrons[i].Pz(), neutralHadrons[i].Energy()));
222	>	}
223	>	ClusterSequence cs(particles, JetDefinition(antikt_algorithm, 0.5));
224		double ptmin = 15.0;
225	<	jets = cs->inclusive_jets(ptmin);
226	<	sortedJets = sorted_by_pt(jets);
116	<	}
117	<
118	<	void UltraFastSim::makeBJets() {
225	>	vector<PseudoJet> allJets = cs.inclusive_jets(ptmin);
226	>	vector<PseudoJet> sortedJets = sorted_by_pt(allJets);
227		for (unsigned int i = 0; i < sortedJets.size(); i++) {
228	<	if(abs(sortedJets[i].rap()) < 2.5) {
229	<	for(unsigned int j = 0; j < selectedBQuarks.size(); j++) {
122	<	float dRap = fabs(sortedJets[i].rap() - selectedBQuarks[j].rap());
123	<	float dPhi = fabs(sortedJets[i].phi() - selectedBQuarks[j].phi());
124	<	float dR = sqrt(dRap*dRap + dPhi*dPhi);
125	<	if(dR < 0.3) {
126	<	bJets.push_back(sortedJets[i]);
127	<	break;
128	<	}
129	<	}
130	<	}
228	>	TLorentzVector jet(sortedJets[i].px(), sortedJets[i].py(), sortedJets[i].pz(), sortedJets[i].e());
229	>	jets.push_back(jet);
230		}
231		}
232
233		void UltraFastSim::makeTaus() {
234	<	for (unsigned int i = 0; i < selectedParticles.size(); i++) {
234	>	for (unsigned int i = 0; i < chargedHadrons.size(); i++) {
235		// Make sure that there is high PT track seed
236	<	if(selectedParticles[i].Et() > 5. && abs(selectedParticles[i].rap()) < 2.5) {
237	<	int nObjectsInInnerCone = 0;
238	<	float isolationEnergy = 0.;
239	<	PseudoJet tau(selectedParticles[i]);
240	<	for(unsigned int j = 0; j < selectedParticles.size(); j++) {
236	>	int nObjectsInInnerCone = 0;
237	>	float isolationEnergy = 0.;
238	>	bool itCouldStillBeATau = true;
239	>	TParticle tau(chargedHadrons[i]);
240	>	int charge = chargedHadrons[i].GetPdgCode()/abs(chargedHadrons[i].GetPdgCode());
241	>	if(chargedHadrons[i].Pt() > 5. && abs(chargedHadrons[i].Eta()) < 2.5) {
242	>	for(unsigned int j = 0; j < chargedHadrons.size(); j++) {
243		if(i != j) {
244	<	if(selectedParticles[i].Et() >= selectedParticles[j].Et()) {
245	<	float dRap = fabs(selectedParticles[i].rap() - selectedParticles[j].rap());
246	<	float dPhi = fabs(selectedParticles[i].phi() - selectedParticles[j].phi());
244	>	if(chargedHadrons[i].Pt() >= chargedHadrons[j].Pt()) {
245	>	float dRap = fabs(chargedHadrons[i].Eta() - chargedHadrons[j].Eta());
246	>	float dPhi = fabs(chargedHadrons[i].Phi() - chargedHadrons[j].Phi());
247	>	if ( dPhi > M_PI ) dPhi = 2. * M_PI - dPhi;
248		float dR = sqrt(dRap*dRap + dPhi*dPhi);
249		// Make sure that there are no more than five objects in the inner 0.3 cone
250		if(dR < 0.3) {
251		nObjectsInInnerCone++;
252	<	if(nObjectsInInnerCone > 5) break;
253	<	tau += selectedParticles[j];
252	>	if(nObjectsInInnerCone > 3) {
253	>	itCouldStillBeATau = false;
254	>	break;
255	>	}
256	>	TLorentzVector tauP;
257	>	tau.Momentum(tauP);
258	>	TLorentzVector chP;
259	>	chargedHadrons[j].Momentum(chP);
260	>	tauP += chP;
261	>	tau.SetMomentum(tauP);
262	>	charge += (chargedHadrons[j].GetPdgCode()/abs(chargedHadrons[j].GetPdgCode()));
263		}
264		else if(dR < 0.5) {
265	<	isolationEnergy += selectedParticles[j].e();
265	>	isolationEnergy += chargedHadrons[j].Energy();
266		}
267		}
268		else {
269	+	itCouldStillBeATau = false;
270		break;
271		}
272		}
273		}
274	+	if(abs(charge) != 1) itCouldStillBeATau = false;
275	+	if(itCouldStillBeATau) {
276	+	for(unsigned int j = 0; j < photons.size(); j++) {
277	+	float dRap = fabs(chargedHadrons[i].Eta() - photons[j].Eta());
278	+	float dPhi = fabs(chargedHadrons[i].Phi() - photons[j].Phi());
279	+	if ( dPhi > M_PI ) dPhi = 2. * M_PI - dPhi;
280	+	float dR = sqrt(dRap*dRap + dPhi*dPhi);
281	+	// Make sure that there are no more than five objects in the inner 0.3 cone
282	+	if(dR < 0.3) {
283	+	nObjectsInInnerCone++;
284	+	if(nObjectsInInnerCone > 5) {
285	+	itCouldStillBeATau = false;
286	+	break;
287	+	}
288	+	TLorentzVector tauP;
289	+	tau.Momentum(tauP);
290	+	TLorentzVector phP;
291	+	photons[j].Momentum(phP);
292	+	tauP += phP;
293	+	tau.SetMomentum(tauP);
294	+	}
295	+	else if(dR < 0.5) {
296	+	isolationEnergy += photons[j].Energy();
297	+	}
298	+	}
299	+	}
300		// Make sure that there is no more than 5% energy in the 0.5
301	<	if(tau.Et() > 15. && (isolationEnergy / tau.e()) < 0.05) {
301	>	if(itCouldStillBeATau && tau.Pt() > 15. && (isolationEnergy / tau.Energy()) < 0.05) {
302	>	tau.SetPdgCode(15 * charge);
303		taus.push_back(tau);
304		}
305		}
306		}
307		}
308	+
309	+	void UltraFastSim::makeETSums() {
310	+	// Use all measured particles
311	+	for (unsigned int i = 0; i < chargedHadrons.size(); i++) {
312	+	MET -= TLorentzVector(chargedHadrons[i].Px(), chargedHadrons[i].Py(), 0, 0);
313	+	ET += chargedHadrons[i].Pt();
314	+	}
315	+	for (unsigned int i = 0; i < neutralHadrons.size(); i++) {
316	+	MET -= TLorentzVector(neutralHadrons[i].Px(), neutralHadrons[i].Py(), 0, 0);
317	+	ET += neutralHadrons[i].Pt();
318	+	}
319	+	for (unsigned int i = 0; i < photons.size(); i++) {
320	+	MET -= TLorentzVector(photons[i].Px(), photons[i].Py(), 0, 0);
321	+	ET += photons[i].Pt();
322	+	}
323	+	for (unsigned int i = 0; i < electrons.size(); i++) {
324	+	MET -= TLorentzVector(electrons[i].Px(), electrons[i].Py(), 0, 0);
325	+	ET += electrons[i].Pt();
326	+	}
327	+	for (unsigned int i = 0; i < muons.size(); i++) {
328	+	MET -= TLorentzVector(muons[i].Px(), muons[i].Py(), 0, 0);
329	+	ET += muons[i].Pt();
330	+	}
331	+	// Use jets [includes photons and electrons] and muons
332	+	for (unsigned int i = 0; i < jets.size(); i++) {
333	+	MHT -= TLorentzVector(jets[i].Px(), jets[i].Py(), 0, 0);
334	+	HT += jets[i].Et();
335	+	}
336	+	}
337	+
338	+	void UltraFastSim::setCommon(Particle& pyParticle, TParticle& smParticle) {
339	+	smParticle.SetPdgCode(pyParticle.id());
340	+	smParticle.SetStatusCode(pyParticle.status());
341	+	smParticle.SetFirstMother(pyParticle.mother1());
342	+	smParticle.SetLastMother(pyParticle.mother2());
343	+	smParticle.SetFirstDaughter(pyParticle.daughter1());
344	+	smParticle.SetLastDaughter(pyParticle.daughter2());
345	+	smParticle.SetCalcMass(pyParticle.m());
346	+	smParticle.SetProductionVertex(pyParticle.xProd(), pyParticle.yProd(), pyParticle.zProd(), pyParticle.tProd());
347	+	}
348	+
349	+	void UltraFastSim::tkSmear(Particle& pyParticle, TParticle& smParticle) {
350	+	if(abs(pyParticle.eta()) < 2.5) {
351	+	double sigma = pow((0.15 * pyParticle.pT() / 1000) * (0.15 * pyParticle.pT() / 1000) + (0.005 * 0.005), 0.5);
352	+	double smear = (1. + rndmPtr->gauss() * sigma);
353	+	smParticle.SetMomentum(pyParticle.px()*smear, pyParticle.py()*smear, pyParticle.pz(), pyParticle.e());
354	+	}
355	+	else {
356	+	hdSmear(pyParticle, smParticle);
357	+	}
358	+	}
359	+
360	+	void UltraFastSim::emSmear(Particle& pyParticle, TParticle& smParticle) {
361	+	if(abs(pyParticle.eta()) < 3.0) {
362	+	double sqrte = pow(pyParticle.e(), 0.5);
363	+	double rsqr = ((0.027 / sqrte) * (0.027 / sqrte)) + 0.005 * 0.005 + (0.150 / pyParticle.e()) * (0.150 / pyParticle.e());
364	+	double sigma = pow((pyParticle.e()*rsqr), 0.5);
365	+	double smear = (1. + rndmPtr->gauss() * sigma);
366	+	smParticle.SetMomentum(pyParticle.px()*smear, pyParticle.py()*smear, pyParticle.pz()*smear, pyParticle.e()*smear);
367	+	}
368	+	else {
369	+	hdSmear(pyParticle, smParticle);
370	+	}
371	+	}
372	+
373	+	void UltraFastSim::hdSmear(Particle& pyParticle, TParticle& smParticle) {
374	+	double sqrte = pow(pyParticle.e(), 0.5);
375	+	double rsqr = ((1.15 / sqrte) * (1.15 / sqrte)) + 0.055 * 0.055;
376	+	double sigma = pow((pyParticle.e()*rsqr), 0.5);
377	+	double smear = (1. + rndmPtr->gauss() * sigma);
378	+	smParticle.SetMomentum(pyParticle.px()*smear, pyParticle.py()*smear, pyParticle.pz()*smear, pyParticle.e()*smear);
379	+	}

Diff Legend

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