213{
215 static const G4double expxl = -expxu;
216
220
221 static const G4int numMul = 1200;
222 static const G4int numMulAn = 400;
223 static const G4int numSec = 60;
224
226
229
230 static G4bool first =
true;
231 static G4double protmul[numMul], protnorm[numSec];
232 static G4double protmulAn[numMulAn],protnormAn[numSec];
233 static G4double neutmul[numMul], neutnorm[numSec];
234 static G4double neutmulAn[numMulAn],neutnormAn[numSec];
235
236
237
238
239 G4int i, counter, nt, npos, nneg, nzero;
240
241 if( first )
242 {
243 first = false;
244 for( i=0; i<numMul ; i++ ) protmul[i] = 0.0;
245 for( i=0; i<numSec ; i++ ) protnorm[i] = 0.0;
246 counter = -1;
247 for( npos=0; npos<(numSec/3); npos++ )
248 {
249 for( nneg=std::max(0,npos-2); nneg<=(npos+1); nneg++ )
250 {
251 for( nzero=0; nzero<numSec/3; nzero++ )
252 {
253 if( ++counter < numMul )
254 {
255 nt = npos+nneg+nzero;
256 if( (nt>0) && (nt<=numSec) )
257 {
258 protmul[counter] =
pmltpc(npos,nneg,nzero,nt,protb,c);
259 protnorm[nt-1] += protmul[counter];
260 }
261 }
262 }
263 }
264 }
265 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
266 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
267 counter = -1;
268 for( npos=0; npos<numSec/3; npos++ )
269 {
270 for( nneg=std::max(0,npos-1); nneg<=(npos+2); nneg++ )
271 {
272 for( nzero=0; nzero<numSec/3; nzero++ )
273 {
274 if( ++counter < numMul )
275 {
276 nt = npos+nneg+nzero;
277 if( (nt>0) && (nt<=numSec) )
278 {
279 neutmul[counter] =
pmltpc(npos,nneg,nzero,nt,neutb,c);
280 neutnorm[nt-1] += neutmul[counter];
281 }
282 }
283 }
284 }
285 }
286 for( i=0; i<numSec; i++ )
287 {
288 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
289 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
290 }
291
292 for( i=0; i<numMulAn ; i++ ) protmulAn[i] = 0.0;
293 for( i=0; i<numSec ; i++ ) protnormAn[i] = 0.0;
294 counter = -1;
295 for( npos=1; npos<(numSec/3); npos++ )
296 {
297 nneg = std::max(0,npos-1);
298 for( nzero=0; nzero<numSec/3; nzero++ )
299 {
300 if( ++counter < numMulAn )
301 {
302 nt = npos+nneg+nzero;
303 if( (nt>1) && (nt<=numSec) )
304 {
305 protmulAn[counter] =
pmltpc(npos,nneg,nzero,nt,protb,c);
306 protnormAn[nt-1] += protmulAn[counter];
307 }
308 }
309 }
310 }
311 for( i=0; i<numMulAn; i++ ) neutmulAn[i] = 0.0;
312 for( i=0; i<numSec; i++ ) neutnormAn[i] = 0.0;
313 counter = -1;
314 for( npos=0; npos<numSec/3; npos++ )
315 {
316 nneg = npos;
317 for( nzero=0; nzero<numSec/3; nzero++ )
318 {
319 if( ++counter < numMulAn )
320 {
321 nt = npos+nneg+nzero;
322 if( (nt>1) && (nt<=numSec) )
323 {
324 neutmulAn[counter] =
pmltpc(npos,nneg,nzero,nt,neutb,c);
325 neutnormAn[nt-1] += neutmulAn[counter];
326 }
327 }
328 }
329 }
330 for( i=0; i<numSec; i++ )
331 {
332 if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i];
333 if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i];
334 }
335 }
336
337
338
339
340 pv[0] = incidentParticle;
341 pv[1] = targetParticle;
342 vecLen = 2;
343
344 if( !inElastic )
345 {
346 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
347
348 G4int iplab = std::min(9,
G4int( incidentTotalMomentum*2.5));
349 if(
G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
350 {
352
353 if ( targetCode == protonCode)
354 {
355 if(ran < 0.2)
356 {
358 }
359 else if (ran < 0.4)
360 {
363 }
364 else if (ran < 0.6)
365 {
368 }
369 else if (ran < 0.8)
370 {
373 }
374 else
375 {
378 }
379 }
380 else
381 {
382 if (ran < 0.2)
383 {
385 }
386 else if (ran < 0.4)
387 {
390 }
391 else if (ran < 0.6)
392 {
395 }
396 else if (ran < 0.8)
397 {
400 }
401 else
402 {
405 }
406 }
407 }
408 return;
409 }
411 return;
412
413
414
415 npos = 0; nneg = 0; nzero = 0;
416 G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.97, 0.88,
417 0.85, 0.81, 0.75, 0.64, 0.64, 0.55, 0.55, 0.45, 0.47, 0.40,
418 0.39, 0.36, 0.33, 0.10, 0.01};
419 G4int iplab =
G4int( incidentTotalMomentum*10.);
420 if ( iplab > 9) iplab = 10 +
G4int( (incidentTotalMomentum -1.)*5. );
421 if ( iplab > 14) iplab = 15 +
G4int( incidentTotalMomentum -2. );
422 if ( iplab > 22) iplab = 23 +
G4int( (incidentTotalMomentum -10.)/10.);
423 iplab = std::min(24, iplab);
424
426
427
428 G4double aleab = std::log(availableEnergy);
429 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
430 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
431
432
433
435
436 for (nt = 1; nt <= numSec; nt++) {
437 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
438 dum =
pi*nt/(2.0*
n*
n);
439 if (std::fabs(dum) < 1.0) {
440 if (test >= 1.0e-10) anpn += dum*test;
441 } else {
442 anpn += dum*test;
443 }
444 }
445
448 if (targetCode == protonCode) {
449 counter = -1;
450 for (npos = 0; npos < numSec/3; npos++) {
451 for (nneg = std::max(0,npos-2); nneg <= (npos+1); nneg++) {
452 for (nzero = 0; nzero < numSec/3; nzero++) {
453 if (++counter < numMul) {
454 nt = npos+nneg+nzero;
455 if ((nt > 0) && (nt <= numSec) ) {
456 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
457 dum = (
pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
458
459 if (std::fabs(dum) < 1.0) {
460 if (test >= 1.0e-10) excs += dum*test;
461 } else {
462 excs += dum*test;
463 }
464
465 if (ran < excs) goto outOfLoop;
466 }
467 }
468 }
469 }
470 }
471
472 inElastic = false;
473 return;
474 } else {
475 counter = -1;
476 for( npos=0; npos<numSec/3; npos++ )
477 {
478 for( nneg=std::max(0,npos-1); nneg<=(npos+2); nneg++ )
479 {
480 for( nzero=0; nzero<numSec/3; nzero++ )
481 {
482 if( ++counter < numMul )
483 {
484 nt = npos+nneg+nzero;
485 if( (nt>0) && (nt<=numSec) )
486 {
487 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
488 dum = (
pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
489 if (std::fabs(dum) < 1.0) {
490 if( test >= 1.0e-10 )excs += dum*test;
491 } else {
492 excs += dum*test;
493 }
494
495 if (ran < excs) goto outOfLoop;
496 }
497 }
498 }
499 }
500 }
501
502 inElastic = false;
503 return;
504 }
505
506 outOfLoop:
507
509
510 if( targetCode == protonCode)
511 {
512 if( npos == nneg)
513 {
514 if (ran < 0.40)
515 {
516 }
517 else if (ran < 0.8)
518 {
520 }
521 else
522 {
525 }
526 }
527 else if (npos == (nneg+1))
528 {
529 if( ran < 0.25)
530 {
532 }
533 else if (ran < 0.5)
534 {
537 }
538 else
539 {
541 }
542 }
543 else if (npos == (nneg-1))
544 {
546 }
547 else
548 {
551 }
552 }
553 else
554 {
555 if( npos == nneg)
556 {
557 if (ran < 0.4)
558 {
559 }
560 else if(ran < 0.8)
561 {
563 }
564 else
565 {
568 }
569 }
570 else if ( npos == (nneg-1))
571 {
572 if (ran < 0.5)
573 {
575 }
576 else if (ran < 0.75)
577 {
579 }
580 else
581 {
584 }
585 }
586 else if (npos == (nneg+1))
587 {
589 }
590 else
591 {
594 }
595 }
596
597 }
598 else
599 {
601 {
602
603 G4double aleab = std::log(availableEnergy);
604 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
605 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
606
607
608
610
611 for (nt=2; nt<=numSec; nt++) {
612 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
613 dum =
pi*nt/(2.0*
n*
n);
614
615 if (std::fabs(dum) < 1.0) {
616 if( test >= 1.0e-10 )anpn += dum*test;
617 } else {
618 anpn += dum*test;
619 }
620 }
621
624 if (targetCode == protonCode) {
625 counter = -1;
626 for (npos=1; npos<numSec/3; npos++) {
627 nneg = npos-1;
628 for( nzero=0; nzero<numSec/3; nzero++ )
629 {
630 if( ++counter < numMulAn )
631 {
632 nt = npos+nneg+nzero;
633 if( (nt>1) && (nt<=numSec) ) {
634 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
635 dum = (
pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n);
636
637 if (std::fabs(dum) < 1.0) {
638 if( test >= 1.0e-10 )excs += dum*test;
639 } else {
640 excs += dum*test;
641 }
642
643 if (ran < excs) goto outOfLoopAn;
644 }
645 }
646 }
647 }
648
649 inElastic = false;
650 return;
651
652 } else {
653 counter = -1;
654 for (npos=0; npos<numSec/3; npos++) {
655 nneg = npos;
656 for( nzero=0; nzero<numSec/3; nzero++ ) {
657 if (++counter < numMulAn) {
658 nt = npos+nneg+nzero;
659 if( (nt>1) && (nt<=numSec) ) {
660 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
661 dum = (
pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n);
662
663 if (std::fabs(dum) < 1.0) {
664 if( test >= 1.0e-10 )excs += dum*test;
665 } else {
666 excs += dum*test;
667 }
668
669 if (ran < excs) goto outOfLoopAn;
670 }
671 }
672 }
673 }
674
675 inElastic = false;
676 return;
677 }
678 outOfLoopAn:
679 vecLen = 0;
680 }
681 }
682
683 nt = npos + nneg + nzero;
684 while ( nt > 0)
685 {
688 {
689 if( npos > 0 )
691 npos--;
692 }
693 }
694 else if ( ran < (
G4double)(npos+nneg)/nt)
695 {
696 if( nneg > 0 )
697 {
699 nneg--;
700 }
701 }
702 else
703 {
704 if( nzero > 0 )
705 {
707 nzero--;
708 }
709 }
710 nt = npos + nneg + nzero;
711 }
713 {
714 G4cout <<
"Particles produced: " ;
717 for (i=2; i < vecLen; i++)
718 {
720 }
722 }
723 return;
724 }
G4double pmltpc(G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
G4HEVector AntiSigmaMinus
G4double getTotalMomentum() const