211 {
212
213
215
217
221 return;
222 }
223
224
226
228
229 int i,more;
231
232 int ndaugjs;
233 int kf[100];
234 EvtId evtnumstable[100],evtnumparton[100];
235 int stableindex[100],partonindex[100];
236 int numstable;
237 int numparton;
238 int km[100];
240
242 if(
mp==0){std::cout<<
"Particle "<<
EvtPDL::name(p->
getId())<<
" has zero mass"<<std::endl;abort();}
244
245 double px[100],py[100],pz[100],e[100];
247
249
250 do{
251
253
254
255 numstable=0;
256 numparton=0;
257
258 for(i=0;i<ndaugjs;i++){
259
261 report(
ERROR,
"EvtGen") <<
"Pythia returned particle:"<<kf[i]<<endl;
262 report(
ERROR,
"EvtGen") <<
"This can not be translated to evt number"<<endl;
263 report(
ERROR,
"EvtGen") <<
"and the decay will be rejected!"<<endl;
264 report(
ERROR,
"EvtGen") <<
"The decay was of particle:"<<ip<<endl;
265 int i=1;
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.0000000000001;
289
290 }
291
292 p4[i].
set(e[i],px[i],py[i],pz[i]);
293
294
295 }
296
298
299 more=(channel!=-1);
300
301
303
304 }
while( more && (
count<10000) );
305
306
308 report(
INFO,
"EvtGen") <<
"Too many loops in EvtPythia!!!"<<endl;
310 for(i=0;i<numstable;i++){
312 }
313
314 }
315
316
317
318 if (numparton==0){
319
321
322 for(i=0;i<numstable;i++){
323 p->
getDaug(i)->
init(evtnumstable[i],p4[stableindex[i]]);
324 }
325
326 fixPolarizations(p);
327
328 return ;
329
330 }
331 else{
332
333
334
336
337 for(i=0;i<numparton;i++){
338 p4string+=p4[partonindex[i]];
339 }
340
341 int nprimary=1;
342 type[0]=STRNG;
343 for(i=0;i<numstable;i++){
344 if (km[stableindex[i]]==0){
345 type[nprimary++]=evtnumstable[i];
346 }
347 }
348
350
352
354
355 for(i=0;i<numparton;i++){
356 p4partons[i]=p4[partonindex[i]];
357 }
358
360
361
362
363 nprimary=1;
364
365 for(i=0;i<numstable;i++){
366
367 if (km[stableindex[i]]==0){
368 p->
getDaug(nprimary++)->
init(evtnumstable[i],p4[stableindex[i]]);
369 }
370 }
371
372
373 int nsecond=0;
374 for(i=0;i<numstable;i++){
375 if (km[stableindex[i]]!=0){
376 type[nsecond++]=evtnumstable[i];
377 }
378 }
379
380
382
383 nsecond=0;
384 for(i=0;i<numstable;i++){
385 if (km[stableindex[i]]!=0){
386 p4[stableindex[i]]=
boostTo(p4[stableindex[i]],p4string);
390 nsecond++;
391 }
392 }
393
394 fixPolarizations(p);
395
396 return ;
397
398 }
399
400}
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
DOUBLE_PRECISION count[3]
void pythiadec_(int *, double *, int *, int *, int *, double *, double *, double *, double *)
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)
static void pythiaInit(int f)
void set(int i, double d)