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 EvtGen/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
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/EvtJetSet.hh"
28 #include "EvtGenBase/EvtReport.hh"
30 #include "EvtGenBase/EvtId.hh"
42 using std::resetiosflags;
43 using std::setiosflags;
47 int EvtJetSet::njetsetdecays=0;
48 EvtDecayBasePtr* EvtJetSet::jetsetdecays=0;
49 int EvtJetSet::ntable=0;
51 int EvtJetSet::ncommand=0;
52 int EvtJetSet::lcommand=0;
53 std::string* EvtJetSet::commands=0;
56 extern void evtjetsetinit_(char* fname, int len);
61 extern void jetset1_(int *,double *,int *,int *,int *,
62 double *,double *,double *,double *);
66 extern void lugive_(const char *cnfgstr,int length);
70 extern int lucomp_(int* kf);
74 EvtJetSet::EvtJetSet(){}
76 EvtJetSet::~EvtJetSet(){
82 //the deletion of commands is really uggly!
84 if (njetsetdecays==0) {
90 for(i=0;i<njetsetdecays;i++){
91 if (jetsetdecays[i]==this){
92 jetsetdecays[i]=jetsetdecays[njetsetdecays-1];
94 if (njetsetdecays==0) {
102 report(ERROR,"EvtGen") << "Error in destroying JetSet model!"<<endl;
107 std::string EvtJetSet::getName(){
113 EvtDecayBase* EvtJetSet::clone(){
115 return new EvtJetSet;
120 void EvtJetSet::initProbMax(){
127 void EvtJetSet::init(){
132 if (getParentId().isAlias()){
134 report(ERROR,"EvtGen") << "EvtJetSet finds that you are decaying the"<<endl
135 << " aliased particle "
136 << EvtPDL::name(getParentId()).c_str()
137 << " with the JetSet model"<<endl
138 << " this does not work, please modify decay table."
140 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
150 std::string EvtJetSet::commandName(){
152 return std::string("JetSetPar");
157 void EvtJetSet::command(std::string cmd){
159 if (ncommand==lcommand){
161 lcommand=10+2*lcommand;
163 std::string* newcommands=new std::string[lcommand];
167 for(i=0;i<ncommand;i++){
168 newcommands[i]=commands[i];
173 commands=newcommands;
177 commands[ncommand]=cmd;
185 void EvtJetSet::decay( EvtParticle *p){
188 //added by Lange Jan4,2000
189 static EvtId STRNG=EvtPDL::getId("string");
191 int istdheppar=EvtPDL::getStdHep(p->getId());
193 if (lucomp_(&istdheppar)==0){
194 report(ERROR,"EvtGen") << "Jetset can not decay:"
195 <<EvtPDL::name(p->getId()).c_str()<<endl;
204 int ip=EvtPDL::getStdHep(p->getId());
207 EvtId evtnumstable[100],evtnumparton[100];
208 int stableindex[100],partonindex[100];
212 EvtId type[MAX_DAUG];
216 double px[100],py[100],pz[100],e[100];
218 if ( p->getNDaug() != 0 ) { p->deleteDaughters(true);}
223 //report(INFO,"EvtGen") << "calling jetset " << ip<< " " << mp <<endl;
224 jetset1_(&ip,&mp,&ndaugjs,kf,km,px,py,pz,e);
229 //report(INFO,"EvtGen") << "found some daughters " << ndaugjs << endl;
230 for(i=0;i<ndaugjs;i++){
232 if (EvtPDL::evtIdFromStdHep(kf[i])==EvtId(-1,-1)) {
233 report(ERROR,"EvtGen") << "JetSet returned particle:"<<kf[i]<<endl;
234 report(ERROR,"EvtGen") << "This can not be translated to evt number"<<endl;
235 report(ERROR,"EvtGen") << "and the decay will be rejected!"<<endl;
236 report(ERROR,"EvtGen") << "The decay was of particle:"<<ip<<endl;
240 //sort out the partons
241 if (abs(kf[i])<=6||kf[i]==21){
242 partonindex[numparton]=i;
243 evtnumparton[numparton]=EvtPDL::evtIdFromStdHep(kf[i]);
247 stableindex[numstable]=i;
248 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(kf[i]);
253 // have to protect against negative mass^2 for massless particles
254 // i.e. neutrinos and photons.
255 // this is uggly but I need to fix it right now....
257 if (px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i]>=e[i]*e[i]){
259 e[i]=sqrt(px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i])+0.0000000000001;
263 p4[i].set(e[i],px[i],py[i],pz[i]);
268 int channel=EvtDecayTable::inChannelList(p->getId(),numstable,evtnumstable);
278 }while( more && (count<10000) );
281 report(INFO,"EvtGen") << "Too many loops in EvtJetSet!!!"<<endl;
282 report(INFO,"EvtGen") << "Parent:"<<EvtPDL::name(getParentId()).c_str()<<endl;
283 for(i=0;i<numstable;i++){
284 report(INFO,"EvtGen") << "Daug("<<i<<")"<<EvtPDL::name(evtnumstable[i]).c_str()<<endl;
293 p->makeDaughters(numstable,evtnumstable);
295 for(i=0;i<numstable;i++){
296 p->getDaug(i)->init(evtnumstable[i],p4[stableindex[i]]);
299 if ( ndaugFound == 0 ) {
300 report(ERROR,"EvtGen") << "Jetset has failed to do a decay ";
301 report(ERROR,"EvtGen") << EvtPDL::name(p->getId()).c_str() << " " << p->mass()<<endl;
312 //have partons in JETSET
314 EvtVector4R p4string(0.0,0.0,0.0,0.0);
316 for(i=0;i<numparton;i++){
317 p4string+=p4[partonindex[i]];
322 for(i=0;i<numstable;i++){
323 if (km[stableindex[i]]==0){
324 type[nprimary++]=evtnumstable[i];
328 p->makeDaughters(nprimary,type);
330 p->getDaug(0)->init(STRNG,p4string);
332 EvtVector4R p4partons[10];
334 for(i=0;i<numparton;i++){
335 p4partons[i]=p4[partonindex[i]];
338 ((EvtStringParticle*)p->getDaug(0))->initPartons(numparton,p4partons,evtnumparton);
344 for(i=0;i<numstable;i++){
346 if (km[stableindex[i]]==0){
347 p->getDaug(nprimary++)->init(evtnumstable[i],p4[stableindex[i]]);
353 for(i=0;i<numstable;i++){
354 if (km[stableindex[i]]!=0){
355 type[nsecond++]=evtnumstable[i];
360 p->getDaug(0)->makeDaughters(nsecond,type);
362 EvtVector4R p4stringboost(p4string.get(0),-p4string.get(1),
363 -p4string.get(2),-p4string.get(3));
366 for(i=0;i<numstable;i++){
367 if (km[stableindex[i]]!=0){
368 p4[stableindex[i]]=boostTo(p4[stableindex[i]],p4stringboost);
369 p->getDaug(0)->getDaug(nsecond)->init(evtnumstable[i],p4[stableindex[i]]);
370 p->getDaug(0)->getDaug(nsecond)->setDiagonalSpinDensity();
371 p->getDaug(0)->getDaug(nsecond)->decay();
376 if ( nsecond == 0 ) {
377 report(ERROR,"EvtGen") << "Jetset has failed to do a decay ";
378 report(ERROR,"EvtGen") << EvtPDL::name(p->getId()).c_str() << " " << p->mass() <<endl;
390 void EvtJetSet::fixPolarizations(EvtParticle *p){
392 //special case for now to handle the J/psi polarization
394 int ndaug=p->getNDaug();
398 static EvtId Jpsi=EvtPDL::getId("J/psi");
400 for(i=0;i<ndaug;i++){
401 if(p->getDaug(i)->getId()==Jpsi){
418 EvtVector4R p4Psi=p->getDaug(i)->getP4();
420 double alpha=atan2(p4Psi.get(2),p4Psi.get(1));
421 double beta=acos(p4Psi.get(3)/p4Psi.d3mag());
424 p->getDaug(i)->setSpinDensityForwardHelicityBasis(rho,alpha,beta,0.0);
425 setDaughterSpinDensity(i);
432 void EvtJetSet::store(EvtDecayBase* jsdecay){
434 if (njetsetdecays==ntable){
436 EvtDecayBasePtr* newjetsetdecays=new EvtDecayBasePtr[2*ntable+10];
438 for(i=0;i<ntable;i++){
439 newjetsetdecays[i]=jetsetdecays[i];
442 delete [] jetsetdecays;
443 jetsetdecays=newjetsetdecays;
446 jetsetdecays[njetsetdecays++]=jsdecay;
453 void EvtJetSet::WriteJetSetEntryHeader(ofstream &outdec, int lundkc,
454 EvtId evtnum,std::string name,
455 int chg, int cchg, int spin2,double mass,
456 double width, double maxwidth,double ctau,
457 int stable,double rawbrfrsum){
468 if (ctau>1000000.0) ctau=0.0;
470 strcpy(sname,name.c_str());
478 // strip up to two + or -
480 if(evtnum.getId()>=0) {
481 if (sname[i-1]=='+'||sname[i-1]=='-'){
485 if (sname[i-1]=='+'||sname[i-1]=='-'){
489 // strip 0 except for _0 and chi...0
490 if (sname[i-1]=='0' && sname[i-2]!='_' && !(sname[0]=='c' && sname[1]=='h')){
497 for(j=1;j<namelength;j++){
498 sname[j]=sname[j+i-namelength];
505 if(evtnum.getId()>=0) {
506 if (abs(EvtPDL::getStdHep(evtnum))==21) cchg=2;
507 if (abs(EvtPDL::getStdHep(evtnum))==90) cchg=-1;
508 if ((abs(EvtPDL::getStdHep(evtnum))<=8)&&
509 (abs(EvtPDL::getStdHep(evtnum))!=0)) cchg=1;
513 outdec << setw(5) << lundkc << " ";
514 outdec.width(namelength);
515 outdec << setiosflags(ios::left) << sname << resetiosflags(ios::left);
516 outdec << setw(3) << chg;
517 outdec << setw(3) << cchg;
519 if (evtnum.getId()>=0) {
520 if (EvtPDL::chargeConj(evtnum)==evtnum) {
530 outdec.setf(ios::fixed);
532 outdec << setw(12) << mass;
533 outdec << setw(12) << width;
535 if (fabs(width)<0.0000000001) {
541 outdec << setw(14) << ctau;
543 if (evtnum.getId()>=0) {
544 if (ctau>1.0 || rawbrfrsum<0.000001) {
554 void EvtJetSet::WriteJetSetParticle(ofstream &outdec,EvtId ipar,
555 EvtId iparname,int &first){
561 for(ijetset=0;ijetset<njetsetdecays;ijetset++){
563 if (jetsetdecays[ijetset]->getParentId()==ipar){
564 br_sum+=jetsetdecays[ijetset]->getBranchingFraction();
566 if (jetsetdecays[ijetset]->getParentId()!=
567 EvtPDL::chargeConj(jetsetdecays[ijetset]->getParentId())&&
568 EvtPDL::chargeConj(jetsetdecays[ijetset]->getParentId())==ipar){
569 br_sum+=jetsetdecays[ijetset]->getBranchingFraction();
575 double br_sum_true=br_sum;
577 if (br_sum<0.000001) br_sum=1.0;
579 for(ijetset=0;ijetset<njetsetdecays;ijetset++){
580 if (jetsetdecays[ijetset]->getParentId()==ipar){
582 double br=jetsetdecays[ijetset]->getBranchingFraction();
589 if(i<jetsetdecays[ijetset]->getNDaug()){
590 daugs[i]=EvtPDL::getStdHep(
591 jetsetdecays[ijetset]->getDaugs()[i]);
592 cdaugs[i]=EvtPDL::chargeConj(jetsetdecays[ijetset]->getDaugs()[i]);
601 channel=EvtDecayTable::findChannel(EvtPDL::chargeConj(ipar),
602 jetsetdecays[ijetset]->getModelName(),
603 jetsetdecays[ijetset]->getNDaug(),
605 jetsetdecays[ijetset]->getNArg(),
606 jetsetdecays[ijetset]->getArgsStr());
608 if (jetsetdecays[ijetset]->getModelName()=="JETSET"){
612 WriteJetSetEntryHeader(outdec,
613 EvtPDL::getLundKC(iparname),
615 EvtPDL::name(iparname),
616 EvtPDL::chg3(iparname),
617 0,0,EvtPDL::getMeanMass(ipar),
618 EvtPDL::getWidth(ipar),
619 EvtPDL::getMeanMass(ipar)-EvtPDL::getMinMass(ipar),
620 EvtPDL::getctau(ipar),1,br_sum_true);
625 if (EvtPDL::getStdHep(ipar)<0) {
627 for(i=0;i<jetsetdecays[ijetset]->getNDaug();i++){
628 daugs[i]=EvtPDL::getStdHep(cdaugs[i]);
633 //now lets check to make sure that jetset, lucomp, knows
634 //about all particles!
636 for(i=0;i<jetsetdecays[ijetset]->getNDaug();i++){
637 if (lucomp_(&daugs[i])==0) {
639 report(ERROR,"EvtGen") << "JetSet (lucomp) does not "
640 << "know the particle:"<<
641 EvtPDL::name(jetsetdecays[ijetset]->getDaugs()[i]).c_str()<<endl;
645 int istdheppar=EvtPDL::getStdHep(ipar);
647 if (lucomp_(&istdheppar)==0){
649 report(ERROR,"EvtGen") << "JetSet (lucomp) does not "
650 << "know the particle:"<<
651 EvtPDL::name(ipar).c_str()<<endl;
657 report(ERROR,"EvtGen") << "Therfore the decay:"<<endl;
658 report(ERROR,"EvtGen") << EvtPDL::name(jetsetdecays[ijetset]->getParentId()).c_str()<<" -> ";
659 for(i=0;i<jetsetdecays[ijetset]->getNDaug();i++){
660 report(ERROR,"") << EvtPDL::name(jetsetdecays[ijetset]->getDaugs()[i]).c_str()<<" ";
662 report(ERROR,"")<<endl;
663 report(ERROR,"EvtGen")<<"Will not be generated."<<endl;
668 if (EvtPDL::chargeConj(ipar)==ipar) {
670 //report(INFO,"EvtGen") << EvtPDL::name(iparname) << " dflag=1 because C(ipar)=ipar!"<<endl;
676 //report(INFO,"EvtGen") << EvtPDL::name(iparname) << " dflag=1 because channel>=0"<<endl;
679 // if (!(EvtPDL::getStdHep(ipar)<0&&channel>=0)){
685 outdec <<(int)jetsetdecays[ijetset]->getArgs()[0];
687 if (fabs(br)<0.000000001) {
712 void EvtJetSet::MakeJetSetFile(char* fname){
717 //int part_list[MAX_PART];
723 //outdec << ";"<<endl;
724 //outdec << ";This decayfile has been automatically created by"<<endl;
725 //outdec << ";EvtGen from the DECAY.DEC file"<<endl;
726 //outdec << ";"<<endl;
730 for(lundkc=1;lundkc<=500;lundkc++){
736 for(iipar=0;iipar< static_cast<int>( EvtPDL::entries()) ;iipar++){
738 ipar=EvtId(iipar,iipar);
739 //no aliased particles!
740 std::string tempStr = EvtPDL::name(ipar);
741 EvtId realId = EvtPDL::getId(tempStr);
742 if ( realId.isAlias() != 0 ) continue;
743 if (lundkc==EvtPDL::getLundKC(ipar)){
749 WriteJetSetParticle(outdec,ipar,ipar,first);
752 EvtId ipar2=EvtPDL::chargeConj(ipar);
756 WriteJetSetParticle(outdec,ipar2,ipar,first);
760 WriteJetSetEntryHeader(outdec,
761 EvtPDL::getLundKC(ipar),
765 0,0,EvtPDL::getMeanMass(ipar),
766 EvtPDL::getWidth(ipar),
767 EvtPDL::getMeanMass(ipar)-EvtPDL::getMinMass(ipar),
768 EvtPDL::getctau(ipar),0,0.0);
775 WriteJetSetEntryHeader(outdec,
776 lundkc,EvtId(-1,-1)," ",
777 0,0,0,EvtPDL::getMeanMass(ipar),0.0,0.0,
778 EvtPDL::getctau(ipar),0,0.0);
785 void EvtJetSet::jetSetInit(){
793 report(INFO,"EvtGen") << "Will initialize JetSet."<<endl;
797 char hostBuffer[100];
799 if ( gethostname( hostBuffer, 100 ) != 0 ){
800 report(ERROR,"EvtGen") << " couldn't get hostname." << endl;
801 strncpy( hostBuffer, "hostnameNotFound", 100 );
808 if ( sprintf( pid, "%d", thePid ) == 0 ){
809 report(ERROR,"EvtGen") << " couldn't get process ID." << endl;
810 strncpy( pid, "666", 100 );
813 strcpy(fname,"jet.d-");
814 strcat(fname,hostBuffer);
818 MakeJetSetFile(fname);
819 evtjetsetinit_(fname,strlen(fname));
821 if (0==getenv("EVTSAVEJETD")){
823 strcpy(delcmd,"rm -f ");
824 strcat(delcmd,fname);
830 for(i=0;i<ncommand;i++){
831 lugive_(commands[i].c_str(),strlen(commands[i].c_str()));
835 report(INFO,"EvtGen") << "Done initializing JetSet."<<endl;