TEvtGen/Tauola/TauolaParticle.cxx

   1 #include "TauolaParticle.h"
   2 #include "Log.h"
   3 #include <stdexcept>
   4
   5 namespace Tauolapp
   6 {
   7
   8 double TauolaParticle::getPolarimetricX(){
   9   return m_pol_x;
  10 }
  11
  12 double TauolaParticle::getPolarimetricY(){
  13   return m_pol_y;
  14 }
  15
  16 double TauolaParticle::getPolarimetricZ(){
  17   return m_pol_z;
  18 }
  19
  20
  21 TauolaParticle * TauolaParticle::clone(){
  22
  23   return createNewParticle(getPdgID(),getStatus(),getMass(),
  24                            getPx(),getPy(),getPz(),getE());
  25
  26 }
  27
  28 // NOTE: Not executed by release examples
  29 double TauolaParticle::getAngle(TauolaParticle * other_particle){
  30
  31   //use the dot product
  32   double x1 = getPx();
  33   double y1 = getPy();
  34   double z1 = getPz();
  35   double x2 = other_particle->getPx();
  36   double y2 = other_particle->getPy();
  37   double z2 = other_particle->getPz();
  38
  39   return acos( (x1*x2+y1*y2+z1*z2) / sqrt((x1*x1+y1*y1+z1*z1)*(x2*x2+y2*y2+z2*z2)) );
  40
  41 }
  42
  43 // NOTE: Not executed by release examples
  44 void TauolaParticle::add(TauolaParticle * other_particle){
  45
  46   setPx(getPx() + other_particle->getPx());
  47   setPy(getPy() + other_particle->getPy());
  48   setPz(getPz() + other_particle->getPz());
  49   setE(getE() + other_particle->getE());
  50   setMass( sqrt( getE()*getE()-getPx()*getPx()-getPy()*getPy()-getPz()*getPz() ));
  51 }
  52
  53 void TauolaParticle::subtract(TauolaParticle * other_particle){
  54
  55   setPx(getPx() - other_particle->getPx());
  56   setPy(getPy() - other_particle->getPy());
  57   setPz(getPz() - other_particle->getPz());
  58   setE(getE() - other_particle->getE());
  59   setMass( sqrt( getE()*getE()-getPx()*getPx()-getPy()*getPy()-getPz()*getPz() ));
  60 }
  61
  62 int TauolaParticle::getSign(){
  63   if(getPdgID()== Tauola::getDecayingParticle())
  64     return SAME_SIGN;
  65   else if(getPdgID()== -1 * Tauola::getDecayingParticle())
  66     return OPPOSITE_SIGN;
  67   else
  68     return NA_SIGN;
  69 }
  70
  71 bool TauolaParticle::hasDaughters(){
  72   if(getDaughters().size()==0)
  73     return 0;
  74   else
  75     return 1;
  76 }
  77
  78 TauolaParticle * TauolaParticle::findLastSelf(){
  79   vector<TauolaParticle*> daughters = getDaughters();
  80   vector<TauolaParticle*>::iterator pcl_itr = daughters.begin();
  81
  82   //get all daughters and look for stable with same pgd id
  83   for(;pcl_itr != daughters.end();pcl_itr++){
  84     if((*pcl_itr)->getPdgID()==this->getPdgID())
  85       return (*pcl_itr)->findLastSelf();
  86   }
  87
  88   return this;
  89 }
  90
  91 std::vector<TauolaParticle*> TauolaParticle::findProductionMothers(){
  92   vector<TauolaParticle*> mothers = getMothers();
  93   vector<TauolaParticle*>::iterator pcl_itr = mothers.begin();
  94
  95     //get all mothers and check none have pdg id of this one
  96     for(;pcl_itr != mothers.end();pcl_itr++){
  97       if((*pcl_itr)->getPdgID()==this->getPdgID())
  98         return (*pcl_itr)->findProductionMothers();
  99     }
 100     return mothers;
 101 }
 102
 103 void TauolaParticle::decay(){
 104
 105   //Do the decay and set the polarimetric vectors
 106   TauolaDecay(getSign(),&m_pol_x, &m_pol_y, &m_pol_z, &m_pol_n);
 107 }
 108
 109 void TauolaParticle::addDecayToEventRecord(){
 110
 111   //Add to decay list used by f_filhep.c
 112   DecayList::addToEnd(this);
 113   TauolaWriteDecayToEventRecord(getSign());
 114
 115   double xmom[4]={0};
 116   double *pp=xmom;
 117
 118   if (Tauola::ion[2])
 119     for(int i=1;1;i++)
 120     {
 121       TauolaParticle *x;
 122       try{ x=DecayList::getParticle(i); }
 123       catch(std::out_of_range d) {break;}
 124       if(x->getPdgID()==221){
 125
 126         // Fix 28.04.2011 The pp vector must have boost for eta undone
 127         TauolaParticle *x_copy = x->clone();
 128         if(getP(TauolaParticle::Z_AXIS)>0)
 129           x_copy->boostAlongZ(-getP(),getE());
 130         else
 131           x_copy->boostAlongZ(getP(),getE());
 132
 133         pp[3]=x_copy->getE();
 134         pp[0]=x_copy->getPx();
 135         pp[1]=x_copy->getPy();
 136         pp[2]=x_copy->getPz();
 137         taueta_(pp,&i);
 138       }
 139     }
 140
 141   if (Tauola::ion[1])
 142     for(int i=1;1;i++)
 143     {
 144       TauolaParticle *x;
 145       try{ x=DecayList::getParticle(i); }
 146       catch(std::out_of_range d) {break;}
 147       if(x->getPdgID()==310){
 148
 149         // Fix 28.04.2011 The pp vector must have boost for k0 undone
 150         TauolaParticle *x_copy = x->clone();
 151         if(getP(TauolaParticle::Z_AXIS)>0)
 152           x_copy->boostAlongZ(-getP(),getE());
 153         else
 154           x_copy->boostAlongZ(getP(),getE());
 155
 156         pp[3]=x_copy->getE();
 157         pp[0]=x_copy->getPx();
 158         pp[1]=x_copy->getPy();
 159         pp[2]=x_copy->getPz();
 160         tauk0s_(pp,&i);
 161       }
 162     }
 163
 164   if (Tauola::ion[0])
 165     for(int i=1;1;i++)
 166     {
 167       TauolaParticle *x;
 168       try{ x=DecayList::getParticle(i); }
 169       catch(std::out_of_range d) {break;}
 170       if(x->getPdgID()==111){
 171
 172         // Fix 28.04.2011 The pp vector must have boost for pi0 undone
 173         TauolaParticle *x_copy = x->clone();
 174         if(getP(TauolaParticle::Z_AXIS)>0)
 175           x_copy->boostAlongZ(-getP(),getE());
 176         else
 177           x_copy->boostAlongZ(getP(),getE());
 178
 179         pp[3]=x_copy->getE();
 180         pp[0]=x_copy->getPx();
 181         pp[1]=x_copy->getPy();
 182         pp[2]=x_copy->getPz();
 183         taupi0_(pp,&i);
 184       }
 185     }
 186   DecayList::clear();
 187
 188   if(!hasDaughters())
 189     Log::Fatal("TAUOLA failed. No decay was created",5);
 190   //  checkMomentumConservation();
 191   //  decayEndgame(); // vertex shift was wrongly calculated,
 192   //                     used 4-momenta should be in lab frame,
 193   //                     thanks to Sho Iwamoto for debug info.
 194
 195 }
 196
 197
 198 void TauolaParticle::boostDaughtersFromRestFrame(TauolaParticle * tau_momentum){
 199
 200   if(!hasDaughters()) //if there are no daughters
 201     return;
 202
 203   vector<TauolaParticle*> daughters = getDaughters();
 204   vector<TauolaParticle*>::iterator pcl_itr = daughters.begin();
 205
 206   //get all daughters then rotate and boost them.
 207   for(;pcl_itr != daughters.end();pcl_itr++){
 208
 209     (*pcl_itr)->boostFromRestFrame(tau_momentum);
 210     (*pcl_itr)->boostDaughtersFromRestFrame(tau_momentum);
 211   }
 212   //checkMomentumConservation();
 213 }
 214
 215 void TauolaParticle::boostDaughtersToRestFrame(TauolaParticle * tau_momentum){
 216
 217   if(!hasDaughters()) //if there are no daughters
 218     return;
 219   // NOTE: Not executed by release examples
 220   //       because !hasDaughters() is always true
 221   vector<TauolaParticle*> daughters = getDaughters();
 222   vector<TauolaParticle*>::iterator pcl_itr = daughters.begin();
 223
 224   //get all daughters then rotate and boost them.
 225   for(;pcl_itr != daughters.end();pcl_itr++){
 226
 227     (*pcl_itr)->boostToRestFrame(tau_momentum);
 228     (*pcl_itr)->boostDaughtersToRestFrame(tau_momentum);
 229   }
 230   //checkMomentumConservation();
 231 }
 232
 233
 234 void TauolaParticle::boostToRestFrame(TauolaParticle * tau_momentum){
 235
 236   double theta = tau_momentum->getRotationAngle(Y_AXIS);
 237   tau_momentum->rotate(Y_AXIS,theta);
 238   double phi = tau_momentum->getRotationAngle(X_AXIS);
 239   tau_momentum->rotate(Y_AXIS,-theta);
 240
 241   //Now rotate coordinates to get boost in Z direction.
 242   rotate(Y_AXIS,theta);
 243   rotate(X_AXIS,phi);
 244   boostAlongZ(-1*tau_momentum->getP(),tau_momentum->getE());
 245   rotate(X_AXIS,-phi);
 246   rotate(Y_AXIS,-theta);
 247
 248 }
 249
 250 void TauolaParticle::boostFromRestFrame(TauolaParticle * tau_momentum){
 251   //get the rotation angles
 252   //and boost z
 253
 254   double theta = tau_momentum->getRotationAngle(Y_AXIS);
 255   tau_momentum->rotate(Y_AXIS,theta);
 256   double phi = tau_momentum->getRotationAngle(X_AXIS);
 257   tau_momentum->rotate(Y_AXIS,-theta);
 258
 259   //Now rotate coordinates to get boost in Z direction.
 260   rotate(Y_AXIS,theta);
 261   rotate(X_AXIS,phi);
 262   boostAlongZ(tau_momentum->getP(),tau_momentum->getE());
 263   rotate(X_AXIS,-phi);
 264   rotate(Y_AXIS,-theta);
 265 }
 266
 267 /** Get the angle needed to rotate the 4 momentum vector so that
 268     the x (y) component disapears. (and the Z component is > 0) */
 269 double TauolaParticle::getRotationAngle(int axis, int second_axis){
 270
 271   /**if(getP(axis)==0){
 272     if(getPz()>0)
 273       return 0; //no rotaion required
 274     else
 275       return M_PI;
 276       }**/
 277   if(getP(second_axis)==0){
 278     if(getP(axis)>0)
 279       return -M_PI/2.0;
 280     else
 281       return M_PI/2.0;
 282   }
 283   if(getP(second_axis)>0)
 284     return -atan(getP(axis)/getP(second_axis));
 285   else
 286     return M_PI-atan(getP(axis)/getP(second_axis));
 287
 288 }
 289
 290 /** Boost this vector along the Z direction.
 291     Assume no momentum components in the X or Y directions. */
 292 void TauolaParticle::boostAlongZ(double boost_pz, double boost_e){
 293
 294   // Boost along the Z axis
 295   double m=sqrt(boost_e*boost_e-boost_pz*boost_pz);
 296
 297   double p=getPz();
 298   double e=getE();
 299
 300   setPz((boost_e*p + boost_pz*e)/m);
 301   setE((boost_pz*p + boost_e*e )/m);
 302 }
 303
 304 /** Rotation around an axis X or Y */
 305 void TauolaParticle::rotate(int axis,double theta, int second_axis){
 306
 307   double temp_px=getP(axis);
 308   double temp_pz=getP(second_axis);
 309   setP(axis,cos(theta)*temp_px + sin(theta)*temp_pz);
 310   setP(second_axis,-sin(theta)*temp_px + cos(theta)*temp_pz);
 311 }
 312
 313 void TauolaParticle::rotateDaughters(int axis,double theta, int second_axis){
 314   if(!hasDaughters()) //if there are no daughters
 315     return;
 316
 317   vector<TauolaParticle*> daughters = getDaughters();
 318   vector<TauolaParticle*>::iterator pcl_itr = daughters.begin();
 319
 320   //get all daughters then rotate and boost them.
 321   for(;pcl_itr != daughters.end();pcl_itr++){
 322
 323     (*pcl_itr)->rotate(axis,theta,second_axis);
 324     (*pcl_itr)->rotateDaughters(axis,theta,second_axis);
 325   }
 326   //checkMomentumConservation();
 327 }
 328
 329 double TauolaParticle::getMass(){
 330   double e_sq=getE()*getE();
 331   double p_sq=getP()*getP();
 332
 333   if(e_sq>p_sq)
 334     return sqrt(e_sq-p_sq);
 335   else
 336     return -1*sqrt(p_sq-e_sq); //if it's negative
 337 }
 338
 339 double TauolaParticle::getP(){
 340   return sqrt(getPx()*getPx()+getPy()*getPy()+getPz()*getPz());
 341 }
 342
 343 double TauolaParticle::getP(int axis){
 344   if(axis==X_AXIS)
 345     return getPx();
 346
 347   if(axis==Y_AXIS)
 348     return getPy();
 349
 350   if(axis==Z_AXIS)
 351     return getPz();
 352
 353   return 0;
 354 }
 355
 356 void TauolaParticle::setP(int axis, double p_component){
 357   if(axis==X_AXIS)
 358     setPx(p_component);
 359   if(axis==Y_AXIS)
 360     setPy(p_component);
 361   if(axis==Z_AXIS)
 362     setPz(p_component);
 363 }
 364
 365 } // namespace Tauolapp