TEvtGen/EvtGenModels/EvtJetSetCDF.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 EvtGen/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 //
  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/EvtJetSetCDF.hh"
  28 #include "EvtGenBase/EvtReport.hh"
  29 #include <string>
  30 #include "EvtGenBase/EvtId.hh"
  31 #include <iostream>
  32 #include <iomanip>
  33 #include <fstream>
  34 #include <string.h>
  35 #include <stdlib.h>
  36 #include <unistd.h>
  37 #include <stdio.h>
  38 using std::endl;
  39 using std::fstream;
  40 using std::ios;
  41 using std::ofstream;
  42 using std::resetiosflags;
  43 using std::setiosflags;
  44 using std::setw;
  45
  46
  47 int EvtJetSetCDF::njetsetdecays=0;
  48   EvtDecayBasePtr* EvtJetSetCDF::jetsetdecays=0;
  49 int EvtJetSetCDF::ntable=0;
  50
  51 int EvtJetSetCDF::ncommand=0;
  52 int EvtJetSetCDF::lcommand=0;
  53 std::string* EvtJetSetCDF::commands=0;
  54
  55 extern "C" {
  56   extern void evtjetsetcdfinit_(char* fname, int len);
  57 }
  58
  59
  60 extern "C" {
  61   extern void jetsetcdf_(int *,double *,int *,int *,int *,
  62                        double *,double *,double *,double *);
  63 }
  64
  65 extern "C" {
  66   extern void lygive_(const char *cnfgstr,int length);
  67 }
  68
  69 extern "C" {
  70   extern int lycomp_(int* kf);
  71 }
  72
  73
  74 EvtJetSetCDF::EvtJetSetCDF(){}
  75
  76 EvtJetSetCDF::~EvtJetSetCDF(){
  77
  78
  79   int i;
  80
  81
  82   //the deletion of commands is really uggly!
  83
  84   if (njetsetdecays==0) {
  85     delete [] commands;
  86     commands=0;
  87     return;
  88   }
  89
  90   for(i=0;i<njetsetdecays;i++){
  91     if (jetsetdecays[i]==this){
  92       jetsetdecays[i]=jetsetdecays[njetsetdecays-1];
  93       njetsetdecays--;
  94       if (njetsetdecays==0) {
  95         delete [] commands;
  96         commands=0;
  97       }
  98       return;
  99     }
 100   }
 101
 102   report(ERROR,"EvtGen") << "Error in destroying JetSet model!"<<endl;
 103
 104 }
 105
 106
 107 std::string EvtJetSetCDF::getName(){
 108
 109   return "JETSETCDF";
 110
 111 }
 112
 113 EvtDecayBase* EvtJetSetCDF::clone(){
 114
 115   return new EvtJetSetCDF;
 116
 117 }
 118
 119
 120 void EvtJetSetCDF::initProbMax(){
 121
 122   noProbMax();
 123
 124 }
 125
 126
 127 void EvtJetSetCDF::init(){
 128
 129   checkNArg(1);
 130
 131
 132   if (getParentId().isAlias()){
 133
 134     report(ERROR,"EvtGen") << "EvtJetSetCDF 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."
 139                            << endl;
 140     report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
 141     ::abort();
 142
 143   }
 144
 145   store(this);
 146
 147 }
 148
 149
 150 std::string EvtJetSetCDF::commandName(){
 151
 152   return std::string("JetSetCDFPar");
 153
 154 }
 155
 156
 157 void EvtJetSetCDF::command(std::string cmd){
 158
 159   if (ncommand==lcommand){
 160
 161     lcommand=10+2*lcommand;
 162
 163     std::string* newcommands=new std::string[lcommand];
 164
 165     int i;
 166
 167     for(i=0;i<ncommand;i++){
 168       newcommands[i]=commands[i];
 169     }
 170
 171     delete [] commands;
 172
 173     commands=newcommands;
 174
 175   }
 176
 177   commands[ncommand]=cmd;
 178
 179   ncommand++;
 180
 181 }
 182
 183
 184
 185 void EvtJetSetCDF::decay( EvtParticle *p){
 186
 187
 188   //added by Lange Jan4,2000
 189   static EvtId STRNG=EvtPDL::getId("string");
 190
 191   int istdheppar=EvtPDL::getStdHep(p->getId());
 192
 193   if (lycomp_(&istdheppar)==0){
 194     report(ERROR,"EvtGen") << "Jetset can not decay:"
 195       <<EvtPDL::name(p->getId()).c_str()<<endl;
 196     return;
 197   }
 198
 199   double mp=p->mass();
 200
 201   EvtVector4R p4[20];
 202
 203   int i,more;
 204   int ip=EvtPDL::getStdHep(p->getId());
 205   int ndaugjs;
 206   int kf[100];
 207   EvtId evtnumstable[100],evtnumparton[100];
 208   int stableindex[100],partonindex[100];
 209   int numstable;
 210   int numparton;
 211   int km[100];
 212   EvtId type[MAX_DAUG];
 213
 214   jetSetInit();
 215
 216   double px[100],py[100],pz[100],e[100];
 217
 218   if ( p->getNDaug() != 0 ) { p->deleteDaughters(true);}
 219
 220   int count=0;
 221
 222   do{
 223     //report(INFO,"EvtGen") << "calling jetset " << ip<< " " << mp <<endl;
 224     jetsetcdf_(&ip,&mp,&ndaugjs,kf,km,px,py,pz,e);
 225
 226
 227     numstable=0;
 228     numparton=0;
 229     //report(INFO,"EvtGen") << "found some daughters " << ndaugjs << endl;
 230     for(i=0;i<ndaugjs;i++){
 231
 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;
 237
 238       }
 239
 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]);
 244         numparton++;
 245       }
 246       else{
 247         stableindex[numstable]=i;
 248         evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(kf[i]);
 249         numstable++;
 250       }
 251
 252
 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....
 256
 257       if (px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i]>=e[i]*e[i]){
 258
 259         e[i]=sqrt(px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i])+0.0000000000001;
 260
 261       }
 262
 263       p4[i].set(e[i],px[i],py[i],pz[i]);
 264
 265
 266     }
 267
 268     int channel=EvtDecayTable::inChannelList(p->getId(),numstable,evtnumstable);
 269
 270
 271    more=(channel!=-1);
 272
 273
 274
 275
 276   count++;
 277
 278   }while( more && (count<10000) );
 279
 280   if (count>9999) {
 281     report(INFO,"EvtGen") << "Too many loops in EvtJetSetCDF!!!"<<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;
 285     }
 286
 287   }
 288
 289
 290
 291   if (numparton==0){
 292
 293     p->makeDaughters(numstable,evtnumstable);
 294     int ndaugFound=0;
 295     for(i=0;i<numstable;i++){
 296       p->getDaug(i)->init(evtnumstable[i],p4[stableindex[i]]);
 297       ndaugFound++;
 298     }
 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;
 302       assert(0);
 303     }
 304
 305     fixPolarizations(p);
 306
 307     return ;
 308
 309   }
 310   else{
 311
 312     //have partons in JETSET
 313
 314     EvtVector4R p4string(0.0,0.0,0.0,0.0);
 315
 316     for(i=0;i<numparton;i++){
 317       p4string+=p4[partonindex[i]];
 318     }
 319
 320     int nprimary=1;
 321     type[0]=STRNG;
 322     for(i=0;i<numstable;i++){
 323       if (km[stableindex[i]]==0){
 324         type[nprimary++]=evtnumstable[i];
 325       }
 326     }
 327
 328     p->makeDaughters(nprimary,type);
 329
 330     p->getDaug(0)->init(STRNG,p4string);
 331
 332     EvtVector4R p4partons[10];
 333
 334     for(i=0;i<numparton;i++){
 335       p4partons[i]=p4[partonindex[i]];
 336     }
 337
 338     ((EvtStringParticle*)p->getDaug(0))->initPartons(numparton,p4partons,evtnumparton);
 339
 340
 341
 342     nprimary=1;
 343
 344     for(i=0;i<numstable;i++){
 345
 346       if (km[stableindex[i]]==0){
 347         p->getDaug(nprimary++)->init(evtnumstable[i],p4[stableindex[i]]);
 348       }
 349     }
 350
 351
 352     int nsecond=0;
 353     for(i=0;i<numstable;i++){
 354       if (km[stableindex[i]]!=0){
 355         type[nsecond++]=evtnumstable[i];
 356       }
 357     }
 358
 359
 360     p->getDaug(0)->makeDaughters(nsecond,type);
 361
 362     EvtVector4R p4stringboost(p4string.get(0),-p4string.get(1),
 363                               -p4string.get(2),-p4string.get(3));
 364
 365     nsecond=0;
 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();
 372         nsecond++;
 373       }
 374     }
 375
 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;
 379       assert(0);
 380     }
 381
 382     fixPolarizations(p);
 383
 384     return ;
 385
 386   }
 387
 388 }
 389
 390 void EvtJetSetCDF::fixPolarizations(EvtParticle *p){
 391
 392   //special case for now to handle the J/psi polarization
 393
 394   int ndaug=p->getNDaug();
 395
 396   int i;
 397
 398   static EvtId Jpsi=EvtPDL::getId("J/psi");
 399
 400   for(i=0;i<ndaug;i++){
 401     if(p->getDaug(i)->getId()==Jpsi){
 402
 403       EvtSpinDensity rho;
 404
 405       rho.setDim(3);
 406       rho.set(0,0,0.5);
 407       rho.set(0,1,0.0);
 408       rho.set(0,2,0.0);
 409
 410       rho.set(1,0,0.0);
 411       rho.set(1,1,1.0);
 412       rho.set(1,2,0.0);
 413
 414       rho.set(2,0,0.0);
 415       rho.set(2,1,0.0);
 416       rho.set(2,2,0.5);
 417
 418       EvtVector4R p4Psi=p->getDaug(i)->getP4();
 419
 420       double alpha=atan2(p4Psi.get(2),p4Psi.get(1));
 421       double beta=acos(p4Psi.get(3)/p4Psi.d3mag());
 422
 423
 424       p->getDaug(i)->setSpinDensityForwardHelicityBasis(rho,alpha,beta,0.0);
 425       setDaughterSpinDensity(i);
 426
 427     }
 428   }
 429
 430 }
 431
 432 void EvtJetSetCDF::store(EvtDecayBase* jsdecay){
 433
 434   if (njetsetdecays==ntable){
 435
 436     EvtDecayBasePtr* newjetsetdecays=new EvtDecayBasePtr[2*ntable+10];
 437     int i;
 438     for(i=0;i<ntable;i++){
 439       newjetsetdecays[i]=jetsetdecays[i];
 440     }
 441     ntable=2*ntable+10;
 442     delete [] jetsetdecays;
 443     jetsetdecays=newjetsetdecays;
 444   }
 445
 446   jetsetdecays[njetsetdecays++]=jsdecay;
 447
 448
 449
 450 }
 451
 452
 453 void EvtJetSetCDF::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){
 458
 459
 460   char sname[100];
 461
 462   int namelength=8;
 463
 464   int i,j;
 465   int temp;
 466   temp = spin2;
 467
 468   if (ctau>1000000.0) ctau=0.0;
 469
 470   strcpy(sname,name.c_str());
 471
 472   i=0;
 473
 474   while (sname[i]!=0){
 475     i++;
 476   }
 477
 478   // strip up to two + or -
 479
 480   if(evtnum.getId()>=0) {
 481     if (sname[i-1]=='+'||sname[i-1]=='-'){
 482       sname[i-1]=0;
 483       i--;
 484     }
 485     if (sname[i-1]=='+'||sname[i-1]=='-'){
 486       sname[i-1]=0;
 487       i--;
 488     }
 489     // strip 0 except for _0 and chi...0
 490     if (sname[i-1]=='0' && sname[i-2]!='_' && !(sname[0]=='c' && sname[1]=='h')){
 491       sname[i-1]=0;
 492       i--;
 493     }
 494   }
 495
 496   if (i>namelength) {
 497     for(j=1;j<namelength;j++){
 498       sname[j]=sname[j+i-namelength];
 499     }
 500     sname[namelength]=0;
 501   }
 502
 503   cchg=0;
 504
 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;
 510
 511   }
 512
 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;
 518   outdec.width(3);
 519   if (evtnum.getId()>=0) {
 520     if (EvtPDL::chargeConj(evtnum)==evtnum) {
 521       outdec << 0;
 522     }
 523     else{
 524       outdec << 1;
 525     }
 526   }
 527   else{
 528     outdec << 0;
 529   }
 530   outdec.setf(ios::fixed);
 531   outdec.precision(5);
 532   outdec << setw(12) << mass;
 533   outdec << setw(12) << width;
 534   outdec.width(12);
 535   if (fabs(width)<0.0000000001) {
 536     outdec << 0.0 ;
 537   }
 538   else{
 539     outdec << maxwidth;
 540   }
 541   outdec << setw(14) << ctau;
 542   outdec.width(3);
 543   if (evtnum.getId()>=0) {
 544     if (ctau>1.0 || rawbrfrsum<0.000001) {
 545       stable=0;
 546     }
 547   }
 548   outdec << stable;
 549   outdec << endl;
 550   outdec.width(0);
 551
 552 }
 553
 554 void EvtJetSetCDF::WriteJetSetParticle(ofstream &outdec,EvtId ipar,
 555                                     EvtId iparname,int &first){
 556
 557   int ijetset;
 558
 559   double br_sum=0.0;
 560
 561   for(ijetset=0;ijetset<njetsetdecays;ijetset++){
 562
 563     if (jetsetdecays[ijetset]->getParentId()==ipar){
 564       br_sum+=jetsetdecays[ijetset]->getBranchingFraction();
 565     }
 566     if (jetsetdecays[ijetset]->getParentId()!=
 567         EvtPDL::chargeConj(jetsetdecays[ijetset]->getParentId())&&
 568         EvtPDL::chargeConj(jetsetdecays[ijetset]->getParentId())==ipar){
 569       br_sum+=jetsetdecays[ijetset]->getBranchingFraction();
 570     }
 571
 572
 573   }
 574
 575   double br_sum_true=br_sum;
 576
 577   if (br_sum<0.000001) br_sum=1.0;
 578
 579   for(ijetset=0;ijetset<njetsetdecays;ijetset++){
 580     if (jetsetdecays[ijetset]->getParentId()==ipar){
 581
 582       double br=jetsetdecays[ijetset]->getBranchingFraction();
 583
 584       int i,daugs[5];
 585       EvtId cdaugs[5];
 586
 587       for(i=0;i<5;i++){
 588
 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]);
 593         }
 594         else{
 595           daugs[i]=0;
 596         }
 597       }
 598
 599       int channel;
 600
 601       channel=EvtDecayTable::findChannel(EvtPDL::chargeConj(ipar),
 602                              jetsetdecays[ijetset]->getModelName(),
 603                              jetsetdecays[ijetset]->getNDaug(),
 604                              cdaugs,
 605                              jetsetdecays[ijetset]->getNArg(),
 606                              jetsetdecays[ijetset]->getArgsStr());
 607
 608       if (jetsetdecays[ijetset]->getModelName()=="JETSET"){
 609
 610         if (first) {
 611           first=0;
 612           WriteJetSetEntryHeader(outdec,
 613                                  EvtPDL::getLundKC(iparname),
 614                                  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);
 621         }
 622
 623         int dflag=2;
 624
 625         if (EvtPDL::getStdHep(ipar)<0) {
 626           dflag=3;
 627           for(i=0;i<jetsetdecays[ijetset]->getNDaug();i++){
 628             daugs[i]=EvtPDL::getStdHep(cdaugs[i]);
 629           }
 630
 631         }
 632
 633         //now lets check to make sure that jetset, lycomp, knows
 634         //about all particles!
 635         int unknown=0;
 636         for(i=0;i<jetsetdecays[ijetset]->getNDaug();i++){
 637           if (lycomp_(&daugs[i])==0) {
 638             unknown=1;
 639             report(ERROR,"EvtGen") << "JetSet (lycomp) does not "
 640                                   << "know the particle:"<<
 641               EvtPDL::name(jetsetdecays[ijetset]->getDaugs()[i]).c_str()<<endl;
 642           }
 643         }
 644
 645         int istdheppar=EvtPDL::getStdHep(ipar);
 646
 647         if (lycomp_(&istdheppar)==0){
 648           unknown=1;
 649           report(ERROR,"EvtGen") << "JetSet (lycomp) does not "
 650                   << "know the particle:"<<
 651               EvtPDL::name(ipar).c_str()<<endl;
 652         }
 653
 654
 655
 656         if (unknown){
 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()<<" ";
 661           }
 662           report(ERROR,"")<<endl;
 663           report(ERROR,"EvtGen")<<"Will not be generated."<<endl;
 664           return;
 665         }
 666
 667
 668         if (EvtPDL::chargeConj(ipar)==ipar) {
 669           dflag=1;
 670           //report(INFO,"EvtGen") << EvtPDL::name(iparname) << " dflag=1 because C(ipar)=ipar!"<<endl;
 671         }
 672
 673
 674         //if (channel>=0) {
 675         //  dflag=1;
 676           //report(INFO,"EvtGen") << EvtPDL::name(iparname) << " dflag=1 because channel>=0"<<endl;
 677         //}
 678
 679         //      if (!(EvtPDL::getStdHep(ipar)<0&&channel>=0)){
 680         if (1){
 681
 682           outdec.width(10);
 683           outdec <<dflag;
 684           outdec.width(5);
 685           outdec <<(int)jetsetdecays[ijetset]->getArgs()[0];
 686           outdec.width(12);
 687           if (fabs(br)<0.000000001) {
 688             outdec <<"0.00000";
 689           }
 690           else{
 691             outdec <<br/br_sum;
 692           }
 693           outdec.width(8);
 694           outdec <<daugs[0];
 695           outdec.width(8);
 696           outdec <<daugs[1];
 697           outdec.width(8);
 698           outdec <<daugs[2];
 699           outdec.width(8);
 700           outdec <<daugs[3];
 701           outdec.width(8);
 702           outdec <<daugs[4];
 703           outdec<<endl;
 704           outdec.width(0);
 705         }
 706       }
 707     }
 708   }
 709 }
 710
 711
 712 void EvtJetSetCDF::MakeJetSetFile(char* fname){
 713
 714   EvtId ipar;
 715   int lundkc;
 716
 717   //int part_list[MAX_PART];
 718
 719   ofstream outdec;
 720
 721   outdec.open(fname);
 722
 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;
 727
 728   int nokcentry;
 729
 730   for(lundkc=1;lundkc<=500;lundkc++){
 731
 732     nokcentry=1;
 733
 734     for(size_t iipar=0;iipar<EvtPDL::entries();iipar++){
 735
 736       ipar=EvtId(iipar,iipar);
 737       //no aliased particles!
 738       std::string tempStr = EvtPDL::name(ipar);
 739       EvtId realId = EvtPDL::getId(tempStr);
 740       if ( realId.isAlias() != 0 ) continue;
 741       if (lundkc==EvtPDL::getLundKC(ipar)){
 742
 743         nokcentry=0;
 744
 745         int first=1;
 746
 747         WriteJetSetParticle(outdec,ipar,ipar,first);
 748
 749
 750         EvtId ipar2=EvtPDL::chargeConj(ipar);
 751
 752
 753         if (ipar2!=ipar){
 754           WriteJetSetParticle(outdec,ipar2,ipar,first);
 755         }
 756
 757         if (first){
 758           WriteJetSetEntryHeader(outdec,
 759                                     EvtPDL::getLundKC(ipar),
 760                                     ipar,
 761                                     EvtPDL::name(ipar),
 762                                     EvtPDL::chg3(ipar),
 763                                     0,0,EvtPDL::getMeanMass(ipar),
 764                                     EvtPDL::getWidth(ipar),
 765                                     EvtPDL::getMeanMass(ipar)-EvtPDL::getMinMass(ipar),
 766                                     EvtPDL::getctau(ipar),0,0.0);
 767
 768         }
 769       }
 770     }
 771     if (nokcentry){
 772
 773       WriteJetSetEntryHeader(outdec,
 774                                 lundkc,EvtId(-1,-1),"  ",
 775                                 0,0,0,EvtPDL::getMeanMass(ipar),0.0,0.0,
 776                                 EvtPDL::getctau(ipar),0,0.0);
 777
 778     }
 779   }
 780     outdec.close();
 781 }
 782
 783 void EvtJetSetCDF::jetSetInit(){
 784
 785   static int first=1;
 786
 787   if (first){
 788
 789     first=0;
 790
 791     report(INFO,"EvtGen") << "Will initialize JetSet."<<endl;
 792
 793     char fname[200];
 794
 795     char hostBuffer[100];
 796
 797     if ( gethostname( hostBuffer, 100 ) != 0 ){
 798       report(ERROR,"EvtGen") << " couldn't get hostname." << endl;
 799       strncpy( hostBuffer, "hostnameNotFound", 100 );
 800     }
 801
 802     char pid[100];
 803
 804     int thePid=getpid();
 805
 806     if ( sprintf( pid, "%d", thePid ) == 0 ){
 807       report(ERROR,"EvtGen") << " couldn't get process ID." << endl;
 808       strncpy( pid, "666", 100 );
 809     }
 810
 811     strcpy(fname,"jet.d-");
 812     strcat(fname,hostBuffer);
 813     strcat(fname,"-");
 814     strcat(fname,pid);
 815
 816     MakeJetSetFile(fname);
 817     evtjetsetcdfinit_(fname,strlen(fname));
 818
 819     if (0==getenv("EVTSAVEJETD")){
 820       char delcmd[300];
 821       strcpy(delcmd,"rm -f ");
 822       strcat(delcmd,fname);
 823       system(delcmd);
 824     }
 825
 826     int i;
 827
 828     for(i=0;i<ncommand;i++){
 829       lygive_(commands[i].c_str(),strlen(commands[i].c_str()));
 830
 831     }
 832
 833     report(INFO,"EvtGen") << "Done initializing JetSet."<<endl;
 834
 835     //int itest = 5122;
 836     //report(ERROR,"EvtGen TEST:") << lycomp_(&itest) << endl;
 837
 838
 839   }
 840
 841 }