00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "EvtGenBase/EvtPatches.hh"
00022 #include "EvtGenBase/EvtPatches.hh"
00023 #include "EvtGenBase/EvtParticle.hh"
00024 #include "EvtGenBase/EvtStringParticle.hh"
00025 #include "EvtGenBase/EvtDecayTable.hh"
00026 #include "EvtGenBase/EvtPDL.hh"
00027 #include "EvtGenModels/EvtPythia.hh"
00028 #include "EvtGenBase/EvtReport.hh"
00029 #include "EvtGenBase/EvtId.hh"
00030 #include <iostream>
00031 #include <iomanip>
00032 #include <fstream>
00033 #include <string.h>
00034 #include <stdlib.h>
00035 #include <unistd.h>
00036 #include <stdio.h>
00037 using std::endl;
00038 using std::fstream;
00039 using std::ios;
00040 using std::ofstream;
00041 using std::resetiosflags;
00042 using std::setiosflags;
00043 using std::setw;
00044
00045 using std::string;
00046
00047 int EvtPythia::njetsetdecays=0;
00048 EvtDecayBasePtr* EvtPythia::jetsetdecays=0;
00049 int EvtPythia::ntable=0;
00050
00051 int EvtPythia::ncommand=0;
00052 int EvtPythia::lcommand=0;
00053 std::string* EvtPythia::commands=0;
00054
00055 extern "C" {
00056 extern void pycontinuum_(double *,int *, int *,
00057 double *,double *,double *,double *);
00058 }
00059
00060 extern "C" {
00061 extern void evtpythiainit_(const char* fname, int len);
00062 }
00063
00064 extern "C" {
00065 extern void init_cont_();
00066 }
00067
00068 extern "C" {
00069 extern void pythiadec_(int *,double *,int *,int *,int *,
00070 double *,double *,double *,double *);
00071 }
00072
00073 extern "C" {
00074 extern void initpythia_(int *);
00075 }
00076
00077 extern "C" {
00078 extern void pygive_(const char *cnfgstr,int length);
00079 }
00080
00081 extern "C" {
00082 extern int pycomp_(int* kf);
00083 }
00084
00085 extern "C" {
00086 extern void pylist_(int &);
00087 }
00088
00089
00090 extern "C" struct {
00091 double maximum;
00092 } cbbeam_;
00093
00094 EvtPythia::EvtPythia(){}
00095
00096 EvtPythia::~EvtPythia(){
00097
00098
00099 int i;
00100
00101
00102
00103
00104 if (njetsetdecays==0) {
00105 delete [] commands;
00106 commands=0;
00107 return;
00108 }
00109
00110 for(i=0;i<njetsetdecays;i++){
00111 if (jetsetdecays[i]==this){
00112 jetsetdecays[i]=jetsetdecays[njetsetdecays-1];
00113 njetsetdecays--;
00114 if (njetsetdecays==0) {
00115 delete [] commands;
00116 commands=0;
00117 }
00118 return;
00119 }
00120 }
00121
00122 report(ERROR,"EvtGen") << "Error in destroying Pythia model!"<<endl;
00123
00124 }
00125
00126
00127 void EvtPythia::getName(std::string& model_name){
00128
00129 model_name="PYTHIA";
00130
00131 }
00132
00133 EvtDecayBase* EvtPythia::clone(){
00134
00135 return new EvtPythia;
00136
00137 }
00138
00139
00140 void EvtPythia::initProbMax(){
00141
00142 noProbMax();
00143
00144 }
00145
00146
00147 void EvtPythia::init(){
00148
00149 checkNArg(1);
00150
00151
00152 if (getParentId().isAlias()){
00153
00154 report(ERROR,"EvtGen") << "EvtPythia finds that you are decaying the"<<endl
00155 << " aliased particle "
00156 << EvtPDL::name(getParentId()).c_str()
00157 << " with the Pythia model"<<endl
00158 << " this does not work, please modify decay table."
00159 << endl;
00160 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00161 ::abort();
00162
00163 }
00164
00165 store(this);
00166
00167 }
00168
00169
00170 std::string EvtPythia::commandName(){
00171
00172 return std::string("JetSetPar");
00173
00174 }
00175
00176
00177 void EvtPythia::command(std::string cmd){
00178
00179 if (ncommand==lcommand){
00180
00181 lcommand=10+2*lcommand;
00182
00183 std::string* newcommands=new std::string[lcommand];
00184
00185 int i;
00186
00187 for(i=0;i<ncommand;i++){
00188 newcommands[i]=commands[i];
00189 }
00190
00191 delete [] commands;
00192
00193 commands=newcommands;
00194
00195 }
00196
00197 commands[ncommand]=cmd;
00198
00199 ncommand++;
00200
00201 }
00202
00203 void EvtPythia::pythiacont(double *energy, int *ndaugjs, int *kf,
00204 double *px, double *py, double *pz, double *e)
00205 {
00206 pycontinuum_(energy,ndaugjs,kf,px,py,pz,e);
00207 }
00208
00209
00210
00211 void EvtPythia::decay( EvtParticle *p){
00212
00213
00214 static EvtId STRNG=EvtPDL::getId("string");
00215
00216 int istdheppar=EvtPDL::getStdHep(p->getId());
00217
00218 if (pycomp_(&istdheppar)==0){
00219 report(ERROR,"EvtGen") << "Pythia can not decay:"
00220 <<EvtPDL::name(p->getId()).c_str()<<endl;
00221 return;
00222 }
00223
00224
00225 double mp=p->mass();
00226
00227 EvtVector4R p4[20];
00228
00229 int i,more;
00230 int ip=EvtPDL::getStdHep(p->getId());
00231
00232 int ndaugjs;
00233 int kf[100];
00234 EvtId evtnumstable[100],evtnumparton[100];
00235 int stableindex[100],partonindex[100];
00236 int numstable;
00237 int numparton;
00238 int km[100];
00239 EvtId type[MAX_DAUG];
00240
00241 cbbeam_.maximum = mp;
00242 if(mp==0){std::cout<<"Particle "<<EvtPDL::name(p->getId())<<" has zero mass"<<std::endl;abort();}
00243 pythiaInit(0);
00244
00245 double px[100],py[100],pz[100],e[100];
00246 if ( p->getNDaug() != 0 ) { p->deleteDaughters(true);}
00247
00248 int count=0;
00249
00250 do{
00251
00252 pythiadec_(&ip,&mp,&ndaugjs,kf,km,px,py,pz,e);
00253
00254
00255 numstable=0;
00256 numparton=0;
00257
00258 for(i=0;i<ndaugjs;i++){
00259
00260 if (EvtPDL::evtIdFromStdHep(kf[i])==EvtId(-1,-1)) {
00261 report(ERROR,"EvtGen") << "Pythia returned particle:"<<kf[i]<<endl;
00262 report(ERROR,"EvtGen") << "This can not be translated to evt number"<<endl;
00263 report(ERROR,"EvtGen") << "and the decay will be rejected!"<<endl;
00264 report(ERROR,"EvtGen") << "The decay was of particle:"<<ip<<endl;
00265 int i=1;
00266 pylist_(i);
00267 }
00268
00269
00270 if (abs(kf[i])<=6||kf[i]==21){
00271 partonindex[numparton]=i;
00272 evtnumparton[numparton]=EvtPDL::evtIdFromStdHep(kf[i]);
00273 numparton++;
00274 }
00275 else{
00276 stableindex[numstable]=i;
00277 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(kf[i]);
00278 numstable++;
00279 }
00280
00281
00282
00283
00284
00285
00286 if (px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i]>=e[i]*e[i]){
00287
00288 e[i]=sqrt(px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i])+0.0000000000001;
00289
00290 }
00291
00292 p4[i].set(e[i],px[i],py[i],pz[i]);
00293
00294
00295 }
00296
00297 int channel= EvtDecayTable::inChannelList(p->getId(),numstable,evtnumstable);
00298
00299 more=(channel!=-1);
00300
00301
00302 count++;
00303
00304 }while( more && (count<10000) );
00305
00306
00307 if (count>9999) {
00308 report(INFO,"EvtGen") << "Too many loops in EvtPythia!!!"<<endl;
00309 report(INFO,"EvtGen") << "Parent:"<<EvtPDL::name(getParentId()).c_str()<<endl;
00310 for(i=0;i<numstable;i++){
00311 report(INFO,"EvtGen") << "Daug("<<i<<")"<<EvtPDL::name(evtnumstable[i]).c_str()<<endl;
00312 }
00313
00314 }
00315
00316
00317
00318 if (numparton==0){
00319
00320 p->makeDaughters(numstable,evtnumstable);
00321
00322 for(i=0;i<numstable;i++){
00323 p->getDaug(i)->init(evtnumstable[i],p4[stableindex[i]]);
00324 }
00325
00326 fixPolarizations(p);
00327
00328 return ;
00329
00330 }
00331 else{
00332
00333
00334
00335 EvtVector4R p4string(0.0,0.0,0.0,0.0);
00336
00337 for(i=0;i<numparton;i++){
00338 p4string+=p4[partonindex[i]];
00339 }
00340
00341 int nprimary=1;
00342 type[0]=STRNG;
00343 for(i=0;i<numstable;i++){
00344 if (km[stableindex[i]]==0){
00345 type[nprimary++]=evtnumstable[i];
00346 }
00347 }
00348
00349 p->makeDaughters(nprimary,type);
00350
00351 p->getDaug(0)->init(STRNG,p4string);
00352
00353 EvtVector4R p4partons[10];
00354
00355 for(i=0;i<numparton;i++){
00356 p4partons[i]=p4[partonindex[i]];
00357 }
00358
00359 ((EvtStringParticle*)p->getDaug(0))->initPartons(numparton,p4partons,evtnumparton);
00360
00361
00362
00363 nprimary=1;
00364
00365 for(i=0;i<numstable;i++){
00366
00367 if (km[stableindex[i]]==0){
00368 p->getDaug(nprimary++)->init(evtnumstable[i],p4[stableindex[i]]);
00369 }
00370 }
00371
00372
00373 int nsecond=0;
00374 for(i=0;i<numstable;i++){
00375 if (km[stableindex[i]]!=0){
00376 type[nsecond++]=evtnumstable[i];
00377 }
00378 }
00379
00380
00381 p->getDaug(0)->makeDaughters(nsecond,type);
00382
00383 nsecond=0;
00384 for(i=0;i<numstable;i++){
00385 if (km[stableindex[i]]!=0){
00386 p4[stableindex[i]]=boostTo(p4[stableindex[i]],p4string);
00387 p->getDaug(0)->getDaug(nsecond)->init(evtnumstable[i],p4[stableindex[i]]);
00388 p->getDaug(0)->getDaug(nsecond)->setDiagonalSpinDensity();
00389 p->getDaug(0)->getDaug(nsecond)->decay();
00390 nsecond++;
00391 }
00392 }
00393
00394 fixPolarizations(p);
00395
00396 return ;
00397
00398 }
00399
00400 }
00401
00402 void EvtPythia::fixPolarizations(EvtParticle *p){
00403
00404
00405
00406 int ndaug=p->getNDaug();
00407
00408 int i;
00409
00410 static EvtId Jpsi=EvtPDL::getId("J/psi");
00411
00412 for(i=0;i<ndaug;i++){
00413 if(p->getDaug(i)->getId()==Jpsi){
00414
00415 EvtSpinDensity rho;
00416
00417 rho.SetDim(3);
00418 rho.Set(0,0,0.5);
00419 rho.Set(0,1,0.0);
00420 rho.Set(0,2,0.0);
00421
00422 rho.Set(1,0,0.0);
00423 rho.Set(1,1,1.0);
00424 rho.Set(1,2,0.0);
00425
00426 rho.Set(2,0,0.0);
00427 rho.Set(2,1,0.0);
00428 rho.Set(2,2,0.5);
00429
00430 EvtVector4R p4Psi=p->getDaug(i)->getP4();
00431
00432 double alpha=atan2(p4Psi.get(2),p4Psi.get(1));
00433 double beta=acos(p4Psi.get(3)/p4Psi.d3mag());
00434
00435
00436 p->getDaug(i)->setSpinDensityForwardHelicityBasis(rho,alpha,beta,0.0);
00437 setDaughterSpinDensity(i);
00438
00439 }
00440 }
00441
00442 }
00443
00444 void EvtPythia::store(EvtDecayBase* jsdecay){
00445
00446 if (njetsetdecays==ntable){
00447
00448 EvtDecayBasePtr* newjetsetdecays=new EvtDecayBasePtr[2*ntable+10];
00449 int i;
00450 for(i=0;i<ntable;i++){
00451 newjetsetdecays[i]=jetsetdecays[i];
00452 }
00453 ntable=2*ntable+10;
00454 delete [] jetsetdecays;
00455 jetsetdecays=newjetsetdecays;
00456 }
00457
00458 jetsetdecays[njetsetdecays++]=jsdecay;
00459
00460
00461
00462 }
00463
00464
00465 void EvtPythia::WritePythiaEntryHeader(ofstream &outdec, int lundkc,
00466 EvtId evtnum,std::string name,
00467 int chg, int cchg, int spin2,double mass,
00468 double width, double maxwidth,double ctau,
00469 int stable,double rawbrfrsum){
00470
00471 char sname[100];
00472 char ccsname[100];
00473
00474
00475 int namelength=16;
00476
00477 int i,j;
00478 int temp;
00479 temp = spin2;
00480
00481 if (ctau>1000000.0) ctau=0.0;
00482
00483 strcpy(sname,name.c_str());
00484
00485 i=0;
00486
00487 while (sname[i]!=0){
00488 i++;
00489 }
00490
00491
00492
00493 if(evtnum.getId()>=0) {
00494 if (sname[i-1]=='+'||sname[i-1]=='-'){
00495 sname[i-1]=0;
00496 i--;
00497 }
00498 if (sname[i-1]=='+'||sname[i-1]=='-'){
00499 sname[i-1]=0;
00500 i--;
00501 }
00502
00503 if (sname[i-1]=='0' && sname[i-2]!='_' && !(sname[0]=='c' && sname[1]=='h')){
00504 sname[i-1]=0;
00505 i--;
00506 }
00507 }
00508
00509 if (i>namelength) {
00510 for(j=1;j<namelength;j++){
00511 sname[j]=sname[j+i-namelength];
00512 }
00513 sname[namelength]=0;
00514 }
00515
00516
00517 for(j=0;j<=namelength;j++)
00518 ccsname[j]=sname[j];
00519 i=0;
00520 while (ccsname[i]!=' '){
00521 i++;
00522 if(ccsname[i]==0) break;
00523 }
00524 if(i<namelength)
00525 {
00526 ccsname[i]='b';
00527 ccsname[i+1]=0;
00528 }
00529
00530
00531
00532 if(evtnum.getId()>=0) {
00533 if (abs(EvtPDL::getStdHep(evtnum))==21) cchg=2;
00534 if (abs(EvtPDL::getStdHep(evtnum))==90) cchg=-1;
00535 if ((abs(EvtPDL::getStdHep(evtnum))<=8)&&
00536 (abs(EvtPDL::getStdHep(evtnum))!=0)) cchg=1;
00537
00538 }
00539
00540
00541 outdec << " " << setw(9) << lundkc;
00542 outdec << " ";
00543 outdec.width(namelength);
00544 outdec << setiosflags(ios::left)
00545 << sname;
00546
00547 if ((evtnum.getId()>=0) && (EvtPDL::chargeConj(evtnum)!=evtnum))
00548 {
00549 outdec << " ";
00550 outdec.width(namelength);
00551 outdec << ccsname;
00552 }else{
00553
00554 outdec << " ";
00555 }
00556
00557 outdec << resetiosflags(ios::left);
00558 outdec << setw(3) << chg;
00559 outdec << setw(3) << cchg;
00560 outdec.width(3);
00561 if (evtnum.getId()>=0) {
00562 if (EvtPDL::chargeConj(evtnum)==evtnum) {
00563 outdec << 0;
00564 }
00565 else{
00566 outdec << 1;
00567 }
00568 }
00569 else{
00570 outdec << 0;
00571 }
00572 outdec.setf(ios::fixed,ios::floatfield);
00573 outdec.precision(5);
00574 outdec << setw(12) << mass;
00575 outdec.setf(ios::fixed,ios::floatfield);
00576 outdec << setw(12) << width;
00577 outdec.setf(ios::fixed,ios::floatfield);
00578 outdec.width(12);
00579 if (fabs(width)<0.0000000001) {
00580 outdec << 0.0 ;
00581 }
00582 else{
00583 outdec << maxwidth;
00584 }
00585
00586 outdec.setf(ios::scientific,ios::floatfield);
00587 outdec << " ";
00588 outdec << ctau;
00589 outdec.width(3);
00590 if (evtnum.getId()>=0) {
00591 if (ctau>1.0 || rawbrfrsum<0.000001) {
00592 stable=0;
00593 }
00594 }
00595
00596 outdec.width(3);
00597 outdec << 0;
00598 outdec.width(3);
00599 outdec << stable;
00600 outdec << endl;
00601 outdec.width(0);
00602
00603
00604 }
00605
00606 void EvtPythia::WritePythiaParticle(ofstream &outdec,EvtId ipar,
00607 EvtId iparname,int &first){
00608
00609 int ijetset;
00610
00611 double br_sum=0.0;
00612
00613 for(ijetset=0;ijetset<njetsetdecays;ijetset++){
00614
00615 if (jetsetdecays[ijetset]->getParentId()==ipar){
00616 br_sum+=jetsetdecays[ijetset]->getBranchingFraction();
00617 }
00618 if (jetsetdecays[ijetset]->getParentId()!=
00619 EvtPDL::chargeConj(jetsetdecays[ijetset]->getParentId())&&
00620 EvtPDL::chargeConj(jetsetdecays[ijetset]->getParentId())==ipar){
00621 br_sum+=jetsetdecays[ijetset]->getBranchingFraction();
00622 }
00623
00624
00625 }
00626
00627 double br_sum_true=br_sum;
00628
00629 if (br_sum<0.000001) br_sum=1.0;
00630
00631 for(ijetset=0;ijetset<njetsetdecays;ijetset++){
00632 if (jetsetdecays[ijetset]->getParentId()==ipar){
00633
00634 double br=jetsetdecays[ijetset]->getBranchingFraction();
00635
00636 int i,daugs[5];
00637 EvtId cdaugs[5];
00638
00639 for(i=0;i<5;i++){
00640
00641 if(i<jetsetdecays[ijetset]->getNDaug()){
00642 daugs[i]=EvtPDL::getStdHep(
00643 jetsetdecays[ijetset]->getDaugs()[i]);
00644 cdaugs[i]=EvtPDL::chargeConj(jetsetdecays[ijetset]->getDaugs()[i]);
00645 }
00646 else{
00647 daugs[i]=0;
00648 }
00649 }
00650
00651 int channel;
00652
00653 channel=EvtDecayTable::findChannel(EvtPDL::chargeConj(ipar),
00654 jetsetdecays[ijetset]->getModelName(),
00655 jetsetdecays[ijetset]->getNDaug(),
00656 cdaugs,
00657 jetsetdecays[ijetset]->getNArg(),
00658 jetsetdecays[ijetset]->getArgsStr());
00659
00660 if (jetsetdecays[ijetset]->getModelName()=="PYTHIA"){
00661
00662 if (first) {
00663 first=0;
00664 WritePythiaEntryHeader(outdec,
00665
00666 EvtPDL::getStdHep(iparname),
00667 iparname,
00668 EvtPDL::name(iparname),
00669 EvtPDL::chg3(iparname),
00670 0,0,EvtPDL::getMeanMass(ipar),
00671 EvtPDL::getWidth(ipar),
00672 EvtPDL::getMeanMass(ipar)-EvtPDL::getMinMass(ipar),
00673 EvtPDL::getctau(ipar),1,br_sum_true);
00674 }
00675
00676 int dflag=2;
00677
00678 if (EvtPDL::getStdHep(ipar)<0) {
00679 dflag=3;
00680 for(i=0;i<jetsetdecays[ijetset]->getNDaug();i++){
00681 daugs[i]=EvtPDL::getStdHep(cdaugs[i]);
00682 }
00683
00684 }
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725 if (EvtPDL::chargeConj(ipar)==ipar) {
00726 dflag=1;
00727
00728 }
00729
00730
00731
00732
00733
00734
00735
00736
00737 if (1){
00738
00739
00740 outdec << " ";
00741 outdec.width(5);
00742 outdec <<dflag;
00743 outdec.width(5);
00744 outdec <<(int)jetsetdecays[ijetset]->getArgs()[0];
00745 outdec.width(12);
00746 if (fabs(br)<0.000000001) {
00747 outdec <<"0.00000";
00748 }
00749 else{
00750 outdec <<br/br_sum;
00751 }
00752 outdec.width(10);
00753 outdec <<daugs[0];
00754 outdec.width(10);
00755 outdec <<daugs[1];
00756 outdec.width(10);
00757 outdec <<daugs[2];
00758 outdec.width(10);
00759 outdec <<daugs[3];
00760 outdec.width(10);
00761 outdec <<daugs[4];
00762 outdec<<endl;
00763 outdec.width(0);
00764 }
00765 }
00766 }
00767 }
00768 }
00769
00770 bool
00771 EvtPythia::diquark(int ID)
00772 {
00773 switch(ID)
00774 {
00775 case 1103:
00776 case 2101:
00777 case 2103:
00778 case 2203:
00779 case 3101:
00780 case 3103:
00781 case 3201:
00782 case 3203:
00783 case 3303:
00784 case 4101:
00785 case 4103:
00786 case 4201:
00787 case 4203:
00788 case 4301:
00789 case 4303:
00790 case 4403:
00791 case 5101:
00792 case 5103:
00793 case 5201:
00794 case 5203:
00795 case 5301:
00796 case 5303:
00797 case 5401:
00798 case 5403:
00799 case 5503:
00800 return true;
00801 break;
00802 default:
00803 return false;
00804 break;
00805 }
00806 }
00807
00808 double
00809 EvtPythia::NominalMass(int ID)
00810 {
00811
00812 switch(ID)
00813 {
00814 case 1103:
00815 return 0.77133;
00816 case 2101:
00817 return 0.57933;
00818 case 2103:
00819 return 0.77133;
00820 case 2203:
00821 return 0.77133;
00822 case 3101:
00823 return 0.80473;
00824 case 3103:
00825 return 0.92953;
00826 case 3201:
00827 return 0.80473;
00828 case 3203:
00829 return 0.92953;
00830 case 3303:
00831 return 1.09361;
00832 case 4101:
00833 return 1.96908;
00834 case 4103:
00835 return 2.00808;
00836 case 4201:
00837 return 1.96908;
00838 case 4203:
00839 return 2.00808;
00840 case 4301:
00841 return 2.15432;
00842 case 4303:
00843 return 2.17967;
00844 case 4403:
00845 return 3.27531;
00846 case 5101:
00847 return 5.38897;
00848 case 5103:
00849 return 5.40145;
00850 case 5201:
00851 return 5.38897;
00852 case 5203:
00853 return 5.40145;
00854 case 5301:
00855 return 5.56725;
00856 case 5303:
00857 return 5.57536;
00858 case 5401:
00859 return 6.67143;
00860 case 5403:
00861 return 6.67397;
00862 case 5503:
00863 return 10.07354;
00864 break;
00865 default:
00866 return 0.0;
00867 break;
00868 }
00869 }
00870
00871 int
00872 NominalCharge(int ID)
00873 {
00874
00875 switch(ID)
00876 {
00877 case 1103:
00878 return -2;
00879 case 2101:
00880 return 1;
00881 case 2103:
00882 return 1;
00883 case 2203:
00884 return 4;
00885 case 3101:
00886 return -2;
00887 case 3103:
00888 return -2;
00889 case 3201:
00890 return 1;
00891 case 3203:
00892 return 1;
00893 case 3303:
00894 return -2;
00895 case 4101:
00896 return 1;
00897 case 4103:
00898 return 1;
00899 case 4201:
00900 return 4;
00901 case 4203:
00902 return 4;
00903 case 4301:
00904 return 1;
00905 case 4303:
00906 return 1;
00907 case 4403:
00908 return 4;
00909 case 5101:
00910 return -2;
00911 case 5103:
00912 return -2;
00913 case 5201:
00914 return 1;
00915 case 5203:
00916 return 1;
00917 case 5301:
00918 return -2;
00919 case 5303:
00920 return -2;
00921 case 5401:
00922 return 1;
00923 case 5403:
00924 return 1;
00925 case 5503:
00926 return -2;
00927 break;
00928 default:
00929 return 0;
00930 break;
00931 }
00932 }
00933
00934 void EvtPythia::MakePythiaFile(char* fname){
00935
00936 EvtId ipar;
00937 int lundkc;
00938
00939
00940
00941 ofstream outdec;
00942
00943 outdec.open(fname);
00944
00945
00946
00947
00948
00949
00950
00951 int nokcentry;
00952
00953 for(lundkc=1;lundkc<500;lundkc++){
00954
00955 nokcentry=1;
00956
00957 int iipar;
00958
00959 for(iipar=0;iipar<EvtPDL::entries();iipar++){
00960
00961 ipar=EvtId(iipar,iipar);
00962
00963 std::string tempStr = EvtPDL::name(ipar);
00964 EvtId realId = EvtPDL::getId(tempStr);
00965 if ( realId.isAlias() != 0 ) continue;
00966
00967 if(!(
00968 EvtPDL::getStdHep(ipar)==21 ||
00969 EvtPDL::getStdHep(ipar)==22 ||
00970 EvtPDL::getStdHep(ipar)==23))
00971 {
00972
00973 if (lundkc==EvtPDL::getLundKC(ipar)){
00974
00975 nokcentry=0;
00976
00977 int first=1;
00978
00979 WritePythiaParticle(outdec,ipar,ipar,first);
00980
00981
00982 EvtId ipar2=EvtPDL::chargeConj(ipar);
00983
00984
00985 if (ipar2!=ipar){
00986 WritePythiaParticle(outdec,ipar2,ipar,first);
00987 }
00988
00989 if (first){
00990 WritePythiaEntryHeader(outdec,
00991
00992 EvtPDL::getStdHep(ipar),
00993 ipar,
00994 EvtPDL::name(ipar),
00995 EvtPDL::chg3(ipar),
00996 0,0,EvtPDL::getMeanMass(ipar),
00997 EvtPDL::getWidth(ipar),
00998 EvtPDL::getMeanMass(ipar)-EvtPDL::getMinMass(ipar),
00999 EvtPDL::getctau(ipar),0,0.0);
01000
01001 }
01002 }
01003 }
01004 }
01005 if (lundkc==99999)
01006 for(iipar=0;iipar<EvtPDL::entries();iipar++){
01007
01008 ipar=EvtId(iipar,iipar);
01009
01010 if (diquark(EvtPDL::getStdHep(ipar))){
01011
01012 nokcentry=0;
01013
01014 int first=1;
01015
01016 WritePythiaParticle(outdec,ipar,ipar,first);
01017
01018
01019 EvtId ipar2=EvtPDL::chargeConj(ipar);
01020
01021
01022 if (ipar2!=ipar){
01023 WritePythiaParticle(outdec,ipar2,ipar,first);
01024 }
01025
01026 if (first){
01027 WritePythiaEntryHeader(outdec,
01028 EvtPDL::getStdHep(ipar),
01029 ipar,
01030 EvtPDL::name(ipar),
01031 NominalCharge(EvtPDL::getStdHep(ipar)),
01032 -1,0,
01033 NominalMass(EvtPDL::getStdHep(ipar)),
01034 0, 0, 0, 0,0.0);
01035
01036 }
01037 }
01038 }
01039
01040
01041
01042
01043
01044
01045
01046
01047 }
01048 outdec.close();
01049 }
01050
01051 void EvtPythia::pythiaInit(int dummy){
01052
01053 static int first=1;
01054 if (first){
01055
01056 first=0;
01057
01058 report(INFO,"EvtGen") << "Will initialize Pythia."<<endl;
01059 for(int i=0;i<ncommand;i++)
01060 pygive_(commands[i].c_str(),strlen(commands[i].c_str()));
01061
01062 char fname[200];
01063
01064 char hostBuffer[100];
01065
01066 if ( gethostname( hostBuffer, 100 ) != 0 ){
01067 report(ERROR,"EvtGen") << " couldn't get hostname." << endl;
01068 strncpy( hostBuffer, "hostnameNotFound", 100 );
01069 }
01070
01071 char pid[100];
01072
01073 int thePid=getpid();
01074
01075 if ( sprintf( pid, "%d", thePid ) == 0 ){
01076 report(ERROR,"EvtGen") << " couldn't get process ID." << endl;
01077 strncpy( pid, "666", 100 );
01078 }
01079
01080 strcpy(fname,"jet.d-");
01081 strcat(fname,hostBuffer);
01082 strcat(fname,"-");
01083 strcat(fname,pid);
01084
01085 MakePythiaFile(fname);
01086 evtpythiainit_(fname,strlen(fname));
01087 initpythia_(&dummy);
01088
01089 if (0==getenv("EVTSAVEJETD")){
01090 char delcmd[300];
01091 strcpy(delcmd,"rm -f ");
01092 strcat(delcmd,fname);
01093 system(delcmd);
01094 }
01095
01096 report(INFO,"EvtGen") << "Done initializing Pythia."<<endl;
01097
01098 }
01099
01100 }