1 //--------------------------------------------------------------------------
4 // This software is part of the EvtGen package developed jointly
5 // for the BaBar and CLEO collaborations. If you use all or part
6 // of it, please give an appropriate acknowledgement.
8 // Copyright Information: See BelEvtGen/COPYRIGHT
9 // Copyright (C) 1998 Caltech, UCSB
11 // Module: EvtJetSet.cc
13 // Description: Routine to use JetSet for decaying particles.
15 // Modification history:
17 // RYD July 24, 1997 Module created
18 // RS October 28, 2002 copied from JETSET module
19 //------------------------------------------------------------------------
21 #include "EvtGenBase/EvtPatches.hh"
22 #include "EvtGenBase/EvtPatches.hh"
23 #include "EvtGenBase/EvtParticle.hh"
24 #include "EvtGenBase/EvtStringParticle.hh"
25 #include "EvtGenBase/EvtDecayTable.hh"
26 #include "EvtGenBase/EvtPDL.hh"
27 #include "EvtGenModels/EvtPythia.hh"
28 #include "EvtGenBase/EvtReport.hh"
29 #include "EvtGenBase/EvtId.hh"
41 using std::resetiosflags;
42 using std::setiosflags;
47 int EvtPythia::njetsetdecays=0;
48 EvtDecayBasePtr* EvtPythia::jetsetdecays=0;
49 int EvtPythia::ntable=0;
51 int EvtPythia::ncommand=0;
52 int EvtPythia::lcommand=0;
53 std::string* EvtPythia::commands=0;
57 extern void pycontinuum_(double *,int *, int *,
58 double *,double *,double *,double *);
62 extern void evtpythiainit_(const char* fname, int len);
66 extern void init_cont_();
70 extern void pythiadec_(int *,double *,int *,int *,int *,
71 double *,double *,double *,double *);
75 extern void initpythia_(int *);
79 extern void pygive_(const char *cnfgstr,int length);
83 extern int pycomp_(int* kf);
87 extern void pylist_(int &);
91 EvtPythia::EvtPythia(){}
93 EvtPythia::~EvtPythia(){
99 //the deletion of commands is really uggly!
101 if (njetsetdecays==0) {
107 for(i=0;i<njetsetdecays;i++){
108 if (jetsetdecays[i]==this){
109 jetsetdecays[i]=jetsetdecays[njetsetdecays-1];
111 if (njetsetdecays==0) {
119 report(ERROR,"EvtGen") << "Error in destroying Pythia model!"<<endl;
124 std::string EvtPythia::getName(){
130 EvtDecayBase* EvtPythia::clone(){
132 return new EvtPythia;
137 void EvtPythia::initProbMax(){
144 void EvtPythia::init(){
149 if (getParentId().isAlias()){
151 report(ERROR,"EvtGen") << "EvtPythia finds that you are decaying the"<<endl
152 << " aliased particle "
153 << EvtPDL::name(getParentId()).c_str()
154 << " with the Pythia model"<<endl
155 << " this does not work, please modify decay table."
157 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
167 std::string EvtPythia::commandName(){
169 return std::string("JetSetPar");
174 void EvtPythia::command(std::string cmd){
176 if (ncommand==lcommand){
178 lcommand=10+2*lcommand;
180 std::string* newcommands=new std::string[lcommand];
184 for(i=0;i<ncommand;i++){
185 newcommands[i]=commands[i];
190 commands=newcommands;
194 commands[ncommand]=cmd;
200 void EvtPythia::pythiacont(double *energy, int *ndaugjs, int *kf,
201 double *px, double *py, double *pz, double *e)
203 pycontinuum_(energy,ndaugjs,kf,px,py,pz,e);
208 void EvtPythia::decay( EvtParticle *p){
211 //added by Lange Jan4,2000
212 static EvtId STRNG=EvtPDL::getId("string");
214 int istdheppar=EvtPDL::getStdHep(p->getId());
216 if (pycomp_(&istdheppar)==0){
217 report(ERROR,"EvtGen") << "Pythia can not decay:"
218 <<EvtPDL::name(p->getId()).c_str()<<endl;
227 int ip=EvtPDL::getStdHep(p->getId());
230 EvtId evtnumstable[100],evtnumparton[100];
231 int stableindex[100],partonindex[100];
235 EvtId type[MAX_DAUG];
239 double px[100],py[100],pz[100],e[100];
240 if ( p->getNDaug() != 0 ) { p->deleteDaughters(true);}
246 pythiadec_(&ip,&mp,&ndaugjs,kf,km,px,py,pz,e);
252 for(i=0;i<ndaugjs;i++){
254 if (EvtPDL::evtIdFromStdHep(kf[i])==EvtId(-1,-1)) {
255 report(ERROR,"EvtGen") << "Pythia returned particle:"<<kf[i]<<endl;
256 report(ERROR,"EvtGen") << "This can not be translated to evt number"<<endl;
257 report(ERROR,"EvtGen") << "and the decay will be rejected!"<<endl;
258 report(ERROR,"EvtGen") << "The decay was of particle:"<<ip<<endl;
263 //sort out the partons
264 if (abs(kf[i])<=6||kf[i]==21){
265 partonindex[numparton]=i;
266 evtnumparton[numparton]=EvtPDL::evtIdFromStdHep(kf[i]);
270 stableindex[numstable]=i;
271 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(kf[i]);
276 // have to protect against negative mass^2 for massless particles
277 // i.e. neutrinos and photons.
278 // this is uggly but I need to fix it right now....
280 if (px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i]>=e[i]*e[i]){
282 e[i]=sqrt(px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i])+0.0000000000001;
286 p4[i].set(e[i],px[i],py[i],pz[i]);
291 int channel=EvtDecayTable::inChannelList(p->getId(),numstable,evtnumstable);
301 }while( more && (count<10000) );
304 report(INFO,"EvtGen") << "Too many loops in EvtPythia!!!"<<endl;
305 report(INFO,"EvtGen") << "Parent:"<<EvtPDL::name(getParentId()).c_str()<<endl;
306 for(i=0;i<numstable;i++){
307 report(INFO,"EvtGen") << "Daug("<<i<<")"<<EvtPDL::name(evtnumstable[i]).c_str()<<endl;
316 p->makeDaughters(numstable,evtnumstable);
318 for(i=0;i<numstable;i++){
319 p->getDaug(i)->init(evtnumstable[i],p4[stableindex[i]]);
329 //have partons in JETSET
331 EvtVector4R p4string(0.0,0.0,0.0,0.0);
333 for(i=0;i<numparton;i++){
334 p4string+=p4[partonindex[i]];
339 for(i=0;i<numstable;i++){
340 if (km[stableindex[i]]==0){
341 type[nprimary++]=evtnumstable[i];
345 p->makeDaughters(nprimary,type);
347 p->getDaug(0)->init(STRNG,p4string);
349 EvtVector4R p4partons[10];
351 for(i=0;i<numparton;i++){
352 p4partons[i]=p4[partonindex[i]];
355 ((EvtStringParticle*)p->getDaug(0))->initPartons(numparton,p4partons,evtnumparton);
361 for(i=0;i<numstable;i++){
363 if (km[stableindex[i]]==0){
364 p->getDaug(nprimary++)->init(evtnumstable[i],p4[stableindex[i]]);
370 for(i=0;i<numstable;i++){
371 if (km[stableindex[i]]!=0){
372 type[nsecond++]=evtnumstable[i];
377 p->getDaug(0)->makeDaughters(nsecond,type);
380 for(i=0;i<numstable;i++){
381 if (km[stableindex[i]]!=0){
382 p4[stableindex[i]]=boostTo(p4[stableindex[i]],p4string);
383 p->getDaug(0)->getDaug(nsecond)->init(evtnumstable[i],p4[stableindex[i]]);
384 p->getDaug(0)->getDaug(nsecond)->setDiagonalSpinDensity();
385 p->getDaug(0)->getDaug(nsecond)->decay();
398 void EvtPythia::fixPolarizations(EvtParticle *p){
400 //special case for now to handle the J/psi polarization
402 int ndaug=p->getNDaug();
406 static EvtId Jpsi=EvtPDL::getId("J/psi");
408 for(i=0;i<ndaug;i++){
409 if(p->getDaug(i)->getId()==Jpsi){
426 EvtVector4R p4Psi=p->getDaug(i)->getP4();
428 double alpha=atan2(p4Psi.get(2),p4Psi.get(1));
429 double beta=acos(p4Psi.get(3)/p4Psi.d3mag());
432 p->getDaug(i)->setSpinDensityForwardHelicityBasis(rho,alpha,beta,0.0);
433 setDaughterSpinDensity(i);
440 void EvtPythia::store(EvtDecayBase* jsdecay){
442 if (njetsetdecays==ntable){
444 EvtDecayBasePtr* newjetsetdecays=new EvtDecayBasePtr[2*ntable+10];
446 for(i=0;i<ntable;i++){
447 newjetsetdecays[i]=jetsetdecays[i];
450 delete [] jetsetdecays;
451 jetsetdecays=newjetsetdecays;
454 jetsetdecays[njetsetdecays++]=jsdecay;
461 void EvtPythia::WritePythiaEntryHeader(ofstream &outdec, int lundkc,
462 EvtId evtnum,std::string name,
463 int chg, int cchg, int spin2,double mass,
464 double width, double maxwidth,double ctau,
465 int stable,double rawbrfrsum){
470 // RS changed to 16, new PYTHIA standard
477 if (ctau>1000000.0) ctau=0.0;
479 strcpy(sname,name.c_str());
487 // strip up to two + or -
489 if(evtnum.getId()>=0) {
490 if (sname[i-1]=='+'||sname[i-1]=='-'){
494 if (sname[i-1]=='+'||sname[i-1]=='-'){
498 // strip 0 except for _0 and chi...0
499 if (sname[i-1]=='0' && sname[i-2]!='_' && !(sname[0]=='c' && sname[1]=='h')){
506 for(j=1;j<namelength;j++){
507 sname[j]=sname[j+i-namelength];
512 // RS: copy name for cc particle
513 for(j=0;j<=namelength;j++)
516 while (ccsname[i]!=' '){
518 if(ccsname[i]==0) break;
528 if(evtnum.getId()>=0) {
529 if (abs(EvtPDL::getStdHep(evtnum))==21) cchg=2;
530 if (abs(EvtPDL::getStdHep(evtnum))==90) cchg=-1;
531 if ((abs(EvtPDL::getStdHep(evtnum))<=8)&&
532 (abs(EvtPDL::getStdHep(evtnum))!=0)) cchg=1;
536 // RS output format changed to new PYTHIA style
537 outdec << " " << setw(9) << lundkc;
539 outdec.width(namelength);
540 outdec << setiosflags(ios::left)
542 // RS: name for cc paricle
543 if ((evtnum.getId()>=0) && (EvtPDL::chargeConj(evtnum)!=evtnum))
546 outdec.width(namelength);
553 outdec << resetiosflags(ios::left);
554 outdec << setw(3) << chg;
555 outdec << setw(3) << cchg;
557 if (evtnum.getId()>=0) {
558 if (EvtPDL::chargeConj(evtnum)==evtnum) {
568 outdec.setf(ios::fixed,ios::floatfield);
570 outdec << setw(12) << mass;
571 outdec.setf(ios::fixed,ios::floatfield);
572 outdec << setw(12) << width;
573 outdec.setf(ios::fixed,ios::floatfield);
575 if (fabs(width)<0.0000000001) {
581 // scientific notation ... outdec << setw(14) << ctau;
582 outdec.setf(ios::scientific,ios::floatfield);
586 if (evtnum.getId()>=0) {
587 if (ctau>1.0 || rawbrfrsum<0.000001) {
591 //resonance width treatment
602 void EvtPythia::WritePythiaParticle(ofstream &outdec,EvtId ipar,
603 EvtId iparname,int &first){
609 for(ijetset=0;ijetset<njetsetdecays;ijetset++){
611 if (jetsetdecays[ijetset]->getParentId()==ipar){
612 br_sum+=jetsetdecays[ijetset]->getBranchingFraction();
614 if (jetsetdecays[ijetset]->getParentId()!=
615 EvtPDL::chargeConj(jetsetdecays[ijetset]->getParentId())&&
616 EvtPDL::chargeConj(jetsetdecays[ijetset]->getParentId())==ipar){
617 br_sum+=jetsetdecays[ijetset]->getBranchingFraction();
623 double br_sum_true=br_sum;
625 if (br_sum<0.000001) br_sum=1.0;
627 for(ijetset=0;ijetset<njetsetdecays;ijetset++){
628 if (jetsetdecays[ijetset]->getParentId()==ipar){
630 double br=jetsetdecays[ijetset]->getBranchingFraction();
637 if(i<jetsetdecays[ijetset]->getNDaug()){
638 daugs[i]=EvtPDL::getStdHep(
639 jetsetdecays[ijetset]->getDaugs()[i]);
640 cdaugs[i]=EvtPDL::chargeConj(jetsetdecays[ijetset]->getDaugs()[i]);
649 channel=EvtDecayTable::findChannel(EvtPDL::chargeConj(ipar),
650 jetsetdecays[ijetset]->getModelName(),
651 jetsetdecays[ijetset]->getNDaug(),
653 jetsetdecays[ijetset]->getNArg(),
654 jetsetdecays[ijetset]->getArgsStr());
656 if (jetsetdecays[ijetset]->getModelName()=="PYTHIA"){
660 WritePythiaEntryHeader(outdec,
661 //EvtPDL::getLundKC(iparname),
662 EvtPDL::getStdHep(iparname),
664 EvtPDL::name(iparname),
665 EvtPDL::chg3(iparname),
666 0,0,EvtPDL::getMeanMass(ipar),
667 EvtPDL::getWidth(ipar),
668 EvtPDL::getMeanMass(ipar)-EvtPDL::getMinMass(ipar),
669 EvtPDL::getctau(ipar),1,br_sum_true);
674 if (EvtPDL::getStdHep(ipar)<0) {
676 for(i=0;i<jetsetdecays[ijetset]->getNDaug();i++){
677 daugs[i]=EvtPDL::getStdHep(cdaugs[i]);
683 PYTHIA allows to introduce new particles via a call to PYUPDA
684 so no need for this check any more
686 //now lets check to make sure that jetset, lucomp, knows
687 //about all particles!
689 for(i=0;i<jetsetdecays[ijetset]->getNDaug();i++){
690 if (pycomp_(&daugs[i])==0) {
692 report(ERROR,"EvtGen") << "Pythia (pycomp) does not "
693 << "know the particle:"<<
694 EvtPDL::name(jetsetdecays[ijetset]->getDaugs()[i])<<endl;
698 int istdheppar=EvtPDL::getStdHep(ipar);
700 if (pycomp_(&istdheppar)==0){
702 report(ERROR,"EvtGen") << "Pythia (pycomp) does not "
703 << "know the particle:"<<
704 EvtPDL::name(ipar)<<endl;
710 report(ERROR,"EvtGen") << "Therfore the decay:"<<endl;
711 report(ERROR,"EvtGen") << EvtPDL::name(jetsetdecays[ijetset]->getParentId())<<" -> ";
712 for(i=0;i<jetsetdecays[ijetset]->getNDaug();i++){
713 report(ERROR,"") << EvtPDL::name(jetsetdecays[ijetset]->getDaugs()[i])<<" ";
715 report(ERROR,"")<<endl;
716 report(ERROR,"EvtGen")<<"Will not be generated."<<endl;
721 if (EvtPDL::chargeConj(ipar)==ipar) {
723 //report(INFO,"EvtGen") << EvtPDL::name(iparname) << " dflag=1 because C(ipar)=ipar!"<<endl;
729 //report(INFO,"EvtGen") << EvtPDL::name(iparname) << " dflag=1 because channel>=0"<<endl;
732 // if (!(EvtPDL::getStdHep(ipar)<0&&channel>=0)){
735 // RS changed format to new PYTHIA one
740 outdec <<(int)jetsetdecays[ijetset]->getArgs()[0];
742 if (fabs(br)<0.000000001) {
767 EvtPythia::diquark(int ID)
805 EvtPythia::NominalMass(int ID)
807 // return default mass in PYTHIA
868 NominalCharge(int ID)
870 // return default mass in PYTHIA
930 void EvtPythia::MakePythiaFile(char* fname){
935 //int part_list[MAX_PART];
941 //outdec << "ERROR;"<<endl;
942 //outdec << ";"<<endl;
943 //outdec << ";This decayfile has been automatically created by"<<endl;
944 //outdec << ";EvtGen from the DECAY.DEC file"<<endl;
945 //outdec << ";"<<endl;
949 for(lundkc=1;lundkc<500;lundkc++){
953 for(size_t iipar=0;iipar<EvtPDL::entries();iipar++){
955 ipar=EvtId(iipar,iipar);
956 //no aliased particles!
957 std::string tempStr = EvtPDL::name(ipar);
958 EvtId realId = EvtPDL::getId(tempStr);
959 if ( realId.isAlias() != 0 ) continue;
962 EvtPDL::getStdHep(ipar)==21 ||
963 EvtPDL::getStdHep(ipar)==22 ||
964 EvtPDL::getStdHep(ipar)==23))
967 if (lundkc==EvtPDL::getLundKC(ipar)){
973 WritePythiaParticle(outdec,ipar,ipar,first);
976 EvtId ipar2=EvtPDL::chargeConj(ipar);
980 WritePythiaParticle(outdec,ipar2,ipar,first);
984 WritePythiaEntryHeader(outdec,
985 //EvtPDL::getLundKC(ipar),
986 EvtPDL::getStdHep(ipar),
990 0,0,EvtPDL::getMeanMass(ipar),
991 EvtPDL::getWidth(ipar),
992 EvtPDL::getMeanMass(ipar)-EvtPDL::getMinMass(ipar),
993 EvtPDL::getctau(ipar),0,0.0);
999 if (lundkc==99999) // Write out diquarks after quarks, but only once
1000 for(size_t iipar=0;iipar<EvtPDL::entries();iipar++){
1002 ipar=EvtId(iipar,iipar);
1004 if (diquark(EvtPDL::getStdHep(ipar))){
1010 WritePythiaParticle(outdec,ipar,ipar,first);
1013 EvtId ipar2=EvtPDL::chargeConj(ipar);
1017 WritePythiaParticle(outdec,ipar2,ipar,first);
1021 WritePythiaEntryHeader(outdec,
1022 EvtPDL::getStdHep(ipar),
1025 NominalCharge(EvtPDL::getStdHep(ipar)),
1027 NominalMass(EvtPDL::getStdHep(ipar)),
1035 WritePythiaEntryHeader(outdec,
1036 lundkc,EvtId(-1,-1)," ",
1037 0,0,0,EvtPDL::getNominalMass(ipar),0.0,0.0,
1038 EvtPDL::getctau(ipar),0,0.0);
1045 void EvtPythia::pythiaInit(int dummy){
1047 //static int first=1;
1048 static int first=0; //if first=0 Pythia is not reinitialize
1053 report(INFO,"EvtGen") << "Will initialize Pythia."<<endl;
1054 for(int i=0;i<ncommand;i++)
1055 pygive_(commands[i].c_str(),strlen(commands[i].c_str()));
1059 char hostBuffer[100];
1061 if ( gethostname( hostBuffer, 100 ) != 0 ){
1062 report(ERROR,"EvtGen") << " couldn't get hostname." << endl;
1063 strncpy( hostBuffer, "hostnameNotFound", 100 );
1068 int thePid=getpid();
1070 if ( sprintf( pid, "%d", thePid ) == 0 ){
1071 report(ERROR,"EvtGen") << " couldn't get process ID." << endl;
1072 strncpy( pid, "666", 100 );
1075 strcpy(fname,"jet.d-");
1076 strcat(fname,hostBuffer);
1080 MakePythiaFile(fname);
1081 evtpythiainit_(fname,strlen(fname));
1082 initpythia_(&dummy);
1084 if (0==getenv("EVTSAVEJETD")){
1086 strcpy(delcmd,"rm -f ");
1087 strcat(delcmd,fname);
1091 report(INFO,"EvtGen") << "Done initializing Pythia."<<endl;