205{
207 static const G4double expxl = -expxu;
208
212
213 static const G4int numMul = 1200;
214 static const G4int numSec = 60;
215
221
224
225 static G4bool first =
true;
226 static G4double protmul[numMul], protnorm[numSec];
227 static G4double neutmul[numMul], neutnorm[numSec];
228
229
230
231
232 G4int i, counter, nt, npos, nneg, nzero;
233
234 if (first) {
235
236 first = false;
237 for (i = 0; i < numMul; i++) protmul[i] = 0.0;
238 for (i = 0; i < numSec; i++) protnorm[i] = 0.0;
239 counter = -1;
240 for (npos = 0; npos < (numSec/3); npos++) {
241 for( nneg=
Imax(0,npos-1); nneg<=npos+1; nneg++ )
242 {
243 for( nzero=0; nzero<numSec/3; nzero++ )
244 {
245 if( ++counter < numMul )
246 {
247 nt = npos+nneg+nzero;
248 if( (nt>0) && (nt<=numSec) )
249 {
250 protmul[counter] =
251 pmltpc(npos,nneg,nzero,nt,protb,c) ;
252 protnorm[nt-1] += protmul[counter];
253 }
254 }
255 }
256 }
257 }
258 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
259 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
260 counter = -1;
261 for (npos = 0; npos < numSec/3; npos++) {
262 for (nneg = npos; nneg <= (npos+2); nneg++) {
263 for( nzero=0; nzero<numSec/3; nzero++ )
264 {
265 if( ++counter < numMul )
266 {
267 nt = npos+nneg+nzero;
268 if( (nt>0) && (nt<=numSec) )
269 {
270 neutmul[counter] =
271 pmltpc(npos,nneg,nzero,nt,neutb,c);
272 neutnorm[nt-1] += neutmul[counter];
273 }
274 }
275 }
276 }
277 }
278
279 for (i = 0; i < numSec; i++) {
280 if (protnorm[i] > 0.0) protnorm[i] = 1.0/protnorm[i];
281 if (neutnorm[i] > 0.0) neutnorm[i] = 1.0/neutnorm[i];
282 }
283 }
284
285 pv[0] = incidentParticle;
286 pv[1] = targetParticle;
287 vecLen = 2;
288
290
291
292 npos = 0, nneg = 0, nzero = 0;
293 G4double cech[] = { 1., 1., 1., 0.70, 0.60, 0.55, 0.35, 0.25, 0.18, 0.15};
294 G4int iplab =
G4int( incidentTotalMomentum*5.);
297 G4double cnk0[] = {0.17, 0.18, 0.17, 0.24, 0.26, 0.20, 0.22, 0.21, 0.34, 0.45,
298 0.58, 0.55, 0.36, 0.29, 0.29, 0.32, 0.32, 0.33, 0.33, 0.33};
300 if (targetCode == protonCode) {
303 return;
304 } else {
305 return;
306 }
307 }
309 if (targetCode == protonCode) {
310 if( ran < 0.25 )
311 {
314 }
315 else if (ran < 0.50)
316 {
319 }
320 else if (ran < 0.75)
321 {
324 }
325 else
326 {
329 }
330 } else {
331 if( ran < 0.25 )
332 {
333 }
334 else if (ran < 0.50)
335 {
338 }
339 else if (ran < 0.75)
340 {
343 }
344 else
345 {
348 }
349 }
350 return;
351 } else {
352
353 G4double aleab = std::log(availableEnergy);
354 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
355 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
356
357
358
360
361 for (nt=1; nt<=numSec; nt++) {
362 test = std::exp(
Amin( expxu,
Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
363 dum =
pi*nt/(2.0*
n*
n);
364 if (std::fabs(dum) < 1.0) {
365 if( test >= 1.0e-10 )anpn += dum*test;
366 } else {
367 anpn += dum*test;
368 }
369 }
370
373 if( targetCode == protonCode )
374 {
375 counter = -1;
376 for( npos=0; npos<numSec/3; npos++ )
377 {
378 for( nneg=
Imax(0,npos-1); nneg<=npos+1; nneg++ )
379 {
380 for (nzero=0; nzero<numSec/3; nzero++) {
381 if (++counter < numMul) {
382 nt = npos+nneg+nzero;
383 if ( (nt>0) && (nt<=numSec) ) {
384 test = std::exp(
Amin( expxu,
Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
385 dum = (
pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
386 if (std::fabs(dum) < 1.0) {
387 if( test >= 1.0e-10 )excs += dum*test;
388 } else {
389 excs += dum*test;
390 }
391
392 if (ran < excs) goto outOfLoop;
393 }
394 }
395 }
396 }
397 }
398
399
400 inElastic = false;
401 return;
402 }
403 else
404 {
405 counter = -1;
406 for( npos=0; npos<numSec/3; npos++ )
407 {
408 for( nneg=npos; nneg<=(npos+2); nneg++ )
409 {
410 for (nzero=0; nzero<numSec/3; nzero++) {
411 if (++counter < numMul) {
412 nt = npos+nneg+nzero;
413 if ( (nt>=1) && (nt<=numSec) ) {
414 test = std::exp(
Amin( expxu,
Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
415 dum = (
pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
416 if (std::fabs(dum) < 1.0) {
417 if( test >= 1.0e-10 )excs += dum*test;
418 } else {
419 excs += dum*test;
420 }
421 if (ran < excs) goto outOfLoop;
422 }
423 }
424 }
425 }
426 }
427
428 inElastic = false;
429 return;
430 }
431 }
432 outOfLoop:
433
434 if( targetCode == protonCode)
435 {
436 if( npos == (1+nneg))
437 {
439 }
440 else if (npos == nneg)
441 {
443 {
444 }
445 else
446 {
449 }
450 }
451 else
452 {
454 }
455 }
456 else
457 {
458 if( npos == (nneg-1))
459 {
461 {
463 }
464 else
465 {
467 }
468 }
469 else if ( npos == nneg)
470 {
471 }
472 else
473 {
476 }
477 }
478
480 {
481 if( ( (pv[0].getCode() == kaonMinusCode)
483 || ( (pv[0].
getCode() == kaonZeroCode)
484 && (pv[1].getCode() == protonCode) )
485 || ( (pv[0].
getCode() == antiKaonZeroCode)
486 && (pv[1].getCode() == protonCode) ) )
487 {
489 if( pv[1].getCode() == protonCode)
490 {
491 if(ran < 0.68)
492 {
495 }
496 else if (ran < 0.84)
497 {
500 }
501 else
502 {
505 }
506 }
507 else
508 {
509 if(ran < 0.68)
510 {
513 }
514 else if (ran < 0.84)
515 {
518 }
519 else
520 {
523 }
524 }
525 }
526 else
527 {
529 if (ran < 0.67)
530 {
533 }
534 else if (ran < 0.78)
535 {
538 }
539 else if (ran < 0.89)
540 {
543 }
544 else
545 {
548 }
549 }
550 }
551
552
553 nt = npos + nneg + nzero;
554 while ( nt > 0)
555 {
558 {
559 if( npos > 0 )
561 npos--;
562 }
563 }
564 else if ( ran < (
G4double)(npos+nneg)/nt)
565 {
566 if( nneg > 0 )
567 {
569 nneg--;
570 }
571 }
572 else
573 {
574 if( nzero > 0 )
575 {
577 nzero--;
578 }
579 }
580 nt = npos + nneg + nzero;
581 }
583 {
584 G4cout <<
"Particles produced: " ;
587 for (i=2; i < vecLen; i++)
588 {
590 }
592 }
593 return;
594 }
G4double pmltpc(G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
G4double Amin(G4double a, G4double b)
G4int Imin(G4int a, G4int b)
G4double Amax(G4double a, G4double b)
G4int Imax(G4int a, G4int b)
G4double getTotalMomentum() const