TEvtGen/EvtGen/EvtGenBase/EvtParticle.cpp

   1 //--------------------------------------------------------------------------
   2 //
   3 // Environment:
   4 //      This software is part of the EvtGen package developed jointly
   5 //      for the BaBar and CLEO collaborations.  If you use all or part
   6 //      of it, please give an appropriate acknowledgement.
   7 //
   8 // Copyright Information: See EvtGen/COPYRIGHT
   9 //      Copyright (C) 1998      Caltech, UCSB
  10 //
  11 // Module: EvtParticle.cc
  12 //
  13 // Description: Class to describe all particles
  14 //
  15 // Modification history:
  16 //
  17 //    DJL/RYD     September 25, 1996         Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include <iostream>
  23 #include <math.h>
  24 #include <stdio.h>
  25 #include <stdlib.h>
  26 #include <sys/stat.h>
  27 #include "EvtGenBase/EvtParticle.hh"
  28 #include "EvtGenBase/EvtId.hh"
  29 #include "EvtGenBase/EvtRandom.hh"
  30 #include "EvtGenBase/EvtRadCorr.hh"
  31 #include "EvtGenBase/EvtPDL.hh"
  32 #include "EvtGenBase/EvtDecayTable.hh"
  33 #include "EvtGenBase/EvtDiracParticle.hh"
  34 #include "EvtGenBase/EvtScalarParticle.hh"
  35 #include "EvtGenBase/EvtVectorParticle.hh"
  36 #include "EvtGenBase/EvtTensorParticle.hh"
  37 #include "EvtGenBase/EvtPhotonParticle.hh"
  38 #include "EvtGenBase/EvtNeutrinoParticle.hh"
  39 #include "EvtGenBase/EvtRaritaSchwingerParticle.hh"
  40 #include "EvtGenBase/EvtStringParticle.hh"
  41 #include "EvtGenBase/EvtStdHep.hh"
  42 #include "EvtGenBase/EvtSecondary.hh"
  43 #include "EvtGenBase/EvtReport.hh"
  44 #include "EvtGenBase/EvtGenKine.hh"
  45 #include "EvtGenBase/EvtCPUtil.hh"
  46 #include "EvtGenBase/EvtParticleFactory.hh"
  47 #include "EvtGenBase/EvtIdSet.hh"
  48 #include "EvtGenBase/EvtStatus.hh"
  49
  50 using std::endl;
  51
  52
  53
  54 EvtParticle::~EvtParticle() {
  55   delete _decayProb;
  56 }
  57
  58 EvtParticle::EvtParticle() {
  59    _ndaug=0;
  60    _parent=0;
  61    _channel=-10;
  62    _t=0.0;
  63    _genlifetime=1;
  64    _first=1;
  65    _isInit=false;
  66    _validP4=false;
  67    _isDecayed=false;
  68    _decayProb=0;
  69    //   _mix=false;
  70 }
  71
  72 void EvtParticle::setFirstOrNot() {
  73   _first=0;
  74 }
  75 void EvtParticle::resetFirstOrNot() {
  76   _first=1;
  77 }
  78
  79 void EvtParticle::setChannel( int i ) {
  80   _channel=i;
  81 }
  82
  83 EvtParticle *EvtParticle::getDaug(int i) { return _daug[i]; }
  84
  85 EvtParticle *EvtParticle::getParent() const { return _parent;}
  86
  87 void EvtParticle::setLifetime(double tau){
  88   _t=tau;
  89 }
  90
  91 void EvtParticle::setLifetime(){
  92   if (_genlifetime){
  93     _t=-log(EvtRandom::Flat())*EvtPDL::getctau(getId());
  94   }
  95 }
  96
  97 double EvtParticle::getLifetime(){
  98
  99   return _t;
 100 }
 101
 102 void EvtParticle::addDaug(EvtParticle *node) {
 103   node->_daug[node->_ndaug++]=this;
 104   _ndaug=0;
 105   _parent=node;
 106 }
 107
 108
 109 int EvtParticle::firstornot() const { return _first;}
 110
 111 EvtId EvtParticle::getId() const { return _id;}
 112
 113 int EvtParticle::getPDGId() const {return EvtPDL::getStdHep(_id);}
 114
 115 EvtSpinType::spintype EvtParticle::getSpinType() const
 116       { return EvtPDL::getSpinType(_id);}
 117
 118 int EvtParticle::getSpinStates() const
 119   { return EvtSpinType::getSpinStates(EvtPDL::getSpinType(_id));}
 120
 121 const EvtVector4R& EvtParticle::getP4() const { return _p;}
 122
 123 int EvtParticle::getChannel() const { return _channel;}
 124
 125 size_t EvtParticle::getNDaug() const { return _ndaug;}
 126
 127 double EvtParticle::mass() const {
 128
 129      return _p.mass();
 130 }
 131
 132
 133 void EvtParticle::setDiagonalSpinDensity(){
 134
 135   _rhoForward.setDiag(getSpinStates());
 136 }
 137
 138 void EvtParticle::setVectorSpinDensity(){
 139
 140   if (getSpinStates()!=3) {
 141     report(ERROR,"EvtGen")<<"Error in EvtParticle::setVectorSpinDensity"<<endl;
 142     report(ERROR,"EvtGen")<<"spin_states:"<<getSpinStates()<<endl;
 143     report(ERROR,"EvtGen")<<"particle:"<<EvtPDL::name(_id).c_str()<<endl;
 144     ::abort();
 145   }
 146
 147   EvtSpinDensity rho;
 148
 149   //Set helicity +1 and -1 to 1.
 150   rho.setDiag(getSpinStates());
 151   rho.set(1,1,EvtComplex(0.0,0.0));
 152
 153   setSpinDensityForwardHelicityBasis(rho);
 154
 155 }
 156
 157
 158 void EvtParticle::setSpinDensityForwardHelicityBasis(const EvtSpinDensity& rho){
 159
 160   EvtSpinDensity R=rotateToHelicityBasis();
 161
 162   assert(R.getDim()==rho.getDim());
 163
 164   int n=rho.getDim();
 165
 166   _rhoForward.setDim(n);
 167
 168   int i,j,k,l;
 169
 170   for(i=0;i<n;i++){
 171     for(j=0;j<n;j++){
 172       EvtComplex tmp=0.0;
 173       for(k=0;k<n;k++){
 174         for(l=0;l<n;l++){
 175           tmp+=R.get(l,i)*rho.get(l,k)*conj(R.get(k,j));
 176         }
 177       }
 178       _rhoForward.set(i,j,tmp);
 179     }
 180   }
 181
 182 }
 183
 184 void EvtParticle::setSpinDensityForwardHelicityBasis(const EvtSpinDensity& rho,
 185                                                      double alpha,
 186                                                      double beta,
 187                                                      double gamma){
 188
 189   EvtSpinDensity R=rotateToHelicityBasis(alpha,beta,gamma);
 190
 191   assert(R.getDim()==rho.getDim());
 192
 193   int n=rho.getDim();
 194
 195   _rhoForward.setDim(n);
 196
 197   int i,j,k,l;
 198
 199   for(i=0;i<n;i++){
 200     for(j=0;j<n;j++){
 201       EvtComplex tmp=0.0;
 202       for(k=0;k<n;k++){
 203         for(l=0;l<n;l++){
 204           tmp+=R.get(l,i)*rho.get(l,k)*conj(R.get(k,j));
 205         }
 206       }
 207       _rhoForward.set(i,j,tmp);
 208     }
 209   }
 210
 211 }
 212
 213 void EvtParticle::initDecay(bool useMinMass) {
 214
 215   EvtParticle* p=this;
 216   // carefull - the parent mass might be fixed in stone..
 217   EvtParticle* par=p->getParent();
 218   double parMass=-1.;
 219   if ( par != 0 ) {
 220     if ( par->hasValidP4() ) parMass=par->mass();
 221     for (size_t i=0;i<par->getNDaug();i++) {
 222       EvtParticle *tDaug=par->getDaug(i);
 223       if ( p != tDaug )
 224         parMass-=EvtPDL::getMinMass(tDaug->getId());
 225     }
 226   }
 227
 228   if ( _isInit ) {
 229     //we have already been here - just reroll the masses!
 230     if ( _ndaug>0) {
 231       for(size_t ii=0;ii<_ndaug;ii++){
 232         if ( _ndaug==1 ||  EvtPDL::getWidth(p->getDaug(ii)->getId()) > 0.0000001)
 233           p->getDaug(ii)->initDecay(useMinMass);
 234         else p->getDaug(ii)->setMass(EvtPDL::getMeanMass(p->getDaug(ii)->getId()));
 235       }
 236     }
 237
 238     EvtId *dauId=0;
 239     double *dauMasses=0;
 240     if ( _ndaug > 0) {
 241       dauId=new EvtId[_ndaug];
 242       dauMasses=new double[_ndaug];
 243       for (size_t j=0;j<_ndaug;j++) {
 244         dauId[j]=p->getDaug(j)->getId();
 245         dauMasses[j]=p->getDaug(j)->mass();
 246       }
 247     }
 248     EvtId *parId=0;
 249     EvtId *othDauId=0;
 250     EvtParticle *tempPar=p->getParent();
 251     if (tempPar) {
 252       parId=new EvtId(tempPar->getId());
 253       if ( tempPar->getNDaug()==2 ) {
 254         if ( tempPar->getDaug(0) == this ) othDauId=new EvtId(tempPar->getDaug(1)->getId());
 255         else othDauId=new EvtId(tempPar->getDaug(0)->getId());
 256       }
 257     }
 258     if ( p->getParent() && _validP4==false ) {
 259       if ( !useMinMass ) {
 260         p->setMass(EvtPDL::getRandMass(p->getId(),parId,_ndaug,dauId,othDauId,parMass,dauMasses));
 261       }
 262       else p->setMass(EvtPDL::getMinMass(p->getId()));
 263     }
 264     if ( parId) delete parId;
 265     if ( othDauId) delete othDauId;
 266     if ( dauId) delete [] dauId;
 267     if ( dauMasses) delete [] dauMasses;
 268     return;
 269   }
 270
 271
 272   //Will include effects of mixing here
 273   //added by Lange Jan4,2000
 274   static EvtId BS0=EvtPDL::getId("B_s0");
 275   static EvtId BSB=EvtPDL::getId("anti-B_s0");
 276   static EvtId BD0=EvtPDL::getId("B0");
 277   static EvtId BDB=EvtPDL::getId("anti-B0");
 278   static EvtId D0=EvtPDL::getId("D0");
 279   static EvtId D0B=EvtPDL::getId("anti-D0");
 280   static EvtId U4S=EvtPDL::getId("Upsilon(4S)");
 281   static EvtIdSet borUps(BS0,BSB,BD0,BDB,U4S);
 282
 283   //only makes sense if there is no parent particle which is a B or an Upsilon
 284   bool hasBorUps=false;
 285   if ( getParent() && borUps.contains(getParent()->getId()) ) hasBorUps=true;
 286   //    if ( (getNDaug()==0)&&(getParent()==0) && (getId()==BS0||getId()==BSB||getId()==BD0||getId()==BDB)){
 287   EvtId thisId=getId();
 288   // remove D0 mixing for now.
 289   //  if ( (getNDaug()==0 && !hasBorUps) && (thisId==BS0||thisId==BSB||thisId==BD0||thisId==BDB||thisId==D0||thisId==D0B)){
 290   if ( (getNDaug()==0 && !hasBorUps) && (thisId==BS0||thisId==BSB||thisId==BD0||thisId==BDB)){
 291     double t;
 292     int mix;
 293     EvtCPUtil::getInstance()->incoherentMix(getId(), t, mix);
 294     setLifetime(t);
 295
 296     if (mix) {
 297
 298       EvtScalarParticle* scalar_part;
 299
 300       scalar_part=new EvtScalarParticle;
 301       if (getId()==BS0) {
 302         EvtVector4R p_init(EvtPDL::getMass(BSB),0.0,0.0,0.0);
 303         scalar_part->init(EvtPDL::chargeConj(getId()),p_init);
 304       }
 305       else if (getId()==BSB) {
 306         EvtVector4R p_init(EvtPDL::getMass(BS0),0.0,0.0,0.0);
 307         scalar_part->init(EvtPDL::chargeConj(getId()),p_init);
 308       }
 309       else if (getId()==BD0) {
 310         EvtVector4R p_init(EvtPDL::getMass(BDB),0.0,0.0,0.0);
 311         scalar_part->init(EvtPDL::chargeConj(getId()),p_init);
 312       }
 313       else if (getId()==BDB) {
 314         EvtVector4R p_init(EvtPDL::getMass(BD0),0.0,0.0,0.0);
 315         scalar_part->init(EvtPDL::chargeConj(getId()),p_init);
 316       }
 317       else if (getId()==D0) {
 318         EvtVector4R p_init(EvtPDL::getMass(D0B),0.0,0.0,0.0);
 319         scalar_part->init(EvtPDL::chargeConj(getId()),p_init);
 320       }
 321       else if (getId()==D0B) {
 322         EvtVector4R p_init(EvtPDL::getMass(D0),0.0,0.0,0.0);
 323         scalar_part->init(EvtPDL::chargeConj(getId()),p_init);
 324       }
 325
 326       scalar_part->setLifetime(0);
 327       scalar_part->setDiagonalSpinDensity();
 328
 329       insertDaugPtr(0,scalar_part);
 330
 331       _ndaug=1;
 332       _isInit=true;
 333       p=scalar_part;
 334       p->initDecay(useMinMass);
 335       return;
 336
 337
 338     }
 339   }
 340
 341   EvtDecayBase *decayer;
 342   decayer = EvtDecayTable::getInstance()->getDecayFunc(p);
 343
 344   if ( decayer ) {
 345     p->makeDaughters(decayer->nRealDaughters(),decayer->getDaugs());
 346     //then loop over the daughters and init their decay
 347     for(size_t i=0;i<p->getNDaug();i++){
 348       //      std::cout << EvtPDL::name(p->getDaug(i)->getId()) << " " << i << " " << p->getDaug(i)->getSpinType() << " " << EvtPDL::name(p->getId()) << std::endl;
 349       if ( EvtPDL::getWidth(p->getDaug(i)->getId()) > 0.0000001)
 350         p->getDaug(i)->initDecay(useMinMass);
 351       else p->getDaug(i)->setMass(EvtPDL::getMeanMass(p->getDaug(i)->getId()));
 352     }
 353   }
 354
 355   int j;
 356   EvtId *dauId=0;
 357   double *dauMasses=0;
 358   int nDaugT=p->getNDaug();
 359   if ( nDaugT > 0) {
 360     dauId=new EvtId[nDaugT];
 361     dauMasses=new double[nDaugT];
 362     for (j=0;j<nDaugT;j++) {
 363       dauId[j]=p->getDaug(j)->getId();
 364       dauMasses[j]=p->getDaug(j)->mass();
 365     }
 366   }
 367
 368   EvtId *parId=0;
 369   EvtId *othDauId=0;
 370   EvtParticle *tempPar=p->getParent();
 371   if (tempPar) {
 372     parId=new EvtId(tempPar->getId());
 373     if ( tempPar->getNDaug()==2 ) {
 374       if ( tempPar->getDaug(0) == this ) othDauId=new EvtId(tempPar->getDaug(1)->getId());
 375       else othDauId=new EvtId(tempPar->getDaug(0)->getId());
 376     }
 377   }
 378   if ( p->getParent() && p->hasValidP4()==false ) {
 379     if ( !useMinMass ) {
 380       p->setMass(EvtPDL::getRandMass(p->getId(),parId,p->getNDaug(),dauId,othDauId,parMass,dauMasses));
 381     }
 382     else {
 383       p->setMass(EvtPDL::getMinMass(p->getId()));
 384     }
 385   }
 386   if ( parId) delete parId;
 387   if ( othDauId) delete othDauId;
 388   if ( dauId) delete [] dauId;
 389   if ( dauMasses) delete [] dauMasses;
 390   _isInit=true;
 391 }
 392
 393
 394 void EvtParticle::decay(){
 395   //P is particle to decay, typically 'this' but sometime
 396   //modified by mixing
 397   EvtParticle* p=this;
 398   //Did it mix?
 399   //if ( p->getMixed() ) {
 400     //should take C(p) - this should only
 401     //happen the first time we call decay for this
 402     //particle
 403   //p->takeCConj();
 404   // p->setUnMixed();
 405   //}
 406
 407   EvtDecayBase *decayer;
 408   decayer = EvtDecayTable::getInstance()->getDecayFunc(p);
 409   //  if ( decayer ) {
 410   //    report(INFO,"EvtGen") << "calling decay for " << EvtPDL::name(p->getId()) << " " << p->mass() << " " << p->getP4() << " " << p->getNDaug() << " " << p << endl;
 411   //    report(INFO,"EvtGen") << "NDaug= " << decayer->getNDaug() << endl;
 412   //    int ti;
 413   //    for ( ti=0; ti<decayer->getNDaug(); ti++)
 414   //      report(INFO,"EvtGen") << "Daug " << ti << " " << EvtPDL::name(decayer->getDaug(ti)) << endl;
 415   //  }
 416   //if (p->_ndaug>0) {
 417   //      report(INFO,"EvtGen") <<"Is decaying particle with daughters!!!!!"<<endl;
 418   //     ::abort();
 419     //return;
 420     //call initdecay first - April 29,2002 - Lange
 421   //}
 422
 423   //if there are already daughters, then this step is already done!
 424   // figure out the masses
 425   bool massTreeOK(true);
 426   if ( _ndaug == 0 ) {
 427     massTreeOK = generateMassTree();
 428   }
 429
 430   if (massTreeOK == false) {
 431     report(INFO,"EvtGen")<<"Could not decay "<<EvtPDL::name(p->getId())
 432                          <<" with mass "<<p->mass()
 433                          <<" to decay channel number "<<_channel<<endl;
 434     _isDecayed = false;
 435     return;
 436   }
 437
 438   static EvtId BS0=EvtPDL::getId("B_s0");
 439   static EvtId BSB=EvtPDL::getId("anti-B_s0");
 440   static EvtId BD0=EvtPDL::getId("B0");
 441   static EvtId BDB=EvtPDL::getId("anti-B0");
 442   // static EvtId D0=EvtPDL::getId("D0");
 443   // static EvtId D0B=EvtPDL::getId("anti-D0");
 444
 445   EvtId thisId=getId();
 446   // remove D0 mixing for now..
 447   //  if ( _ndaug==1 &&  (thisId==BS0||thisId==BSB||thisId==BD0||thisId==BDB||thisId==D0||thisId==D0B) ) {
 448   if ( _ndaug==1 &&  (thisId==BS0||thisId==BSB||thisId==BD0||thisId==BDB) ) {
 449     p=p->getDaug(0);
 450     decayer = EvtDecayTable::getInstance()->getDecayFunc(p);
 451   }
 452   //now we have accepted a set of masses - time
 453   if ( decayer != 0) {
 454     decayer->makeDecay(p);
 455   }
 456   else{
 457     p->_rhoBackward.setDiag(p->getSpinStates());
 458   }
 459
 460   _isDecayed=true;
 461   return;
 462 }
 463
 464 bool EvtParticle::generateMassTree() {
 465
 466   bool isOK(true);
 467
 468   double massProb=1.;
 469   double ranNum=2.;
 470   int counter=0;
 471   EvtParticle *p=this;
 472   while (massProb<ranNum) {
 473     //check it out the first time.
 474     p->initDecay();
 475     massProb=p->compMassProb();
 476     ranNum=EvtRandom::Flat();
 477     counter++;
 478
 479     if ( counter > 10000 ) {
 480       if ( counter == 10001 ) {
 481         report(INFO,"EvtGen") << "Too many iterations to determine the mass tree. Parent mass= "<< p->mass() << " " << massProb <<endl;
 482         p->printTree();
 483         report(INFO,"EvtGen") << "will take next combo with non-zero likelihood\n";
 484       }
 485       if ( massProb>0. ) massProb=2.0;
 486       if ( counter > 20000 ) {
 487         // one last try - take the minimum masses
 488         p->initDecay(true);
 489         p->printTree();
 490         massProb=p->compMassProb();
 491         if ( massProb>0. ) {
 492           massProb=2.0;
 493           report(INFO,"EvtGen") << "Taking the minimum mass of all particles in the chain\n";
 494         }
 495         else {
 496           report(INFO,"EvtGen") << "Sorry, no luck finding a valid set of masses.  This may be a pathological combo\n";
 497           isOK = false;
 498           break;
 499         }
 500       }
 501     }
 502   }
 503
 504   return isOK;
 505
 506 }
 507
 508 double EvtParticle::compMassProb() {
 509
 510   EvtParticle *p=this;
 511   double mass=p->mass();
 512   double parMass=0.;
 513   if ( p->getParent()) {
 514     parMass=p->getParent()->mass();
 515   }
 516
 517   int nDaug=p->getNDaug();
 518   double *dMasses=0;
 519
 520   int i;
 521   if ( nDaug>0 ) {
 522     dMasses=new double[nDaug];
 523     for (i=0; i<nDaug; i++) dMasses[i]=p->getDaug(i)->mass();
 524   }
 525
 526   double temp=1.0;
 527   temp=EvtPDL::getMassProb(p->getId(), mass, parMass, nDaug, dMasses);
 528
 529   //If the particle already has a mass, we dont need to include
 530   //it in the probability calculation
 531   if ( (!p->getParent() || _validP4 ) && temp>0.0 ) temp=1.;
 532
 533   delete [] dMasses;
 534   for (i=0; i<nDaug; i++) {
 535     temp*=p->getDaug(i)->compMassProb();
 536   }
 537   return temp;
 538 }
 539
 540 void EvtParticle::deleteDaughters(bool keepChannel){
 541
 542   for(size_t i=0;i<_ndaug;i++){
 543     _daug[i]->deleteTree();
 544   }
 545
 546   _ndaug=0;
 547   if ( !keepChannel) _channel=-10;
 548   _first=1;
 549   _isInit=false;
 550 }
 551
 552 void EvtParticle::deleteTree(){
 553
 554   this->deleteDaughters();
 555
 556   delete this;
 557
 558 }
 559
 560 EvtVector4C EvtParticle::epsParent(int i) const {
 561   EvtVector4C temp;
 562   printParticle();
 563   report(ERROR,"EvtGen") << "and you have asked for the:"<<i
 564                          <<"th polarization vector."
 565                          <<" I.e. you thought it was a"
 566                          <<" vector particle!" << endl;
 567   ::abort();
 568   return temp;
 569 }
 570
 571 EvtVector4C EvtParticle::eps(int i) const {
 572   EvtVector4C temp;
 573   printParticle();
 574   report(ERROR,"EvtGen") << "and you have asked for the:"<<i
 575                          <<"th polarization vector."
 576                          <<" I.e. you thought it was a"
 577                          <<" vector particle!" << endl;
 578   ::abort();
 579   return temp;
 580 }
 581
 582 EvtVector4C EvtParticle::epsParentPhoton(int i){
 583   EvtVector4C temp;
 584   printParticle();
 585   report(ERROR,"EvtGen") << "and you have asked for the:"<<i
 586                          <<"th polarization vector of photon."
 587                          <<" I.e. you thought it was a"
 588                          <<" photon particle!" << endl;
 589   ::abort();
 590   return temp;
 591 }
 592
 593 EvtVector4C EvtParticle::epsPhoton(int i){
 594   EvtVector4C temp;
 595   printParticle();
 596   report(ERROR,"EvtGen") << "and you have asked for the:"<<i
 597                          <<"th polarization vector of a photon."
 598                          <<" I.e. you thought it was a"
 599                          <<" photon particle!" << endl;
 600   ::abort();
 601   return temp;
 602 }
 603
 604 EvtDiracSpinor EvtParticle::spParent(int i) const {
 605   EvtDiracSpinor tempD;
 606   printParticle();
 607   report(ERROR,"EvtGen") << "and you have asked for the:"<<i
 608                          <<"th dirac spinor."
 609                          <<" I.e. you thought it was a"
 610                          <<" Dirac particle!" << endl;
 611   ::abort();
 612   return tempD;
 613 }
 614
 615 EvtDiracSpinor EvtParticle::sp(int i) const {
 616   EvtDiracSpinor tempD;
 617   printParticle();
 618   report(ERROR,"EvtGen") << "and you have asked for the:"<<i
 619                          <<"th dirac spinor."
 620                          <<" I.e. you thought it was a"
 621                          <<" Dirac particle!" << endl;
 622   ::abort();
 623   return tempD;
 624 }
 625
 626 EvtDiracSpinor EvtParticle::spParentNeutrino() const {
 627   EvtDiracSpinor tempD;
 628   printParticle();
 629   report(ERROR,"EvtGen") << "and you have asked for the "
 630                          <<"dirac spinor."
 631                          <<" I.e. you thought it was a"
 632                          <<" neutrino particle!" << endl;
 633   ::abort();
 634   return tempD;
 635 }
 636
 637 EvtDiracSpinor EvtParticle::spNeutrino() const {
 638   EvtDiracSpinor tempD;
 639   printParticle();
 640   report(ERROR,"EvtGen") << "and you have asked for the "
 641                          <<"dirac spinor."
 642                          <<" I.e. you thought it was a"
 643                          <<" neutrino particle!" << endl;
 644   ::abort();
 645   return tempD;
 646 }
 647
 648 EvtTensor4C EvtParticle::epsTensorParent(int i) const {
 649   EvtTensor4C tempC;
 650   printParticle();
 651   report(ERROR,"EvtGen") << "and you have asked for the:"<<i
 652                          <<"th tensor."
 653                          <<" I.e. you thought it was a"
 654                          <<" Tensor particle!" << endl;
 655   ::abort();
 656   return tempC;
 657 }
 658
 659 EvtTensor4C EvtParticle::epsTensor(int i) const {
 660   EvtTensor4C tempC;
 661   printParticle();
 662   report(ERROR,"EvtGen") << "and you have asked for the:"<<i
 663                          <<"th tensor."
 664                          <<" I.e. you thought it was a"
 665                          <<" Tensor particle!" << endl;
 666   ::abort();
 667   return tempC;
 668 }
 669
 670
 671 EvtRaritaSchwinger EvtParticle::spRSParent(int i) const {
 672   EvtRaritaSchwinger tempD;
 673   printParticle();
 674   report(ERROR,"EvtGen") << "and you have asked for the:"<<i
 675                          <<"th Rarita-Schwinger spinor."
 676                          <<" I.e. you thought it was a"
 677                          <<" RaritaSchwinger particle!" << std::endl;
 678   ::abort();
 679   return tempD;
 680 }
 681
 682 EvtRaritaSchwinger EvtParticle::spRS(int i) const {
 683   EvtRaritaSchwinger tempD;
 684   printParticle();
 685   report(ERROR,"EvtGen") << "and you have asked for the:"<<i
 686                          <<"th Rarita-Schwinger spinor."
 687                          <<" I.e. you thought it was a"
 688                          <<" RaritaSchwinger particle!" << std::endl;
 689   ::abort();
 690   return tempD;
 691 }
 692
 693
 694
 695 EvtVector4R EvtParticle::getP4Lab() const {
 696   EvtVector4R temp,mom;
 697   const EvtParticle *ptemp;
 698
 699   temp=this->getP4();
 700   ptemp=this;
 701
 702   while (ptemp->getParent()!=0) {
 703     ptemp=ptemp->getParent();
 704     mom=ptemp->getP4();
 705     temp=boostTo(temp,mom);
 706   }
 707   return temp;
 708 }
 709
 710 EvtVector4R EvtParticle::getP4LabBeforeFSR() {
 711   EvtVector4R temp,mom;
 712   EvtParticle *ptemp;
 713
 714   temp=this->_pBeforeFSR;
 715   ptemp=this;
 716
 717   while (ptemp->getParent()!=0) {
 718     ptemp=ptemp->getParent();
 719     mom=ptemp->getP4();
 720     temp=boostTo(temp,mom);
 721   }
 722   return temp;
 723 }
 724
 725
 726
 727 EvtVector4R EvtParticle::getP4Restframe() const {
 728
 729   return EvtVector4R(mass(),0.0,0.0,0.0);
 730
 731 }
 732
 733 EvtVector4R EvtParticle::get4Pos() const {
 734
 735   EvtVector4R temp,mom;
 736   EvtParticle *ptemp;
 737
 738   temp.set(0.0,0.0,0.0,0.0);
 739   ptemp=getParent();
 740
 741   if (ptemp==0) return temp;
 742
 743   temp=(ptemp->_t/ptemp->mass())*(ptemp->getP4());
 744
 745   while (ptemp->getParent()!=0) {
 746     ptemp=ptemp->getParent();
 747     mom=ptemp->getP4();
 748     temp=boostTo(temp,mom);
 749     temp=temp+(ptemp->_t/ptemp->mass())*(ptemp->getP4());
 750   }
 751
 752   return temp;
 753 }
 754
 755
 756 EvtParticle * EvtParticle::nextIter(EvtParticle *rootOfTree) {
 757
 758   EvtParticle *bpart;
 759   EvtParticle *current;
 760
 761   current=this;
 762   size_t i;
 763
 764   if (_ndaug!=0) return _daug[0];
 765
 766   do{
 767     bpart=current->_parent;
 768     if (bpart==0) return 0;
 769     i=0;
 770     while (bpart->_daug[i]!=current) {i++;}
 771
 772     if ( bpart==rootOfTree ) {
 773       if ( i+1 == bpart->_ndaug ) return 0;
 774     }
 775
 776     i++;
 777     current=bpart;
 778
 779   }while(i>=bpart->_ndaug);
 780
 781   return bpart->_daug[i];
 782
 783 }
 784
 785
 786 void EvtParticle::makeStdHep(EvtStdHep& stdhep,EvtSecondary& secondary,
 787                              EvtId *list_of_stable){
 788
 789   //first add particle to the stdhep list;
 790   stdhep.createParticle(getP4Lab(),get4Pos(),-1,-1,
 791                         EvtPDL::getStdHep(getId()));
 792
 793   int ii=0;
 794
 795   //lets see if this is a longlived particle and terminate the
 796   //list building!
 797
 798   while (list_of_stable[ii]!=EvtId(-1,-1)) {
 799     if (getId()==list_of_stable[ii]){
 800       secondary.createSecondary(0,this);
 801       return;
 802     }
 803     ii++;
 804   }
 805
 806
 807
 808   for(size_t i=0;i<_ndaug;i++){
 809     stdhep.createParticle(_daug[i]->getP4Lab(),_daug[i]->get4Pos(),0,0,
 810                           EvtPDL::getStdHep(_daug[i]->getId()));
 811   }
 812
 813   for(size_t i=0;i<_ndaug;i++){
 814     _daug[i]->makeStdHepRec(1+i,1+i,stdhep,secondary,list_of_stable);
 815   }
 816   return;
 817
 818 }
 819
 820 void EvtParticle::makeStdHep(EvtStdHep& stdhep){
 821
 822   //first add particle to the stdhep list;
 823   stdhep.createParticle(getP4Lab(),get4Pos(),-1,-1,
 824                         EvtPDL::getStdHep(getId()));
 825
 826   for(size_t i=0;i<_ndaug;i++){
 827     stdhep.createParticle(_daug[i]->getP4Lab(),_daug[i]->get4Pos(),0,0,
 828                           EvtPDL::getStdHep(_daug[i]->getId()));
 829   }
 830
 831   for(size_t i=0;i<_ndaug;i++){
 832     _daug[i]->makeStdHepRec(1+i,1+i,stdhep);
 833   }
 834   return;
 835
 836 }
 837
 838
 839 void EvtParticle::makeStdHepRec(int firstparent,int lastparent,
 840                                 EvtStdHep& stdhep,
 841                                 EvtSecondary& secondary,
 842                                 EvtId *list_of_stable){
 843
 844
 845   int ii=0;
 846
 847   //lets see if this is a longlived particle and terminate the
 848   //list building!
 849
 850   while (list_of_stable[ii]!=EvtId(-1,-1)) {
 851     if (getId()==list_of_stable[ii]){
 852       secondary.createSecondary(firstparent,this);
 853       return;
 854     }
 855     ii++;
 856   }
 857
 858
 859
 860   int parent_num=stdhep.getNPart();
 861   for(size_t i=0;i<_ndaug;i++){
 862     stdhep.createParticle(_daug[i]->getP4Lab(),_daug[i]->get4Pos(),
 863                           firstparent,lastparent,
 864                           EvtPDL::getStdHep(_daug[i]->getId()));
 865   }
 866
 867   for(size_t i=0;i<_ndaug;i++){
 868     _daug[i]->makeStdHepRec(parent_num+i,parent_num+i,stdhep,
 869                            secondary,list_of_stable);
 870   }
 871   return;
 872
 873 }
 874
 875 void EvtParticle::makeStdHepRec(int firstparent,int lastparent,
 876                                 EvtStdHep& stdhep){
 877
 878   int parent_num=stdhep.getNPart();
 879   for(size_t i=0;i<_ndaug;i++){
 880     stdhep.createParticle(_daug[i]->getP4Lab(),_daug[i]->get4Pos(),
 881                           firstparent,lastparent,
 882                           EvtPDL::getStdHep(_daug[i]->getId()));
 883   }
 884
 885   for(size_t i=0;i<_ndaug;i++){
 886     _daug[i]->makeStdHepRec(parent_num+i,parent_num+i,stdhep);
 887   }
 888   return;
 889
 890 }
 891
 892 void EvtParticle::printTreeRec(unsigned int level) const {
 893
 894   size_t newlevel,i;
 895   newlevel = level +1;
 896
 897
 898   if (_ndaug!=0) {
 899     if ( level > 0 ) {
 900       for (i=0;i<(5*level);i++) {
 901         report(INFO,"") <<" ";
 902       }
 903     }
 904     report(INFO,"") << EvtPDL::name(_id).c_str();
 905     report(INFO,"") << " -> ";
 906     for(i=0;i<_ndaug;i++){
 907       report(INFO,"") << EvtPDL::name(_daug[i]->getId()).c_str()<<" ";
 908     }
 909     for(i=0;i<_ndaug;i++){
 910       report(INFO,"") << _daug[i]->mass()<< " " << _daug[i]->getP4() << " " <<_daug[i]->getSpinStates() << "; ";
 911     }
 912     report(INFO,"")<<endl;
 913     for(i=0;i<_ndaug;i++){
 914       _daug[i]->printTreeRec(newlevel);
 915     }
 916   }
 917 }
 918
 919 void EvtParticle::printTree() const {
 920
 921   report(INFO,"EvtGen") << "This is the current decay chain"<<endl;
 922   report(INFO,"") << "This top particle is "<<
 923     EvtPDL::name(_id).c_str()<<" " << this->mass() << " " << this->getP4() << endl;
 924
 925   this->printTreeRec(0);
 926   report(INFO,"EvtGen") << "End of decay chain."<<endl;
 927
 928 }
 929
 930 std::string EvtParticle::treeStrRec(unsigned int level) const {
 931
 932   size_t newlevel,i;
 933   newlevel = level +1;
 934
 935   std::string retval="";
 936
 937   for(i=0;i<_ndaug;i++){
 938     retval+=EvtPDL::name(_daug[i]->getId());
 939     if ( _daug[i]->getNDaug() > 0 ) {
 940       retval+= " (";
 941       retval+= _daug[i]->treeStrRec(newlevel);
 942       retval+= ") ";
 943     }
 944     else{
 945       if ( i+1 !=_ndaug) retval+=" ";
 946     }
 947   }
 948
 949   return retval;
 950 }
 951
 952
 953 std::string EvtParticle::treeStr() const {
 954
 955   std::string retval=EvtPDL::name(_id);
 956   retval+=" -> ";
 957
 958   retval+=treeStrRec(0);
 959
 960   return retval;
 961 }
 962
 963 void EvtParticle::printParticle() const {
 964
 965   switch (EvtPDL::getSpinType(_id)){
 966   case EvtSpinType::SCALAR:
 967     report(INFO,"EvtGen") << "This is a scalar particle:"<<EvtPDL::name(_id).c_str()<<"\n";
 968     break;
 969   case EvtSpinType::VECTOR:
 970     report(INFO,"EvtGen") << "This is a vector particle:"<<EvtPDL::name(_id).c_str()<<"\n";
 971     break;
 972   case EvtSpinType::TENSOR:
 973     report(INFO,"EvtGen") << "This is a tensor particle:"<<EvtPDL::name(_id).c_str()<<"\n";
 974     break;
 975   case EvtSpinType::DIRAC:
 976     report(INFO,"EvtGen") << "This is a dirac particle:"<<EvtPDL::name(_id).c_str()<<"\n";
 977     break;
 978   case EvtSpinType::PHOTON:
 979     report(INFO,"EvtGen") << "This is a photon:"<<EvtPDL::name(_id).c_str()<<"\n";
 980     break;
 981   case EvtSpinType::NEUTRINO:
 982     report(INFO,"EvtGen") << "This is a neutrino:"<<EvtPDL::name(_id).c_str()<<"\n";
 983     break;
 984   case EvtSpinType::STRING:
 985     report(INFO,"EvtGen") << "This is a string:"<<EvtPDL::name(_id).c_str()<<"\n";
 986     break;
 987   default:
 988     report(INFO,"EvtGen") <<"Unknown particle type in EvtParticle::printParticle()"<<endl;
 989     break;
 990   }
 991   report(INFO,"EvtGen") << "Number of daughters:"<<_ndaug<<"\n";
 992
 993
 994 }
 995
 996
 997
 998 void init_vector( EvtParticle **part ){
 999   *part = (EvtParticle *) new EvtVectorParticle;
1000 }
1001
1002
1003 void init_scalar( EvtParticle **part ){
1004   *part = (EvtParticle *) new EvtScalarParticle;
1005 }
1006
1007 void init_tensor( EvtParticle **part ){
1008   *part = (EvtParticle *) new EvtTensorParticle;
1009 }
1010
1011 void init_dirac( EvtParticle **part ){
1012   *part = (EvtParticle *) new EvtDiracParticle;
1013 }
1014
1015 void init_photon( EvtParticle **part ){
1016   *part = (EvtParticle *) new EvtPhotonParticle;
1017 }
1018
1019 void init_neutrino( EvtParticle **part ){
1020   *part = (EvtParticle *) new EvtNeutrinoParticle;
1021 }
1022
1023 void init_string( EvtParticle **part ){
1024   *part = (EvtParticle *) new EvtStringParticle;
1025 }
1026
1027 double EvtParticle::initializePhaseSpace(
1028                    unsigned int numdaughter,EvtId *daughters,
1029                    bool forceDaugMassReset, double poleSize,
1030                    int whichTwo1, int whichTwo2) {
1031
1032   double m_b;
1033   unsigned int i;
1034   //lange
1035   //  this->makeDaughters(numdaughter,daughters);
1036
1037   static EvtVector4R p4[100];
1038   static double mass[100];
1039
1040   m_b = this->mass();
1041
1042   //lange - Jan2,2002 - Need to check to see if the daughters of the parent
1043   // have changed. If so, delete them and start over.
1044   //report(INFO,"EvtGen") << "the parent is\n";
1045   //if ( this->getParent() ) {
1046   //  if ( this->getParent()->getParent() ) this->getParent()->getParent()->printTree();
1047     //    this->getParent()->printTree();
1048   //}
1049   //report(INFO,"EvtGen") << "and this is\n";
1050   //if ( this) this->printTree();
1051   bool resetDaughters=false;
1052
1053   if ( numdaughter != this->getNDaug() && this->getNDaug() > 0 ) resetDaughters=true;
1054   if ( numdaughter == this->getNDaug() )
1055     for (i=0; i<numdaughter;i++) {
1056       if ( this->getDaug(i)->getId() != daughters[i] ) resetDaughters=true;
1057       //report(INFO,"EvtGen") << EvtPDL::name(this->getDaug(i)->getId())
1058       //                    << " " << EvtPDL::name(daughters[i]) << endl;
1059     }
1060
1061   if ( resetDaughters || forceDaugMassReset) {
1062     bool t1=true;
1063     //but keep the decay channel of the parent.
1064     this->deleteDaughters(t1);
1065     this->makeDaughters(numdaughter,daughters);
1066     bool massTreeOK = this->generateMassTree();
1067     if (massTreeOK == false) {return 0.0;}
1068   }
1069
1070   double weight=0.;
1071   for (i=0; i<numdaughter;i++) {
1072     mass[i]=this->getDaug(i)->mass();
1073   }
1074
1075   if ( poleSize<-0.1) {
1076     //special case to enforce 4-momentum conservation in 1->1 decays
1077     if (numdaughter==1) {
1078       this->getDaug(0)->init(daughters[0],EvtVector4R(m_b,0.0,0.0,0.0));
1079     }
1080     else{
1081       EvtGenKine::PhaseSpace( numdaughter, mass, p4, m_b );
1082       for(i=0;i<numdaughter;i++){
1083         this->getDaug(i)->init(daughters[i],p4[i]);
1084       }
1085     }
1086   }
1087   else  {
1088     if ( numdaughter != 3 ) {
1089       report(ERROR,"EvtGen") << "Only can generate pole phase space "
1090                              << "distributions for 3 body final states"
1091                              << endl<<"Will terminate."<<endl;
1092       ::abort();
1093     }
1094     bool ok=false;
1095     if ( (whichTwo1 == 1 && whichTwo2 == 0 ) ||
1096          (whichTwo1 == 0 && whichTwo2 == 1 ) ) {
1097       weight=EvtGenKine::PhaseSpacePole( m_b, mass[0], mass[1], mass[2],
1098                                           poleSize, p4);
1099       this->getDaug(0)->init(daughters[0],p4[0]);
1100       this->getDaug(1)->init(daughters[1],p4[1]);
1101       this->getDaug(2)->init(daughters[2],p4[2]);
1102       ok=true;
1103     }
1104     if ( (whichTwo1 == 1 && whichTwo2 == 2 ) ||
1105          (whichTwo1 == 2 && whichTwo2 == 1 ) ) {
1106       weight=EvtGenKine::PhaseSpacePole( m_b, mass[2], mass[1], mass[0],
1107                                           poleSize, p4);
1108       this->getDaug(0)->init(daughters[0],p4[2]);
1109       this->getDaug(1)->init(daughters[1],p4[1]);
1110       this->getDaug(2)->init(daughters[2],p4[0]);
1111       ok=true;
1112     }
1113     if ( (whichTwo1 == 0 && whichTwo2 == 2 ) ||
1114          (whichTwo1 == 2 && whichTwo2 == 0 ) ) {
1115       weight=EvtGenKine::PhaseSpacePole( m_b, mass[1], mass[0], mass[2],
1116                                           poleSize, p4);
1117       this->getDaug(0)->init(daughters[0],p4[1]);
1118       this->getDaug(1)->init(daughters[1],p4[0]);
1119       this->getDaug(2)->init(daughters[2],p4[2]);
1120       ok=true;
1121     }
1122     if ( !ok) {
1123       report(ERROR,"EvtGen") << "Invalid pair of particle to generate a pole dist "
1124                              << whichTwo1 << " " << whichTwo2
1125                              << endl<<"Will terminate."<<endl;
1126       ::abort();
1127     }
1128   }
1129
1130   return weight;
1131 }
1132
1133 void EvtParticle::makeDaughters(unsigned int ndaugstore, std::vector<EvtId> idVector) {
1134
1135   // Convert the STL vector method to use the array method for now, since the
1136   // array method pervades most of the EvtGen code...
1137
1138   unsigned int nVector = idVector.size();
1139   if (nVector < ndaugstore) {
1140     report(ERROR,"EvtGen") << "Asking to make "<<ndaugstore<<" daughters when there "
1141                            << "are only "<<nVector<<" EvtId values available"<<endl;
1142     return;
1143   }
1144
1145   EvtId *idArray=new EvtId[ndaugstore];
1146   unsigned int i;
1147   for (i = 0; i < ndaugstore; i++) {
1148     idArray[i] = idVector[i];
1149   }
1150
1151   this->makeDaughters(ndaugstore, idArray);
1152
1153   delete[] idArray;
1154 }
1155
1156 void EvtParticle::makeDaughters( unsigned int ndaugstore, EvtId *id){
1157
1158   unsigned int i;
1159   if ( _channel < 0 ) {
1160     setChannel(0);
1161   }
1162   EvtParticle* pdaug;
1163   if (_ndaug!=0 ){
1164     if (_ndaug!=ndaugstore){
1165       report(ERROR,"EvtGen") << "Asking to make a different number of "
1166                              << "daughters than what was previously created."<<endl;
1167       report(ERROR,"EvtGen") << "Original parent:"<<EvtPDL::name(_id)<<endl;
1168       for (size_t i=0;i<_ndaug;i++){
1169           report(ERROR,"EvtGen") << "Original daugther:"<<EvtPDL::name(getDaug(i)->getId())<<endl;
1170       }
1171       for (size_t i=0;i<ndaugstore;i++){
1172           report(ERROR,"EvtGen") << "New Daug:"<<EvtPDL::name(id[i])<<endl;
1173       }
1174       report(ERROR,"EvtGen") << "Will terminate."<<endl;
1175       ::abort();
1176     }
1177   }
1178   else{
1179     for(i=0;i<ndaugstore;i++){
1180       pdaug=EvtParticleFactory::particleFactory(EvtPDL::getSpinType(id[i]));
1181       pdaug->setId(id[i]);
1182       pdaug->addDaug(this);
1183     }
1184
1185   } //else
1186 } //makeDaughters
1187
1188
1189 void EvtParticle::setDecayProb(double prob) {
1190
1191   if ( _decayProb == 0 )  _decayProb=new double;
1192   *_decayProb=prob;
1193 }
1194
1195 std::string EvtParticle::getName() {
1196
1197   std::string theName = _id.getName();
1198   return theName;
1199
1200 }