152 {
153
154 static int iniflag = 0;
155
157
159
160
161
162
163
164
165
166
167
168
169
170
171
172 if ( istdheppar != 443 && istdheppar != 100443 && istdheppar != 10441 &&
173 istdheppar != 20443 && istdheppar != 445 && istdheppar != 10443 && istdheppar != 441 &&
174 istdheppar != 30443 )
175 {
176 std::cout << "EvtGen: EvtLundCharm cann't not decay the particle pid= " << istdheppar
177 << endl;
178 ::abort();
179 }
180
183 double totEn = 0;
184
185
186 EvtVector4R p4[50];
187
188 int i, more, pflag;
189 ;
191 int ndaugjs;
192 static int kf[100];
193 EvtId evtnumstable[100], evtnumparton[100];
194 int stableindex[100], partonindex[100];
195 int numstable;
196 int numparton;
197 static int km[100];
199
200 static double px[100], py[100], pz[100], e[100];
201 static int myflag;
202 if ( iniflag == 0 )
203 lundcrm_( &iniflag, &istdheppar, &xmp, &ndaugjs, kf, km, px, py, pz, e, &myflag );
205
207
211 double totalM = 0;
212 do {
213
214 iniflag = iniflag + 1;
215 lundcrm_( &iniflag, &istdheppar, &xmp, &ndaugjs, kf, km, px, py, pz, e, &myflag );
216
217
219
220 numstable = 0;
221 numparton = 0;
222
223 totEn = 0;
224 double parityf = 1;
225 for ( i = 0; i < ndaugjs; i++ )
226 {
227
228
229 totEn += e[i];
234 {
235 report(
ERROR,
"EvtGen" ) <<
"LundCharm returned particle:" << kf[i] << endl;
236 report(
ERROR,
"EvtGen" ) <<
"This can not be translated to evt number" << endl;
237 report(
ERROR,
"EvtGen" ) <<
"and the decay will be rejected!" << endl;
238 report(
ERROR,
"EvtGen" ) <<
"The decay was of particle:" << ip << endl;
239 }
240
241
242 if (
abs( kf[i] ) <= 6 || kf[i] == 21 )
243 {
244 partonindex[numparton] = i;
246 numparton++;
247 }
248 else
249 {
250 stableindex[numstable] = i;
252 numstable++;
253 }
254
255
256
257
258
259 if ( px[i] * px[i] + py[i] * py[i] + pz[i] * pz[i] >= e[i] * e[i] )
260 { e[i] = sqrt( px[i] * px[i] + py[i] * py[i] + pz[i] * pz[i] ) + 0.0000000001; }
261
262 p4[i].
set( e[i], px[i], py[i], pz[i] );
263 }
264
266
267
268
269
270
271
272
273
274
275
276
277 if ( parityi != 0 && parityf != 0 ) { pflag = ( parityi != parityf ); }
278 else { pflag = 2; }
279
280 bool eck = fabs( xmp - totEn ) > 0.01;
281
282 more = ( channel != -1 || pflag == 1 || eck );
283
284
285
286
287
288
290
291 }
while ( more && (
count < 10000 ) );
292
293
294
295
296
297
298
299
301 {
302 report(
INFO,
"EvtGen" ) <<
"Too many loops in EvtLundCharm!!!" << endl;
304 for ( i = 0; i < numstable; i++ )
305 {
306 report(
INFO,
"EvtGen" ) <<
"Daug(" << i <<
")"
308 }
309 }
310
311 if ( numparton == 0 )
312 {
313
315 int ndaugFound = 0;
316 for ( i = 0; i < numstable; i++ )
317 {
318 p->
getDaug( i )->
init( evtnumstable[i], p4[stableindex[i]] );
319 ndaugFound++;
320 }
321 if ( ndaugFound == 0 )
322 {
323 report(
ERROR,
"EvtGen" ) <<
"Lundcharm has failed to do a decay ";
325 << endl;
326 assert( 0 );
327 }
328
329 fixPolarizations( p );
330
331 return;
332 }
333 else
334 {
335
336
337
338 EvtVector4R p4string( 0.0, 0.0, 0.0, 0.0 );
339
340 for ( i = 0; i < numparton; i++ ) { p4string += p4[partonindex[i]]; }
341
342 int nprimary = 1;
343 type[0] = STRNG;
344 for ( i = 0; i < numstable; i++ )
345 {
346 if ( km[stableindex[i]] == 0 ) { type[nprimary++] = evtnumstable[i]; }
347 }
348
350
352
353 EvtVector4R p4partons[10];
354
355 for ( i = 0; i < numparton; i++ ) { p4partons[i] = p4[partonindex[i]]; }
356
357 ( (EvtStringParticle*)p->
getDaug( 0 ) )->initPartons( numparton, p4partons, evtnumparton );
358
359 nprimary = 1;
360
361 for ( i = 0; i < numstable; i++ )
362 {
363
364 if ( km[stableindex[i]] == 0 )
365 { p->
getDaug( nprimary++ )->
init( evtnumstable[i], p4[stableindex[i]] ); }
366 }
367
368 int nsecond = 0;
369 for ( i = 0; i < numstable; i++ )
370 {
371 if ( km[stableindex[i]] != 0 ) { type[nsecond++] = evtnumstable[i]; }
372 }
373
375
376 EvtVector4R p4stringboost( p4string.get( 0 ), -p4string.get( 1 ), -p4string.get( 2 ),
377 -p4string.get( 3 ) );
378
379 nsecond = 0;
380 for ( i = 0; i < numstable; i++ )
381 {
382 if ( km[stableindex[i]] != 0 )
383 {
384 p4[stableindex[i]] =
boostTo( p4[stableindex[i]], p4stringboost );
388 nsecond++;
389 }
390 }
391
392 if ( nsecond == 0 )
393 {
394 report(
ERROR,
"EvtGen" ) <<
"Jetset has failed to do a decay ";
396 << endl;
397 assert( 0 );
398 }
399
400 fixPolarizations( p );
401
402 return;
403 }
404}
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)
void set(int i, double d)
static double getC(string parname)