TEvtGen/EvtGenBase/EvtParticle.cxx

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