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 double totalM=0;
236 do{
237
238 iniflag=iniflag+1;
239 lundcrm_(&iniflag,&istdheppar,&xmp,&ndaugjs,kf,km,px,py,pz,e, &myflag);
240
241
243
244 numstable=0;
245 numparton=0;
246
247 totEn=0;
248 double parityf=1;
249 for(i=0;i<ndaugjs;i++){
250
251 totEn +=e[i];
256 report(
ERROR,
"EvtGen") <<
"LundCharm 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
261 }
262
263
264 if (
abs(kf[i])<=6||kf[i]==21){
265 partonindex[numparton]=i;
267 numparton++;
268 }
269 else{
270 stableindex[numstable]=i;
272 numstable++;
273 }
274
275
276
277
278
279
280 if (px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i]>=e[i]*e[i]){
281
282 e[i]=sqrt(px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i])+0.0000000001;
283
284 }
285
286 p4[i].set(e[i],px[i],py[i],pz[i]);
287
288
289 }
290
292
293
294
295
296
297
298
299
300
301
302
303 if(parityi!=0 && parityf!=0){
304 pflag=(parityi!=parityf);
305 }else{pflag=2;}
306
307 bool eck = fabs(xmp-totEn)>0.01;
308
309 more=(channel!=-1 || pflag ==1 || eck );
310
311
312
313
314
315
316
318
319 }
while( more && (
count<10000) );
320
321
322
323
324
325
326
327
329 report(
INFO,
"EvtGen") <<
"Too many loops in EvtLundCharm!!!"<<endl;
331 for(i=0;i<numstable;i++){
333 }
334
335 }
336
337
338
339 if (numparton==0){
340
342 int ndaugFound=0;
343 for(i=0;i<numstable;i++){
345 ndaugFound++;
346 }
347 if ( ndaugFound == 0 ) {
348 report(
ERROR,
"EvtGen") <<
"Lundcharm has failed to do a decay ";
350 assert(0);
351 }
352
353 fixPolarizations(p);
354
355 return ;
356
357 }
358 else{
359
360
361
363
364 for(i=0;i<numparton;i++){
365 p4string+=
p4[partonindex[i]];
366 }
367
368 int nprimary=1;
369 type[0]=STRNG;
370 for(i=0;i<numstable;i++){
371 if (km[stableindex[i]]==0){
372 type[nprimary++]=evtnumstable[i];
373 }
374 }
375
377
379
381
382 for(i=0;i<numparton;i++){
383 p4partons[i]=
p4[partonindex[i]];
384 }
385
387
388
389
390 nprimary=1;
391
392 for(i=0;i<numstable;i++){
393
394 if (km[stableindex[i]]==0){
395 p->
getDaug(nprimary++)->
init(evtnumstable[i],
p4[stableindex[i]]);
396 }
397 }
398
399
400 int nsecond=0;
401 for(i=0;i<numstable;i++){
402 if (km[stableindex[i]]!=0){
403 type[nsecond++]=evtnumstable[i];
404 }
405 }
406
407
409
410 EvtVector4R p4stringboost(p4string.get(0),-p4string.get(1),
411 -p4string.get(2),-p4string.get(3));
412
413 nsecond=0;
414 for(i=0;i<numstable;i++){
415 if (km[stableindex[i]]!=0){
416 p4[stableindex[i]]=
boostTo(
p4[stableindex[i]],p4stringboost);
420 nsecond++;
421 }
422 }
423
424 if ( nsecond == 0 ) {
425 report(
ERROR,
"EvtGen") <<
"Jetset has failed to do a decay ";
427 assert(0);
428 }
429
430 fixPolarizations(p);
431
432 return ;
433
434 }
435
436}
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
void lundcrm_(int *, int *, float *, int *, int *, int *, double *, double *, double *, double *, int *)
DOUBLE_PRECISION count[3]
static int inChannelList(EvtId parent, int ndaug, EvtId *daugs)
static void LundcrmInit(int f)
static int getStdHep(EvtId id)
static double getMeanMass(EvtId i)
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)
static double getC(string parname)
double double double * p4