208 {
209
210
212
214
218 return;
219 }
220
221
223
225
226 int i,more;
228
229 int ndaugjs;
230 int kf[100];
231 EvtId evtnumstable[100],evtnumparton[100];
232 int stableindex[100],partonindex[100];
233 int numstable;
234 int numparton;
235 int km[100];
237
239
240 double px[100],py[100],pz[100],e[100];
242
244
245 do{
246
248
249
250 numstable=0;
251 numparton=0;
252
253 for(i=0;i<ndaugjs;i++){
254
256 report(
ERROR,
"EvtGen") <<
"Pythia returned particle:"<<kf[i]<<endl;
257 report(
ERROR,
"EvtGen") <<
"This can not be translated to evt number"<<endl;
258 report(
ERROR,
"EvtGen") <<
"and the decay will be rejected!"<<endl;
259 report(
ERROR,
"EvtGen") <<
"The decay was of particle:"<<ip<<endl;
260 int i=1;
262 }
263
264
265 if (
abs(kf[i])<=6||kf[i]==21){
266 partonindex[numparton]=i;
268 numparton++;
269 }
270 else{
271 stableindex[numstable]=i;
273 numstable++;
274 }
275
276
277
278
279
280
281 if (px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i]>=e[i]*e[i]){
282
283 e[i]=sqrt(px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i])+0.0000000000001;
284
285 }
286
287 p4[i].
set(e[i],px[i],py[i],pz[i]);
288
289
290 }
291
293
294 more=(channel!=-1);
295
296
298
299 }while( more && (count<10000) );
300
301
302 if (count>9999) {
303 report(
INFO,
"EvtGen") <<
"Too many loops in EvtPythia!!!"<<endl;
305 for(i=0;i<numstable;i++){
307 }
308
309 }
310
311
312
313 if (numparton==0){
314
316
317 for(i=0;i<numstable;i++){
318 p->
getDaug(i)->
init(evtnumstable[i],p4[stableindex[i]]);
319 }
320
321 fixPolarizations(p);
322
323 return ;
324
325 }
326 else{
327
328
329
331
332 for(i=0;i<numparton;i++){
333 p4string+=p4[partonindex[i]];
334 }
335
336 int nprimary=1;
337 type[0]=STRNG;
338 for(i=0;i<numstable;i++){
339 if (km[stableindex[i]]==0){
340 type[nprimary++]=evtnumstable[i];
341 }
342 }
343
345
347
349
350 for(i=0;i<numparton;i++){
351 p4partons[i]=p4[partonindex[i]];
352 }
353
355
356
357
358 nprimary=1;
359
360 for(i=0;i<numstable;i++){
361
362 if (km[stableindex[i]]==0){
363 p->
getDaug(nprimary++)->
init(evtnumstable[i],p4[stableindex[i]]);
364 }
365 }
366
367
368 int nsecond=0;
369 for(i=0;i<numstable;i++){
370 if (km[stableindex[i]]!=0){
371 type[nsecond++]=evtnumstable[i];
372 }
373 }
374
375
377
378 nsecond=0;
379 for(i=0;i<numstable;i++){
380 if (km[stableindex[i]]!=0){
381 p4[stableindex[i]]=
boostTo(p4[stableindex[i]],p4string);
385 nsecond++;
386 }
387 }
388
389 fixPolarizations(p);
390
391 return ;
392
393 }
394
395}
double abs(const EvtComplex &c)
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
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)