]>
Commit | Line | Data |
---|---|---|
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 EvtGen/COPYRIGHT | |
9 | // Copyright (C) 1998 Caltech, UCSB | |
10 | // | |
11 | // Module: EvtPycont.cc | |
12 | // | |
13 | // Description: Routine to generate e+e- --> q\barq via Jetset | |
14 | // | |
15 | // Modification history: | |
16 | // | |
17 | // PCK August 4, 1997 Module created | |
18 | // RS October 28, 2002 copied from EvtJscont.cc | |
19 | // | |
20 | //------------------------------------------------------------------------ | |
21 | // | |
22 | #include "EvtGenBase/EvtPatches.hh" | |
23 | #include <stdlib.h> | |
24 | #include "EvtGenBase/EvtParticle.hh" | |
25 | #include "EvtGenBase/EvtDecayTable.hh" | |
26 | #include "EvtGenBase/EvtPDL.hh" | |
27 | #include "EvtGenModels/EvtPycont.hh" | |
28 | #include "EvtGenModels/EvtPythia.hh" | |
29 | #include "EvtGenBase/EvtId.hh" | |
30 | #include "EvtGenBase/EvtReport.hh" | |
31 | #include <string> | |
32 | #include <iostream> | |
33 | using std::endl; | |
34 | ||
35 | extern "C" { | |
36 | extern void pystat_(int &); | |
c2aad3ae | 37 | extern struct |
38 | { | |
39 | int dc[18]; | |
40 | } decaych_; | |
da0e9ce3 | 41 | } |
42 | ||
da0e9ce3 | 43 | |
44 | ||
45 | EvtPycont::~EvtPycont() { | |
46 | // int i=1; | |
47 | // pystat_(i); | |
48 | } | |
49 | ||
50 | std::string EvtPycont::getName() | |
51 | { | |
52 | return "PYCONT"; | |
53 | } | |
54 | ||
55 | EvtDecayBase* EvtPycont::clone() | |
56 | { | |
57 | return new EvtPycont; | |
58 | } | |
59 | ||
60 | void EvtPycont::init() | |
61 | { | |
62 | // check that there are 1 argument | |
63 | if ( getNArg() != 12 && getNArg() != 0 ) { | |
64 | report(ERROR,"EvtGen") << "EvtPYCONT expects " | |
65 | << " 12 arguments (d u s c b t e nu_e mu nu_mu tau nu_tau) but found: " | |
66 | << getNArg() <<endl; | |
67 | ||
68 | } | |
69 | checkNArg(0,12); | |
70 | ||
c2aad3ae | 71 | for( int i=0; i<18; i++) |
72 | decaych_.dc[i]=0; | |
73 | if ( getNArg() == 12 ) { | |
74 | decaych_.dc[0]=(int)getArg(0); | |
75 | decaych_.dc[1]=(int)getArg(1); | |
76 | decaych_.dc[2]=(int)getArg(2); | |
77 | decaych_.dc[3]=(int)getArg(3); | |
78 | decaych_.dc[4]=(int)getArg(4); | |
79 | decaych_.dc[5]=(int)getArg(5); | |
80 | decaych_.dc[10]=(int)getArg(6); | |
81 | decaych_.dc[11]=(int)getArg(7); | |
82 | decaych_.dc[12]=(int)getArg(8); | |
83 | decaych_.dc[13]=(int)getArg(9); | |
84 | decaych_.dc[14]=(int)getArg(10); | |
85 | decaych_.dc[15]=(int)getArg(11); | |
86 | } | |
87 | else{ | |
88 | decaych_.dc[0]=1; | |
89 | decaych_.dc[1]=1; | |
90 | decaych_.dc[2]=1; | |
91 | decaych_.dc[3]=1; | |
92 | } | |
da0e9ce3 | 93 | |
94 | } | |
95 | ||
96 | void EvtPycont::initProbMax() | |
97 | { | |
98 | noProbMax(); | |
99 | } | |
100 | ||
101 | void EvtPycont::decay( EvtParticle *p) | |
102 | { | |
103 | EvtPythia::pythiaInit(0); | |
104 | EvtVector4R p4[100]; | |
105 | ||
106 | double energy=p->mass(); | |
107 | ||
108 | int i,more; | |
109 | int ndaugjs; | |
110 | int kf[100]; | |
111 | EvtId id[100]; | |
112 | int type[MAX_DAUG]; | |
113 | ||
114 | double px[100],py[100],pz[100],e[100]; | |
115 | ||
116 | if ( p->getNDaug() != 0 ) { return;} | |
117 | do{ | |
118 | EvtPythia::pythiacont(&energy,&ndaugjs,kf,px,py,pz,e); | |
119 | ||
120 | for(i=0;i<ndaugjs;i++) | |
121 | { | |
122 | ||
123 | id[i]=EvtPDL::evtIdFromStdHep(kf[i]); | |
124 | ||
125 | type[i]=EvtPDL::getSpinType(id[i]); | |
126 | ||
127 | // have to protect against negative mass^2 for massless particles | |
128 | // i.e. neutrinos and photons. | |
129 | // this is uggly but I need to fix it right now.... | |
130 | ||
131 | if (px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i]>=e[i]*e[i]) | |
132 | e[i]=sqrt(px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i])+0.0000000000001; | |
133 | ||
134 | p4[i].set(e[i],px[i],py[i],pz[i]); | |
135 | ||
136 | } | |
137 | ||
138 | int channel=EvtDecayTable::inChannelList(p->getId(),ndaugjs,id); | |
139 | ||
140 | more=((channel!=-1)&&(channel!=p->getChannel())); | |
141 | ||
142 | }while(more); | |
143 | ||
144 | p->makeDaughters(ndaugjs,id); | |
145 | ||
146 | for(i=0;i<ndaugjs;i++) | |
147 | p->getDaug(i)->init( id[i], p4[i] ); | |
148 | ||
149 | return ; | |
150 | } | |
151 |