TEvtGen/Tauola/TauolaHEPEVTParticle.cxx

   1 #include "TauolaHEPEVTParticle.h"
   2
   3 #include "Log.h"
   4
   5 namespace Tauolapp
   6 {
   7
   8 TauolaHEPEVTParticle::~TauolaHEPEVTParticle()
   9 {
  10   // Cleanup particles that do not belong to event
  11   for(unsigned int i=0;i<cache.size();i++)
  12     if(cache[i]->m_barcode<0)
  13       delete cache[i];
  14 }
  15
  16 TauolaHEPEVTParticle::TauolaHEPEVTParticle(int pdgid, int status, double px, double py, double pz, double e, double m, int ms, int me, int ds, int de){
  17   m_px = px;
  18   m_py = py;
  19   m_pz = pz;
  20   m_e  = e;
  21   m_generated_mass = m;
  22
  23   m_pdgid  = pdgid;
  24   m_status = status;
  25
  26   m_first_mother   = ms;
  27   m_second_mother  = me;
  28   m_daughter_start = ds;
  29   m_daughter_end   = de;
  30
  31   m_barcode = -1;
  32   m_event = NULL;
  33 }
  34
  35 void TauolaHEPEVTParticle::undecay(){
  36
  37   Log::Info()<<"TauolaHEPEVTParticle::undecay not implemented for HEPEVT"<<endl;
  38 }
  39
  40 void TauolaHEPEVTParticle::setMothers(vector<TauolaParticle*> mothers){
  41
  42   // If this particle has not yet been added to the event record
  43   // then add it to the mothers' event record
  44   if(m_barcode<0 && mothers.size()>0)
  45   {
  46     TauolaHEPEVTEvent *evt = ((TauolaHEPEVTParticle*)mothers[0])->m_event;
  47     evt->addParticle(this);
  48   }
  49
  50   if(mothers.size()>2) Log::Fatal("TauolaHEPEVTParticle::setMothers: HEPEVT does not allow more than two mothers!");
  51
  52   if(mothers.size()>0) m_first_mother  = mothers[0]->getBarcode();
  53   if(mothers.size()>1) m_second_mother = mothers[1]->getBarcode();
  54 }
  55
  56 void TauolaHEPEVTParticle::setDaughters(vector<TauolaParticle*> daughters){
  57
  58   // This particle must be inside some event record to be able to add daughters
  59   if(m_event==NULL) Log::Fatal("TauolaHEPEVTParticle::setDaughters: particle not inside event record.");
  60
  61   int beg = 65535, end = -1;
  62
  63   for(unsigned int i=0;i<daughters.size();i++)
  64   {
  65     int bc = daughters[i]->getBarcode();
  66     if(bc<0) Log::Fatal("TauolaHEPEVTParticle::setDaughters: all daughters has to be in event record first");
  67     if(bc<beg) beg = bc;
  68     if(bc>end) end = bc;
  69   }
  70   if(end == -1) beg = -1;
  71
  72   m_daughter_start = beg;
  73   m_daughter_end   = end;
  74 }
  75
  76 std::vector<TauolaParticle*> TauolaHEPEVTParticle::getMothers(){
  77
  78   std::vector<TauolaParticle*> mothers;
  79
  80   TauolaParticle *p1 = NULL;
  81   TauolaParticle *p2 = NULL;
  82
  83   if(m_first_mother>=0)  p1 = m_event->getParticle(m_first_mother);
  84   if(m_second_mother>=0) p2 = m_event->getParticle(m_second_mother);
  85
  86   if(p1) mothers.push_back(p1);
  87   if(p2) mothers.push_back(p2);
  88
  89   return mothers;
  90 }
  91
  92 // WARNING: this method also corrects daughter indices
  93 //          if such were not defined
  94 std::vector<TauolaParticle*> TauolaHEPEVTParticle::getDaughters(){
  95
  96   std::vector<TauolaParticle*> daughters;
  97
  98   if(!m_event) return daughters;
  99
 100   // Check if m_daughter_start and m_daughter_end are set
 101   // If not - try to get list of daughters from event
 102   if(m_daughter_end<0)
 103   {
 104     int min_d=65535, max_d=-1;
 105     for(int i=0;i<m_event->getParticleCount();i++)
 106     {
 107       if(m_event->getParticle(i)->isDaughterOf(this))
 108       {
 109         if(i<min_d) min_d = i;
 110         if(i>max_d) max_d = i;
 111       }
 112     }
 113     if(max_d>=0)
 114     {
 115       m_daughter_start = min_d;
 116       m_daughter_end   = max_d;
 117       m_status         = 2;
 118     }
 119   }
 120
 121   // If m_daughter_end is still not set - there are no daughters
 122   // Otherwsie - get daughters
 123   if(m_daughter_end>=0)
 124   {
 125     for(int i=m_daughter_start;i<=m_daughter_end;i++)
 126     {
 127       TauolaParticle *p = m_event->getParticle(i);
 128       if(p==NULL)
 129       {
 130         Log::Warning()<<"TauolaHEPEVTParticle::getDaughters(): No particle with index "<<i<<endl;
 131         return daughters;
 132       }
 133
 134       daughters.push_back(p);
 135     }
 136   }
 137
 138   return daughters;
 139 }
 140
 141 void TauolaHEPEVTParticle::checkMomentumConservation(){
 142
 143   if(!m_event)           return;
 144   if(m_daughter_end < 0) return;
 145
 146   TauolaHEPEVTParticle *buf = m_event->getParticle(m_daughter_start);
 147
 148   int first_mother_idx  = buf->getFirstMotherIndex();
 149   int second_mother_idx = buf->getSecondMotherIndex();
 150
 151   double px =0.0, py =0.0, pz =0.0, e =0.0;
 152   double px2=0.0, py2=0.0, pz2=0.0, e2=0.0;
 153
 154   for(int i=m_daughter_start;i<=m_daughter_end;i++)
 155   {
 156     buf = m_event->getParticle(i);
 157     px += buf->getPx();
 158     py += buf->getPy();
 159     pz += buf->getPz();
 160     e  += buf->getE ();
 161   }
 162
 163   if(first_mother_idx>=0)
 164   {
 165     buf = m_event->getParticle(first_mother_idx);
 166     px2 += buf->getPx();
 167     py2 += buf->getPy();
 168     pz2 += buf->getPz();
 169     e2  += buf->getE();
 170   }
 171
 172   if(second_mother_idx>=0)
 173   {
 174     buf = m_event->getParticle(second_mother_idx);
 175     px2 += buf->getPx();
 176     py2 += buf->getPy();
 177     pz2 += buf->getPz();
 178     e2  += buf->getE();
 179   }
 180   // 3-momentum  // test HepMC style
 181   double dp = sqrt( (px-px2)*(px-px2) + (py-py2)*(py-py2) + (pz-pz2)*(pz-pz2) );
 182   // virtuality test as well.
 183   double m1 = sqrt( fabs( e*e   - px*px   - py*py   - pz*pz   ) );
 184   double m2 = sqrt( fabs( e2*e2 - px2*px2 - py2*py2 - pz2*pz2 ) );
 185
 186   if( fabs(m1-m2) > 0.0001 || dp > 0.0001*(e+e2))
 187   {
 188     Log::RedirectOutput( Log::Warning()<<"Momentum not conserved in vertex: " );
 189     if(first_mother_idx >=0) m_event->getParticle(first_mother_idx) ->print();
 190     if(second_mother_idx>=0) m_event->getParticle(second_mother_idx)->print();
 191     for(int i=m_daughter_start;i<=m_daughter_end;i++) m_event->getParticle(i)->print();
 192     Log::RevertOutput();
 193   }
 194 }
 195
 196 TauolaHEPEVTParticle * TauolaHEPEVTParticle::createNewParticle(
 197                         int pdg_id, int status, double mass,
 198                         double px, double py, double pz, double e){
 199
 200   // New particles created using this method are added to cache
 201   // They will be deleted when this particle will be deleted
 202
 203   cache.push_back(new TauolaHEPEVTParticle(pdg_id,status,px,py,pz,e,mass,-1,-1,-1,-1));
 204   return cache.back();
 205 }
 206
 207 bool TauolaHEPEVTParticle::isDaughterOf(TauolaHEPEVTParticle *p)
 208 {
 209   int bc = p->getBarcode();
 210   if(bc==m_first_mother || bc==m_second_mother) return true;
 211
 212   return false;
 213 }
 214
 215 bool TauolaHEPEVTParticle::isMotherOf  (TauolaHEPEVTParticle *p)
 216 {
 217   int bc = p->getBarcode();
 218   if(bc>=m_daughter_start && bc<=m_daughter_end) return true;
 219
 220   return false;
 221 }
 222
 223 void TauolaHEPEVTParticle::print(){
 224   char buf[256];
 225   sprintf(buf,"P: (%2i) %6i %2i | %11.4e %11.4e %11.4e %11.4e | %11.4e | M: %2i %2i | D: %2i %2i\n",
 226           m_barcode, m_pdgid, m_status, m_px, m_py, m_pz, m_e, m_generated_mass,
 227           m_first_mother, m_second_mother,   m_daughter_start, m_daughter_end);
 228
 229   cout<<buf;
 230 }
 231
 232 /******** Getter and Setter methods: ***********************/
 233
 234 void TauolaHEPEVTParticle::setPdgID(int pdg_id){
 235   m_pdgid = pdg_id;
 236 }
 237
 238 void TauolaHEPEVTParticle::setStatus(int status){
 239   m_status = status;
 240 }
 241
 242 void TauolaHEPEVTParticle::setMass(double mass){
 243   m_generated_mass = mass;
 244 }
 245
 246 int TauolaHEPEVTParticle::getPdgID(){
 247   return m_pdgid;
 248 }
 249
 250 int TauolaHEPEVTParticle::getStatus(){
 251   return m_status;
 252 }
 253
 254 double TauolaHEPEVTParticle::getMass(){
 255   return m_generated_mass;
 256 }
 257
 258 inline double TauolaHEPEVTParticle::getPx(){
 259   return m_px;
 260 }
 261
 262 inline double TauolaHEPEVTParticle::getPy(){
 263   return m_py;
 264 }
 265
 266 double TauolaHEPEVTParticle::getPz(){
 267   return m_pz;
 268 }
 269
 270 double TauolaHEPEVTParticle::getE(){
 271   return m_e;
 272 }
 273
 274 void TauolaHEPEVTParticle::setPx(double px){
 275   m_px = px;
 276 }
 277
 278 void TauolaHEPEVTParticle::setPy(double py){
 279   m_py = py;
 280 }
 281
 282
 283 void TauolaHEPEVTParticle::setPz(double pz){
 284   m_pz = pz;
 285 }
 286
 287 void TauolaHEPEVTParticle::setE(double e){
 288   m_e = e;
 289 }
 290
 291 int TauolaHEPEVTParticle::getBarcode(){
 292   return m_barcode;
 293 }
 294
 295 void TauolaHEPEVTParticle::setBarcode(int barcode){
 296   m_barcode = barcode;
 297 }
 298
 299 void TauolaHEPEVTParticle::setEvent(TauolaHEPEVTEvent *event){
 300   m_event = event;
 301 }
 302
 303 int TauolaHEPEVTParticle::getFirstMotherIndex(){
 304   return m_first_mother;
 305 }
 306
 307 int TauolaHEPEVTParticle::getSecondMotherIndex(){
 308   return m_second_mother;
 309 }
 310
 311 int TauolaHEPEVTParticle::getDaughterRangeStart(){
 312   return m_daughter_start;
 313 }
 314
 315 int TauolaHEPEVTParticle::getDaughterRangeEnd(){
 316   return m_daughter_end;
 317 }
 318
 319 } // namespace Tauolapp