C++ Interface to Tauola
TauolaHepMC3Particle.cxx
1#include "TauolaHepMC3Particle.h"
2#include "Log.h"
3
4#include "HepMC3/GenVertex.h"
5#include "HepMC3/Print.h"
6
7
8namespace Tauolapp
9{
10
12 m_particle = make_shared<GenParticle>();
13}
14
15TauolaHepMC3Particle::~TauolaHepMC3Particle(){
16
17 //delete the mother and daughter pointers
18 while(m_mothers.size()!=0){
19 TauolaParticle * temp = m_mothers.back();
20 m_mothers.pop_back();
21 delete temp;
22 }
23 while(m_daughters.size()!=0){
24 TauolaParticle * temp = m_daughters.back();
25 m_daughters.pop_back();
26 delete temp;
27 }
28
29 while(m_created_particles.size()!=0){
31 m_created_particles.pop_back();
32 //if(temp->getHepMC3()->id()==0) delete temp->getHepMC3();
33 delete temp;
34 }
35
36}
37
38// NOTE: Not executed by release examples
39TauolaHepMC3Particle::TauolaHepMC3Particle(int pdg_id, int status, double mass){
40 m_particle = make_shared<GenParticle>();
41 m_particle->set_pid(pdg_id);
42 m_particle->set_status(status);
43 m_particle->set_generated_mass(mass);
44}
45
47 m_particle = particle;
48}
49
51 return m_particle;
52}
53
55
56 std::vector<TauolaParticle*> daughters = getDaughters();
57 std::vector<TauolaParticle*>::iterator dIter = daughters.begin();
58
59 for(; dIter != daughters.end(); dIter++)
60 (*dIter)->undecay();
61
62 if(m_particle->end_vertex())
63 {
64 while(m_particle->end_vertex()->particles_out().size())
65 {
66 m_particle->end_vertex()->remove_particle_out(m_particle->end_vertex()->particles_out().front());
67 }
68 }
69
70 m_daughters.clear();
72
73 for(unsigned int i=0;i<daughters.size();i++)
74 delete daughters[i];
75}
76
77void TauolaHepMC3Particle::setMothers(vector<TauolaParticle*> mothers){
78
79 /******** Deal with mothers ***********/
80
81 //If there are mothers
82 if(mothers.size()>0){
83
84 GenParticlePtr part;
85 part=dynamic_cast<TauolaHepMC3Particle*>(mothers.at(0))->getHepMC3();
86
87 //Use end vertex of first mother as production vertex for particle
88 GenVertexPtr production_vertex = part->end_vertex();
89 GenVertexPtr orig_production_vertex = production_vertex;
90
91 //If production_vertex does not exist - create it
92 //If it's tau decay - set the time and position including the tau lifetime correction
93 //otherwise - copy the time and position of decaying particle
94 if(!production_vertex){
95 production_vertex = make_shared<GenVertex>();
96 FourVector point = part->production_vertex()->position();
97 production_vertex->set_position(point);
98 part->parent_event()->add_vertex(production_vertex);
99 }
100
101 //Loop over all mothers to check that the end points to the right place
102 vector<TauolaParticle*>::iterator mother_itr;
103 for(mother_itr = mothers.begin(); mother_itr != mothers.end();
104 mother_itr++){
105
106 GenParticlePtr moth;
107 moth = dynamic_cast<TauolaHepMC3Particle*>(*mother_itr)->getHepMC3();
108
109 if(moth->end_vertex()!=orig_production_vertex)
110 Log::Fatal("Mother production_vertices point to difference places. Can not override. Please delete vertices first.",1);
111 else
112 production_vertex->add_particle_in(moth);
113
114 //update status info
115 if(moth->status()==TauolaParticle::STABLE)
116 moth->set_status(TauolaParticle::DECAYED);
117 }
118 production_vertex->add_particle_out(m_particle);
119 }
120}
121
122void TauolaHepMC3Particle::setDaughters(vector<TauolaParticle*> daughters){
123
124 if(!m_particle->parent_event())
125 Log::Fatal("New particle needs the event set before it's daughters can be added",2);
126
127 //If there are daughters
128 if(daughters.size()>0){
129 // NOTE: Not executed by release examples
130 // because daughters.size() is always 0
131
132 //Use production vertex of first daughter as end vertex for particle
133 GenParticlePtr first_daughter;
134 first_daughter = (dynamic_cast<TauolaHepMC3Particle*>(daughters.at(0)))->getHepMC3();
135
136 GenVertexPtr end_vertex;
137 end_vertex=first_daughter->production_vertex();
138 GenVertexPtr orig_end_vertex = end_vertex;
139
140 if(!end_vertex){ //if it does not exist create it
141 end_vertex = make_shared<GenVertex>();
142 m_particle->parent_event()->add_vertex(end_vertex);
143 }
144
145 //Loop over all daughters to check that the end points to the right place
146 vector<TauolaParticle*>::iterator daughter_itr;
147 for(daughter_itr = daughters.begin(); daughter_itr != daughters.end();
148 daughter_itr++){
149
150 GenParticlePtr daug;
151 daug = dynamic_cast<TauolaHepMC3Particle*>(*daughter_itr)->getHepMC3();
152
153
154 if(daug->production_vertex()!=orig_end_vertex)
155 Log::Fatal("Daughter production_vertices point to difference places. Can not override. Please delete vertices first.",3);
156 else
157 end_vertex->add_particle_out(daug);
158 }
159 end_vertex->add_particle_in(m_particle);
160 }
161}
162
163std::vector<TauolaParticle*> TauolaHepMC3Particle::getMothers(){
164
165 if(m_mothers.size()==0&&m_particle->production_vertex()){
166 for(auto p: m_particle->production_vertex()->particles_in() ) {
167 m_mothers.push_back(new TauolaHepMC3Particle(p));
168 }
169 }
170 return m_mothers;
171}
172
173std::vector<TauolaParticle*> TauolaHepMC3Particle::getDaughters(){
174
175 if(m_daughters.size()==0&&m_particle->end_vertex()){
176 for(auto p: m_particle->end_vertex()->particles_out() ) {
177 m_daughters.push_back(new TauolaHepMC3Particle(p));
178 }
179 }
180 return m_daughters;
181}
182
184
185 if(!m_particle->end_vertex()) return;
186
187 // HepMC version of check_momentum_conservation
188 // with added energy check
189
190 FourVector sum;
191 for(ConstGenParticlePtr p: m_particle->end_vertex()->particles_in() ) {
192 sum += p->momentum();
193 }
194
195 for(ConstGenParticlePtr p: m_particle->end_vertex()->particles_out() ) {
196 sum -= p->momentum();
197 }
198
199 if( sum.length() > Tauola::momentum_conservation_threshold ) {
200 Log::Warning()<<"Momentum not conserved in the vertex:"<<endl;
201 Log::RedirectOutput(Log::Warning(false));
202 Print::line(m_particle->end_vertex());
204 return;
205 }
206
207 return;
208}
209
210// NOTE: Not executed by release examples
212 m_particle->set_pid(pdg_id);
213}
214
216 m_particle->set_generated_mass(mass);
217}
218
219// NOTE: Not executed by release examples
221 m_particle->set_status(status);
222}
223
225 return m_particle->pid();
226}
227
229 return m_particle->status();
230}
231
233 return m_particle->pid();
234}
235
236// Set (X,T) Position of tau decay trees
238
239 double lifetime = Tauola::tau_lifetime * (-log( Tauola::randomDouble() ));
240 FourVector tau_momentum = m_particle->momentum();
241
242 double mass = sqrt(abs( tau_momentum.e()*tau_momentum.e()
243 - tau_momentum.px()*tau_momentum.px()
244 - tau_momentum.py()*tau_momentum.py()
245 - tau_momentum.pz()*tau_momentum.pz()
246 ) );
247
248 // Get previous position
249 FourVector previous_position = m_particle->production_vertex()->position();
250
251 // Calculate new position
252 FourVector new_position(previous_position.x()+tau_momentum.px()/mass*lifetime,
253 previous_position.y()+tau_momentum.py()/mass*lifetime,
254 previous_position.z()+tau_momentum.pz()/mass*lifetime,
255 previous_position.t()+tau_momentum.e() /mass*lifetime);
256
257 // Set new position
258 m_particle->end_vertex()->set_position(new_position);
259 recursiveSetPosition(m_particle,new_position);
260}
261
262void TauolaHepMC3Particle::recursiveSetPosition(GenParticlePtr p, FourVector pos){
263
264 if(!p->end_vertex()) return;
265
266 // Iterate over all outgoing particles
267 for(auto pp: p->end_vertex()->particles_out() ) {
268 if( !pp->end_vertex() ) continue;
269
270 // Set position
271 pp->end_vertex()->set_position(pos);
272 recursiveSetPosition(pp,pos);
273 }
274}
275
277 int pdg_id, int status, double mass,
278 double px, double py, double pz, double e){
279
280 TauolaHepMC3Particle * new_particle = new TauolaHepMC3Particle();
281 new_particle->getHepMC3()->set_pid(pdg_id);
282 new_particle->getHepMC3()->set_status(status);
283 new_particle->getHepMC3()->set_generated_mass(mass);
284
285 FourVector momentum(px,py,pz,e);
286 new_particle->getHepMC3()->set_momentum(momentum);
287
288 m_created_particles.push_back(new_particle);
289
290 return new_particle;
291}
292
294 Print::line(m_particle);
295}
296
297
298/******** Getter and Setter methods: ***********************/
299
301 return m_particle->momentum().px();
302}
303
305 return m_particle->momentum().py();
306}
307
309 return m_particle->momentum().pz();
310}
311
313 return m_particle->momentum().e();
314}
315
317 //make new momentum as something is wrong with
318 //the HepMC momentum setters
319
320 FourVector momentum(m_particle->momentum());
321 momentum.setPx(px);
322 m_particle->set_momentum(momentum);
323}
324
326 FourVector momentum(m_particle->momentum());
327 momentum.setPy(py);
328 m_particle->set_momentum(momentum);
329}
330
331
333 FourVector momentum(m_particle->momentum());
334 momentum.setPz(pz);
335 m_particle->set_momentum(momentum);
336}
337
339 FourVector momentum(m_particle->momentum());
340 momentum.setE(e);
341 m_particle->set_momentum(momentum);
342}
343
344} // namespace Tauolapp
Interface to GenParticle objects.
Abstract base class for particle in the event. This class also handles boosting.
static void Fatal(string text, unsigned short int code=0)
static void RevertOutput()
static void RedirectOutput(void(*func)(), ostream &where=*out)
Definition: Log.cxx:93
std::vector< TauolaParticle * > getDaughters()
void setMothers(std::vector< TauolaParticle * > mothers)
std::vector< TauolaParticle * > m_created_particles
void setDaughters(std::vector< TauolaParticle * > daughters)
void recursiveSetPosition(GenParticlePtr p, FourVector pos)
std::vector< TauolaParticle * > getMothers()
TauolaHepMC3Particle * createNewParticle(int pdg_id, int status, double mass, double px, double py, double pz, double e)