#include <EvtRaritaSchwingerParticle.hh>
Inheritance diagram for EvtRaritaSchwingerParticle:
|
|
00038 {
00039
00040 return;
00041
00042 }
|
|
|
00034 {
00035 }
|
|
|
|
|
|
Add another daughter to the particle
|
|
|
00480 {
00481
00482 EvtParticle *p=this;
00483 //report(INFO,"EvtGen") << "compMassProb " << endl;
00484 //p->printTree();
00485 double mass=p->mass();
00486 double parMass=0.;
00487 if ( p->getParent()) {
00488 parMass=p->getParent()->mass();
00489 }
00490
00491 int nDaug=p->getNDaug();
00492 double *dMasses=0;
00493
00494 int i;
00495 if ( nDaug>0 ) {
00496 dMasses=new double[nDaug];
00497 for (i=0; i<nDaug; i++) dMasses[i]=p->getDaug(i)->mass();
00498 }
00499
00500 double temp=1.0;
00501 temp=EvtPDL::getMassProb(p->getId(), mass, parMass, nDaug, dMasses);
00502
00503 //report(INFO,"EvtGen") << temp << " " << EvtPDL::name(p->getId()) << endl;
00504 //If the particle already has a mass, we dont need to include
00505 //it in the probability calculation
00506 if ( (!p->getParent() || _validP4 ) && temp>0.0 ) temp=1.;
00507
00508 delete [] dMasses;
00509 // if ( temp < 0.9999999 )
00510 for (i=0; i<nDaug; i++) {
00511 temp*=p->getDaug(i)->compMassProb();
00512 }
00513 return temp;
00514 }
|
|
|
Decay particle 00383 {
00384 //P is particle to decay, typically 'this' but sometime
00385 //modified by mixing
00386 EvtParticle* p=this;
00387 //Did it mix?
00388 //if ( p->getMixed() ) {
00389 //should take C(p) - this should only
00390 //happen the first time we call decay for this
00391 //particle
00392 //p->takeCConj();
00393 // p->setUnMixed();
00394 //}
00395 EvtSpinDensity myRho; //pingrg test code
00396 EvtDecayBase *decayer;
00397 decayer = EvtDecayTable::getDecayFunc(p);
00398 // if ( decayer ) {
00399 // report(INFO,"EvtGen") << "calling decay for " << EvtPDL::name(p->getId()) << " " << p->mass() << " " << p->getP4() << " " << p->getNDaug() << " " << p << endl;
00400 // report(INFO,"EvtGen") << "NDaug= " << decayer->getNDaug() << endl;
00401 // int ti;
00402 // for ( ti=0; ti<decayer->getNDaug(); ti++)
00403 // report(INFO,"EvtGen") << "Daug " << ti << " " << EvtPDL::name(decayer->getDaug(ti)) << endl;
00404 // }
00405 //if (p->_ndaug>0) {
00406 // report(INFO,"EvtGen") <<"Is decaying particle with daughters!!!!!"<<endl;
00407 // ::abort();
00408 //return;
00409 //call initdecay first - April 29,2002 - Lange
00410 //}
00411
00412 //if there are already daughters, then this step is already done!
00413 // figure out the masses
00414 if ( _ndaug == 0 ) generateMassTree();
00415 static EvtId BS0=EvtPDL::getId("B_s0");
00416 static EvtId BSB=EvtPDL::getId("anti-B_s0");
00417 static EvtId BD0=EvtPDL::getId("B0");
00418 static EvtId BDB=EvtPDL::getId("anti-B0");
00419 if ( _ndaug==1 && (getId()==BS0||getId()==BSB||getId()==BD0||getId()==BDB) )
00420 {getDaug(0)->decay();
00421 std::cout<<"if ( _ndaug==1 && (getId()==BS0||getId()==BSB||getId()==BD0||getId()==BDB) )"<<endl;
00422 }
00423
00424 else{
00425 //now we have accepted a set of masses - time
00426 if ( decayer ) {
00427 decayer->makeDecay(p);
00428
00429 }
00430 else{
00431 p->_rhoBackward.SetDiag(p->getSpinStates());
00432
00433 }
00434 }
00435 _isDecayed=true;
00436 return;
00437 }
|
|
|
00375 {return _decayProb;}
|
|
|
00516 {
00517 int i;
00518
00519 for(i=0;i<_ndaug;i++){
00520 _daug[i]->deleteTree();
00521 }
00522
00523 _ndaug=0;
00524 //if ( keepChannel ) report(INFO,"EvtGen") << "keeping \n";
00525 if ( !keepChannel) _channel=-10;
00526 //_t=0.0;
00527 //_genlifetime=1;
00528 _first=1;
00529 _isInit=false;
00530 //report(INFO,"EvtGen") << "calling deletedaughters " << EvtPDL::name(this->getId()) <<endl;
00531 }
|
|
|
Delete a decay chain 00533 {
00534
00535 this->deleteDaughters();
00536
00537 delete this;
00538
00539 }
|
|
|
00954 { //pingrg, Mar. 25,2008
00955
00956 report(INFO,"EvtGen") << "This is the current decay chain to dump"<<endl;
00957 report(INFO,"") << "This top particle is "<<
00958 EvtPDL::name(_id).c_str()<<endl;
00959
00960 this->dumpTreeRec(0,0);
00961 report(INFO,"EvtGen") << "End of decay chain."<<endl;
00962
00963 }
|
|
||||||||||||
|
00924 { //pingrg, Mar. 25,2008
00925
00926 int newlevel,i;
00927 newlevel = level +1;
00928
00929
00930 if (_ndaug!=0) {
00931
00932 int Ids= EvtPDL::getStdHep(_id);
00933 std::string c1,cid;
00934 if(Ids>0) c1="p";
00935 if(Ids<0) {
00936 c1="m";Ids=-1*Ids;
00937 }
00938
00939 cid=c1+IntToStr(Ids);
00940
00941 report(INFO,"") <<newlevel<<" "<< cid<<" "<<dj;
00942 report(INFO,"") << " = ";
00943
00944 int Nchannel=this->getChannel()+1;
00945 report(INFO,"") <<Nchannel<<endl;
00946
00947 for(i=0;i<_ndaug;i++){
00948 _daug[i]->dumpTreeRec(newlevel,i);
00949 }
00950 }
00951 }
|
|
|
Returns polarization vector in the particles own restframe. Reimplemented in EvtVectorParticle. 00552 {
00553 EvtVector4C temp;
00554 printParticle();
00555 report(ERROR,"EvtGen") << "and you have asked for the:"<<i
00556 <<"th polarization vector."
00557 <<" I.e. you thought it was a"
00558 <<" vector particle!" << endl;
00559 ::abort();
00560 return temp;
00561 }
|
|
|
Returns polarization vector in the parents restframe. Reimplemented in EvtVectorParticle. 00541 {
00542 EvtVector4C temp;
00543 printParticle();
00544 report(ERROR,"EvtGen") << "and you have asked for the:"<<i
00545 <<"th polarization vector."
00546 <<" I.e. you thought it was a"
00547 <<" vector particle!" << endl;
00548 ::abort();
00549 return temp;
00550 }
|
|
|
Returns polarization vector in the parents restframe for a photon. Reimplemented in EvtPhotonParticle. 00563 {
00564 EvtVector4C temp;
00565 printParticle();
00566 report(ERROR,"EvtGen") << "and you have asked for the:"<<i
00567 <<"th polarization vector of photon."
00568 <<" I.e. you thought it was a"
00569 <<" photon particle!" << endl;
00570 ::abort();
00571 return temp;
00572 }
|
|
|
Returns polarization vector in the particles own restframe for a photon. Reimplemented in EvtPhotonParticle. 00574 {
00575 EvtVector4C temp;
00576 printParticle();
00577 report(ERROR,"EvtGen") << "and you have asked for the:"<<i
00578 <<"th polarization vector of a photon."
00579 <<" I.e. you thought it was a"
00580 <<" photon particle!" << endl;
00581 ::abort();
00582 return temp;
00583 }
|
|
|
Returns tensor in the particles own restframe for a spin 2 particle. Reimplemented in EvtTensorParticle. 00646 {
00647 int temp;
00648 temp = i;
00649 EvtTensor4C tempC;
00650 printParticle();
00651 report(ERROR,"EvtGen") << "and you have asked for the:"<<i
00652 <<"th tensor."
00653 <<" I.e. you thought it was a"
00654 <<" Tensor particle!" << endl;
00655 ::abort();
00656 return tempC;
00657 }
|
|
|
Returns tensor in the parents restframe for a spin 2 particle. Reimplemented in EvtTensorParticle. 00633 {
00634 int temp;
00635 temp = i;
00636 EvtTensor4C tempC;
00637 printParticle();
00638 report(ERROR,"EvtGen") << "and you have asked for the:"<<i
00639 <<"th tensor."
00640 <<" I.e. you thought it was a"
00641 <<" Tensor particle!" << endl;
00642 ::abort();
00643 return tempC;
00644 }
|
|
|
Used internally to decide if first time particle is decayed. 00111 { return _first;}
|
|
|
00439 {
00440 double massProb=1.;
00441 double ranNum=2.;
00442 int counter=0;
00443 EvtParticle *p=this;
00444 while (massProb<ranNum) {
00445 //check it out the first time.
00446 p->initDecay();
00447 //report(INFO,"EvtGen") << "calling massProb \n";
00448 massProb=p->compMassProb();
00449 ranNum=EvtRandom::Flat();
00450 // report(INFO,"EvtGen") << "end of iter " << massProb <<endl;
00451 counter++;
00452
00453 if ( counter > 10000 ) {
00454 if ( counter == 10001 ) {
00455 report(INFO,"EvtGen") << "Too many iterations to determine the mass tree. Parent mass= "<< p->mass() << " " << massProb <<endl;
00456 p->printTree();
00457 report(INFO,"EvtGen") << "will take next combo with non-zero likelihood\n";
00458 }
00459 if ( massProb>0. ) massProb=2.0;
00460 if ( counter > 20000 ) {
00461 // one last try - take the minimum masses
00462 p->initDecay(true);
00463 p->printTree();
00464 massProb=p->compMassProb();
00465 if ( massProb>0. ) {
00466 massProb=2.0;
00467 report(INFO,"EvtGen") << "Taking the minimum mass of all particles in the chain\n";
00468 }
00469 else {
00470 report(INFO,"EvtGen") << "Sorry, no luck finding a valid set of masses. This may be a pathological combo\n";
00471 assert(0);
00472 }
00473 }
00474 }
00475 }
00476 // report(INFO,"EvtGen") << counter << endl;
00477 // p->printTree();
00478 }
|
|
|
Returns the 4position of the particle in the lab frame. 00682 {
00683
00684 EvtVector4R temp,mom;
00685 EvtParticle *ptemp;
00686
00687 temp.set(0.0,0.0,0.0,0.0);
00688 ptemp=getParent();
00689
00690 if (ptemp==0) return temp;
00691
00692 temp=(ptemp->_t/ptemp->mass())*(ptemp->getP4());
00693
00694 while (ptemp->getParent()!=0) {
00695 ptemp=ptemp->getParent();
00696 mom=ptemp->getP4();
00697 temp=boostTo(temp,mom);
00698 temp=temp+(ptemp->_t/ptemp->mass())*(ptemp->getP4());
00699 }
00700
00701 return temp;
00702 }
|
|
|
Retunrs the decay channel. 00123 { return _channel;}
|
|
|
Get pointer the the i:th daugther. 00085 { return _daug[i]; }
|
|
|
get generator information; pingrg-2011-1-6 00146 { return _generatorFlag ;}
|
|
|
Returns Id of particle. 00113 { return _id;}
|
|
|
Returns the lifetime. 00099 {
00100
00101 return _t;
00102 }
|
|
|
Returns number of daugthers. 00125 { return _ndaug;}
|
|
|
Returns 4momentum in parents restframe. 00121 { return _p;}
|
|
|
Gets 4vector in the labframe, i.e., the frame in which the root particles momentum is measured. 00661 {
00662 EvtVector4R temp,mom;
00663 EvtParticle *ptemp;
00664
00665 temp=this->getP4();
00666 ptemp=this;
00667
00668 while (ptemp->getParent()!=0) {
00669 ptemp=ptemp->getParent();
00670 mom=ptemp->getP4();
00671 temp=boostTo(temp,mom);
00672 }
00673 return temp;
00674 }
|
|
|
Gets 4vector in the particles restframe, i.e. this functiont will return (m,0,0,0) 00676 {
00677
00678 return EvtVector4R(mass(),0.0,0.0,0.0);
00679
00680 }
|
|
|
Returns pointer to parent particle. 00087 { return _parent;}
|
|
|
Get backward spin density matrix. 00342 {return _rhoBackward;}
|
|
|
Get forward spin density matrix. 00332 {return _rhoForward;}
|
|
|
Returns number of spin states of the particle. 00119 { return EvtSpinType::getSpinStates(EvtPDL::getSpinType(_id));}
|
|
|
Returns particle type. 00116 { return EvtPDL::getSpinType(_id);}
|
|
|
00368 {return _validP4;}
|
|
||||||||||||
|
Initialiaze particle with id and 4momentum. Implements EvtParticle. 00044 {
00045
00046 _validP4=true;
00047 setp(p4);
00048 setpart_num(id);
00049
00050 if (EvtPDL::getStdHep(id)==0){
00051 report(ERROR,"EvtGen") << "Error in EvtRaritaSchwingerParticle::init, part_n="
00052 << id.getId()<<endl;
00053 ::abort();
00054 }
00055
00056 static EvtVector4R e1(0.0,1.0,0.0,0.0);
00057 static EvtVector4R e2(0.0,0.0,1.0,0.0);
00058 static EvtVector4R e3(0.0,0.0,0.0,1.0);
00059
00060
00061 if (EvtPDL::getStdHep(id)>0){
00062
00063 EvtDiracSpinor u1,u2;
00064
00065 u1.set(EvtComplex(sqrt(2.0*mass()),0.0),EvtComplex(0.0,0.0),
00066 EvtComplex(0.0,0.0),EvtComplex(0.0,0.0));
00067 u2.set(EvtComplex(0.0,0.0),EvtComplex(sqrt(2.0*mass()),0.0),
00068 EvtComplex(0.0,0.0),EvtComplex(0.0,0.0));
00069
00070
00071
00072 _spinorRest[0]=dirProd(e3,u1+u2);
00073 _spinorRest[1]=dirProd(e1+EvtComplex(0.0,1.0)*e2,u1);
00074 _spinorRest[2]=dirProd(e1-EvtComplex(0.0,1.0)*e2,u2);
00075 _spinorRest[3]=dirProd(e1,(u1+u2))+dirProd(EvtComplex(0.0,1.0)*e2,(u1-u2));
00076
00077 _spinor[0]=boostTo(_spinorRest[0],p4);
00078 _spinor[1]=boostTo(_spinorRest[1],p4);
00079 _spinor[2]=boostTo(_spinorRest[2],p4);
00080 _spinor[3]=boostTo(_spinorRest[3],p4);
00081
00082 }
00083 else{
00084
00085 EvtDiracSpinor u1,u2;
00086
00087 u1.set(EvtComplex(0.0,0.0),EvtComplex(0.0,0.0),
00088 EvtComplex(sqrt(2.0*mass()),0.0),EvtComplex(0.0,0.0));
00089 u2.set(EvtComplex(0.0,0.0),EvtComplex(0.0,0.0),
00090 EvtComplex(0.0,0.0),EvtComplex(sqrt(2.0*mass()),0.0));
00091
00092
00093
00094 _spinorRest[0]=dirProd(e3,(u1+u2));
00095 _spinorRest[1]=dirProd(e1+EvtComplex(0.0,1.0)*e2,u1);
00096 _spinorRest[2]=dirProd(e1-EvtComplex(0.0,1.0)*e2,u2);
00097 _spinorRest[3]=dirProd(e1,(u1+u2))+dirProd(EvtComplex(0.0,1.0)*e2,(u1-u2));
00098
00099 _spinor[0]=boostTo(_spinorRest[0],p4);
00100 _spinor[1]=boostTo(_spinorRest[1],p4);
00101 _spinor[2]=boostTo(_spinorRest[2],p4);
00102 _spinor[3]=boostTo(_spinorRest[3],p4);
00103
00104 }
00105
00106 setLifetime();
00107 }
|
|
|
00216 {
00217
00218 EvtParticle* p=this;
00219 // carefull - the parent mass might be fixed in stone..
00220 EvtParticle* par=p->getParent();
00221 double parMass=-1.;
00222 if ( par != 0 ) {
00223 if ( par->hasValidP4() ) parMass=par->mass();
00224 int i;
00225 for ( i=0;i<par->getNDaug();i++) {
00226 EvtParticle *tDaug=par->getDaug(i);
00227 if ( p != tDaug )
00228 parMass-=EvtPDL::getMinMass(tDaug->getId());
00229 }
00230 }
00231
00232 if ( _isInit ) {
00233 //we have already been here - just reroll the masses!
00234 if ( _ndaug>0) {
00235 int ii;
00236 for(ii=0;ii<_ndaug;ii++){
00237 if ( _ndaug==1 || EvtPDL::getWidth(p->getDaug(ii)->getId()) > 0.0000001)
00238 p->getDaug(ii)->initDecay(useMinMass);
00239 else p->getDaug(ii)->setMass(EvtPDL::getMeanMass(p->getDaug(ii)->getId()));
00240 }
00241 }
00242
00243 int j;
00244 EvtId *dauId=0;
00245 double *dauMasses=0;
00246 if ( _ndaug > 0) {
00247 dauId=new EvtId[_ndaug];
00248 dauMasses=new double[_ndaug];
00249 for (j=0;j<_ndaug;j++) {
00250 dauId[j]=p->getDaug(j)->getId();
00251 dauMasses[j]=p->getDaug(j)->mass();
00252 }
00253 }
00254 EvtId *parId=0;
00255 EvtId *othDauId=0;
00256 EvtParticle *tempPar=p->getParent();
00257 if (tempPar) {
00258 parId=new EvtId(tempPar->getId());
00259 if ( tempPar->getNDaug()==2 ) {
00260 if ( tempPar->getDaug(0) == this ) othDauId=new EvtId(tempPar->getDaug(1)->getId());
00261 else othDauId=new EvtId(tempPar->getDaug(0)->getId());
00262 }
00263 }
00264 if ( p->getParent() && _validP4==false ) {
00265 if ( !useMinMass ) p->setMass(EvtPDL::getRandMass(p->getId(),parId,_ndaug,dauId,othDauId,parMass,dauMasses));
00266 else p->setMass(EvtPDL::getMinMass(p->getId()));
00267 }
00268 if ( parId) delete parId;
00269 if ( othDauId) delete othDauId;
00270 if ( dauId) delete [] dauId;
00271 if ( dauMasses) delete [] dauMasses;
00272 return;
00273 }
00274
00275
00276 //Will include effects of mixing here
00277 //added by Lange Jan4,2000
00278 static EvtId BS0=EvtPDL::getId("B_s0");
00279 static EvtId BSB=EvtPDL::getId("anti-B_s0");
00280 static EvtId BD0=EvtPDL::getId("B0");
00281 static EvtId BDB=EvtPDL::getId("anti-B0");
00282
00283 //only makes sense if there is no parent particle
00284 if ( (getNDaug()==0 && getParent()==0) && (getId()==BS0||getId()==BSB||getId()==BD0||getId()==BDB)){
00285 double t;
00286 int mix;
00287 EvtCPUtil::incoherentMix(getId(), t, mix);
00288 setLifetime(t);
00289
00290 if (mix) {
00291
00292 EvtScalarParticle* scalar_part;
00293
00294 scalar_part=new EvtScalarParticle;
00295 if (getId()==BS0) {
00296 EvtVector4R p_init(EvtPDL::getMass(BSB),0.0,0.0,0.0);
00297 scalar_part->init(BSB,p_init);
00298 }
00299 else if (getId()==BSB) {
00300 EvtVector4R p_init(EvtPDL::getMass(BS0),0.0,0.0,0.0);
00301 scalar_part->init(BS0,p_init);
00302 }
00303 else if (getId()==BD0) {
00304 EvtVector4R p_init(EvtPDL::getMass(BDB),0.0,0.0,0.0);
00305 scalar_part->init(BDB,p_init);
00306 }
00307 else if (getId()==BDB) {
00308 EvtVector4R p_init(EvtPDL::getMass(BD0),0.0,0.0,0.0);
00309 scalar_part->init(BD0,p_init);
00310 }
00311
00312 scalar_part->setLifetime(0);
00313
00314 scalar_part->setDiagonalSpinDensity();
00315
00316 insertDaugPtr(0,scalar_part);
00317
00318 _ndaug=1;
00319 _isInit=true;
00320 p=scalar_part;
00321 p->initDecay(useMinMass);
00322 return;
00323
00324
00325 }
00326 }
00327 if ( _ndaug==1 ) std::cout << "hi " << EvtPDL::name(this->getId()) << std::endl;
00328
00329 EvtDecayBase *decayer;
00330 decayer = EvtDecayTable::getDecayFunc(p);
00331
00332 if ( decayer ) {
00333 p->makeDaughters(decayer->nRealDaughters(),decayer->getDaugs());
00334 // report(INFO,"EvtGen") << "has found decay " << decayer->nRealDaughters() << endl;
00335 //then loop over the daughters and init their decay
00336 int i;
00337 for(i=0;i<p->getNDaug();i++){
00338 if ( EvtPDL::getWidth(p->getDaug(i)->getId()) > 0.0000001)
00339 p->getDaug(i)->initDecay(useMinMass);
00340 else p->getDaug(i)->setMass(EvtPDL::getMeanMass(p->getDaug(i)->getId()));
00341 // report(INFO,"EvtGen") << "has inited " << EvtPDL::name(p->getDaug(i)->getId()) <<
00342 // " " << EvtPDL::name(p->getId()) << endl;
00343 }
00344 }
00345 //figure masses in separate step...
00346 // if ( p->getParent() && _validP4==false ) EvtDecayBase::findMass(p);
00347
00348 int j;
00349 EvtId *dauId=0;
00350 double *dauMasses=0;
00351 int nDaugT=p->getNDaug();
00352 if ( nDaugT > 0) {
00353 dauId=new EvtId[nDaugT];
00354 dauMasses=new double[nDaugT];
00355 for (j=0;j<nDaugT;j++) {
00356 dauId[j]=p->getDaug(j)->getId();
00357 dauMasses[j]=p->getDaug(j)->mass();
00358 }
00359 }
00360
00361 EvtId *parId=0;
00362 EvtId *othDauId=0;
00363 EvtParticle *tempPar=p->getParent();
00364 if (tempPar) {
00365 parId=new EvtId(tempPar->getId());
00366 if ( tempPar->getNDaug()==2 ) {
00367 if ( tempPar->getDaug(0) == this ) othDauId=new EvtId(tempPar->getDaug(1)->getId());
00368 else othDauId=new EvtId(tempPar->getDaug(0)->getId());
00369 }
00370 }
00371 if ( p->getParent() && p->hasValidP4()==false ) {
00372 if ( !useMinMass ) p->setMass(EvtPDL::getRandMass(p->getId(),parId,p->getNDaug(),dauId,othDauId,parMass,dauMasses));
00373 else p->setMass(EvtPDL::getMinMass(p->getId()));
00374 }
00375 if ( parId) delete parId;
00376 if ( othDauId) delete othDauId;
00377 if ( dauId) delete [] dauId;
00378 if ( dauMasses) delete [] dauMasses;
00379 _isInit=true;
00380 }
|
|
||||||||||||||||||||||||
|
Similar to the routine above except that here momentum is generated according to phase space daughters are filled with this momentum. 01049 {
01050
01051 double m_b;
01052 int i;
01053 //lange
01054 // this->makeDaughters(numdaughter,daughters);
01055
01056 static EvtVector4R p4[100];
01057 static double mass[100];
01058
01059 m_b = this->mass();
01060
01061 //lange - Jan2,2002 - Need to check to see if the daughters of the parent
01062 // have changed. If so, delete them and start over.
01063 //report(INFO,"EvtGen") << "the parent is\n";
01064 //if ( this->getParent() ) {
01065 // if ( this->getParent()->getParent() ) this->getParent()->getParent()->printTree();
01066 // this->getParent()->printTree();
01067 //}
01068 //report(INFO,"EvtGen") << "and this is\n";
01069 //if ( this) this->printTree();
01070 bool resetDaughters=false;
01071 if ( numdaughter != this->getNDaug() && this->getNDaug() > 0 ) resetDaughters=true;
01072 if ( numdaughter == this->getNDaug() )
01073 for (i=0; i<numdaughter;i++) {
01074 if ( this->getDaug(i)->getId() != daughters[i] ) resetDaughters=true;
01075 //report(INFO,"EvtGen") << this->getDaug(i)->getId() << " " << daughters[i] << endl;
01076 }
01077
01078 if ( resetDaughters ) {
01079 // report(INFO,"EvtGen") << "reseting daughters\n";
01080 //for (i=0; i<numdaughter;i++) {
01081 // report(INFO,"EvtGen") << "reset " <<i<< " "<< EvtPDL::name(this->getDaug(i)->getId()) << " " << EvtPDL::name(daughters[i]) << endl;
01082 //}
01083 bool t1=true;
01084 //but keep the decay channel of the parent.
01085 this->deleteDaughters(t1);
01086 this->makeDaughters(numdaughter,daughters);
01087 this->generateMassTree();
01088 }
01089
01090 double weight=0.;
01091 // EvtDecayBase::findMasses( this, numdaughter, daughters, mass );
01092 //get the list of masses
01093 //report(INFO,"EvtGen") << "mpar= " << m_b << " " << this <<endl;
01094 for (i=0; i<numdaughter;i++) {
01095 mass[i]=this->getDaug(i)->mass();
01096 // report(INFO,"EvtGen") << "mass " << i << " " << mass[i] << " " << this->getDaug(i) << endl;
01097 }
01098
01099 if ( poleSize<-0.1) {
01100 EvtGenKine::PhaseSpace( numdaughter, mass, p4, m_b );
01101 for(i=0;i<numdaughter;i++){
01102 this->getDaug(i)->init(daughters[i],p4[i]);
01103 }
01104
01105 }
01106 else {
01107 if ( numdaughter != 3 ) {
01108 report(ERROR,"EvtGen") << "Only can generate pole phase space "
01109 << "distributions for 3 body final states"
01110 << endl<<"Will terminate."<<endl;
01111 ::abort();
01112 }
01113 bool ok=false;
01114 if ( (whichTwo1 == 1 && whichTwo2 == 0 ) ||
01115 (whichTwo1 == 0 && whichTwo2 == 1 ) ) {
01116 weight=EvtGenKine::PhaseSpacePole( m_b, mass[0], mass[1], mass[2],
01117 poleSize, p4);
01118 //report(INFO,"EvtGen") << "here " << weight << " " << poleSize << endl;
01119 this->getDaug(0)->init(daughters[0],p4[0]);
01120 this->getDaug(1)->init(daughters[1],p4[1]);
01121 this->getDaug(2)->init(daughters[2],p4[2]);
01122 ok=true;
01123 }
01124 if ( (whichTwo1 == 1 && whichTwo2 == 2 ) ||
01125 (whichTwo1 == 2 && whichTwo2 == 1 ) ) {
01126 weight=EvtGenKine::PhaseSpacePole( m_b, mass[2], mass[1], mass[0],
01127 poleSize, p4);
01128 this->getDaug(0)->init(daughters[0],p4[2]);
01129 this->getDaug(1)->init(daughters[1],p4[1]);
01130 this->getDaug(2)->init(daughters[2],p4[0]);
01131 ok=true;
01132 }
01133 if ( (whichTwo1 == 0 && whichTwo2 == 2 ) ||
01134 (whichTwo1 == 2 && whichTwo2 == 0 ) ) {
01135 weight=EvtGenKine::PhaseSpacePole( m_b, mass[1], mass[0], mass[2],
01136 poleSize, p4);
01137 this->getDaug(0)->init(daughters[0],p4[1]);
01138 this->getDaug(1)->init(daughters[1],p4[0]);
01139 this->getDaug(2)->init(daughters[2],p4[2]);
01140 ok=true;
01141 }
01142 if ( !ok) {
01143 report(ERROR,"EvtGen") << "Invalid pair of particle to generate a pole dist"
01144 << whichTwo1 << " " << whichTwo2
01145 << endl<<"Will terminate."<<endl;
01146 ::abort();
01147 }
01148 }
01149
01150 return weight;
01151 }
|
|
||||||||||||
|
Makes partptr the idaug:th daugther. 00206 { _daug[idaug]=partptr;
00207 partptr->_parent=this; }
|
|
|
00369 {return _isDecayed;}
|
|
|
00367 {return _isInit;}
|
|
||||||||||||
|
Creates the daughters in the list of ids and adds them to the parent. Note that momentum is left uninitialized, this is _only_ creation. 01153 {
01154
01155 int i;
01156 if ( _channel < 0 ) {
01157 //report(INFO,"EvtGen") << "setting channel " << EvtPDL::name(this->getId()) << endl;
01158 setChannel(0);
01159 }
01160 EvtParticle* pdaug;
01161 if (_ndaug!=0 ){
01162 if (_ndaug!=ndaugstore){
01163 report(ERROR,"EvtGen") << "Asking to make a different number of "
01164 << "daughters than what was previously created."
01165 << endl<<"Will terminate."<<endl;
01166 ::abort();
01167 }
01168 }
01169 else{
01170 for(i=0;i<ndaugstore;i++){
01171 pdaug=EvtParticleFactory::particleFactory(EvtPDL::getSpinType(id[i]));
01172 pdaug->setId(id[i]);
01173 pdaug->addDaug(this);
01174 }
01175
01176 } //else
01177 } //makeDaughters
|
|
|
00771 {
00772
00773 //first add particle to the stdhep list;
00774 stdhep.createParticle(getP4Lab(),get4Pos(),-1,-1,
00775 EvtPDL::getStdHep(getId()));
00776
00777 int i;
00778 for(i=0;i<_ndaug;i++){
00779 stdhep.createParticle(_daug[i]->getP4Lab(),_daug[i]->get4Pos(),0,0,
00780 EvtPDL::getStdHep(_daug[i]->getId()));
00781 }
00782
00783 for(i=0;i<_ndaug;i++){
00784 _daug[i]->makeStdHepRec(1+i,1+i,stdhep);
00785 }
00786 return;
00787
00788 }
|
|
||||||||||||||||
|
Makes stdhep list 00736 {
00737
00738 //first add particle to the stdhep list;
00739 stdhep.createParticle(getP4Lab(),get4Pos(),-1,-1,
00740 EvtPDL::getStdHep(getId()));
00741
00742 int ii=0;
00743
00744 //lets see if this is a longlived particle and terminate the
00745 //list building!
00746
00747 while (list_of_stable[ii]!=EvtId(-1,-1)) {
00748 if (getId()==list_of_stable[ii]){
00749 secondary.createSecondary(0,this);
00750 return;
00751 }
00752 ii++;
00753 }
00754
00755
00756
00757
00758 int i;
00759 for(i=0;i<_ndaug;i++){
00760 stdhep.createParticle(_daug[i]->getP4Lab(),_daug[i]->get4Pos(),0,0,
00761 EvtPDL::getStdHep(_daug[i]->getId()));
00762 }
00763
00764 for(i=0;i<_ndaug;i++){
00765 _daug[i]->makeStdHepRec(1+i,1+i,stdhep,secondary,list_of_stable);
00766 }
00767 return;
00768
00769 }
|
|
|
Returns mass of particle. 00127 {
00128
00129 return _p.mass();
00130 }
|
|
|
Iterates over the particles in a decay chain. 00705 {
00706
00707 EvtParticle *bpart;
00708 EvtParticle *current;
00709
00710 current=this;
00711 int i;
00712
00713 if (_ndaug!=0) return _daug[0];
00714
00715 do{
00716 bpart=current->_parent;
00717 if (bpart==0) return 0;
00718 i=0;
00719 while (bpart->_daug[i]!=current) {i++;}
00720
00721 if ( bpart==rootOfTree ) {
00722 if ( i== bpart->_ndaug-1 ) return 0;
00723 }
00724
00725 i++;
00726 current=bpart;
00727
00728 }while(i>=bpart->_ndaug);
00729
00730 return bpart->_daug[i];
00731
00732 }
|
|
|
00347 { _genlifetime=0; }
|
|
|
|
|
|
Prints information for the particle. 00976 {
00977
00978 switch (EvtPDL::getSpinType(_id)){
00979 case EvtSpinType::SCALAR:
00980 report(INFO,"EvtGen") << "This is a scalar particle:"<<EvtPDL::name(_id).c_str()<<"\n";
00981 break;
00982 case EvtSpinType::VECTOR:
00983 report(INFO,"EvtGen") << "This is a vector particle:"<<EvtPDL::name(_id).c_str()<<"\n";
00984 break;
00985 case EvtSpinType::TENSOR:
00986 report(INFO,"EvtGen") << "This is a tensor particle:"<<EvtPDL::name(_id).c_str()<<"\n";
00987 break;
00988 case EvtSpinType::DIRAC:
00989 report(INFO,"EvtGen") << "This is a dirac particle:"<<EvtPDL::name(_id).c_str()<<"\n";
00990 break;
00991 case EvtSpinType::PHOTON:
00992 report(INFO,"EvtGen") << "This is a photon:"<<EvtPDL::name(_id).c_str()<<"\n";
00993 break;
00994 case EvtSpinType::NEUTRINO:
00995 report(INFO,"EvtGen") << "This is a neutrino:"<<EvtPDL::name(_id).c_str()<<"\n";
00996 break;
00997 case EvtSpinType::RARITASCHWINGER:
00998 report(INFO,"EvtGen") << "This is a raritaschwinger:"<<EvtPDL::name(_id).c_str()<<"\n";
00999 break;
01000 case EvtSpinType::STRING:
01001 report(INFO,"EvtGen") << "This is a string:"<<EvtPDL::name(_id).c_str()<<"\n";
01002 break;
01003 default:
01004 report(INFO,"EvtGen") <<"Unknown particle type in EvtParticle::printParticle()"<<endl;
01005 break;
01006 }
01007 report(INFO,"EvtGen") << "Number of daughters:"<<_ndaug<<"\n";
01008
01009
01010 }
|
|
|
Prints out the particle "tree" of a given particle. The tree consists of all daughters (and their daughters, etc) and their properties. 00873 {
00874
00875 report(INFO,"EvtGen") << "This is the current decay chain"<<endl;
00876 report(INFO,"") << "This top particle is "<<
00877 EvtPDL::name(_id).c_str()<<endl;
00878
00879 this->printTreeRec(0);
00880 report(INFO,"EvtGen") << "End of decay chain."<<endl;
00881
00882 }
|
|
|
00846 {
00847
00848 int newlevel,i;
00849 newlevel = level +1;
00850
00851
00852 if (_ndaug!=0) {
00853 if ( level > 0 ) {
00854 for (i=0;i<(5*level);i++) {
00855 report(INFO,"") <<" ";
00856 }
00857 }
00858 report(INFO,"") << EvtPDL::name(_id).c_str();
00859 report(INFO,"") << " -> ";
00860 for(i=0;i<_ndaug;i++){
00861 report(INFO,"") << EvtPDL::name(_daug[i]->getId()).c_str()<<" ";
00862 }
00863 for(i=0;i<_ndaug;i++){
00864 report(INFO,"") << _daug[i]->mass()<<" ";
00865 }
00866 report(INFO,"")<<endl;
00867 for(i=0;i<_ndaug;i++){
00868 _daug[i]->printTreeRec(newlevel);
00869 }
00870 }
00871 }
|
|
|
00077 {
00078 _first=1;
00079 }
|
|
|
00255 {_ndaug=0; return;}
|
|
||||||||||||||||
|
Implements EvtParticle. 00170 {
00171
00172 EvtDiracSpinor spplus;
00173 EvtDiracSpinor spminus;
00174
00175 if (EvtPDL::getStdHep(getId())>0){
00176 spplus.set(1.0,0.0,0.0,0.0);
00177 spminus.set(0.0,1.0,0.0,0.0);
00178 } else {
00179 spplus.set(0.0,0.0,1.0,0.0);
00180 spminus.set(0.0,0.0,0.0,1.0);
00181 }
00182
00183 EvtVector4C eplus(0.0,-1.0/sqrt(2.0),EvtComplex(0.0,-1.0/sqrt(2.0)),0.0);
00184 EvtVector4C ezero(0.0,0.0,0.0,1.0);
00185 EvtVector4C eminus(0.0,1.0/sqrt(2.0),EvtComplex(0.0,-1.0/sqrt(2.0)),0.0);
00186
00187 EvtRaritaSchwinger sppp=dirProd(eplus,spplus);
00188 EvtRaritaSchwinger spp=dirProd(sqrt(2.0/3.0)*ezero,spplus)+
00189 dirProd(sqrt(1.0/3.0)*eplus,spminus);
00190 EvtRaritaSchwinger spm=dirProd(sqrt(2.0/3.0)*ezero,spminus)+
00191 dirProd(sqrt(1.0/3.0)*eminus,spplus);
00192 EvtRaritaSchwinger spmm=dirProd(eminus,spminus);
00193 //
00194
00195 sppp.applyRotateEuler(alpha,beta,gamma);
00196 spp.applyRotateEuler(alpha,beta,gamma);
00197 spm.applyRotateEuler(alpha,beta,gamma);
00198 spmm.applyRotateEuler(alpha,beta,gamma);
00199
00200 EvtSpinDensity R;
00201 R.SetDim(4);
00202
00203 double sqmt2=sqrt(2.0*(this->getP4().mass()));
00204
00205 for ( int i=0; i<4; i++) {
00206 R.Set(0,i,(sppp*_spinorRest[i])/sqmt2);
00207 R.Set(1,i,(spp*_spinorRest[i])/sqmt2);
00208 R.Set(2,i,(spm*_spinorRest[i])/sqmt2);
00209 R.Set(3,i,(spmm*_spinorRest[i])/sqmt2);
00210 }
00211
00212 return R;
00213
00214 }
|
|
|
Returns a rotation matrix need to rotate the basis state to the helicity basis. The EvtSpinDensity matrix is just use as a matrix here. This function is to be implemented in each derived class. Implements EvtParticle. 00124 {
00125
00126 double sqmt2=sqrt(2.0*(this->getP4().mass()));
00127
00128 EvtDiracSpinor spplus;
00129 EvtDiracSpinor spminus;
00130
00131 if (EvtPDL::getStdHep(getId())>0){
00132 spplus.set(1.0,0.0,0.0,0.0);
00133 spminus.set(0.0,1.0,0.0,0.0);
00134 } else {
00135 spplus.set(0.0,0.0,1.0,0.0);
00136 spminus.set(0.0,0.0,0.0,1.0);
00137 }
00138
00139 EvtVector4C eplus(0.0,-1.0/sqrt(2.0),EvtComplex(0.0,-1.0/sqrt(2.0)),0.0);
00140 EvtVector4C ezero(0.0,0.0,0.0,1.0);
00141 EvtVector4C eminus(0.0,1.0/sqrt(2.0),EvtComplex(0.0,-1.0/sqrt(2.0)),0.0);
00142
00143 EvtRaritaSchwinger sppp=dirProd(eplus,spplus);
00144 EvtRaritaSchwinger spp=dirProd(sqrt(2.0/3.0)*ezero,spplus)+
00145 dirProd(sqrt(1.0/3.0)*eplus,spminus);
00146 EvtRaritaSchwinger spm=dirProd(sqrt(2.0/3.0)*ezero,spminus)+
00147 dirProd(sqrt(1.0/3.0)*eminus,spplus);
00148 EvtRaritaSchwinger spmm=dirProd(eminus,spminus);
00149
00150 //
00151
00152 EvtSpinDensity R;
00153 R.SetDim(4);
00154
00155
00156 for ( int i=0; i<4; i++) {
00157 R.Set(0,i,(sppp*_spinorRest[i])/sqmt2);
00158 R.Set(1,i,(spp*_spinorRest[i])/sqmt2);
00159 R.Set(2,i,(spm*_spinorRest[i])/sqmt2);
00160 R.Set(3,i,(spmm*_spinorRest[i])/sqmt2);
00161 }
00162
00163 return R;
00164
00165 }
|
|
|
Should only be used internally. 00081 {
00082 _channel=i;
00083 }
|
|
|
01180 {
01181
01182 if ( _decayProb == 0 ) _decayProb=new double;
01183 *_decayProb=prob;
01184 }
|
|
|
Set diagonal spindensity matrix. 00133 {
00134
00135 _rhoForward.SetDiag(getSpinStates());
00136 }
|
|
|
00074 {
00075 _first=0;
00076 }
|
|
|
set generator information; pingrg-2011-1-6 00141 {_generatorFlag = flag;}
|
|
|
00350 { _id=id;}
|
|
|
Generate lifetime according to pure exponential. 00093 {
00094 if (_genlifetime){
00095 _t=-log(EvtRandom::Flat())*EvtPDL::getctau(getId());
00096 }
00097 }
|
|
|
Set lifetime of the particle in parents restframe. 00089 {
00090 _t=tau;
00091 }
|
|
|
00357 { _p=EvtVector4R(m,0.0,0.0,0.0);}
|
|
|
00382 { _p =p4; }
|
|
||||||||||||||||||||
|
00381 { _p.set(e,px,py,pz); }
|
|
|
Sets the 4momentum in the parents restframe. 00244 {_p=p4;}
|
|
|
00384 {
00385 assert(_channel==-10||
00386 _id.getId()==particle_number.getId()||
00387 _id.getId()==-1);
00388 _id = particle_number;
00389 }
|
|
|
Set backward spin density matrix. 00337 {_rhoBackward=rho;}
|
|
|
Set forward spin density matrix. 00306 {_rhoForward=rho;}
|
|
||||||||||||||||||||
|
00188 {
00189
00190 EvtSpinDensity R=rotateToHelicityBasis(alpha,beta,gamma);
00191
00192 assert(R.GetDim()==rho.GetDim());
00193
00194 int n=rho.GetDim();
00195
00196 _rhoForward.SetDim(n);
00197
00198 int i,j,k,l;
00199
00200 for(i=0;i<n;i++){
00201 for(j=0;j<n;j++){
00202 EvtComplex tmp=0.0;
00203 for(k=0;k<n;k++){
00204 for(l=0;l<n;l++){
00205 tmp+=R.Get(l,i)*rho.Get(l,k)*conj(R.Get(k,j));
00206 }
00207 }
00208 _rhoForward.Set(i,j,tmp);
00209 }
00210 }
00211
00212 //report(INFO,"EvtGen") << "_rhoForward:"<<_rhoForward<<endl;
00213
00214 }
|
|
|
Set forward spin density matrix according to the density matrix rho in the helicity amplitude basis. 00157 {
00158
00159 EvtSpinDensity R=rotateToHelicityBasis();
00160
00161 assert(R.GetDim()==rho.GetDim());
00162
00163 int n=rho.GetDim();
00164
00165 _rhoForward.SetDim(n);
00166
00167 int i,j,k,l;
00168
00169 for(i=0;i<n;i++){
00170 for(j=0;j<n;j++){
00171 EvtComplex tmp=0.0;
00172 for(k=0;k<n;k++){
00173 for(l=0;l<n;l++){
00174 tmp+=R.Get(l,i)*rho.Get(l,k)*conj(R.Get(k,j));
00175 }
00176 }
00177 _rhoForward.Set(i,j,tmp);
00178 }
00179 }
00180
00181 //report(INFO,"EvtGen") << "_rhoForward:"<<_rhoForward<<endl;
00182
00183 }
|
|
|
Set spindensity matrix for e+e- -> V 00138 {
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 //Set helicity +1 and -1 to 1.
00150 rho.SetDiag(getSpinStates());
00151 rho.Set(1,1,EvtComplex(0.0,0.0));
00152 setSpinDensityForwardHelicityBasis(rho);
00153
00154 }
|
|
|
Returns Dirac spinor in the particles own restframe for a Dirac particle. Reimplemented in EvtDiracParticle. 00598 {
00599 EvtDiracSpinor tempD;
00600 int temp;
00601 temp = i;
00602 printParticle();
00603 report(ERROR,"EvtGen") << "and you have asked for the:"<<i
00604 <<"th dirac spinor."
00605 <<" I.e. you thought it was a"
00606 <<" Dirac particle!" << endl;
00607 ::abort();
00608 return tempD;
00609 }
|
|
|
Returns Dirac spinor in the particles own restframe for a Neutrino particle. Reimplemented in EvtNeutrinoParticle. 00622 {
00623 EvtDiracSpinor tempD;
00624 printParticle();
00625 report(ERROR,"EvtGen") << "and you have asked for the "
00626 <<"dirac spinor."
00627 <<" I.e. you thought it was a"
00628 <<" neutrino particle!" << endl;
00629 ::abort();
00630 return tempD;
00631 }
|
|
|
Returns Dirac spinor in the parents restframe for a Dirac particle. Reimplemented in EvtDiracParticle. 00585 {
00586 EvtDiracSpinor tempD;
00587 int temp;
00588 temp = i;
00589 printParticle();
00590 report(ERROR,"EvtGen") << "and you have asked for the:"<<i
00591 <<"th dirac spinor."
00592 <<" I.e. you thought it was a"
00593 <<" Dirac particle!" << endl;
00594 ::abort();
00595 return tempD;
00596 }
|
|
|
Returns Dirac spinor in the parents restframe for a Neutrino particle. Reimplemented in EvtNeutrinoParticle. 00611 {
00612 EvtDiracSpinor tempD;
00613 printParticle();
00614 report(ERROR,"EvtGen") << "and you have asked for the "
00615 <<"dirac spinor."
00616 <<" I.e. you thought it was a"
00617 <<" neutrino particle!" << endl;
00618 ::abort();
00619 return tempD;
00620 }
|
|
|
00116 {
00117
00118 return _spinorRest[i];
00119
00120 }
|
|
|
00110 {
00111
00112 return _spinor[i];
00113
00114 }
|
|
|
00966 {
00967
00968 std::string retval=EvtPDL::name(_id);
00969 retval+=" -> ";
00970
00971 retval+=treeStrRec(0);
00972
00973 return retval;
00974 }
|
|
|
00884 {
00885
00886 int newlevel,i;
00887 newlevel = level +1;
00888
00889 std::string retval="";
00890
00891 for(i=0;i<_ndaug;i++){
00892 retval+=EvtPDL::name(_daug[i]->getId());
00893 if ( _daug[i]->getNDaug() > 0 ) {
00894 retval+= " (";
00895 retval+= _daug[i]->treeStrRec(newlevel);
00896 retval+= ") ";
00897 }
00898 else{
00899 if ( i!=_ndaug-1) retval+=" ";
00900 }
00901 }
00902
00903 return retval;
00904 }
|
|
|
00906 { //pingrg, Jun. 6, 2008
00907 std::string retval="";
00908
00909 if (resonance == EvtPDL::name(_id).c_str() && _ndaug!=0) {
00910 retval=resonance+": "+IntToStr(_ndaug)+" ";
00911 for(int i=0;i<_ndaug;i++){
00912 retval += EvtPDL::name(_daug[i]->getId()).c_str();
00913 retval += " ";
00914 }
00915 }
00916
00917 for(int i=0;i<_ndaug;i++){
00918 _daug[i]->writeTreeRec(resonance);
00919 }
00920 std::cout<<retval;
00921 return retval;
00922 }
|
|
|
|
|
|
|
|
|
|
1.3.9.1