151 {
152
153
155
157
158 if (
lucomp_( &istdheppar ) == 0 )
159 {
161 << endl;
162 return;
163 }
164
166
167 EvtVector4R p4[20];
168
169 int i, more;
171 int ndaugjs;
172 int kf[100];
173 EvtId evtnumstable[100], evtnumparton[100];
174 int stableindex[100], partonindex[100];
175 int numstable;
176 int numparton;
177 int km[100];
179
181
182 double px[100], py[100], pz[100], e[100];
183
185
187
188 do {
189
190 jetset1_( &ip, &
mp, &ndaugjs, kf, km, px, py, pz, e );
191
192 numstable = 0;
193 numparton = 0;
194
195 for ( i = 0; i < ndaugjs; i++ )
196 {
197
199 {
200 report(
ERROR,
"EvtGen" ) <<
"JetSet returned particle:" << kf[i] << endl;
201 report(
ERROR,
"EvtGen" ) <<
"This can not be translated to evt number" << endl;
202 report(
ERROR,
"EvtGen" ) <<
"and the decay will be rejected!" << endl;
203 report(
ERROR,
"EvtGen" ) <<
"The decay was of particle:" << ip << endl;
204 }
205
206
207 if (
abs( kf[i] ) <= 6 || kf[i] == 21 )
208 {
209 partonindex[numparton] = i;
211 numparton++;
212 }
213 else
214 {
215 stableindex[numstable] = i;
217 numstable++;
218 }
219
220
221
222
223
224 if ( px[i] * px[i] + py[i] * py[i] + pz[i] * pz[i] >= e[i] * e[i] )
225 { e[i] = sqrt( px[i] * px[i] + py[i] * py[i] + pz[i] * pz[i] ) + 0.0000000000001; }
226
227 p4[i].
set( e[i], px[i], py[i], pz[i] );
228 }
229
231
232 more = ( channel != -1 );
233
235
236 }
while ( more && (
count < 10000 ) );
237
239 {
240 report(
INFO,
"EvtGen" ) <<
"Too many loops in EvtJetSet!!!" << endl;
242 for ( i = 0; i < numstable; i++ )
243 {
244 report(
INFO,
"EvtGen" ) <<
"Daug(" << i <<
")"
246 }
247 }
248
249 if ( numparton == 0 )
250 {
251
253 int ndaugFound = 0;
254 for ( i = 0; i < numstable; i++ )
255 {
256 p->
getDaug( i )->
init( evtnumstable[i], p4[stableindex[i]] );
257 ndaugFound++;
258 }
259 if ( ndaugFound == 0 )
260 {
261 report(
ERROR,
"EvtGen" ) <<
"Jetset has failed to do a decay ";
263 << endl;
264 assert( 0 );
265 }
266
267 fixPolarizations( p );
268
269 return;
270 }
271 else
272 {
273
274
275
276 EvtVector4R p4string( 0.0, 0.0, 0.0, 0.0 );
277
278 for ( i = 0; i < numparton; i++ ) { p4string += p4[partonindex[i]]; }
279
280 int nprimary = 1;
281 type[0] = STRNG;
282 for ( i = 0; i < numstable; i++ )
283 {
284 if ( km[stableindex[i]] == 0 ) { type[nprimary++] = evtnumstable[i]; }
285 }
286
288
290
291 EvtVector4R p4partons[10];
292
293 for ( i = 0; i < numparton; i++ ) { p4partons[i] = p4[partonindex[i]]; }
294
295 ( (EvtStringParticle*)p->
getDaug( 0 ) )->initPartons( numparton, p4partons, evtnumparton );
296
297 nprimary = 1;
298
299 for ( i = 0; i < numstable; i++ )
300 {
301
302 if ( km[stableindex[i]] == 0 )
303 { p->
getDaug( nprimary++ )->
init( evtnumstable[i], p4[stableindex[i]] ); }
304 }
305
306 int nsecond = 0;
307 for ( i = 0; i < numstable; i++ )
308 {
309 if ( km[stableindex[i]] != 0 ) { type[nsecond++] = evtnumstable[i]; }
310 }
311
313
314 EvtVector4R p4stringboost( p4string.get( 0 ), -p4string.get( 1 ), -p4string.get( 2 ),
315 -p4string.get( 3 ) );
316
317 nsecond = 0;
318 for ( i = 0; i < numstable; i++ )
319 {
320 if ( km[stableindex[i]] != 0 )
321 {
322 p4[stableindex[i]] =
boostTo( p4[stableindex[i]], p4stringboost );
326 nsecond++;
327 }
328 }
329
330 if ( nsecond == 0 )
331 {
332 report(
ERROR,
"EvtGen" ) <<
"Jetset has failed to do a decay ";
334 << endl;
335 assert( 0 );
336 }
337
338 fixPolarizations( p );
339
340 return;
341 }
342}
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
void jetset1_(int *, double *, int *, int *, int *, double *, double *, double *, double *)
DOUBLE_PRECISION count[3]
static int inChannelList(EvtId parent, int ndaug, EvtId *daugs)
static int getStdHep(EvtId id)
static EvtId evtIdFromStdHep(int stdhep)
static std::string name(EvtId i)
static EvtId getId(const std::string &name)
void makeDaughters(int ndaug, EvtId *id)
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
void setDiagonalSpinDensity()
EvtParticle * getDaug(int i)
void deleteDaughters(bool keepChannel=false)
void set(int i, double d)