TEvtGen/EvtGenModels/EvtPythia.cxx

   1 //--------------------------------------------------------------------------
   2 //
   3 // Environment:
   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.
   7 //
   8 // Copyright Information: See BelEvtGen/COPYRIGHT
   9 //      Copyright (C) 1998      Caltech, UCSB
  10 //
  11 // Module: EvtJetSet.cc
  12 //
  13 // Description: Routine to use JetSet for decaying particles.
  14 //
  15 // Modification history:
  16 //
  17 //    RYD     July 24, 1997        Module created
  18 //    RS      October 28, 2002        copied from JETSET module
  19 //------------------------------------------------------------------------
  20 //
  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"
  30 #include <iostream>
  31 #include <iomanip>
  32 #include <fstream>
  33 #include <string.h>
  34 #include <stdlib.h>
  35 #include <unistd.h>
  36 #include <stdio.h>
  37 using std::endl;
  38 using std::fstream;
  39 using std::ios;
  40 using std::ofstream;
  41 using std::resetiosflags;
  42 using std::setiosflags;
  43 using std::setw;
  44
  45 using std::string;
  46
  47 int EvtPythia::njetsetdecays=0;
  48   EvtDecayBasePtr* EvtPythia::jetsetdecays=0;
  49 int EvtPythia::ntable=0;
  50
  51 int EvtPythia::ncommand=0;
  52 int EvtPythia::lcommand=0;
  53 std::string* EvtPythia::commands=0;
  54
  55
  56 extern "C" {
  57   extern void pycontinuum_(double *,int *, int *,
  58                             double *,double *,double *,double *);
  59 }
  60
  61 extern "C" {
  62   extern void evtpythiainit_(const char* fname, int len);
  63 }
  64
  65 extern "C" {
  66   extern void init_cont_();
  67 }
  68
  69 extern "C" {
  70   extern void pythiadec_(int *,double *,int *,int *,int *,
  71                           double *,double *,double *,double *);
  72 }
  73
  74 extern "C" {
  75   extern void initpythia_(int *);
  76 }
  77
  78 extern "C" {
  79   extern void pygive_(const char *cnfgstr,int length);
  80 }
  81
  82 extern "C" {
  83   extern int pycomp_(int* kf);
  84 }
  85
  86 extern "C" {
  87   extern void pylist_(int &);
  88 }
  89
  90
  91 EvtPythia::EvtPythia(){}
  92
  93 EvtPythia::~EvtPythia(){
  94
  95
  96   int i;
  97
  98
  99   //the deletion of commands is really uggly!
 100
 101   if (njetsetdecays==0) {
 102     delete [] commands;
 103     commands=0;
 104     return;
 105   }
 106
 107   for(i=0;i<njetsetdecays;i++){
 108     if (jetsetdecays[i]==this){
 109       jetsetdecays[i]=jetsetdecays[njetsetdecays-1];
 110       njetsetdecays--;
 111       if (njetsetdecays==0) {
 112         delete [] commands;
 113         commands=0;
 114       }
 115       return;
 116     }
 117   }
 118
 119   report(ERROR,"EvtGen") << "Error in destroying Pythia model!"<<endl;
 120
 121 }
 122
 123
 124 std::string EvtPythia::getName(){
 125
 126   return "PYTHIA";
 127
 128 }
 129
 130 EvtDecayBase* EvtPythia::clone(){
 131
 132   return new EvtPythia;
 133
 134 }
 135
 136
 137 void EvtPythia::initProbMax(){
 138
 139   noProbMax();
 140
 141 }
 142
 143
 144 void EvtPythia::init(){
 145
 146   checkNArg(1);
 147
 148
 149   if (getParentId().isAlias()){
 150
 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."
 156                            << endl;
 157     report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
 158     ::abort();
 159
 160   }
 161
 162   store(this);
 163
 164 }
 165
 166
 167 std::string EvtPythia::commandName(){
 168
 169   return std::string("JetSetPar");
 170
 171 }
 172
 173
 174 void EvtPythia::command(std::string cmd){
 175
 176   if (ncommand==lcommand){
 177
 178     lcommand=10+2*lcommand;
 179
 180     std::string* newcommands=new std::string[lcommand];
 181
 182     int i;
 183
 184     for(i=0;i<ncommand;i++){
 185       newcommands[i]=commands[i];
 186     }
 187
 188     delete [] commands;
 189
 190     commands=newcommands;
 191
 192   }
 193
 194   commands[ncommand]=cmd;
 195
 196   ncommand++;
 197
 198 }
 199
 200 void EvtPythia::pythiacont(double *energy, int *ndaugjs, int *kf,
 201                            double *px, double *py, double *pz, double *e)
 202 {
 203   pycontinuum_(energy,ndaugjs,kf,px,py,pz,e);
 204 }
 205
 206
 207
 208 void EvtPythia::decay( EvtParticle *p){
 209
 210
 211   //added by Lange Jan4,2000
 212   static EvtId STRNG=EvtPDL::getId("string");
 213
 214   int istdheppar=EvtPDL::getStdHep(p->getId());
 215
 216   if (pycomp_(&istdheppar)==0){
 217     report(ERROR,"EvtGen") << "Pythia can not decay:"
 218       <<EvtPDL::name(p->getId()).c_str()<<endl;
 219     return;
 220   }
 221
 222   double mp=p->mass();
 223
 224   EvtVector4R p4[20];
 225
 226   int i,more;
 227   int ip=EvtPDL::getStdHep(p->getId());
 228   int ndaugjs;
 229   int kf[100];
 230   EvtId evtnumstable[100],evtnumparton[100];
 231   int stableindex[100],partonindex[100];
 232   int numstable;
 233   int numparton;
 234   int km[100];
 235   EvtId type[MAX_DAUG];
 236
 237   pythiaInit(0);
 238
 239   double px[100],py[100],pz[100],e[100];
 240   if ( p->getNDaug() != 0 ) { p->deleteDaughters(true);}
 241
 242   int count=0;
 243
 244   do{
 245
 246     pythiadec_(&ip,&mp,&ndaugjs,kf,km,px,py,pz,e);
 247
 248
 249     numstable=0;
 250     numparton=0;
 251
 252     for(i=0;i<ndaugjs;i++){
 253
 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;
 259         int i=1;
 260         pylist_(i);
 261       }
 262
 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]);
 267         numparton++;
 268       }
 269       else{
 270         stableindex[numstable]=i;
 271         evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(kf[i]);
 272         numstable++;
 273       }
 274
 275
 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....
 279
 280       if (px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i]>=e[i]*e[i]){
 281
 282         e[i]=sqrt(px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i])+0.0000000000001;
 283
 284       }
 285
 286       p4[i].set(e[i],px[i],py[i],pz[i]);
 287
 288
 289     }
 290
 291     int channel=EvtDecayTable::inChannelList(p->getId(),numstable,evtnumstable);
 292
 293
 294    more=(channel!=-1);
 295
 296
 297
 298
 299    count++;
 300
 301   }while( more && (count<10000) );
 302
 303   if (count>9999) {
 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;
 308     }
 309
 310   }
 311
 312
 313
 314   if (numparton==0){
 315
 316     p->makeDaughters(numstable,evtnumstable);
 317
 318     for(i=0;i<numstable;i++){
 319       p->getDaug(i)->init(evtnumstable[i],p4[stableindex[i]]);
 320     }
 321
 322     fixPolarizations(p);
 323
 324     return ;
 325
 326   }
 327   else{
 328
 329     //have partons in JETSET
 330
 331     EvtVector4R p4string(0.0,0.0,0.0,0.0);
 332
 333     for(i=0;i<numparton;i++){
 334       p4string+=p4[partonindex[i]];
 335     }
 336
 337     int nprimary=1;
 338     type[0]=STRNG;
 339     for(i=0;i<numstable;i++){
 340       if (km[stableindex[i]]==0){
 341         type[nprimary++]=evtnumstable[i];
 342       }
 343     }
 344
 345     p->makeDaughters(nprimary,type);
 346
 347     p->getDaug(0)->init(STRNG,p4string);
 348
 349     EvtVector4R p4partons[10];
 350
 351     for(i=0;i<numparton;i++){
 352       p4partons[i]=p4[partonindex[i]];
 353     }
 354
 355     ((EvtStringParticle*)p->getDaug(0))->initPartons(numparton,p4partons,evtnumparton);
 356
 357
 358
 359     nprimary=1;
 360
 361     for(i=0;i<numstable;i++){
 362
 363       if (km[stableindex[i]]==0){
 364         p->getDaug(nprimary++)->init(evtnumstable[i],p4[stableindex[i]]);
 365       }
 366     }
 367
 368
 369     int nsecond=0;
 370     for(i=0;i<numstable;i++){
 371       if (km[stableindex[i]]!=0){
 372         type[nsecond++]=evtnumstable[i];
 373       }
 374     }
 375
 376
 377     p->getDaug(0)->makeDaughters(nsecond,type);
 378
 379     nsecond=0;
 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();
 386         nsecond++;
 387       }
 388     }
 389
 390     fixPolarizations(p);
 391
 392     return ;
 393
 394   }
 395
 396 }
 397
 398 void EvtPythia::fixPolarizations(EvtParticle *p){
 399
 400   //special case for now to handle the J/psi polarization
 401
 402   int ndaug=p->getNDaug();
 403
 404   int i;
 405
 406   static EvtId Jpsi=EvtPDL::getId("J/psi");
 407
 408   for(i=0;i<ndaug;i++){
 409     if(p->getDaug(i)->getId()==Jpsi){
 410
 411       EvtSpinDensity rho;
 412
 413       rho.setDim(3);
 414       rho.set(0,0,0.5);
 415       rho.set(0,1,0.0);
 416       rho.set(0,2,0.0);
 417
 418       rho.set(1,0,0.0);
 419       rho.set(1,1,1.0);
 420       rho.set(1,2,0.0);
 421
 422       rho.set(2,0,0.0);
 423       rho.set(2,1,0.0);
 424       rho.set(2,2,0.5);
 425
 426       EvtVector4R p4Psi=p->getDaug(i)->getP4();
 427
 428       double alpha=atan2(p4Psi.get(2),p4Psi.get(1));
 429       double beta=acos(p4Psi.get(3)/p4Psi.d3mag());
 430
 431
 432       p->getDaug(i)->setSpinDensityForwardHelicityBasis(rho,alpha,beta,0.0);
 433       setDaughterSpinDensity(i);
 434
 435     }
 436   }
 437
 438 }
 439
 440 void EvtPythia::store(EvtDecayBase* jsdecay){
 441
 442   if (njetsetdecays==ntable){
 443
 444     EvtDecayBasePtr* newjetsetdecays=new EvtDecayBasePtr[2*ntable+10];
 445     int i;
 446     for(i=0;i<ntable;i++){
 447       newjetsetdecays[i]=jetsetdecays[i];
 448     }
 449     ntable=2*ntable+10;
 450     delete [] jetsetdecays;
 451     jetsetdecays=newjetsetdecays;
 452   }
 453
 454   jetsetdecays[njetsetdecays++]=jsdecay;
 455
 456
 457
 458 }
 459
 460
 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){
 466
 467   char sname[100];
 468   char ccsname[100];
 469
 470   // RS changed to 16, new PYTHIA standard
 471   int namelength=16;
 472
 473   int i,j;
 474   int temp;
 475   temp = spin2;
 476
 477   if (ctau>1000000.0) ctau=0.0;
 478
 479   strcpy(sname,name.c_str());
 480
 481   i=0;
 482
 483   while (sname[i]!=0){
 484     i++;
 485   }
 486
 487   // strip up to two + or -
 488
 489   if(evtnum.getId()>=0) {
 490     if (sname[i-1]=='+'||sname[i-1]=='-'){
 491       sname[i-1]=0;
 492       i--;
 493     }
 494     if (sname[i-1]=='+'||sname[i-1]=='-'){
 495       sname[i-1]=0;
 496       i--;
 497     }
 498     // strip 0 except for _0 and chi...0
 499     if (sname[i-1]=='0' && sname[i-2]!='_' && !(sname[0]=='c' && sname[1]=='h')){
 500       sname[i-1]=0;
 501       i--;
 502     }
 503   }
 504
 505   if (i>namelength) {
 506     for(j=1;j<namelength;j++){
 507       sname[j]=sname[j+i-namelength];
 508     }
 509     sname[namelength]=0;
 510   }
 511
 512   // RS: copy name for cc particle
 513   for(j=0;j<=namelength;j++)
 514     ccsname[j]=sname[j];
 515   i=0;
 516   while (ccsname[i]!=' '){
 517     i++;
 518     if(ccsname[i]==0) break;
 519   }
 520   if(i<namelength)
 521     {
 522       ccsname[i]='b';
 523       ccsname[i+1]=0;
 524     }
 525
 526   //cchg=0;
 527
 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;
 533
 534   }
 535
 536   // RS output format changed to new PYTHIA style
 537   outdec << " " << setw(9) << lundkc;
 538   outdec << "  ";
 539   outdec.width(namelength);
 540   outdec << setiosflags(ios::left)
 541          << sname;
 542   // RS: name for cc paricle
 543   if ((evtnum.getId()>=0) && (EvtPDL::chargeConj(evtnum)!=evtnum))
 544     {
 545       outdec << "  ";
 546       outdec.width(namelength);
 547       outdec << ccsname;
 548     }else{
 549       // 2+16 spaces
 550       outdec << "                  ";
 551     }
 552
 553   outdec << resetiosflags(ios::left);
 554   outdec << setw(3) << chg;
 555   outdec << setw(3) << cchg;
 556   outdec.width(3);
 557   if (evtnum.getId()>=0) {
 558     if (EvtPDL::chargeConj(evtnum)==evtnum) {
 559       outdec << 0;
 560     }
 561     else{
 562       outdec << 1;
 563     }
 564   }
 565   else{
 566     outdec << 0;
 567   }
 568   outdec.setf(ios::fixed,ios::floatfield);
 569   outdec.precision(5);
 570   outdec << setw(12) << mass;
 571   outdec.setf(ios::fixed,ios::floatfield);
 572   outdec << setw(12) << width;
 573   outdec.setf(ios::fixed,ios::floatfield);
 574   outdec.width(12);
 575   if (fabs(width)<0.0000000001) {
 576     outdec << 0.0 ;
 577   }
 578   else{
 579     outdec << maxwidth;
 580   }
 581   // scientific notation ...  outdec << setw(14) << ctau;
 582   outdec.setf(ios::scientific,ios::floatfield);
 583   outdec << "  ";
 584   outdec << ctau;
 585   outdec.width(3);
 586   if (evtnum.getId()>=0) {
 587     if (ctau>1.0 || rawbrfrsum<0.000001) {
 588       stable=0;
 589     }
 590   }
 591   //resonance width treatment
 592   outdec.width(3);
 593   outdec << 0;
 594   outdec.width(3);
 595   outdec << stable;
 596   outdec << endl;
 597   outdec.width(0);
 598   //outdec.setf(0,0);
 599
 600 }
 601
 602 void EvtPythia::WritePythiaParticle(ofstream &outdec,EvtId ipar,
 603                                     EvtId iparname,int &first){
 604
 605   int ijetset;
 606
 607   double br_sum=0.0;
 608
 609   for(ijetset=0;ijetset<njetsetdecays;ijetset++){
 610
 611     if (jetsetdecays[ijetset]->getParentId()==ipar){
 612       br_sum+=jetsetdecays[ijetset]->getBranchingFraction();
 613     }
 614     if (jetsetdecays[ijetset]->getParentId()!=
 615         EvtPDL::chargeConj(jetsetdecays[ijetset]->getParentId())&&
 616         EvtPDL::chargeConj(jetsetdecays[ijetset]->getParentId())==ipar){
 617       br_sum+=jetsetdecays[ijetset]->getBranchingFraction();
 618     }
 619
 620
 621   }
 622
 623   double br_sum_true=br_sum;
 624
 625   if (br_sum<0.000001) br_sum=1.0;
 626
 627   for(ijetset=0;ijetset<njetsetdecays;ijetset++){
 628     if (jetsetdecays[ijetset]->getParentId()==ipar){
 629
 630       double br=jetsetdecays[ijetset]->getBranchingFraction();
 631
 632       int i,daugs[5];
 633       EvtId cdaugs[5];
 634
 635       for(i=0;i<5;i++){
 636
 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]);
 641         }
 642         else{
 643           daugs[i]=0;
 644         }
 645       }
 646
 647       int channel;
 648
 649       channel=EvtDecayTable::findChannel(EvtPDL::chargeConj(ipar),
 650                              jetsetdecays[ijetset]->getModelName(),
 651                              jetsetdecays[ijetset]->getNDaug(),
 652                              cdaugs,
 653                              jetsetdecays[ijetset]->getNArg(),
 654                              jetsetdecays[ijetset]->getArgsStr());
 655
 656       if (jetsetdecays[ijetset]->getModelName()=="PYTHIA"){
 657
 658         if (first) {
 659           first=0;
 660           WritePythiaEntryHeader(outdec,
 661                                  //EvtPDL::getLundKC(iparname),
 662                                  EvtPDL::getStdHep(iparname),
 663                                  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);
 670         }
 671
 672         int dflag=2;
 673
 674         if (EvtPDL::getStdHep(ipar)<0) {
 675           dflag=3;
 676           for(i=0;i<jetsetdecays[ijetset]->getNDaug();i++){
 677             daugs[i]=EvtPDL::getStdHep(cdaugs[i]);
 678           }
 679
 680         }
 681
 682         /* RS
 683           PYTHIA allows to introduce new particles via a call to PYUPDA
 684           so no need for this check any more
 685
 686           //now lets check to make sure that jetset, lucomp, knows
 687           //about all particles!
 688           int unknown=0;
 689           for(i=0;i<jetsetdecays[ijetset]->getNDaug();i++){
 690           if (pycomp_(&daugs[i])==0) {
 691           unknown=1;
 692           report(ERROR,"EvtGen") << "Pythia (pycomp) does not "
 693           << "know the particle:"<<
 694           EvtPDL::name(jetsetdecays[ijetset]->getDaugs()[i])<<endl;
 695           }
 696           }
 697
 698           int istdheppar=EvtPDL::getStdHep(ipar);
 699
 700           if (pycomp_(&istdheppar)==0){
 701           unknown=1;
 702           report(ERROR,"EvtGen") << "Pythia (pycomp) does not "
 703           << "know the particle:"<<
 704           EvtPDL::name(ipar)<<endl;
 705           }
 706
 707
 708
 709           if (unknown){
 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])<<" ";
 714           }
 715           report(ERROR,"")<<endl;
 716           report(ERROR,"EvtGen")<<"Will not be generated."<<endl;
 717           return;
 718           }
 719           */
 720
 721         if (EvtPDL::chargeConj(ipar)==ipar) {
 722           dflag=1;
 723           //report(INFO,"EvtGen") << EvtPDL::name(iparname) << " dflag=1 because C(ipar)=ipar!"<<endl;
 724         }
 725
 726
 727         //if (channel>=0) {
 728         //  dflag=1;
 729         //report(INFO,"EvtGen") << EvtPDL::name(iparname) << " dflag=1 because channel>=0"<<endl;
 730         //}
 731
 732         //      if (!(EvtPDL::getStdHep(ipar)<0&&channel>=0)){
 733         if (1){
 734
 735           // RS changed format to new PYTHIA one
 736           outdec << "          ";
 737           outdec.width(5);
 738           outdec <<dflag;
 739           outdec.width(5);
 740           outdec <<(int)jetsetdecays[ijetset]->getArgs()[0];
 741           outdec.width(12);
 742           if (fabs(br)<0.000000001) {
 743             outdec <<"0.00000";
 744           }
 745           else{
 746             outdec <<br/br_sum;
 747           }
 748           outdec.width(10);
 749           outdec <<daugs[0];
 750           outdec.width(10);
 751           outdec <<daugs[1];
 752           outdec.width(10);
 753           outdec <<daugs[2];
 754           outdec.width(10);
 755           outdec <<daugs[3];
 756           outdec.width(10);
 757           outdec <<daugs[4];
 758           outdec<<endl;
 759           outdec.width(0);
 760         }
 761       }
 762     }
 763   }
 764 }
 765
 766 bool
 767 EvtPythia::diquark(int ID)
 768 {
 769   switch(ID)
 770     {
 771     case 1103:
 772     case 2101:
 773     case 2103:
 774     case 2203:
 775     case 3101:
 776     case 3103:
 777     case 3201:
 778     case 3203:
 779     case 3303:
 780     case 4101:
 781     case 4103:
 782     case 4201:
 783     case 4203:
 784     case 4301:
 785     case 4303:
 786     case 4403:
 787     case 5101:
 788     case 5103:
 789     case 5201:
 790     case 5203:
 791     case 5301:
 792     case 5303:
 793     case 5401:
 794     case 5403:
 795     case 5503:
 796       return true;
 797       break;
 798     default:
 799       return false;
 800       break;
 801     }
 802 }
 803
 804 double
 805 EvtPythia::NominalMass(int ID)
 806 {
 807   // return default mass in PYTHIA
 808   switch(ID)
 809     {
 810     case 1103:
 811       return 0.77133;
 812     case 2101:
 813       return 0.57933;
 814     case 2103:
 815       return 0.77133;
 816     case 2203:
 817       return 0.77133;
 818     case 3101:
 819       return 0.80473;
 820     case 3103:
 821       return 0.92953;
 822     case 3201:
 823       return 0.80473;
 824     case 3203:
 825       return 0.92953;
 826     case 3303:
 827       return 1.09361;
 828     case 4101:
 829       return 1.96908;
 830     case 4103:
 831       return 2.00808;
 832     case 4201:
 833       return 1.96908;
 834     case 4203:
 835       return 2.00808;
 836     case 4301:
 837       return 2.15432;
 838     case 4303:
 839       return 2.17967;
 840     case 4403:
 841       return 3.27531;
 842     case 5101:
 843       return 5.38897;
 844     case 5103:
 845       return 5.40145;
 846     case 5201:
 847       return 5.38897;
 848     case 5203:
 849       return 5.40145;
 850     case 5301:
 851       return 5.56725;
 852     case 5303:
 853       return 5.57536;
 854     case 5401:
 855       return 6.67143;
 856     case 5403:
 857       return 6.67397;
 858     case 5503:
 859       return 10.07354;
 860       break;
 861     default:
 862       return 0.0;
 863       break;
 864     }
 865 }
 866
 867 int
 868 NominalCharge(int ID)
 869 {
 870   // return default mass in PYTHIA
 871   switch(ID)
 872     {
 873     case 1103:
 874       return -2;
 875     case 2101:
 876       return  1;
 877     case 2103:
 878       return  1;
 879     case 2203:
 880       return  4;
 881     case 3101:
 882       return -2;
 883     case 3103:
 884       return -2;
 885     case 3201:
 886       return  1;
 887     case 3203:
 888       return  1;
 889     case 3303:
 890       return -2;
 891     case 4101:
 892       return  1;
 893     case 4103:
 894       return  1;
 895     case 4201:
 896       return  4;
 897     case 4203:
 898       return  4;
 899     case 4301:
 900       return  1;
 901     case 4303:
 902       return  1;
 903     case 4403:
 904       return  4;
 905     case 5101:
 906       return -2;
 907     case 5103:
 908       return -2;
 909     case 5201:
 910       return  1;
 911     case 5203:
 912       return  1;
 913     case 5301:
 914       return -2;
 915     case 5303:
 916       return -2;
 917     case 5401:
 918       return  1;
 919     case 5403:
 920       return  1;
 921     case 5503:
 922       return -2;
 923       break;
 924     default:
 925       return 0;
 926       break;
 927     }
 928 }
 929
 930 void EvtPythia::MakePythiaFile(char* fname){
 931
 932   EvtId ipar;
 933   int lundkc;
 934
 935   //int part_list[MAX_PART];
 936
 937   ofstream outdec;
 938
 939   outdec.open(fname);
 940
 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;
 946
 947   int nokcentry;
 948
 949   for(lundkc=1;lundkc<500;lundkc++){
 950
 951     nokcentry=1;
 952
 953     for(size_t iipar=0;iipar<EvtPDL::entries();iipar++){
 954
 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;
 960
 961       if(!(
 962                EvtPDL::getStdHep(ipar)==21 ||
 963                EvtPDL::getStdHep(ipar)==22 ||
 964                EvtPDL::getStdHep(ipar)==23))
 965         {
 966
 967           if (lundkc==EvtPDL::getLundKC(ipar)){
 968
 969             nokcentry=0;
 970
 971             int first=1;
 972
 973             WritePythiaParticle(outdec,ipar,ipar,first);
 974
 975
 976             EvtId ipar2=EvtPDL::chargeConj(ipar);
 977
 978
 979             if (ipar2!=ipar){
 980               WritePythiaParticle(outdec,ipar2,ipar,first);
 981             }
 982
 983             if (first){
 984               WritePythiaEntryHeader(outdec,
 985                                      //EvtPDL::getLundKC(ipar),
 986                                      EvtPDL::getStdHep(ipar),
 987                                      ipar,
 988                                      EvtPDL::name(ipar),
 989                                      EvtPDL::chg3(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);
 994
 995             }
 996           }
 997         }
 998     }
 999     if (lundkc==99999) // Write out diquarks after quarks, but only once
1000       for(size_t iipar=0;iipar<EvtPDL::entries();iipar++){
1001
1002         ipar=EvtId(iipar,iipar);
1003
1004         if (diquark(EvtPDL::getStdHep(ipar))){
1005
1006           nokcentry=0;
1007
1008           int first=1;
1009
1010           WritePythiaParticle(outdec,ipar,ipar,first);
1011
1012
1013           EvtId ipar2=EvtPDL::chargeConj(ipar);
1014
1015
1016           if (ipar2!=ipar){
1017             WritePythiaParticle(outdec,ipar2,ipar,first);
1018           }
1019
1020           if (first){
1021             WritePythiaEntryHeader(outdec,
1022                                    EvtPDL::getStdHep(ipar),
1023                                    ipar,
1024                                    EvtPDL::name(ipar),
1025                                    NominalCharge(EvtPDL::getStdHep(ipar)),
1026                                    -1,0,
1027                                    NominalMass(EvtPDL::getStdHep(ipar)),
1028                                    0, 0, 0, 0,0.0);
1029
1030           }
1031         }
1032       }
1033     /* if (nokcentry){
1034
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);
1039
1040        } */
1041   }
1042   outdec.close();
1043 }
1044
1045 void EvtPythia::pythiaInit(int dummy){
1046
1047   //static int first=1;
1048   static int first=0; //if first=0 Pythia is not reinitialize
1049   if (first){
1050
1051     first=0;
1052
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()));
1056
1057     char fname[200];
1058
1059     char hostBuffer[100];
1060
1061     if ( gethostname( hostBuffer, 100 ) != 0 ){
1062       report(ERROR,"EvtGen") << " couldn't get hostname." << endl;
1063       strncpy( hostBuffer, "hostnameNotFound", 100 );
1064     }
1065
1066     char pid[100];
1067
1068     int thePid=getpid();
1069
1070     if ( sprintf( pid, "%d", thePid ) == 0 ){
1071       report(ERROR,"EvtGen") << " couldn't get process ID." << endl;
1072       strncpy( pid, "666", 100 );
1073     }
1074
1075     strcpy(fname,"jet.d-");
1076     strcat(fname,hostBuffer);
1077     strcat(fname,"-");
1078     strcat(fname,pid);
1079
1080     MakePythiaFile(fname);
1081     evtpythiainit_(fname,strlen(fname));
1082     initpythia_(&dummy);
1083
1084     if (0==getenv("EVTSAVEJETD")){
1085       char delcmd[300];
1086       strcpy(delcmd,"rm -f ");
1087       strcat(delcmd,fname);
1088       system(delcmd);
1089     }
1090
1091     report(INFO,"EvtGen") << "Done initializing Pythia."<<endl;
1092
1093   }
1094
1095 }