00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "EvtGenBase/EvtPatches.hh"
00022 #include <iostream>
00023 #include <math.h>
00024 #include <fstream>
00025 #include <stdio.h>
00026 #include <stdlib.h>
00027 #include <sys/stat.h>
00028 #include <strstream>
00029 #include "EvtGenBase/EvtParticle.hh"
00030 #include "EvtGenBase/EvtId.hh"
00031 #include "EvtGenBase/EvtRandom.hh"
00032 #include "EvtGenBase/EvtRadCorr.hh"
00033 #include "EvtGenBase/EvtPDL.hh"
00034 #include "EvtGenBase/EvtDecayTable.hh"
00035 #include "EvtGenBase/EvtDiracParticle.hh"
00036 #include "EvtGenBase/EvtScalarParticle.hh"
00037 #include "EvtGenBase/EvtVectorParticle.hh"
00038 #include "EvtGenBase/EvtTensorParticle.hh"
00039 #include "EvtGenBase/EvtPhotonParticle.hh"
00040 #include "EvtGenBase/EvtRaritaSchwingerParticle.hh"
00041 #include "EvtGenBase/EvtNeutrinoParticle.hh"
00042 #include "EvtGenBase/EvtStringParticle.hh"
00043 #include "EvtGenBase/EvtStdHep.hh"
00044 #include "EvtGenBase/EvtSecondary.hh"
00045 #include "EvtGenBase/EvtReport.hh"
00046 #include "EvtGenBase/EvtGenKine.hh"
00047 #include "EvtGenBase/EvtCPUtil.hh"
00048 #include "EvtGenBase/EvtParticleFactory.hh"
00049 #include "EvtGenBase/EvtSpinDensity.hh"
00050
00051 using std::endl;
00052 using std::fstream;
00053 using std::strstream;
00054
00055
00056 EvtParticle::~EvtParticle() {
00057 delete _decayProb;
00058 }
00059
00060 EvtParticle::EvtParticle() {
00061 _ndaug=0;
00062 _parent=0;
00063 _channel=-10;
00064 _t=0.0;
00065 _genlifetime=1;
00066 _first=1;
00067 _isInit=false;
00068 _validP4=false;
00069 _isDecayed=false;
00070 _decayProb=0;
00071
00072 }
00073
00074 void EvtParticle::setFirstOrNot() {
00075 _first=0;
00076 }
00077 void EvtParticle::resetFirstOrNot() {
00078 _first=1;
00079 }
00080
00081 void EvtParticle::setChannel( int i ) {
00082 _channel=i;
00083 }
00084
00085 EvtParticle *EvtParticle::getDaug(int i) { return _daug[i]; }
00086
00087 EvtParticle *EvtParticle::getParent() { return _parent;}
00088
00089 void EvtParticle::setLifetime(double tau){
00090 _t=tau;
00091 }
00092
00093 void EvtParticle::setLifetime(){
00094 if (_genlifetime){
00095 _t=-log(EvtRandom::Flat())*EvtPDL::getctau(getId());
00096 }
00097 }
00098
00099 double EvtParticle::getLifetime(){
00100
00101 return _t;
00102 }
00103
00104 void EvtParticle::addDaug(EvtParticle *node) {
00105 node->_daug[node->_ndaug++]=this;
00106 _ndaug=0;
00107 _parent=node;
00108 }
00109
00110
00111 int EvtParticle::firstornot() const { return _first;}
00112
00113 EvtId EvtParticle::getId() const { return _id;}
00114
00115 EvtSpinType::spintype EvtParticle::getSpinType() const
00116 { return EvtPDL::getSpinType(_id);}
00117
00118 int EvtParticle::getSpinStates() const
00119 { return EvtSpinType::getSpinStates(EvtPDL::getSpinType(_id));}
00120
00121 const EvtVector4R& EvtParticle::getP4() const { return _p;}
00122
00123 int EvtParticle::getChannel() const { return _channel;}
00124
00125 int EvtParticle::getNDaug() const { return _ndaug;}
00126
00127 double EvtParticle::mass() const {
00128
00129 return _p.mass();
00130 }
00131
00132
00133 void EvtParticle::setDiagonalSpinDensity(){
00134
00135 _rhoForward.SetDiag(getSpinStates());
00136 }
00137
00138 void EvtParticle::setVectorSpinDensity(){
00139
00140 if (getSpinStates()!=3) {
00141 report(ERROR,"EvtGen")<<"Error in EvtParticle::setVectorSpinDensity"<<endl;
00142 report(ERROR,"EvtGen")<<"spin_states:"<<getSpinStates()<<endl;
00143 report(ERROR,"EvtGen")<<"particle:"<<EvtPDL::name(_id).c_str()<<endl;
00144 ::abort();
00145 }
00146
00147 EvtSpinDensity rho;
00148
00149
00150 rho.SetDiag(getSpinStates());
00151 rho.Set(1,1,EvtComplex(0.0,0.0));
00152 setSpinDensityForwardHelicityBasis(rho);
00153
00154 }
00155
00156
00157 void EvtParticle::setPolarizedSpinDensity(double r00,double r11,double r22){
00158
00159 if (getSpinStates()!=3) {
00160 report(ERROR,"EvtGen")<<"Error in EvtParticle::setVectorSpinDensity"<<endl;
00161 report(ERROR,"EvtGen")<<"spin_states:"<<getSpinStates()<<endl;
00162 report(ERROR,"EvtGen")<<"particle:"<<EvtPDL::name(_id).c_str()<<endl;
00163 ::abort();
00164 }
00165
00166 EvtSpinDensity rho;
00167
00168
00169 rho.SetDiag(getSpinStates());
00170 rho.Set(0,0,EvtComplex(r00,0.));
00171 rho.Set(1,1,EvtComplex(r11,0.));
00172 rho.Set(2,2,EvtComplex(r22,0.));
00173 setSpinDensityForwardHelicityBasis(rho);
00174
00175 }
00176
00177
00178 void EvtParticle::setSpinDensityForwardHelicityBasis(const EvtSpinDensity& rho){
00179
00180 EvtSpinDensity R=rotateToHelicityBasis();
00181
00182 assert(R.GetDim()==rho.GetDim());
00183
00184 int n=rho.GetDim();
00185
00186 _rhoForward.SetDim(n);
00187
00188 int i,j,k,l;
00189
00190 for(i=0;i<n;i++){
00191 for(j=0;j<n;j++){
00192 EvtComplex tmp=0.0;
00193 for(k=0;k<n;k++){
00194 for(l=0;l<n;l++){
00195 tmp+=R.Get(l,i)*rho.Get(l,k)*conj(R.Get(k,j));
00196 }
00197 }
00198 _rhoForward.Set(i,j,tmp);
00199 }
00200 }
00201
00202
00203
00204 }
00205
00206 void EvtParticle::setSpinDensityForwardHelicityBasis(const EvtSpinDensity& rho,
00207 double alpha,
00208 double beta,
00209 double gamma){
00210
00211 EvtSpinDensity R=rotateToHelicityBasis(alpha,beta,gamma);
00212
00213 assert(R.GetDim()==rho.GetDim());
00214
00215 int n=rho.GetDim();
00216
00217 _rhoForward.SetDim(n);
00218
00219 int i,j,k,l;
00220
00221 for(i=0;i<n;i++){
00222 for(j=0;j<n;j++){
00223 EvtComplex tmp=0.0;
00224 for(k=0;k<n;k++){
00225 for(l=0;l<n;l++){
00226 tmp+=R.Get(l,i)*rho.Get(l,k)*conj(R.Get(k,j));
00227 }
00228 }
00229 _rhoForward.Set(i,j,tmp);
00230 }
00231 }
00232
00233
00234
00235 }
00236
00237 void EvtParticle::initDecay(bool useMinMass) {
00238
00239 EvtParticle* p=this;
00240
00241 EvtParticle* par=p->getParent();
00242 double parMass=-1.;
00243 if ( par != 0 ) {
00244 if ( par->hasValidP4() ) parMass=par->mass();
00245 int i;
00246 for ( i=0;i<par->getNDaug();i++) {
00247 EvtParticle *tDaug=par->getDaug(i);
00248 if ( p != tDaug )
00249 parMass-=EvtPDL::getMinMass(tDaug->getId());
00250 }
00251 }
00252
00253 if ( _isInit ) {
00254
00255 if ( _ndaug>0) {
00256 int ii;
00257 for(ii=0;ii<_ndaug;ii++){
00258 if ( _ndaug==1 || EvtPDL::getWidth(p->getDaug(ii)->getId()) > 0.0000001)
00259 p->getDaug(ii)->initDecay(useMinMass);
00260 else p->getDaug(ii)->setMass(EvtPDL::getMeanMass(p->getDaug(ii)->getId()));
00261 }
00262 }
00263
00264 int j;
00265 EvtId *dauId=0;
00266 double *dauMasses=0;
00267 if ( _ndaug > 0) {
00268 dauId=new EvtId[_ndaug];
00269 dauMasses=new double[_ndaug];
00270 for (j=0;j<_ndaug;j++) {
00271 dauId[j]=p->getDaug(j)->getId();
00272 dauMasses[j]=p->getDaug(j)->mass();
00273 }
00274 }
00275 EvtId *parId=0;
00276 EvtId *othDauId=0;
00277 EvtParticle *tempPar=p->getParent();
00278 if (tempPar) {
00279 parId=new EvtId(tempPar->getId());
00280 if ( tempPar->getNDaug()==2 ) {
00281 if ( tempPar->getDaug(0) == this ) othDauId=new EvtId(tempPar->getDaug(1)->getId());
00282 else othDauId=new EvtId(tempPar->getDaug(0)->getId());
00283 }
00284 }
00285 if ( p->getParent() && _validP4==false ) {
00286 if ( !useMinMass ) p->setMass(EvtPDL::getRandMass(p->getId(),parId,_ndaug,dauId,othDauId,parMass,dauMasses));
00287 else p->setMass(EvtPDL::getMinMass(p->getId()));
00288 }
00289 if ( parId) delete parId;
00290 if ( othDauId) delete othDauId;
00291 if ( dauId) delete [] dauId;
00292 if ( dauMasses) delete [] dauMasses;
00293 return;
00294 }
00295
00296
00297
00298
00299 static EvtId BS0=EvtPDL::getId("B_s0");
00300 static EvtId BSB=EvtPDL::getId("anti-B_s0");
00301 static EvtId BD0=EvtPDL::getId("B0");
00302 static EvtId BDB=EvtPDL::getId("anti-B0");
00303
00304
00305 if ( (getNDaug()==0 && getParent()==0) && (getId()==BS0||getId()==BSB||getId()==BD0||getId()==BDB)){
00306 double t;
00307 int mix;
00308 EvtCPUtil::incoherentMix(getId(), t, mix);
00309 setLifetime(t);
00310
00311 if (mix) {
00312
00313 EvtScalarParticle* scalar_part;
00314
00315 scalar_part=new EvtScalarParticle;
00316 if (getId()==BS0) {
00317 EvtVector4R p_init(EvtPDL::getMass(BSB),0.0,0.0,0.0);
00318 scalar_part->init(BSB,p_init);
00319 }
00320 else if (getId()==BSB) {
00321 EvtVector4R p_init(EvtPDL::getMass(BS0),0.0,0.0,0.0);
00322 scalar_part->init(BS0,p_init);
00323 }
00324 else if (getId()==BD0) {
00325 EvtVector4R p_init(EvtPDL::getMass(BDB),0.0,0.0,0.0);
00326 scalar_part->init(BDB,p_init);
00327 }
00328 else if (getId()==BDB) {
00329 EvtVector4R p_init(EvtPDL::getMass(BD0),0.0,0.0,0.0);
00330 scalar_part->init(BD0,p_init);
00331 }
00332
00333 scalar_part->setLifetime(0);
00334
00335 scalar_part->setDiagonalSpinDensity();
00336
00337 insertDaugPtr(0,scalar_part);
00338
00339 _ndaug=1;
00340 _isInit=true;
00341 p=scalar_part;
00342 p->initDecay(useMinMass);
00343 return;
00344
00345
00346 }
00347 }
00348 if ( _ndaug==1 ) std::cout << "hi " << EvtPDL::name(this->getId()) << std::endl;
00349
00350 EvtDecayBase *decayer;
00351 decayer = EvtDecayTable::getDecayFunc(p);
00352
00353 if ( decayer ) {
00354 p->makeDaughters(decayer->nRealDaughters(),decayer->getDaugs());
00355
00356
00357 int i;
00358 for(i=0;i<p->getNDaug();i++){
00359 if ( EvtPDL::getWidth(p->getDaug(i)->getId()) > 0.0000001)
00360 p->getDaug(i)->initDecay(useMinMass);
00361 else p->getDaug(i)->setMass(EvtPDL::getMeanMass(p->getDaug(i)->getId()));
00362
00363
00364 }
00365 }
00366
00367
00368
00369 int j;
00370 EvtId *dauId=0;
00371 double *dauMasses=0;
00372 int nDaugT=p->getNDaug();
00373 if ( nDaugT > 0) {
00374 dauId=new EvtId[nDaugT];
00375 dauMasses=new double[nDaugT];
00376 for (j=0;j<nDaugT;j++) {
00377 dauId[j]=p->getDaug(j)->getId();
00378 dauMasses[j]=p->getDaug(j)->mass();
00379 }
00380 }
00381
00382 EvtId *parId=0;
00383 EvtId *othDauId=0;
00384 EvtParticle *tempPar=p->getParent();
00385 if (tempPar) {
00386 parId=new EvtId(tempPar->getId());
00387 if ( tempPar->getNDaug()==2 ) {
00388 if ( tempPar->getDaug(0) == this ) othDauId=new EvtId(tempPar->getDaug(1)->getId());
00389 else othDauId=new EvtId(tempPar->getDaug(0)->getId());
00390 }
00391 }
00392 if ( p->getParent() && p->hasValidP4()==false ) {
00393 if ( !useMinMass ) p->setMass(EvtPDL::getRandMass(p->getId(),parId,p->getNDaug(),dauId,othDauId,parMass,dauMasses));
00394 else p->setMass(EvtPDL::getMinMass(p->getId()));
00395 }
00396 if ( parId) delete parId;
00397 if ( othDauId) delete othDauId;
00398 if ( dauId) delete [] dauId;
00399 if ( dauMasses) delete [] dauMasses;
00400 _isInit=true;
00401 }
00402
00403
00404 void EvtParticle::decay(){
00405
00406
00407
00408 EvtParticle* p=this;
00409
00410
00411
00412
00413
00414
00415
00416
00417 EvtSpinDensity myRho;
00418 EvtDecayBase *decayer;
00419 decayer = EvtDecayTable::getDecayFunc(p);
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436 if ( _ndaug == 0 ) generateMassTree();
00437 static EvtId BS0=EvtPDL::getId("B_s0");
00438 static EvtId BSB=EvtPDL::getId("anti-B_s0");
00439 static EvtId BD0=EvtPDL::getId("B0");
00440 static EvtId BDB=EvtPDL::getId("anti-B0");
00441 if ( _ndaug==1 && (getId()==BS0||getId()==BSB||getId()==BD0||getId()==BDB) )
00442 {getDaug(0)->decay();
00443 std::cout<<"if ( _ndaug==1 && (getId()==BS0||getId()==BSB||getId()==BD0||getId()==BDB) )"<<endl;
00444 }
00445
00446 else{
00447
00448 if ( decayer ) {
00449 decayer->makeDecay(p);
00450
00451 }
00452 else{
00453 p->_rhoBackward.SetDiag(p->getSpinStates());
00454
00455 }
00456 }
00457 _isDecayed=true;
00458 return;
00459 }
00460
00461 void EvtParticle::generateMassTree() {
00462 double massProb=1.;
00463 double ranNum=2.;
00464 int counter=0;
00465 EvtParticle *p=this;
00466 while (massProb<ranNum) {
00467
00468 p->initDecay();
00469
00470 massProb=p->compMassProb();
00471 ranNum=EvtRandom::Flat();
00472
00473 counter++;
00474
00475 if ( counter > 10000 ) {
00476 if ( counter == 10001 ) {
00477 report(INFO,"EvtGen") << "Too many iterations to determine the mass tree. Parent mass= "<< p->mass() << _p<<" " << massProb <<endl;
00478 p->printTree();
00479 report(INFO,"EvtGen") << "will take next combo with non-zero likelihood\n";
00480 }
00481 if ( massProb>0. ) massProb=2.0;
00482 if ( counter > 20000 ) {
00483
00484 p->initDecay(true);
00485 p->printTree();
00486 massProb=p->compMassProb();
00487 if ( massProb>0. ) {
00488 massProb=2.0;
00489 report(INFO,"EvtGen") << "Taking the minimum mass of all particles in the chain\n";
00490 }
00491 else {
00492 report(INFO,"EvtGen") << "Sorry, no luck finding a valid set of masses. This may be a pathological combo\n";
00493 std::cout<<EvtPDL::name(p->getId())<<": Parent mass "<<p->getP4().mass()<<" with momentum "<<p->getP4()<<std::endl;
00494 if(EvtGlobalSet::ConExcPythia){EvtGlobalSet::ConExcPythia=0;return;}else{abort();}
00495
00496 }
00497 }
00498 }
00499 }
00500
00501
00502 }
00503
00504 double EvtParticle::compMassProb() {
00505
00506 EvtParticle *p=this;
00507
00508
00509 double mass=p->mass();
00510 double parMass=0.;
00511 if ( p->getParent()) {
00512 parMass=p->getParent()->mass();
00513 }
00514
00515 int nDaug=p->getNDaug();
00516 double *dMasses=0;
00517
00518 int i;
00519 if ( nDaug>0 ) {
00520 dMasses=new double[nDaug];
00521 for (i=0; i<nDaug; i++) dMasses[i]=p->getDaug(i)->mass();
00522 }
00523
00524 double temp=1.0;
00525 temp=EvtPDL::getMassProb(p->getId(), mass, parMass, nDaug, dMasses);
00526
00527
00528
00529
00530 if ( (!p->getParent() || _validP4 ) && temp>0.0 ) temp=1.;
00531
00532 delete [] dMasses;
00533
00534 for (i=0; i<nDaug; i++) {
00535 temp*=p->getDaug(i)->compMassProb();
00536 }
00537 return temp;
00538 }
00539
00540 void EvtParticle::deleteDaughters(bool keepChannel){
00541 int i;
00542
00543 for(i=0;i<_ndaug;i++){
00544 _daug[i]->deleteTree();
00545 }
00546
00547 _ndaug=0;
00548
00549 if ( !keepChannel) _channel=-10;
00550
00551
00552 _first=1;
00553 _isInit=false;
00554
00555 }
00556
00557 void EvtParticle::deleteTree(){
00558
00559 this->deleteDaughters();
00560
00561 delete this;
00562
00563 }
00564
00565 EvtVector4C EvtParticle::epsParent(int i) const {
00566 EvtVector4C temp;
00567 printParticle();
00568 report(ERROR,"EvtGen") << "and you have asked for the:"<<i
00569 <<"th polarization vector."
00570 <<" I.e. you thought it was a"
00571 <<" vector particle!" << endl;
00572 ::abort();
00573 return temp;
00574 }
00575
00576 EvtVector4C EvtParticle::eps(int i) const {
00577 EvtVector4C temp;
00578 printParticle();
00579 report(ERROR,"EvtGen") << "and you have asked for the:"<<i
00580 <<"th polarization vector."
00581 <<" I.e. you thought it was a"
00582 <<" vector particle!" << endl;
00583 ::abort();
00584 return temp;
00585 }
00586
00587 EvtVector4C EvtParticle::epsParentPhoton(int i){
00588 EvtVector4C temp;
00589 printParticle();
00590 report(ERROR,"EvtGen") << "and you have asked for the:"<<i
00591 <<"th polarization vector of photon."
00592 <<" I.e. you thought it was a"
00593 <<" photon particle!" << endl;
00594 ::abort();
00595 return temp;
00596 }
00597
00598 EvtVector4C EvtParticle::epsPhoton(int i){
00599 EvtVector4C temp;
00600 printParticle();
00601 report(ERROR,"EvtGen") << "and you have asked for the:"<<i
00602 <<"th polarization vector of a photon."
00603 <<" I.e. you thought it was a"
00604 <<" photon particle!" << endl;
00605 ::abort();
00606 return temp;
00607 }
00608
00609 EvtDiracSpinor EvtParticle::spParent(int i) const {
00610 EvtDiracSpinor tempD;
00611 int temp;
00612 temp = i;
00613 printParticle();
00614 report(ERROR,"EvtGen") << "and you have asked for the:"<<i
00615 <<"th dirac spinor."
00616 <<" I.e. you thought it was a"
00617 <<" Dirac particle!" << endl;
00618 ::abort();
00619 return tempD;
00620 }
00621
00622 EvtDiracSpinor EvtParticle::sp(int i) const {
00623 EvtDiracSpinor tempD;
00624 int temp;
00625 temp = i;
00626 printParticle();
00627 report(ERROR,"EvtGen") << "and you have asked for the:"<<i
00628 <<"th dirac spinor."
00629 <<" I.e. you thought it was a"
00630 <<" Dirac particle!" << endl;
00631 ::abort();
00632 return tempD;
00633 }
00634
00635 EvtDiracSpinor EvtParticle::spParentNeutrino() const {
00636 EvtDiracSpinor tempD;
00637 printParticle();
00638 report(ERROR,"EvtGen") << "and you have asked for the "
00639 <<"dirac spinor."
00640 <<" I.e. you thought it was a"
00641 <<" neutrino particle!" << endl;
00642 ::abort();
00643 return tempD;
00644 }
00645
00646 EvtDiracSpinor EvtParticle::spNeutrino() const {
00647 EvtDiracSpinor tempD;
00648 printParticle();
00649 report(ERROR,"EvtGen") << "and you have asked for the "
00650 <<"dirac spinor."
00651 <<" I.e. you thought it was a"
00652 <<" neutrino particle!" << endl;
00653 ::abort();
00654 return tempD;
00655 }
00656
00657 EvtTensor4C EvtParticle::epsTensorParent(int i) const {
00658 int temp;
00659 temp = i;
00660 EvtTensor4C tempC;
00661 printParticle();
00662 report(ERROR,"EvtGen") << "and you have asked for the:"<<i
00663 <<"th tensor."
00664 <<" I.e. you thought it was a"
00665 <<" Tensor particle!" << endl;
00666 ::abort();
00667 return tempC;
00668 }
00669
00670 EvtTensor4C EvtParticle::epsTensor(int i) const {
00671 int temp;
00672 temp = i;
00673 EvtTensor4C tempC;
00674 printParticle();
00675 report(ERROR,"EvtGen") << "and you have asked for the:"<<i
00676 <<"th tensor."
00677 <<" I.e. you thought it was a"
00678 <<" Tensor particle!" << endl;
00679 ::abort();
00680 return tempC;
00681 }
00682
00683
00684
00685 EvtVector4R EvtParticle::getP4Lab() {
00686 EvtVector4R temp,mom;
00687 EvtParticle *ptemp;
00688
00689 temp=this->getP4();
00690 ptemp=this;
00691
00692 while (ptemp->getParent()!=0) {
00693 ptemp=ptemp->getParent();
00694 mom=ptemp->getP4();
00695 temp=boostTo(temp,mom);
00696 }
00697 return temp;
00698 }
00699
00700 EvtVector4R EvtParticle::getP4Restframe() {
00701
00702 return EvtVector4R(mass(),0.0,0.0,0.0);
00703
00704 }
00705
00706 EvtVector4R EvtParticle::get4Pos() {
00707
00708 EvtVector4R temp,mom;
00709 EvtParticle *ptemp;
00710
00711 temp.set(0.0,0.0,0.0,0.0);
00712 ptemp=getParent();
00713
00714 if (ptemp==0) return temp;
00715
00716 temp=(ptemp->_t/ptemp->mass())*(ptemp->getP4());
00717
00718 while (ptemp->getParent()!=0) {
00719 ptemp=ptemp->getParent();
00720 mom=ptemp->getP4();
00721 temp=boostTo(temp,mom);
00722 temp=temp+(ptemp->_t/ptemp->mass())*(ptemp->getP4());
00723 }
00724
00725 return temp;
00726 }
00727
00728
00729 EvtParticle * EvtParticle::nextIter(EvtParticle *rootOfTree) {
00730
00731 EvtParticle *bpart;
00732 EvtParticle *current;
00733
00734 current=this;
00735 int i;
00736
00737 if (_ndaug!=0) return _daug[0];
00738
00739 do{
00740 bpart=current->_parent;
00741 if (bpart==0) return 0;
00742 i=0;
00743 while (bpart->_daug[i]!=current) {i++;}
00744
00745 if ( bpart==rootOfTree ) {
00746 if ( i== bpart->_ndaug-1 ) return 0;
00747 }
00748
00749 i++;
00750 current=bpart;
00751
00752 }while(i>=bpart->_ndaug);
00753
00754 return bpart->_daug[i];
00755
00756 }
00757
00758
00759 void EvtParticle::makeStdHep(EvtStdHep& stdhep,EvtSecondary& secondary,
00760 EvtId *list_of_stable){
00761
00762
00763 stdhep.createParticle(getP4Lab(),get4Pos(),-1,-1,
00764 EvtPDL::getStdHep(getId()));
00765
00766 int ii=0;
00767
00768
00769
00770
00771 while (list_of_stable[ii]!=EvtId(-1,-1)) {
00772 if (getId()==list_of_stable[ii]){
00773 secondary.createSecondary(0,this);
00774 return;
00775 }
00776 ii++;
00777 }
00778
00779
00780
00781
00782 int i;
00783 for(i=0;i<_ndaug;i++){
00784 stdhep.createParticle(_daug[i]->getP4Lab(),_daug[i]->get4Pos(),0,0,
00785 EvtPDL::getStdHep(_daug[i]->getId()));
00786 }
00787
00788 for(i=0;i<_ndaug;i++){
00789 _daug[i]->makeStdHepRec(1+i,1+i,stdhep,secondary,list_of_stable);
00790 }
00791 return;
00792
00793 }
00794
00795 void EvtParticle::makeStdHep(EvtStdHep& stdhep){
00796
00797
00798 stdhep.createParticle(getP4Lab(),get4Pos(),-1,-1,
00799 EvtPDL::getStdHep(getId()));
00800
00801 int i;
00802 for(i=0;i<_ndaug;i++){
00803 stdhep.createParticle(_daug[i]->getP4Lab(),_daug[i]->get4Pos(),0,0,
00804 EvtPDL::getStdHep(_daug[i]->getId()));
00805 }
00806
00807 for(i=0;i<_ndaug;i++){
00808 _daug[i]->makeStdHepRec(1+i,1+i,stdhep);
00809 }
00810 return;
00811
00812 }
00813
00814
00815 void EvtParticle::makeStdHepRec(int firstparent,int lastparent,
00816 EvtStdHep& stdhep,
00817 EvtSecondary& secondary,
00818 EvtId *list_of_stable){
00819
00820
00821 int ii=0;
00822
00823
00824
00825
00826 while (list_of_stable[ii]!=EvtId(-1,-1)) {
00827 if (getId()==list_of_stable[ii]){
00828 secondary.createSecondary(firstparent,this);
00829 return;
00830 }
00831 ii++;
00832 }
00833
00834
00835
00836 int i;
00837 int parent_num=stdhep.getNPart();
00838 for(i=0;i<_ndaug;i++){
00839 stdhep.createParticle(_daug[i]->getP4Lab(),_daug[i]->get4Pos(),
00840 firstparent,lastparent,
00841 EvtPDL::getStdHep(_daug[i]->getId()));
00842 }
00843
00844 for(i=0;i<_ndaug;i++){
00845 _daug[i]->makeStdHepRec(parent_num+i,parent_num+i,stdhep,
00846 secondary,list_of_stable);
00847 }
00848 return;
00849
00850 }
00851
00852 void EvtParticle::makeStdHepRec(int firstparent,int lastparent,
00853 EvtStdHep& stdhep){
00854
00855 int i;
00856 int parent_num=stdhep.getNPart();
00857 for(i=0;i<_ndaug;i++){
00858 stdhep.createParticle(_daug[i]->getP4Lab(),_daug[i]->get4Pos(),
00859 firstparent,lastparent,
00860 EvtPDL::getStdHep(_daug[i]->getId()));
00861 }
00862
00863 for(i=0;i<_ndaug;i++){
00864 _daug[i]->makeStdHepRec(parent_num+i,parent_num+i,stdhep);
00865 }
00866 return;
00867
00868 }
00869
00870 void EvtParticle::printTreeRec(int level) const {
00871
00872 int newlevel,i;
00873 newlevel = level +1;
00874
00875
00876 if (_ndaug!=0) {
00877 if ( level > 0 ) {
00878 for (i=0;i<(5*level);i++) {
00879 report(INFO,"") <<" ";
00880 }
00881 }
00882 report(INFO,"") << EvtPDL::name(_id).c_str();
00883 report(INFO,"") << " -> ";
00884 for(i=0;i<_ndaug;i++){
00885 report(INFO,"") << EvtPDL::name(_daug[i]->getId()).c_str()<<" ";
00886 }
00887 for(i=0;i<_ndaug;i++){
00888 report(INFO,"") << _daug[i]->mass()<<" ";
00889 }
00890 report(INFO,"")<<endl;
00891 for(i=0;i<_ndaug;i++){
00892 _daug[i]->printTreeRec(newlevel);
00893 }
00894 }
00895 }
00896
00897 void EvtParticle::printTree() const {
00898
00899 report(INFO,"EvtGen") << "This is the current decay chain"<<endl;
00900 report(INFO,"") << "This top particle is "<<
00901 EvtPDL::name(_id).c_str()<<endl;
00902
00903 this->printTreeRec(0);
00904 report(INFO,"EvtGen") << "End of decay chain."<<endl;
00905
00906 }
00907
00908 std::string EvtParticle::treeStrRec(int level) const {
00909
00910 int newlevel,i;
00911 newlevel = level +1;
00912
00913 std::string retval="";
00914
00915 for(i=0;i<_ndaug;i++){
00916 retval+=EvtPDL::name(_daug[i]->getId());
00917 if ( _daug[i]->getNDaug() > 0 ) {
00918 retval+= " (";
00919 retval+= _daug[i]->treeStrRec(newlevel);
00920 retval+= ") ";
00921 }
00922 else{
00923 if ( i!=_ndaug-1) retval+=" ";
00924 }
00925 }
00926
00927 return retval;
00928 }
00929
00930 std::string EvtParticle::writeTreeRec(std::string resonance) const {
00931 std::string retval="";
00932
00933 if (resonance == EvtPDL::name(_id).c_str() && _ndaug!=0) {
00934 retval=resonance+": "+IntToStr(_ndaug)+" ";
00935 for(int i=0;i<_ndaug;i++){
00936 retval += EvtPDL::name(_daug[i]->getId()).c_str();
00937 retval += " ";
00938 }
00939 }
00940
00941 for(int i=0;i<_ndaug;i++){
00942 _daug[i]->writeTreeRec(resonance);
00943 }
00944 std::cout<<retval;
00945 return retval;
00946 }
00947
00948 void EvtParticle::dumpTreeRec(int level,int dj) const {
00949
00950 int newlevel,i;
00951 newlevel = level +1;
00952
00953
00954 if (_ndaug!=0) {
00955
00956 int Ids= EvtPDL::getStdHep(_id);
00957 std::string c1,cid;
00958 if(Ids>0) c1="p";
00959 if(Ids<0) {
00960 c1="m";Ids=-1*Ids;
00961 }
00962
00963 cid=c1+IntToStr(Ids);
00964
00965 report(INFO,"") <<newlevel<<" "<< cid<<" "<<dj;
00966 report(INFO,"") << " = ";
00967
00968 int Nchannel=this->getChannel()+1;
00969 report(INFO,"") <<Nchannel<<endl;
00970
00971 for(i=0;i<_ndaug;i++){
00972 _daug[i]->dumpTreeRec(newlevel,i);
00973 }
00974 }
00975 }
00976
00977
00978 void EvtParticle::dumpTree() const {
00979
00980 report(INFO,"EvtGen") << "This is the current decay chain to dump"<<endl;
00981 report(INFO,"") << "This top particle is "<<
00982 EvtPDL::name(_id).c_str()<<endl;
00983
00984 this->dumpTreeRec(0,0);
00985 report(INFO,"EvtGen") << "End of decay chain."<<endl;
00986
00987 }
00988
00989
00990 std::string EvtParticle::treeStr() const {
00991
00992 std::string retval=EvtPDL::name(_id);
00993 retval+=" -> ";
00994
00995 retval+=treeStrRec(0);
00996
00997 return retval;
00998 }
00999
01000 void EvtParticle::printParticle() const {
01001
01002 switch (EvtPDL::getSpinType(_id)){
01003 case EvtSpinType::SCALAR:
01004 report(INFO,"EvtGen") << "This is a scalar particle:"<<EvtPDL::name(_id).c_str()<<"\n";
01005 break;
01006 case EvtSpinType::VECTOR:
01007 report(INFO,"EvtGen") << "This is a vector particle:"<<EvtPDL::name(_id).c_str()<<"\n";
01008 break;
01009 case EvtSpinType::TENSOR:
01010 report(INFO,"EvtGen") << "This is a tensor particle:"<<EvtPDL::name(_id).c_str()<<"\n";
01011 break;
01012 case EvtSpinType::DIRAC:
01013 report(INFO,"EvtGen") << "This is a dirac particle:"<<EvtPDL::name(_id).c_str()<<"\n";
01014 break;
01015 case EvtSpinType::PHOTON:
01016 report(INFO,"EvtGen") << "This is a photon:"<<EvtPDL::name(_id).c_str()<<"\n";
01017 break;
01018 case EvtSpinType::NEUTRINO:
01019 report(INFO,"EvtGen") << "This is a neutrino:"<<EvtPDL::name(_id).c_str()<<"\n";
01020 break;
01021 case EvtSpinType::RARITASCHWINGER:
01022 report(INFO,"EvtGen") << "This is a raritaschwinger:"<<EvtPDL::name(_id).c_str()<<"\n";
01023 break;
01024 case EvtSpinType::STRING:
01025 report(INFO,"EvtGen") << "This is a string:"<<EvtPDL::name(_id).c_str()<<"\n";
01026 break;
01027 default:
01028 report(INFO,"EvtGen") <<"Unknown particle type in EvtParticle::printParticle()"<<endl;
01029 break;
01030 }
01031 report(INFO,"EvtGen") << "Number of daughters:"<<_ndaug<<"\n";
01032
01033
01034 }
01035
01036
01037
01038 void init_vector( EvtParticle **part ){
01039 *part = (EvtParticle *) new EvtVectorParticle;
01040 }
01041
01042
01043 void init_scalar( EvtParticle **part ){
01044 *part = (EvtParticle *) new EvtScalarParticle;
01045 }
01046
01047 void init_tensor( EvtParticle **part ){
01048 *part = (EvtParticle *) new EvtTensorParticle;
01049 }
01050
01051 void init_dirac( EvtParticle **part ){
01052 *part = (EvtParticle *) new EvtDiracParticle;
01053 }
01054
01055 void init_photon( EvtParticle **part ){
01056 *part = (EvtParticle *) new EvtPhotonParticle;
01057 }
01058
01059 void init_raritaschinger( EvtParticle **part ){
01060 *part = (EvtParticle *) new EvtRaritaSchwingerParticle;
01061 }
01062
01063 void init_neutrino( EvtParticle **part ){
01064 *part = (EvtParticle *) new EvtNeutrinoParticle;
01065 }
01066
01067 void init_string( EvtParticle **part ){
01068 *part = (EvtParticle *) new EvtStringParticle;
01069 }
01070
01071 double EvtParticle::initializePhaseSpace(
01072 int numdaughter,EvtId *daughters, double poleSize,
01073 int whichTwo1, int whichTwo2) {
01074
01075 double m_b;
01076 int i;
01077
01078
01079
01080 static EvtVector4R p4[100];
01081 static double mass[100];
01082
01083 m_b = this->mass();
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094 bool resetDaughters=false;
01095 if ( numdaughter != this->getNDaug() && this->getNDaug() > 0 ) resetDaughters=true;
01096 if ( numdaughter == this->getNDaug() )
01097 for (i=0; i<numdaughter;i++) {
01098 if ( this->getDaug(i)->getId() != daughters[i] ) resetDaughters=true;
01099
01100 }
01101
01102 if ( resetDaughters ) {
01103
01104
01105
01106
01107 bool t1=true;
01108
01109 this->deleteDaughters(t1);
01110 this->makeDaughters(numdaughter,daughters);
01111 this->generateMassTree();
01112 }
01113
01114 double weight=0.;
01115
01116
01117
01118 for (i=0; i<numdaughter;i++) {
01119 mass[i]=this->getDaug(i)->mass();
01120
01121 }
01122
01123 if ( poleSize<-0.1) {
01124 EvtGenKine::PhaseSpace( numdaughter, mass, p4, m_b );
01125 for(i=0;i<numdaughter;i++){
01126 this->getDaug(i)->init(daughters[i],p4[i]);
01127 }
01128
01129 }
01130 else {
01131 if ( numdaughter != 3 ) {
01132 report(ERROR,"EvtGen") << "Only can generate pole phase space "
01133 << "distributions for 3 body final states"
01134 << endl<<"Will terminate."<<endl;
01135 ::abort();
01136 }
01137 bool ok=false;
01138 if ( (whichTwo1 == 1 && whichTwo2 == 0 ) ||
01139 (whichTwo1 == 0 && whichTwo2 == 1 ) ) {
01140 weight=EvtGenKine::PhaseSpacePole( m_b, mass[0], mass[1], mass[2],
01141 poleSize, p4);
01142
01143 this->getDaug(0)->init(daughters[0],p4[0]);
01144 this->getDaug(1)->init(daughters[1],p4[1]);
01145 this->getDaug(2)->init(daughters[2],p4[2]);
01146 ok=true;
01147 }
01148 if ( (whichTwo1 == 1 && whichTwo2 == 2 ) ||
01149 (whichTwo1 == 2 && whichTwo2 == 1 ) ) {
01150 weight=EvtGenKine::PhaseSpacePole( m_b, mass[2], mass[1], mass[0],
01151 poleSize, p4);
01152 this->getDaug(0)->init(daughters[0],p4[2]);
01153 this->getDaug(1)->init(daughters[1],p4[1]);
01154 this->getDaug(2)->init(daughters[2],p4[0]);
01155 ok=true;
01156 }
01157 if ( (whichTwo1 == 0 && whichTwo2 == 2 ) ||
01158 (whichTwo1 == 2 && whichTwo2 == 0 ) ) {
01159 weight=EvtGenKine::PhaseSpacePole( m_b, mass[1], mass[0], mass[2],
01160 poleSize, p4);
01161 this->getDaug(0)->init(daughters[0],p4[1]);
01162 this->getDaug(1)->init(daughters[1],p4[0]);
01163 this->getDaug(2)->init(daughters[2],p4[2]);
01164 ok=true;
01165 }
01166 if ( !ok) {
01167 report(ERROR,"EvtGen") << "Invalid pair of particle to generate a pole dist"
01168 << whichTwo1 << " " << whichTwo2
01169 << endl<<"Will terminate."<<endl;
01170 ::abort();
01171 }
01172 }
01173
01174 return weight;
01175 }
01176
01177 void EvtParticle::makeDaughters( int ndaugstore, EvtId *id){
01178
01179 int i;
01180 if ( _channel < 0 ) {
01181
01182 setChannel(0);
01183 }
01184 EvtParticle* pdaug;
01185 if (_ndaug!=0 ){
01186 if (_ndaug!=ndaugstore){
01187 report(ERROR,"EvtGen") << "Asking to make a different number of "
01188 << "daughters than what was previously created."
01189 << endl<<"Will terminate."<<endl;
01190 ::abort();
01191 }
01192 }
01193 else{
01194 for(i=0;i<ndaugstore;i++){
01195 pdaug=EvtParticleFactory::particleFactory(EvtPDL::getSpinType(id[i]));
01196 pdaug->setId(id[i]);
01197 pdaug->addDaug(this);
01198 }
01199
01200 }
01201 }
01202
01203
01204 void EvtParticle::setDecayProb(double prob) {
01205
01206 if ( _decayProb == 0 ) _decayProb=new double;
01207 *_decayProb=prob;
01208 }
01209
01210
01211 std::string IntToStr( int a)
01212 {
01213 std::string ans;
01214 std::string ans1;
01215 int k = 10 ;
01216 while (a > 0 )
01217 {
01218 ans += char (a % 10 + 48 );
01219 a /= 10 ;
01220 }
01221 for ( int i = ans.size() - 1 ; i >= 0 ; i -- )
01222 {
01223 ans1 += ans[i];
01224 }
01225 return ans1;
01226 }