TEvtGen/Tauola/TauolaHepMCParticle.cxx

   1 #include "TauolaHepMCParticle.h"
   2 #include "TauolaLog.h"
   3
   4 namespace Tauolapp
   5 {
   6
   7 TauolaHepMCParticle::TauolaHepMCParticle(){
   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
  42 TauolaHepMCParticle::TauolaHepMCParticle(HepMC::GenParticle * particle){
  43   m_particle = particle;
  44 }
  45
  46 HepMC::GenParticle * TauolaHepMCParticle::getHepMC(){
  47   return m_particle;
  48 }
  49
  50 void TauolaHepMCParticle::undecay(){
  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();
  68   m_particle->set_status(TauolaParticle::STABLE);
  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)
 113         moth->set_status(TauolaParticle::DECAYED);
 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()){
 163     HepMC::GenVertex::particles_in_const_iterator pcle_itr;
 164     pcle_itr=m_particle->production_vertex()->particles_in_const_begin();
 165
 166     HepMC::GenVertex::particles_in_const_iterator pcle_itr_end;
 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()){
 179     HepMC::GenVertex::particles_out_const_iterator pcle_itr;
 180     pcle_itr=m_particle->end_vertex()->particles_out_const_begin();
 181
 182     HepMC::GenVertex::particles_out_const_iterator pcle_itr_end;
 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
 192 void TauolaHepMCParticle::checkMomentumConservation(){
 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;
 200   for( HepMC::GenVertex::particles_in_const_iterator part1 = m_particle->end_vertex()->particles_in_const_begin();
 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
 209   for( HepMC::GenVertex::particles_out_const_iterator part2 = m_particle->end_vertex()->particles_out_const_begin();
 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();
 222     Log::RevertOutput();
 223     return;
 224   }
 225
 226   return;
 227 }
 228
 229 // NOTE: Not executed by release examples
 230 void TauolaHepMCParticle::setPdgID(int pdg_id){
 231   m_particle->set_pdg_id(pdg_id);
 232 }
 233
 234 void TauolaHepMCParticle::setMass(double mass){
 235   m_particle->set_generated_mass(mass);
 236 }
 237
 238 // NOTE: Not executed by release examples
 239 void TauolaHepMCParticle::setStatus(int status){
 240   m_particle->set_status(status);
 241 }
 242
 243 int TauolaHepMCParticle::getPdgID(){
 244   return m_particle->pdg_id();
 245 }
 246
 247 int TauolaHepMCParticle::getStatus(){
 248   return m_particle->status();
 249 }
 250
 251 int TauolaHepMCParticle::getBarcode(){
 252   return m_particle->barcode();
 253 }
 254
 255 // Set (X,T) Position of tau decay trees
 256 void TauolaHepMCParticle::decayEndgame(){
 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
 286   for(HepMC::GenVertex::particles_out_const_iterator pp = p->end_vertex()->particles_out_const_begin();
 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
 297 TauolaHepMCParticle * TauolaHepMCParticle::createNewParticle(
 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
 314 void TauolaHepMCParticle::print(){
 315   m_particle->print();
 316 }
 317
 318
 319 /******** Getter and Setter methods: ***********************/
 320
 321 inline double TauolaHepMCParticle::getPx(){
 322   return m_particle->momentum().px();
 323 }
 324
 325 inline double TauolaHepMCParticle::getPy(){
 326   return m_particle->momentum().py();
 327 }
 328
 329 double TauolaHepMCParticle::getPz(){
 330   return m_particle->momentum().pz();
 331 }
 332
 333 double TauolaHepMCParticle::getE(){
 334   return m_particle->momentum().e();
 335 }
 336
 337 void TauolaHepMCParticle::setPx(double px){
 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
 346 void TauolaHepMCParticle::setPy(double py){
 347   HepMC::FourVector momentum(m_particle->momentum());
 348   momentum.setPy(py);
 349   m_particle->set_momentum(momentum);
 350 }
 351
 352
 353 void TauolaHepMCParticle::setPz(double pz){
 354   HepMC::FourVector momentum(m_particle->momentum());
 355   momentum.setPz(pz);
 356   m_particle->set_momentum(momentum);
 357 }
 358
 359 void TauolaHepMCParticle::setE(double e){
 360   HepMC::FourVector momentum(m_particle->momentum());
 361   momentum.setE(e);
 362   m_particle->set_momentum(momentum);
 363 }
 364
 365 } // namespace Tauolapp