EvtPythia Class Reference

#include <EvtPythia.hh>

Inheritance diagram for EvtPythia:

List of all members.

Public Member Functions
void	checkNArg (int a1, int a2=-1, int a3=-1, int a4=-1)
void	checkNDaug (int d1, int d2=-1)
void	checkQ ()
void	checkSpinDaughter (int d1, EvtSpinType::spintype sp)
void	checkSpinParent (EvtSpinType::spintype sp)
EvtDecayBase *	clone ()
void	command (std::string cmd)
std::string	commandName ()
void	decay (EvtParticle *p)
void	disableCheckQ ()
	EvtPythia ()
double	getArg (int j)
double *	getArgs ()
std::string *	getArgsStr ()
std::string	getArgStr (int j)
double	getBranchingFraction ()
EvtId	getDaug (int i)
EvtId *	getDaugs ()
int	getDSum ()
std::string	getModelName ()
void	getName (std::string &name)
int	getNArg ()
int	getNDaug ()
EvtId	getParentId ()
int	getPHOTOS ()
double	getProbMax (double prob)
void	init ()
void	initProbMax ()
int	isDaughterSpinDensitySet (int daughter)
void	makeDecay (EvtParticle *p)
virtual bool	matchingDecay (const EvtDecayBase &other) const
void	noProbMax ()
virtual int	nRealDaughters ()
void	printSummary ()
double	resetProbMax (double prob)
void	saveDecayInfo (EvtId ipar, int ndaug, EvtId *daug, int narg, std::vector< std::string > &args, std::string name, double brfr)
void	setDaughterSpinDensity (int daughter)
void	setPHOTOS ()
void	setProbMax (double prbmx)
void	setSummary ()
void	setVerbose ()
int	summary ()
int	verbose ()
virtual	~EvtPythia ()
Static Public Member Functions
void	findMass (EvtParticle *p)
void	findMasses (EvtParticle *p, int ndaugs, EvtId daugs[10], double masses[10])
double	findMaxMass (EvtParticle *p)
void	pythiacont (double *, int *, int *, double *, double *, double *, double *)
void	pythiaInit (int f)
Protected Member Functions
bool	daugsDecayedByParentModel ()
Protected Attributes
bool	_daugsDecayedByParentModel
Private Member Functions
void	fixPolarizations (EvtParticle *p)
void	store (EvtDecayBase *jsdecay)
Static Private Member Functions
bool	diquark (int)
void	MakePythiaFile (char *fname)
double	NominalMass (int)
void	WritePythiaEntryHeader (std::ofstream &outdec, int lundkc, EvtId evtnum, std::string name, int chg, int cchg, int spin2, double mass, double width, double maxwidth, double ctau, int stable, double rawbrfrsum)
void	WritePythiaParticle (std::ofstream &outdec, EvtId ipar, EvtId iparname, int &first)
Static Private Attributes
std::string *	commands = 0
EvtDecayBasePtr *	jetsetdecays = 0
int	lcommand = 0
int	ncommand = 0
int	njetsetdecays = 0
int	ntable = 0

---

Constructor & Destructor Documentation

EvtPythia::EvtPythia (   )

{}

EvtPythia::~EvtPythia (   )  [virtual]

{ int i; //the deletion of commands is really uggly! if (njetsetdecays==0) { delete [] commands; commands=0; return; } for(i=0;i<njetsetdecays;i++){ if (jetsetdecays[i]==this){ jetsetdecays[i]=jetsetdecays[njetsetdecays-1]; njetsetdecays--; if (njetsetdecays==0) { delete [] commands; commands=0; } return; } } report(ERROR,"EvtGen") << "Error in destroying Pythia model!"<<endl; }

---

Member Function Documentation

void EvtDecayBase::checkNArg (  int  a1, int  a2 = -1, int  a3 = -1, int  a4 = -1 )  [inherited]

{ if ( _narg != a1 && _narg != a2 && _narg != a3 && _narg != a4 ) { report(ERROR,"EvtGen") << _modelname.c_str() << " generator expected "<<endl; report(ERROR,"EvtGen") << a1<<endl;; if ( a2>-1) { report(ERROR,"EvtGen") << " or " << a2<<endl; } if ( a3>-1) { report(ERROR,"EvtGen") << " or " << a3<<endl; } if ( a4>-1) { report(ERROR,"EvtGen") << " or " << a4<<endl; } report(ERROR,"EvtGen") << " arguments but found:"<< _narg << endl; printSummary(); report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; ::abort(); } }

void EvtDecayBase::checkNDaug (  int  d1, int  d2 = -1 )  [inherited]

{ if ( _ndaug != d1 && _ndaug != d2 ) { report(ERROR,"EvtGen") << _modelname.c_str() << " generator expected "; report(ERROR,"EvtGen") << d1; if ( d2>-1) { report(ERROR,"EvtGen") << " or " << d2; } report(ERROR,"EvtGen") << " daughters but found:"<< _ndaug << endl; printSummary(); report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; ::abort(); } }

void EvtDecayBase::checkQ (   )  [inherited]

{ int i; int q=0; int qpar; //If there are no daughters (jetset etc) then we do not //want to do this test. Why? Because sometimes the parent //will have a nonzero charge. if ( _ndaug != 0) { for(i=0; i<_ndaug; i++ ) { q += EvtPDL::chg3(_daug[i]); } qpar = EvtPDL::chg3(_parent); if ( q != qpar ) { report(ERROR,"EvtGen") <<_modelname.c_str()<< " generator expected " << " charge to be conserved, found:"<<endl; report(ERROR,"EvtGen") << "Parent charge of "<<(qpar/3)<<endl; report(ERROR,"EvtGen") << "Sum of daughter charge of "<<(q/3)<<endl; report(ERROR,"EvtGen") << "The parent is "<< EvtPDL::name(_parent).c_str()<<endl; for(i=0; i<_ndaug; i++ ) { report(ERROR,"EvtGen") << "Daughter "<< EvtPDL::name(_daug[i]).c_str()<<endl; } report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; ::abort(); } } }

void EvtDecayBase::checkSpinDaughter (  int  d1, EvtSpinType::spintype  sp )  [inherited]

{ EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getDaug(d1)); if ( parenttype != sp ) { report(ERROR,"EvtGen") << _modelname.c_str() << " did not get the correct daughter spin d=" << d1 << endl; printSummary(); report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; ::abort(); } }

void EvtDecayBase::checkSpinParent (  EvtSpinType::spintype  sp  )  [inherited]

{ EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getParentId()); if ( parenttype != sp ) { report(ERROR,"EvtGen") << _modelname.c_str() << " did not get the correct parent spin\n"; printSummary(); report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; ::abort(); } }

