StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
TauolaHepMCParticle.cxx
1 #include "TauolaHepMCParticle.h"
2 #include "Log.h"
3 
4 namespace Tauolapp
5 {
6 
8  m_particle = new HepMC::GenParticle();
9 }
10 
11 TauolaHepMCParticle::~TauolaHepMCParticle(){
12 
13  //delete the mother and daughter pointers
14  while(m_mothers.size()!=0){
15  TauolaParticle * temp = m_mothers.back();
16  m_mothers.pop_back();
17  delete temp;
18  }
19  while(m_daughters.size()!=0){
20  TauolaParticle * temp = m_daughters.back();
21  m_daughters.pop_back();
22  delete temp;
23  }
24 
25  while(m_created_particles.size()!=0){
26  TauolaHepMCParticle * temp = (TauolaHepMCParticle*) m_created_particles.back();
27  m_created_particles.pop_back();
28  if(temp->getHepMC()->barcode()==0) delete temp->getHepMC();
29  delete temp;
30  }
31 
32 }
33 
34 // NOTE: Not executed by release examples
35 TauolaHepMCParticle::TauolaHepMCParticle(int pdg_id, int status, double mass){
36  m_particle = new HepMC::GenParticle();
37  m_particle->set_pdg_id(pdg_id);
38  m_particle->set_status(status);
39  m_particle->set_generated_mass(mass);
40 }
41 
43  m_particle = particle;
44 }
45 
47  return m_particle;
48 }
49 
51  std::vector<TauolaParticle*> daughters = getDaughters();
52  std::vector<TauolaParticle*>::iterator dIter = daughters.begin();
53 
54  for(; dIter != daughters.end(); dIter++)
55  (*dIter)->undecay();
56 
57  if(m_particle->end_vertex())
58  {
59  while(m_particle->end_vertex()->particles_out_size())
60  {
61  HepMC::GenParticle *p = m_particle->end_vertex()->remove_particle(*(m_particle->end_vertex()->particles_out_const_begin()));
62  delete p;
63  }
64  delete m_particle->end_vertex();
65  }
66 
67  m_daughters.clear();
69 
70  for(unsigned int i=0;i<daughters.size();i++)
71  delete daughters[i];
72 }
73 
74 void TauolaHepMCParticle::setMothers(vector<TauolaParticle*> mothers){
75 
76  /******** Deal with mothers ***********/
77 
78  //If there are mothers
79  if(mothers.size()>0){
80 
81  HepMC::GenParticle * part;
82  part=dynamic_cast<TauolaHepMCParticle*>(mothers.at(0))->getHepMC();
83 
84  //Use end vertex of first mother as production vertex for particle
85  HepMC::GenVertex * production_vertex = part->end_vertex();
86  HepMC::GenVertex * orig_production_vertex = production_vertex;
87 
88  //If production_vertex does not exist - create it
89  //If it's tau decay - set the time and position including the tau lifetime correction
90  //otherwise - copy the time and position of decaying particle
91  if(!production_vertex){
92  production_vertex = new HepMC::GenVertex();
93  HepMC::FourVector point = part->production_vertex()->position();
94  production_vertex->set_position(point);
95  part->parent_event()->add_vertex(production_vertex);
96  }
97 
98  //Loop over all mothers to check that the end points to the right place
99  vector<TauolaParticle*>::iterator mother_itr;
100  for(mother_itr = mothers.begin(); mother_itr != mothers.end();
101  mother_itr++){
102 
103  HepMC::GenParticle * moth;
104  moth = dynamic_cast<TauolaHepMCParticle*>(*mother_itr)->getHepMC();
105 
106  if(moth->end_vertex()!=orig_production_vertex)
107  Log::Fatal("Mother production_vertices point to difference places. Can not override. Please delete vertices first.",1);
108  else
109  production_vertex->add_particle_in(moth);
110 
111  //update status info
112  if(moth->status()==TauolaParticle::STABLE)
114  }
115  production_vertex->add_particle_out(m_particle);
116  }
117 }
118 
119 void TauolaHepMCParticle::setDaughters(vector<TauolaParticle*> daughters){
120 
121  if(!m_particle->parent_event())
122  Log::Fatal("New particle needs the event set before it's daughters can be added",2);
123 
124  //If there are daughters
125  if(daughters.size()>0){
126  // NOTE: Not executed by release examples
127  // because daughters.size() is always 0
128 
129  //Use production vertex of first daughter as end vertex for particle
130  HepMC::GenParticle * first_daughter;
131  first_daughter = (dynamic_cast<TauolaHepMCParticle*>(daughters.at(0)))->getHepMC();
132 
133  HepMC::GenVertex * end_vertex;
134  end_vertex=first_daughter->production_vertex();
135  HepMC::GenVertex * orig_end_vertex = end_vertex;
136 
137  if(!end_vertex){ //if it does not exist create it
138  end_vertex = new HepMC::GenVertex();
139  m_particle->parent_event()->add_vertex(end_vertex);
140  }
141 
142  //Loop over all daughters to check that the end points to the right place
143  vector<TauolaParticle*>::iterator daughter_itr;
144  for(daughter_itr = daughters.begin(); daughter_itr != daughters.end();
145  daughter_itr++){
146 
147  HepMC::GenParticle * daug;
148  daug = dynamic_cast<TauolaHepMCParticle*>(*daughter_itr)->getHepMC();
149 
150 
151  if(daug->production_vertex()!=orig_end_vertex)
152  Log::Fatal("Daughter production_vertices point to difference places. Can not override. Please delete vertices first.",3);
153  else
154  end_vertex->add_particle_out(daug);
155  }
156  end_vertex->add_particle_in(m_particle);
157  }
158 }
159 
160 std::vector<TauolaParticle*> TauolaHepMCParticle::getMothers(){
161 
162  if(m_mothers.size()==0&&m_particle->production_vertex()){
164  pcle_itr=m_particle->production_vertex()->particles_in_const_begin();
165 
167  pcle_itr_end=m_particle->production_vertex()->particles_in_const_end();
168 
169  for(;pcle_itr != pcle_itr_end; pcle_itr++){
170  m_mothers.push_back(new TauolaHepMCParticle(*pcle_itr));
171  }
172  }
173  return m_mothers;
174 }
175 
176 std::vector<TauolaParticle*> TauolaHepMCParticle::getDaughters(){
177 
178  if(m_daughters.size()==0&&m_particle->end_vertex()){
180  pcle_itr=m_particle->end_vertex()->particles_out_const_begin();
181 
183  pcle_itr_end=m_particle->end_vertex()->particles_out_const_end();
184 
185  for(;pcle_itr != pcle_itr_end; pcle_itr++){
186  m_daughters.push_back(new TauolaHepMCParticle(*pcle_itr));
187  }
188  }
189  return m_daughters;
190 }
191 
193 
194  if(!m_particle->end_vertex()) return;
195 
196  // HepMC version of check_momentum_conservation
197  // with added energy check
198 
199  double sumpx = 0, sumpy = 0, sumpz = 0, sume = 0;
201  part1 != m_particle->end_vertex()->particles_in_const_end(); part1++ ){
202 
203  sumpx += (*part1)->momentum().px();
204  sumpy += (*part1)->momentum().py();
205  sumpz += (*part1)->momentum().pz();
206  sume += (*part1)->momentum().e();
207  }
208 
210  part2 != m_particle->end_vertex()->particles_out_const_end(); part2++ ){
211 
212  sumpx -= (*part2)->momentum().px();
213  sumpy -= (*part2)->momentum().py();
214  sumpz -= (*part2)->momentum().pz();
215  sume -= (*part2)->momentum().e();
216  }
217 
218  if( sqrt( sumpx*sumpx + sumpy*sumpy + sumpz*sumpz + sume*sume) > Tauola::momentum_conservation_threshold ) {
219  Log::Warning()<<"Momentum not conserved in the vertex:"<<endl;
220  Log::RedirectOutput(Log::Warning(false));
221  m_particle->end_vertex()->print();
223  return;
224  }
225 
226  return;
227 }
228 
229 // NOTE: Not executed by release examples
231  m_particle->set_pdg_id(pdg_id);
232 }
233 
235  m_particle->set_generated_mass(mass);
236 }
237 
238 // NOTE: Not executed by release examples
240  m_particle->set_status(status);
241 }
242 
244  return m_particle->pdg_id();
245 }
246 
248  return m_particle->status();
249 }
250 
252  return m_particle->barcode();
253 }
254 
255 // Set (X,T) Position of tau decay trees
257 
258  double lifetime = Tauola::tau_lifetime * (-log( Tauola::randomDouble() ));
259  HepMC::FourVector tau_momentum = m_particle->momentum();
260 
261  double mass = sqrt(abs( tau_momentum.e()*tau_momentum.e()
262  - tau_momentum.px()*tau_momentum.px()
263  - tau_momentum.py()*tau_momentum.py()
264  - tau_momentum.pz()*tau_momentum.pz()
265  ) );
266 
267  // Get previous position
268  HepMC::FourVector previous_position = m_particle->production_vertex()->position();
269 
270  // Calculate new position
271  HepMC::FourVector new_position(previous_position.x()+tau_momentum.px()/mass*lifetime,
272  previous_position.y()+tau_momentum.py()/mass*lifetime,
273  previous_position.z()+tau_momentum.pz()/mass*lifetime,
274  previous_position.t()+tau_momentum.e() /mass*lifetime);
275 
276  // Set new position
277  m_particle->end_vertex()->set_position(new_position);
278  recursiveSetPosition(m_particle,new_position);
279 }
280 
281 void TauolaHepMCParticle::recursiveSetPosition(HepMC::GenParticle *p, HepMC::FourVector pos){
282 
283  if(!p->end_vertex()) return;
284 
285  // Iterate over all outgoing particles
287  pp != p->end_vertex()->particles_out_const_end();
288  ++pp){
289  if( !(*pp)->end_vertex() ) continue;
290 
291  // Set position
292  (*pp)->end_vertex()->set_position(pos);
293  recursiveSetPosition(*pp,pos);
294  }
295 }
296 
298  int pdg_id, int status, double mass,
299  double px, double py, double pz, double e){
300 
301  TauolaHepMCParticle * new_particle = new TauolaHepMCParticle();
302  new_particle->getHepMC()->set_pdg_id(pdg_id);
303  new_particle->getHepMC()->set_status(status);
304  new_particle->getHepMC()->set_generated_mass(mass);
305 
306  HepMC::FourVector momentum(px,py,pz,e);
307  new_particle->getHepMC()->set_momentum(momentum);
308 
309  m_created_particles.push_back(new_particle);
310 
311  return new_particle;
312 }
313 
315  m_particle->print();
316 }
317 
318 
319 /******** Getter and Setter methods: ***********************/
320 
322  return m_particle->momentum().px();
323 }
324 
326  return m_particle->momentum().py();
327 }
328 
330  return m_particle->momentum().pz();
331 }
332 
334  return m_particle->momentum().e();
335 }
336 
338  //make new momentum as something is wrong with
339  //the HepMC momentum setters
340 
341  HepMC::FourVector momentum(m_particle->momentum());
342  momentum.setPx(px);
343  m_particle->set_momentum(momentum);
344 }
345 
347  HepMC::FourVector momentum(m_particle->momentum());
348  momentum.setPy(py);
349  m_particle->set_momentum(momentum);
350 }
351 
352 
354  HepMC::FourVector momentum(m_particle->momentum());
355  momentum.setPz(pz);
356  m_particle->set_momentum(momentum);
357 }
358 
360  HepMC::FourVector momentum(m_particle->momentum());
361  momentum.setE(e);
362  m_particle->set_momentum(momentum);
363 }
364 
365 } // namespace Tauolapp
GenEvent * parent_event() const
pointer to the event that owns this particle
Definition: GenParticle.cc:123
void set_generated_mass(const double &m)
define the actual generated mass
Definition: GenParticle.cc:240
int pdg_id() const
particle ID
Definition: GenParticle.h:214
double z() const
return z
Definition: SimpleVector.h:77
static const int DECAYED
std::vector< HepMC::GenParticle * >::const_iterator particles_in_const_iterator
const iterator for incoming particles
Definition: GenVertex.h:152
void set_status(int status=0)
set decay status
Definition: GenParticle.h:236
GenParticle * remove_particle(GenParticle *particle)
remove a particle
Definition: GenVertex.cc:295
static const int STABLE
Interface to HepMC::GenParticle objects.
void print(std::ostream &ostr=std::cout) const
print vertex information
Definition: GenVertex.cc:145
void add_particle_in(GenParticle *inparticle)
add incoming particle
Definition: GenVertex.cc:273
int barcode() const
particle barcode
Definition: GenParticle.h:252
particles_out_const_iterator particles_out_const_end() const
end iteration of outgoing particles
Definition: GenVertex.h:450
std::vector< HepMC::GenParticle * >::const_iterator particles_out_const_iterator
const iterator for outgoing particles
Definition: GenVertex.h:155
void setDaughters(std::vector< TauolaParticle * > daughters)
double e() const
return E
Definition: SimpleVector.h:73
Abstract base class for particle in the event. This class also handles boosting.
static void RedirectOutput(void(*func)(), ostream &where=*out)
Definition: Log.cxx:93
double y() const
return y
Definition: SimpleVector.h:76
void set_position(const FourVector &position=FourVector(0, 0, 0, 0))
set vertex position and time
Definition: GenVertex.h:424
double px() const
return px
Definition: SimpleVector.h:70
GenVertex * production_vertex() const
pointer to the production vertex
Definition: GenParticle.h:218
GenVertex contains information about decay vertices.
Definition: GenVertex.h:52
std::vector< TauolaParticle * > getMothers()
double pz() const
return pz
Definition: SimpleVector.h:72
void add_particle_out(GenParticle *outparticle)
add outgoing particle
Definition: GenVertex.cc:284
particles_in_const_iterator particles_in_const_end() const
end iteration of incoming particles
Definition: GenVertex.h:440
bool add_vertex(GenVertex *vtx)
adds to evt and adopts
Definition: GenEvent.cc:334
int status() const
HEPEVT decay status.
Definition: GenParticle.h:216
double x() const
return x
Definition: SimpleVector.h:75
static void Fatal(string text, unsigned short int code=0)
int particles_out_size() const
number of outgoing particles
Definition: GenVertex.h:458
particles_out_const_iterator particles_out_const_begin() const
begin iteration of outgoing particles
Definition: GenVertex.h:445
void set_momentum(const FourVector &vec4)
set standard 4 momentum
Definition: GenParticle.h:231
const FourVector & position() const
vertex position and time
Definition: GenVertex.h:406
void print(std::ostream &ostr=std::cout) const
dump this particle&#39;s full info to ostr
Definition: GenParticle.cc:106
const FourVector & momentum() const
standard 4 momentum
Definition: GenParticle.h:211
FourVector is a simple representation of a physics 4 vector.
Definition: SimpleVector.h:42
GenVertex * end_vertex() const
pointer to the decay vertex
Definition: GenParticle.h:221
TauolaHepMCParticle * createNewParticle(int pdg_id, int status, double mass, double px, double py, double pz, double e)
void set_pdg_id(int id)
set particle ID
Definition: GenParticle.h:234
double py() const
return py
Definition: SimpleVector.h:71
particles_in_const_iterator particles_in_const_begin() const
begin iteration of incoming particles
Definition: GenVertex.h:435
std::vector< TauolaParticle * > getDaughters()
double t() const
return t
Definition: SimpleVector.h:78
void setMothers(std::vector< TauolaParticle * > mothers)
static void RevertOutput()
Definition: Log.h:85
The GenParticle class contains information about generated particles.
Definition: GenParticle.h:60