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 
8 namespace Tauolapp
9 {
10 
12  m_particle = make_shared<GenParticle>();
13 }
14 
15 TauolaHepMC3Particle::~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
39 TauolaHepMC3Particle::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 
77 void 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 
122 void 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 
163 std::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 
173 std::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 
262 void 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)