EvtDecayBase * EvtPythia::clone (   )  [virtual]

Implements EvtDecayBase. { return new EvtPythia; }

void EvtPythia::command (  std::string  cmd  )  [virtual]

Reimplemented from EvtDecayBase. { if (ncommand==lcommand){ lcommand=10+2*lcommand; std::string* newcommands=new std::string[lcommand]; int i; for(i=0;i<ncommand;i++){ newcommands[i]=commands[i]; } delete [] commands; commands=newcommands; } commands[ncommand]=cmd; ncommand++; }

std::string EvtPythia::commandName (   )  [virtual]

Reimplemented from EvtDecayBase. { return std::string("JetSetPar"); }

bool EvtDecayBase::daugsDecayedByParentModel (   )  [inline, protected, inherited]

{return _daugsDecayedByParentModel;}

void EvtPythia::decay (  EvtParticle *  p  )  [virtual]

Implements EvtDecayBase. { //added by Lange Jan4,2000 static EvtId STRNG=EvtPDL::getId("string"); int istdheppar=EvtPDL::getStdHep(p->getId()); if (pycomp_(&istdheppar)==0){ report(ERROR,"EvtGen") << "Pythia can not decay:" <<EvtPDL::name(p->getId()).c_str()<<endl; // return; //comment by Ping RG, due not to decay eta_c(2S) } double mp=p->mass(); EvtVector4R p4[20]; int i,more; int ip=EvtPDL::getStdHep(p->getId()); int ndaugjs; int kf[100]; EvtId evtnumstable[100],evtnumparton[100]; int stableindex[100],partonindex[100]; int numstable; int numparton; int km[100]; EvtId type[MAX_DAUG]; pythiaInit(0); double px[100],py[100],pz[100],e[100]; if ( p->getNDaug() != 0 ) { p->deleteDaughters(true);} int count=0; do{ pythiadec_(&ip,&mp,&ndaugjs,kf,km,px,py,pz,e); numstable=0; numparton=0; for(i=0;i<ndaugjs;i++){ if (EvtPDL::evtIdFromStdHep(kf[i])==EvtId(-1,-1)) { report(ERROR,"EvtGen") << "Pythia returned particle:"<<kf[i]<<endl; report(ERROR,"EvtGen") << "This can not be translated to evt number"<<endl; report(ERROR,"EvtGen") << "and the decay will be rejected!"<<endl; report(ERROR,"EvtGen") << "The decay was of particle:"<<ip<<endl; int i=1; pylist_(i); } //sort out the partons if (abs(kf[i])<=6||kf[i]==21){ partonindex[numparton]=i; evtnumparton[numparton]=EvtPDL::evtIdFromStdHep(kf[i]); numparton++; } else{ stableindex[numstable]=i; evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(kf[i]); numstable++; } // have to protect against negative mass^2 for massless particles // i.e. neutrinos and photons. // this is uggly but I need to fix it right now.... if (px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i]>=e[i]*e[i]){ e[i]=sqrt(px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i])+0.0000000000001; } p4[i].set(e[i],px[i],py[i],pz[i]); } int channel=EvtDecayTable::inChannelList(p->getId(),numstable,evtnumstable); more=(channel!=-1); count++; }while( more && (count<10000) ); if (count>9999) { report(INFO,"EvtGen") << "Too many loops in EvtPythia!!!"<<endl; report(INFO,"EvtGen") << "Parent:"<<EvtPDL::name(getParentId()).c_str()<<endl; for(i=0;i<numstable;i++){ report(INFO,"EvtGen") << "Daug("<<i<<")"<<EvtPDL::name(evtnumstable[i]).c_str()<<endl; } } if (numparton==0){ p->makeDaughters(numstable,evtnumstable); for(i=0;i<numstable;i++){ p->getDaug(i)->init(evtnumstable[i],p4[stableindex[i]]); } fixPolarizations(p); return ; } else{ //have partons in JETSET EvtVector4R p4string(0.0,0.0,0.0,0.0); for(i=0;i<numparton;i++){ p4string+=p4[partonindex[i]]; } int nprimary=1; type[0]=STRNG; for(i=0;i<numstable;i++){ if (km[stableindex[i]]==0){ type[nprimary++]=evtnumstable[i]; } } p->makeDaughters(nprimary,type); p->getDaug(0)->init(STRNG,p4string); EvtVector4R p4partons[10]; for(i=0;i<numparton;i++){ p4partons[i]=p4[partonindex[i]]; } ((EvtStringParticle*)p->getDaug(0))->initPartons(numparton,p4partons,evtnumparton); nprimary=1; for(i=0;i<numstable;i++){ if (km[stableindex[i]]==0){ p->getDaug(nprimary++)->init(evtnumstable[i],p4[stableindex[i]]); } } int nsecond=0; for(i=0;i<numstable;i++){ if (km[stableindex[i]]!=0){ type[nsecond++]=evtnumstable[i]; } } p->getDaug(0)->makeDaughters(nsecond,type); nsecond=0; for(i=0;i<numstable;i++){ if (km[stableindex[i]]!=0){ p4[stableindex[i]]=boostTo(p4[stableindex[i]],p4string); p->getDaug(0)->getDaug(nsecond)->init(evtnumstable[i],p4[stableindex[i]]); p->getDaug(0)->getDaug(nsecond)->setDiagonalSpinDensity(); p->getDaug(0)->getDaug(nsecond)->decay(); nsecond++; } } fixPolarizations(p); return ; } }

bool EvtPythia::diquark (  int   )  [static, private]

{ switch(ID) { case 1103: case 2101: case 2103: case 2203: case 3101: case 3103: case 3201: case 3203: case 3303: case 4101: case 4103: case 4201: case 4203: case 4301: case 4303: case 4403: case 5101: case 5103: case 5201: case 5203: case 5301: case 5303: case 5401: case 5403: case 5503: return true; break; default: return false; break; } }

void EvtDecayBase::disableCheckQ (   )  [inline, inherited]

{_chkCharge=0;};

void EvtDecayBase::findMass (  EvtParticle *  p  )  [static, inherited]

{ //Need to also check that this mass does not screw //up the parent //This code assumes that for the ith daughter, 0..i-1 //already have a mass double maxOkMass=findMaxMass(p); int count=0; double mass; bool massOk=false; int i; while (!massOk) { count++; if ( count > 10000 ) { report(INFO,"EvtGen") << "Can not find a valid mass for: " << EvtPDL::name(p->getId()).c_str() <<endl; report(INFO,"EvtGen") << "Now printing parent and/or grandparent tree\n"; if ( p->getParent() ) { if ( p->getParent()->getParent() ) { p->getParent()->getParent()->printTree(); report(INFO,"EvtGen") << p->getParent()->getParent()->mass() <<endl; report(INFO,"EvtGen") << p->getParent()->mass() <<endl; } else{ p->getParent()->printTree(); report(INFO,"EvtGen") << p->getParent()->mass() <<endl; } } else p->printTree(); report(INFO,"EvtGen") << "maxokmass=" << maxOkMass << " " << EvtPDL::getMinMass(p->getId()) << " " << EvtPDL::getMaxMass(p->getId())<<endl; if ( p->getNDaug() ) { for (i=0; i<p->getNDaug(); i++) { report(INFO,"EvtGen") << p->getDaug(i)->mass()<<" "; } report(INFO,"EvtGen") << endl; } if ( maxOkMass >= EvtPDL::getMinMass(p->getId()) ) { report(INFO,"EvtGen") << "taking a default value\n"; p->setMass(maxOkMass); return; } assert(0); } mass = EvtPDL::getMass(p->getId()); //Just need to check that this mass is > than //the mass of all daughters double massSum=0.; if ( p->getNDaug() ) { for (i=0; i<p->getNDaug(); i++) { massSum+= p->getDaug(i)->mass(); } } //some special cases are handled with 0 (stable) or 1 (k0->ks/kl) daughters if (p->getNDaug()<2) massOk=true; if ( p->getParent() ) { if ( p->getParent()->getNDaug()==1 ) massOk=true; } if ( !massOk ) { if (massSum < mass) massOk=true; if ( mass> maxOkMass) massOk=false; } } p->setMass(mass); }

void EvtDecayBase::findMasses (  EvtParticle *  p, int  ndaugs, EvtId  daugs[10], double  masses[10] )  [static, inherited]

{ int i; double mass_sum; int count=0; if (!( p->firstornot() )) { for (i = 0; i < ndaugs; i++ ) { masses[i] = p->getDaug(i)->mass(); } //for } //if else { p->setFirstOrNot(); // if only one daughter do it if (ndaugs==1) { masses[0]=p->mass(); return; } //until we get a combo whose masses are less than _parent mass. do { mass_sum = 0.0; for (i = 0; i < ndaugs; i++ ) { masses[i] = EvtPDL::getMass(daugs[i]); mass_sum = mass_sum + masses[i]; } count++; if(count==10000) { report(ERROR,"EvtGen") <<"Decaying particle:"<< EvtPDL::name(p->getId()).c_str()<<" (m="<<p->mass()<<")"<<endl; report(ERROR,"EvtGen") <<"To the following daugthers"<<endl; for (i = 0; i < ndaugs; i++ ) { report(ERROR,"EvtGen") << EvtPDL::name(daugs[i]).c_str() << endl; } report(ERROR,"EvtGen") << "Has been rejected "<<count << " times, will now take minimal masses " << " of daugthers"<<endl; mass_sum=0.; for (i = 0; i < ndaugs; i++ ) { masses[i] = EvtPDL::getMinMass(daugs[i]); mass_sum = mass_sum + masses[i]; } if (mass_sum > p->mass()){ report(ERROR,"EvtGen") << "Parent mass="<<p->mass() << "to light for daugthers."<<endl << "Will throw the event away."<<endl; //dont terminate - start over on the event. EvtStatus::setRejectFlag(); mass_sum=0.; // ::abort(); } } } while ( mass_sum > p->mass()); } //else return; }

double EvtDecayBase::findMaxMass (  EvtParticle *  p  )  [static, inherited]

{ double maxOkMass=EvtPDL::getMaxMass(p->getId()); //protect against vphotons if ( maxOkMass < 0.0000000001 ) return 10000000.; //and against already determined masses if ( p->hasValidP4() ) maxOkMass=p->mass(); EvtParticle *par=p->getParent(); if ( par ) { double maxParMass=findMaxMass(par); int i; double minDaugMass=0.; for(i=0;i<par->getNDaug();i++){ EvtParticle *dau=par->getDaug(i); if ( dau!=p) { // it might already have a mass if ( dau->isInitialized() || dau->hasValidP4() ) minDaugMass+=dau->mass(); else //give it a bit of phase space minDaugMass+=1.000001*EvtPDL::getMinMass(dau->getId()); } } if ( maxOkMass>(maxParMass-minDaugMass)) maxOkMass=maxParMass-minDaugMass; } return maxOkMass; }

void EvtPythia::fixPolarizations (  EvtParticle *  p  )  [private]

{ //special case for now to handle the J/psi polarization int ndaug=p->getNDaug(); int i; static EvtId Jpsi=EvtPDL::getId("J/psi"); for(i=0;i<ndaug;i++){ if(p->getDaug(i)->getId()==Jpsi){ EvtSpinDensity rho; rho.SetDim(3); rho.Set(0,0,0.5); rho.Set(0,1,0.0); rho.Set(0,2,0.0); rho.Set(1,0,0.0); rho.Set(1,1,1.0); rho.Set(1,2,0.0); rho.Set(2,0,0.0); rho.Set(2,1,0.0); rho.Set(2,2,0.5); EvtVector4R p4Psi=p->getDaug(i)->getP4(); double alpha=atan2(p4Psi.get(2),p4Psi.get(1)); double beta=acos(p4Psi.get(3)/p4Psi.d3mag()); p->getDaug(i)->setSpinDensityForwardHelicityBasis(rho,alpha,beta,0.0); setDaughterSpinDensity(i); } } }

double EvtDecayBase::getArg (  int  j  )  [inherited]

{ // Verify string const char* str = _args[j].c_str(); int i = 0; while(str[i]!=0){ if (isalpha(str[i]) && str[i]!='e') { report(INFO,"EvtGen") << "String " << str << " is not a number" << endl; assert(0); } i++; } char** tc=0; return strtod(_args[j].c_str(),tc); }

double * EvtDecayBase::getArgs (   )  [inherited]

{ if ( _argsD ) return _argsD; //The user has asked for a list of doubles - the arguments //better all be doubles... if ( _narg==0 ) return _argsD; _argsD = new double[_narg]; int i; char * tc; for(i=0;i<_narg;i++) { _argsD[i] = strtod(_args[i].c_str(),&tc); } return _argsD; }

std::string* EvtDecayBase::getArgsStr (   )  [inline, inherited]

{return _args;}

std::string EvtDecayBase::getArgStr (  int  j  )  [inline, inherited]

{return _args[j];}

double EvtDecayBase::getBranchingFraction (   )  [inline, inherited]

{return _brfr;}

EvtId EvtDecayBase::getDaug (  int  i  )  [inline, inherited]

{return _daug[i];}

EvtId* EvtDecayBase::getDaugs (   )  [inline, inherited]

{return _daug;}

int EvtDecayBase::getDSum (   )  [inline, inherited]

{return _dsum; }

std::string EvtDecayBase::getModelName (   )  [inline, inherited]

{return _modelname; }

void EvtPythia::getName (  std::string &  name  )  [virtual]

Implements EvtDecayBase. { model_name="PYTHIA"; }

int EvtDecayBase::getNArg (   )  [inline, inherited]

{return _narg;}

int EvtDecayBase::getNDaug (   )  [inline, inherited]

{return _ndaug;}

EvtId EvtDecayBase::getParentId (   )  [inline, inherited]

{return _parent;}

int EvtDecayBase::getPHOTOS (   )  [inline, inherited]

{return _photos;}

double EvtDecayBase::getProbMax (  double  prob  )  [inherited]

{ int i; //diagnostics sum_prob+=prob; if (prob>max_prob) max_prob=prob; if ( defaultprobmax && ntimes_prob<=500 ) { //We are building up probmax with this iteration ntimes_prob += 1; if ( prob > probmax ) { probmax = prob;} if (ntimes_prob==500) { probmax*=1.2; } return 1000000.0*prob; } if ( prob> probmax*1.0001) { report(INFO,"EvtGen") << "prob > probmax:("<<prob<<">"<<probmax<<")"; report(INFO,"") << "("<<_modelname.c_str()<<") "; report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> "; for(i=0;i<_ndaug;i++){ report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " "; } report(INFO,"") << endl; if (defaultprobmax) probmax = prob; } ntimes_prob += 1; return probmax; } //getProbMax

void EvtPythia::init (   )  [virtual]

Reimplemented from EvtDecayBase. { checkNArg(1); if (getParentId().isAlias()){ report(ERROR,"EvtGen") << "EvtPythia finds that you are decaying the"<<endl << " aliased particle " << EvtPDL::name(getParentId()).c_str() << " with the Pythia model"<<endl << " this does not work, please modify decay table." << endl; report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; ::abort(); } store(this); }

void EvtPythia::initProbMax (   )  [virtual]

Reimplemented from EvtDecayBase. { noProbMax(); }

int EvtDecayIncoherent::isDaughterSpinDensitySet (  int  daughter  )  [inline, inherited]

{return spinDensitySet[daughter];}

void EvtDecayIncoherent::makeDecay (  EvtParticle *  p  )  [virtual, inherited]

Implements EvtDecayBase. { int i; //initialize this the hard way.. //Lange June 26, 2000 for (i=0; i<MAX_DAUG; i++ ) { spinDensitySet[i]=0;} _daugsDecayedByParentModel=false; decay(p); p->setDecayProb(1.); EvtSpinDensity rho; rho.SetDiag(p->getSpinStates()); p->setSpinDensityBackward(rho); if (getPHOTOS() || EvtRadCorr::alwaysRadCorr()) { EvtRadCorr::doRadCorr(p); } //Now decay the daughters. if ( !daugsDecayedByParentModel()) { for(i=0;i<p->getNDaug();i++){ //Need to set the spin density of the daughters to be //diagonal. rho.SetDiag(p->getDaug(i)->getSpinStates()); //if (p->getDaug(i)->getNDaug()==0){ //only do this if the user has not already set the //spin density matrix herself. //Lange June 26, 2000 if ( isDaughterSpinDensitySet(i)==0 ) { p->getDaug(i)->setSpinDensityForward(rho); } else{ //report(INFO,"EvtGen") << "spinDensitymatrix already set!!!\n"; EvtSpinDensity temp=p->getDaug(i)->getSpinDensityForward(); // report(INFO,"EvtGen") <<temp<<endl; } //Now decay the daughter. Really! p->getDaug(i)->decay(); //} } } }

void EvtPythia::MakePythiaFile (  char *  fname  )  [static, private]

{ EvtId ipar; int lundkc; //int part_list[MAX_PART]; ofstream outdec; outdec.open(fname); //outdec << "ERROR;"<<endl; //outdec << ";"<<endl; //outdec << ";This decayfile has been automatically created by"<<endl; //outdec << ";EvtGen from the DECAY.DEC file"<<endl; //outdec << ";"<<endl; int nokcentry; for(lundkc=1;lundkc<500;lundkc++){ nokcentry=1; int iipar; for(iipar=0;iipar<EvtPDL::entries();iipar++){ ipar=EvtId(iipar,iipar); //no aliased particles! std::string tempStr = EvtPDL::name(ipar); EvtId realId = EvtPDL::getId(tempStr); if ( realId.isAlias() != 0 ) continue; if(!( EvtPDL::getStdHep(ipar)==21 || EvtPDL::getStdHep(ipar)==22 || EvtPDL::getStdHep(ipar)==23)) { if (lundkc==EvtPDL::getLundKC(ipar)){ nokcentry=0; int first=1; WritePythiaParticle(outdec,ipar,ipar,first); EvtId ipar2=EvtPDL::chargeConj(ipar); if (ipar2!=ipar){ WritePythiaParticle(outdec,ipar2,ipar,first); } if (first){ WritePythiaEntryHeader(outdec, //EvtPDL::getLundKC(ipar), EvtPDL::getStdHep(ipar), ipar, EvtPDL::name(ipar), EvtPDL::chg3(ipar), 0,0,EvtPDL::getMeanMass(ipar), EvtPDL::getWidth(ipar), EvtPDL::getMeanMass(ipar)-EvtPDL::getMinMass(ipar), EvtPDL::getctau(ipar),0,0.0); } } } } if (lundkc==99999) // Write out diquarks after quarks, but only once for(iipar=0;iipar<EvtPDL::entries();iipar++){ ipar=EvtId(iipar,iipar); if (diquark(EvtPDL::getStdHep(ipar))){ nokcentry=0; int first=1; WritePythiaParticle(outdec,ipar,ipar,first); EvtId ipar2=EvtPDL::chargeConj(ipar); if (ipar2!=ipar){ WritePythiaParticle(outdec,ipar2,ipar,first); } if (first){ WritePythiaEntryHeader(outdec, EvtPDL::getStdHep(ipar), ipar, EvtPDL::name(ipar), NominalCharge(EvtPDL::getStdHep(ipar)), -1,0, NominalMass(EvtPDL::getStdHep(ipar)), 0, 0, 0, 0,0.0); } } } /* if (nokcentry){ WritePythiaEntryHeader(outdec, lundkc,EvtId(-1,-1)," ", 0,0,0,EvtPDL::getNominalMass(ipar),0.0,0.0, EvtPDL::getctau(ipar),0,0.0); } */ } outdec.close(); }

bool EvtDecayBase::matchingDecay (  const EvtDecayBase &  other  )  const [virtual, inherited]

{ if ( _ndaug != other._ndaug) return false; if ( _parent != other._parent) return false; std::vector<int> useDs; for ( unsigned int i=0; i<_ndaug; i++) useDs.push_back(0); for ( unsigned int i=0; i<_ndaug; i++) { bool foundIt=false; for ( unsigned int j=0; j<_ndaug; j++) { if ( useDs[j] == 1 ) continue; if ( _daug[i] == other._daug[j] && _daug[i].getAlias() == other._daug[j].getAlias()) { foundIt=true; useDs[j]=1; break; } } if ( foundIt==false) return false; } for ( unsigned int i=0; i<_ndaug; i++) if ( useDs[i]==0) return false; return true; }

double EvtPythia::NominalMass (  int   )  [static, private]

{ // return default mass in PYTHIA switch(ID) { case 1103: return 0.77133; case 2101: return 0.57933; case 2103: return 0.77133; case 2203: return 0.77133; case 3101: return 0.80473; case 3103: return 0.92953; case 3201: return 0.80473; case 3203: return 0.92953; case 3303: return 1.09361; case 4101: return 1.96908; case 4103: return 2.00808; case 4201: return 1.96908; case 4203: return 2.00808; case 4301: return 2.15432; case 4303: return 2.17967; case 4403: return 3.27531; case 5101: return 5.38897; case 5103: return 5.40145; case 5201: return 5.38897; case 5203: return 5.40145; case 5301: return 5.56725; case 5303: return 5.57536; case 5401: return 6.67143; case 5403: return 6.67397; case 5503: return 10.07354; break; default: return 0.0; break; } }

void EvtDecayBase::noProbMax (   )  [inherited]

{ defaultprobmax=0; }

virtual int EvtDecayBase::nRealDaughters (   )  [inline, virtual, inherited]

Reimplemented in EvtBtoKD3P, and EvtVSSBMixCPT. { return _ndaug;}

void EvtDecayBase::printSummary (   )  [inherited]

{ int i; if (ntimes_prob>0) { report(INFO,"EvtGen") << "Calls="<<ntimes_prob<<" eff:"<< sum_prob/(probmax*ntimes_prob)<<" frac. max:"<<max_prob/probmax; report(INFO,"") <<" probmax:"<<probmax<<" max:"<<max_prob<<" : "; } report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> "; for(i=0;i<_ndaug;i++){ report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " "; } report(INFO,"") << " ("<<_modelname.c_str()<<"):"<< endl; }

void EvtPythia::pythiacont (  double *  , int *  , int *  , double *  , double *  , double *  , double *  )  [static]

{ pycontinuum_(energy,ndaugjs,kf,px,py,pz,e); }

void EvtPythia::pythiaInit (  int  f  )  [static]

{ static int first=1; if (first){ first=0; report(INFO,"EvtGen") << "Will initialize Pythia."<<endl; for(int i=0;i<ncommand;i++) pygive_(commands[i].c_str(),strlen(commands[i].c_str())); char fname[200]; char hostBuffer[100]; if ( gethostname( hostBuffer, 100 ) != 0 ){ report(ERROR,"EvtGen") << " couldn't get hostname." << endl; strncpy( hostBuffer, "hostnameNotFound", 100 ); } char pid[100]; int thePid=getpid(); if ( sprintf( pid, "%d", thePid ) == 0 ){ report(ERROR,"EvtGen") << " couldn't get process ID." << endl; strncpy( pid, "666", 100 ); } strcpy(fname,"jet.d-"); strcat(fname,hostBuffer); strcat(fname,"-"); strcat(fname,pid); MakePythiaFile(fname); evtpythiainit_(fname,strlen(fname)); initpythia_(&dummy); if (0==getenv("EVTSAVEJETD")){ char delcmd[300]; strcpy(delcmd,"rm -f "); strcat(delcmd,fname); system(delcmd); } report(INFO,"EvtGen") << "Done initializing Pythia."<<endl; } }

double EvtDecayBase::resetProbMax (  double  prob  )  [inherited]

{ report(INFO,"EvtGen") << "Reseting prob max\n"; report(INFO,"EvtGen") << "prob > probmax:("<<prob<<">"<<probmax<<")"; report(INFO,"") << "("<<_modelname.c_str()<<")"; report(INFO,"") << EvtPDL::getStdHep(_parent)<<"->"; for( int i=0;i<_ndaug;i++){ report(INFO,"") << EvtPDL::getStdHep(_daug[i]) << " "; } report(INFO,"") << endl; probmax = 0.0; defaultprobmax = 0; ntimes_prob = 0; return prob; }

void EvtDecayBase::saveDecayInfo (  EvtId  ipar, int  ndaug, EvtId *  daug, int  narg, std::vector< std::string > &  args, std::string  name, double  brfr )  [inherited]

{ int i; _brfr=brfr; _ndaug=ndaug; _narg=narg; _parent=ipar; _dsum=0; if (_ndaug>0) { _daug=new EvtId [_ndaug]; for(i=0;i<_ndaug;i++){ _daug[i]=daug[i]; _dsum+=daug[i].getAlias(); } } else{ _daug=0; } if (_narg>0) { _args=new std::string[_narg+1]; for(i=0;i<_narg;i++){ _args[i]=args[i]; } } else{ _args = 0; } _modelname=name; this->init(); this->initProbMax(); if (_chkCharge){ this->checkQ(); } if (defaultprobmax){ report(INFO,"EvtGen") << "No default probmax for "; report(INFO,"") << "("<<_modelname.c_str()<<") "; report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> "; for(i=0;i<_ndaug;i++){ report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " "; } report(INFO,"") << endl; report(INFO,"") << "This is fine for development, but must be provided for production."<<endl; report(INFO,"EvtGen") << "Never fear though - the decay will use the \n"; report(INFO,"EvtGen") << "500 iterations to build up a good probmax \n"; report(INFO,"EvtGen") << "before accepting a decay. "<<endl; } }

void EvtDecayIncoherent::setDaughterSpinDensity (  int  daughter  )  [inline, inherited]

{ spinDensitySet[daughter]=1; return;}

void EvtDecayBase::setPHOTOS (   )  [inline, inherited]

{_photos=1;}

void EvtDecayBase::setProbMax (  double  prbmx  )  [inherited]

{ defaultprobmax=0; probmax=prbmx; }

void EvtDecayBase::setSummary (   )  [inline, inherited]

{_summary=1;}

void EvtDecayBase::setVerbose (   )  [inline, inherited]

{_verbose=1;}

void EvtPythia::store (  EvtDecayBase *  jsdecay  )  [private]

{ if (njetsetdecays==ntable){ EvtDecayBasePtr* newjetsetdecays=new EvtDecayBasePtr[2*ntable+10]; int i; for(i=0;i<ntable;i++){ newjetsetdecays[i]=jetsetdecays[i]; } ntable=2*ntable+10; delete [] jetsetdecays; jetsetdecays=newjetsetdecays; } jetsetdecays[njetsetdecays++]=jsdecay; }

int EvtDecayBase::summary (   )  [inline, inherited]

{return _summary; }

int EvtDecayBase::verbose (   )  [inline, inherited]

{return _verbose; }

void EvtPythia::WritePythiaEntryHeader (  std::ofstream &  outdec, int  lundkc, EvtId  evtnum, std::string  name, int  chg, int  cchg, int  spin2, double  mass, double  width, double  maxwidth, double  ctau, int  stable, double  rawbrfrsum )  [static, private]

{ char sname[100]; char ccsname[100]; // RS changed to 16, new PYTHIA standard int namelength=16; int i,j; int temp; temp = spin2; if (ctau>1000000.0) ctau=0.0; strcpy(sname,name.c_str()); i=0; while (sname[i]!=0){ i++; } // strip up to two + or - if(evtnum.getId()>=0) { if (sname[i-1]=='+'||sname[i-1]=='-'){ sname[i-1]=0; i--; } if (sname[i-1]=='+'||sname[i-1]=='-'){ sname[i-1]=0; i--; } // strip 0 except for _0 and chi...0 if (sname[i-1]=='0' && sname[i-2]!='_' && !(sname[0]=='c' && sname[1]=='h')){ sname[i-1]=0; i--; } } if (i>namelength) { for(j=1;j<namelength;j++){ sname[j]=sname[j+i-namelength]; } sname[namelength]=0; } // RS: copy name for cc particle for(j=0;j<=namelength;j++) ccsname[j]=sname[j]; i=0; while (ccsname[i]!=' '){ i++; if(ccsname[i]==0) break; } if(i<namelength) { ccsname[i]='b'; ccsname[i+1]=0; } //cchg=0; if(evtnum.getId()>=0) { if (abs(EvtPDL::getStdHep(evtnum))==21) cchg=2; if (abs(EvtPDL::getStdHep(evtnum))==90) cchg=-1; if ((abs(EvtPDL::getStdHep(evtnum))<=8)&& (abs(EvtPDL::getStdHep(evtnum))!=0)) cchg=1; } // RS output format changed to new PYTHIA style outdec << " " << setw(9) << lundkc; outdec << " "; outdec.width(namelength); outdec << setiosflags(ios::left) << sname; // RS: name for cc paricle if ((evtnum.getId()>=0) && (EvtPDL::chargeConj(evtnum)!=evtnum)) { outdec << " "; outdec.width(namelength); outdec << ccsname; }else{ // 2+16 spaces outdec << " "; } outdec << resetiosflags(ios::left); outdec << setw(3) << chg; outdec << setw(3) << cchg; outdec.width(3); if (evtnum.getId()>=0) { if (EvtPDL::chargeConj(evtnum)==evtnum) { outdec << 0; } else{ outdec << 1; } } else{ outdec << 0; } outdec.setf(ios::fixed,ios::floatfield); outdec.precision(5); outdec << setw(12) << mass; outdec.setf(ios::fixed,ios::floatfield); outdec << setw(12) << width; outdec.setf(ios::fixed,ios::floatfield); outdec.width(12); if (fabs(width)<0.0000000001) { outdec << 0.0 ; } else{ outdec << maxwidth; } // scientific notation ... outdec << setw(14) << ctau; outdec.setf(ios::scientific,ios::floatfield); outdec << " "; outdec << ctau; outdec.width(3); if (evtnum.getId()>=0) { if (ctau>1.0 || rawbrfrsum<0.000001) { stable=0; } } //resonance width treatment outdec.width(3); outdec << 0; outdec.width(3); outdec << stable; outdec << endl; outdec.width(0); //outdec.setf(0,0); }

void EvtPythia::WritePythiaParticle (  std::ofstream &  outdec, EvtId  ipar, EvtId  iparname, int &  first )  [static, private]

{ int ijetset; double br_sum=0.0; for(ijetset=0;ijetset<njetsetdecays;ijetset++){ if (jetsetdecays[ijetset]->getParentId()==ipar){ br_sum+=jetsetdecays[ijetset]->getBranchingFraction(); } if (jetsetdecays[ijetset]->getParentId()!= EvtPDL::chargeConj(jetsetdecays[ijetset]->getParentId())&& EvtPDL::chargeConj(jetsetdecays[ijetset]->getParentId())==ipar){ br_sum+=jetsetdecays[ijetset]->getBranchingFraction(); } } double br_sum_true=br_sum; if (br_sum<0.000001) br_sum=1.0; for(ijetset=0;ijetset<njetsetdecays;ijetset++){ if (jetsetdecays[ijetset]->getParentId()==ipar){ double br=jetsetdecays[ijetset]->getBranchingFraction(); int i,daugs[5]; EvtId cdaugs[5]; for(i=0;i<5;i++){ if(i<jetsetdecays[ijetset]->getNDaug()){ daugs[i]=EvtPDL::getStdHep( jetsetdecays[ijetset]->getDaugs()[i]); cdaugs[i]=EvtPDL::chargeConj(jetsetdecays[ijetset]->getDaugs()[i]); } else{ daugs[i]=0; } } int channel; channel=EvtDecayTable::findChannel(EvtPDL::chargeConj(ipar), jetsetdecays[ijetset]->getModelName(), jetsetdecays[ijetset]->getNDaug(), cdaugs, jetsetdecays[ijetset]->getNArg(), jetsetdecays[ijetset]->getArgsStr()); if (jetsetdecays[ijetset]->getModelName()=="PYTHIA"){ if (first) { first=0; WritePythiaEntryHeader(outdec, //EvtPDL::getLundKC(iparname), EvtPDL::getStdHep(iparname), iparname, EvtPDL::name(iparname), EvtPDL::chg3(iparname), 0,0,EvtPDL::getMeanMass(ipar), EvtPDL::getWidth(ipar), EvtPDL::getMeanMass(ipar)-EvtPDL::getMinMass(ipar), EvtPDL::getctau(ipar),1,br_sum_true); } int dflag=2; if (EvtPDL::getStdHep(ipar)<0) { dflag=3; for(i=0;i<jetsetdecays[ijetset]->getNDaug();i++){ daugs[i]=EvtPDL::getStdHep(cdaugs[i]); } } /* RS PYTHIA allows to introduce new particles via a call to PYUPDA so no need for this check any more //now lets check to make sure that jetset, lucomp, knows //about all particles! int unknown=0; for(i=0;i<jetsetdecays[ijetset]->getNDaug();i++){ if (pycomp_(&daugs[i])==0) { unknown=1; report(ERROR,"EvtGen") << "Pythia (pycomp) does not " << "know the particle:"<< EvtPDL::name(jetsetdecays[ijetset]->getDaugs()[i])<<endl; } } int istdheppar=EvtPDL::getStdHep(ipar); if (pycomp_(&istdheppar)==0){ unknown=1; report(ERROR,"EvtGen") << "Pythia (pycomp) does not " << "know the particle:"<< EvtPDL::name(ipar)<<endl; } if (unknown){ report(ERROR,"EvtGen") << "Therfore the decay:"<<endl; report(ERROR,"EvtGen") << EvtPDL::name(jetsetdecays[ijetset]->getParentId())<<" -> "; for(i=0;i<jetsetdecays[ijetset]->getNDaug();i++){ report(ERROR,"") << EvtPDL::name(jetsetdecays[ijetset]->getDaugs()[i])<<" "; } report(ERROR,"")<<endl; report(ERROR,"EvtGen")<<"Will not be generated."<<endl; return; } */ if (EvtPDL::chargeConj(ipar)==ipar) { dflag=1; //report(INFO,"EvtGen") << EvtPDL::name(iparname) << " dflag=1 because C(ipar)=ipar!"<<endl; } //if (channel>=0) { // dflag=1; //report(INFO,"EvtGen") << EvtPDL::name(iparname) << " dflag=1 because channel>=0"<<endl; //} // if (!(EvtPDL::getStdHep(ipar)<0&&channel>=0)){ if (1){ // RS changed format to new PYTHIA one outdec << " "; outdec.width(5); outdec <<dflag; outdec.width(5); outdec <<(int)jetsetdecays[ijetset]->getArgs()[0]; outdec.width(12); if (fabs(br)<0.000000001) { outdec <<"0.00000"; } else{ outdec <<br/br_sum; } outdec.width(10); outdec <<daugs[0]; outdec.width(10); outdec <<daugs[1]; outdec.width(10); outdec <<daugs[2]; outdec.width(10); outdec <<daugs[3]; outdec.width(10); outdec <<daugs[4]; outdec<<endl; outdec.width(0); } } } } }

---

Member Data Documentation

bool EvtDecayBase::_daugsDecayedByParentModel [protected, inherited]

std::string * EvtPythia::commands = 0 [static, private]

EvtDecayBasePtr * EvtPythia::jetsetdecays = 0 [static, private]

int EvtPythia::lcommand = 0 [static, private]

int EvtPythia::ncommand = 0 [static, private]

int EvtPythia::njetsetdecays = 0 [static, private]

int EvtPythia::ntable = 0 [static, private]

---

The documentation for this class was generated from the following files:

• /data0/mdestefa/bes3soft/Boss_6.5.5/dist/6.5.5/Generator/BesEvtGen/BesEvtGen-00-01-11/src/EvtGen/EvtGenModels/EvtPythia.hh
• /data0/mdestefa/bes3soft/Boss_6.5.5/dist/6.5.5/Generator/BesEvtGen/BesEvtGen-00-01-11/src/EvtGen/EvtGenModels/EvtPythia.cc

---