183 {
184
185 static int iniflag=0;
186
188
190
191
192
193
194
195
196
197
198
199
200
201
202 if(istdheppar != 443 && istdheppar != 100443 && istdheppar != 10441 &&istdheppar != 20443 &&istdheppar != 445 &&istdheppar != 10443 &&istdheppar != 441 &&istdheppar!= 30443){
203 std::cout<<"EvtGen: EvtLundCharm cann't not decay the particle pid= "<<istdheppar<<endl;
204 ::abort();
205 }
206
209 double totEn=0;
210
211
213
214 int i,more, pflag;;
216 int ndaugjs;
217 static int kf[100];
218 EvtId evtnumstable[100],evtnumparton[100];
219 int stableindex[100],partonindex[100];
220 int numstable;
221 int numparton;
222 static int km[100];
224
225 static double px[100],py[100],pz[100],e[100];
226 static int myflag;
227 if (iniflag==0)
lundcrm_(&iniflag,&istdheppar,&xmp,&ndaugjs,kf,km,px,py,pz,e, &myflag);
229
231
235 do{
236
237 iniflag=iniflag+1;
238 lundcrm_(&iniflag,&istdheppar,&xmp,&ndaugjs,kf,km,px,py,pz,e, &myflag);
239
240
242
243 numstable=0;
244 numparton=0;
245
246 totEn=0;
247 double parityf=1;
248 for(i=0;i<ndaugjs;i++){
249
250 totEn +=e[i];
253
255 report(
ERROR,
"EvtGen") <<
"LundCharm returned particle:"<<kf[i]<<endl;
256 report(
ERROR,
"EvtGen") <<
"This can not be translated to evt number"<<endl;
257 report(
ERROR,
"EvtGen") <<
"and the decay will be rejected!"<<endl;
258 report(
ERROR,
"EvtGen") <<
"The decay was of particle:"<<ip<<endl;
259
260 }
261
262
263 if (
abs(kf[i])<=6||kf[i]==21){
264 partonindex[numparton]=i;
266 numparton++;
267 }
268 else{
269 stableindex[numstable]=i;
271 numstable++;
272 }
273
274
275
276
277
278
279 if (px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i]>=e[i]*e[i]){
280
281 e[i]=sqrt(px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i])+0.0000000001;
282
283 }
284
285 p4[i].
set(e[i],px[i],py[i],pz[i]);
286
287
288 }
289
291
292
293
294
295
296
297
298
299
300
301
302 if(parityi!=0 && parityf!=0){
303 pflag=(parityi!=parityf);
304 }else{pflag=2;}
305
306 bool eck = fabs(xmp-totEn)>0.01;
307
308 more=(channel!=-1 || pflag ==1 || eck);
309
310
311
312
313
314
315
317
318 }while( more && (count<10000) );
319
320
321
322
323
324
325
326
327 if (count>9999) {
328 report(
INFO,
"EvtGen") <<
"Too many loops in EvtLundCharm!!!"<<endl;
330 for(i=0;i<numstable;i++){
332 }
333
334 }
335
336
337
338 if (numparton==0){
339
341 int ndaugFound=0;
342 for(i=0;i<numstable;i++){
343 p->
getDaug(i)->
init(evtnumstable[i],p4[stableindex[i]]);
344 ndaugFound++;
345 }
346 if ( ndaugFound == 0 ) {
347 report(
ERROR,
"EvtGen") <<
"Lundcharm has failed to do a decay ";
349 assert(0);
350 }
351
352 fixPolarizations(p);
353
354 return ;
355
356 }
357 else{
358
359
360
362
363 for(i=0;i<numparton;i++){
364 p4string+=p4[partonindex[i]];
365 }
366
367 int nprimary=1;
368 type[0]=STRNG;
369 for(i=0;i<numstable;i++){
370 if (km[stableindex[i]]==0){
371 type[nprimary++]=evtnumstable[i];
372 }
373 }
374
376
378
380
381 for(i=0;i<numparton;i++){
382 p4partons[i]=p4[partonindex[i]];
383 }
384
386
387
388
389 nprimary=1;
390
391 for(i=0;i<numstable;i++){
392
393 if (km[stableindex[i]]==0){
394 p->
getDaug(nprimary++)->
init(evtnumstable[i],p4[stableindex[i]]);
395 }
396 }
397
398
399 int nsecond=0;
400 for(i=0;i<numstable;i++){
401 if (km[stableindex[i]]!=0){
402 type[nsecond++]=evtnumstable[i];
403 }
404 }
405
406
408
409 EvtVector4R p4stringboost(p4string.get(0),-p4string.get(1),
410 -p4string.get(2),-p4string.get(3));
411
412 nsecond=0;
413 for(i=0;i<numstable;i++){
414 if (km[stableindex[i]]!=0){
415 p4[stableindex[i]]=
boostTo(p4[stableindex[i]],p4stringboost);
419 nsecond++;
420 }
421 }
422
423 if ( nsecond == 0 ) {
424 report(
ERROR,
"EvtGen") <<
"Jetset has failed to do a decay ";
426 assert(0);
427 }
428
429 fixPolarizations(p);
430
431 return ;
432
433 }
434
435}
double abs(const EvtComplex &c)
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
void lundcrm_(int *, int *, float *, int *, int *, int *, double *, double *, double *, double *, int *)
static int inChannelList(EvtId parent, int ndaug, EvtId *daugs)
static void LundcrmInit(int f)
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 setGeneratorFlag(int flag)
void setDiagonalSpinDensity()
EvtParticle * getDaug(int i)
void deleteDaughters(bool keepChannel=false)
void set(int i, double d)
static double getC(string parname)