]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TEvtGen/EvtGenModels/EvtPythia.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Removing the flat makefiles
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtPythia.cxx
CommitLineData
da0e9ce3 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>
37using std::endl;
38using std::fstream;
39using std::ios;
40using std::ofstream;
41using std::resetiosflags;
42using std::setiosflags;
43using std::setw;
44
45using std::string;
46
47int EvtPythia::njetsetdecays=0;
48 EvtDecayBasePtr* EvtPythia::jetsetdecays=0;
49int EvtPythia::ntable=0;
50
51int EvtPythia::ncommand=0;
52int EvtPythia::lcommand=0;
53std::string* EvtPythia::commands=0;
54
55
56extern "C" {
57 extern void pycontinuum_(double *,int *, int *,
58 double *,double *,double *,double *);
59}
60
61extern "C" {
62 extern void evtpythiainit_(const char* fname, int len);
63}
64
65extern "C" {
66 extern void init_cont_();
67}
68
69extern "C" {
70 extern void pythiadec_(int *,double *,int *,int *,int *,
71 double *,double *,double *,double *);
72}
73
74extern "C" {
75 extern void initpythia_(int *);
76}
77
78extern "C" {
79 extern void pygive_(const char *cnfgstr,int length);
80}
81
82extern "C" {
83 extern int pycomp_(int* kf);
84}
85
86extern "C" {
87 extern void pylist_(int &);
88}
89
90
91EvtPythia::EvtPythia(){}
92
93EvtPythia::~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
124std::string EvtPythia::getName(){
125
126 return "PYTHIA";
127
128}
129
130EvtDecayBase* EvtPythia::clone(){
131
132 return new EvtPythia;
133
134}
135
136
137void EvtPythia::initProbMax(){
138
139 noProbMax();
140
141}
142
143
144void 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
167std::string EvtPythia::commandName(){
168
169 return std::string("JetSetPar");
170
171}
172
173
174void 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
200void 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
208void 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
398void 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
440void 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
461void 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
602void 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
766bool
767EvtPythia::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
804double
805EvtPythia::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
867int
868NominalCharge(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
930void 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
1045void 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}