199 {
200
201 static int iniflag=0;
202
204
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
225
226
228
229 int i,more;
231 int ndaugjs;
232 static int kf[100];
233 EvtId evtnumstable[100],evtnumparton[100];
234 int stableindex[100],partonindex[100];
235 int numstable;
236 int numparton;
237 static int km[100];
239
240 int isr=1;
242
243 static double px[100],py[100],pz[100],e[100];
244 if (iniflag==0)
lundar_(&isr,&iniflag,&xmp,&ndaugjs,kf,km,px,py,pz,e);
245
247
249
250 do{
251
252 iniflag=iniflag+1;
253 lundar_(&isr,&iniflag,&xmp,&ndaugjs,kf,km,px,py,pz,e);
255 numstable=0;
256 numparton=0;
257
258 for(i=0;i<ndaugjs;i++){
259
260
262 report(
ERROR,
"EvtGen") <<
"Lunda returned particle:"<<kf[i]<<endl;
263 report(
ERROR,
"EvtGen") <<
"This can not be translated to evt number"<<endl;
264 report(
ERROR,
"EvtGen") <<
"and the decay will be rejected!"<<endl;
265 report(
ERROR,
"EvtGen") <<
"The decay was of particle:"<<ip<<endl;
266
267 }
268
269
270 if (
abs(kf[i])<=6||kf[i]==21){
271 partonindex[numparton]=i;
273 numparton++;
274 }
275 else{
276 stableindex[numstable]=i;
278 numstable++;
279 }
280
281
282
283
284
285
286 if (px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i]>=e[i]*e[i]){
287
288 e[i]=sqrt(px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i])+0.0000000001;
289
290 }
291
292 p4[i].
set(e[i],px[i],py[i],pz[i]);
293
294
295 }
296
298
299
300
301
302
303
304
305
306
307
308
309 more=(channel!=-1);
310
311
312
313
315
316 }while( more && (count<10000) );
317
318 if (count>9999) {
319 report(
INFO,
"EvtGen") <<
"Too many loops in EvtLunda!!!"<<endl;
321 for(i=0;i<numstable;i++){
323 }
324
325 }
326
327
328
329 if (numparton==0){
330
332 int ndaugFound=0;
333 for(i=0;i<numstable;i++){
334 p->
getDaug(i)->
init(evtnumstable[i],p4[stableindex[i]]);
335 ndaugFound++;
336 }
337 if ( ndaugFound == 0 ) {
338 report(
ERROR,
"EvtGen") <<
"Lunda has failed to do a decay ";
340 assert(0);
341 }
342
343 fixPolarizations(p);
344
345 return ;
346
347 }
348 else{
349
350
351
353
354 for(i=0;i<numparton;i++){
355 p4string+=p4[partonindex[i]];
356 }
357
358 int nprimary=1;
359 type[0]=STRNG;
360 for(i=0;i<numstable;i++){
361 if (km[stableindex[i]]==0){
362 type[nprimary++]=evtnumstable[i];
363 }
364 }
365
367
369
371
372 for(i=0;i<numparton;i++){
373 p4partons[i]=p4[partonindex[i]];
374 }
375
377
378
379
380 nprimary=1;
381
382 for(i=0;i<numstable;i++){
383
384 if (km[stableindex[i]]==0){
385 p->
getDaug(nprimary++)->
init(evtnumstable[i],p4[stableindex[i]]);
386 }
387 }
388
389
390 int nsecond=0;
391 for(i=0;i<numstable;i++){
392 if (km[stableindex[i]]!=0){
393 type[nsecond++]=evtnumstable[i];
394 }
395 }
396
397
399
400 EvtVector4R p4stringboost(p4string.get(0),-p4string.get(1),
401 -p4string.get(2),-p4string.get(3));
402
403 nsecond=0;
404 for(i=0;i<numstable;i++){
405 if (km[stableindex[i]]!=0){
406 p4[stableindex[i]]=
boostTo(p4[stableindex[i]],p4stringboost);
410 nsecond++;
411 }
412 }
413
414 if ( nsecond == 0 ) {
415 report(
ERROR,
"EvtGen") <<
"Jetset has failed to do a decay ";
417 assert(0);
418 }
419
420 fixPolarizations(p);
421
422 return ;
423
424 }
425
426}
double abs(const EvtComplex &c)
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
void lundar_(int *, int *, float *, int *, int *, int *, double *, double *, double *, double *)
static int inChannelList(EvtId parent, int ndaug, EvtId *daugs)
void LundaInit(int dummy)
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)