61 G4cout <<
" **************************************************** " <<
G4endl;
62 G4cout <<
" * The RPG model is currently under development and * " <<
G4endl;
64 G4cout <<
" **************************************************** " <<
G4endl;
77 for( i=2; i<=np; i++ )npf +=
G4Log((
double)i);
78 for( i=2; i<=nneg; i++ )nmf +=
G4Log((
double)i);
79 for( i=2; i<=nz; i++ )nzf +=
G4Log((
double)i);
81 r =
std::min( expxu,
std::max( expxl, -(np-nneg+nz+b)*(np-nneg+nz+b)/(2*c*c*n*n)-npf-nmf-nzf ) );
90 if (j <= 1)
return result;
91 for (
G4int i = 2; i <= j; ++i) result *= i;
113 leadParticle = currentParticle;
120 leadParticle = targetParticle;
131 if( np+nneg+nz == 0 )
return;
134 for( i=0; i<np; ++i )
141 for( i=np; i<np+nneg; ++i )
148 for( i=np+nneg; i<np+nneg+nz; ++i )
165 const G4int numSec = 60;
181 n = 3.62567 + aleab*(0.665843 + aleab*(0.336514 + aleab*(0.117712 + 0.0136912*aleab)));
188 for(
G4int i=iBegin; i<=numSec; ++i )
190 temp =
pi*i/(2.0*n*
n);
194 if( test >= 1.0
e-10 )anpn += temp*
test;
210 G4bool& incidentHasChanged,
229 incidentHasChanged, originalTarget,
230 targetParticle, targetHasChanged,
231 targetNucleus, currentParticle,
233 false, leadingStrangeParticle);
239 leadingStrangeParticle );
244 G4bool finishedGenXPt =
false;
245 G4bool annihilation =
false;
247 currentParticle.
GetMass() == 0.0 && targetParticle.
GetMass() == 0.0 )
256 if( ek > 1.0*
GeV )ekcor = 1./(ek/
GeV);
258 ek = 2*tarmas + ek*(1.+ekcor/atomicWeight);
280 if( ekOrg <= 0.0001 )
289 const G4double twsup[] = { 1.0, 0.7, 0.5, 0.3, 0.2, 0.1 };
296 std::vector<G4ReactionProduct> savevec;
297 for (
G4int i = 0; i < vecLen; i++) savevec.push_back(*vec[i]);
306 if( annihilation || vecLen > 5 ||
314 || rand2 > twsup[vecLen]) ) )
318 incidentHasChanged, originalTarget,
319 targetParticle, targetHasChanged,
320 targetNucleus, currentParticle,
322 leadFlag, leadingStrangeParticle);
324 if (finishedGenXPt)
return;
326 G4bool finishedTwoClu =
false;
329 for (
G4int i = 0; i < vecLen; i++)
delete vec[i];
337 if (!finishedGenXPt && annihilation) {
338 currentParticle = saveCurrent;
339 targetParticle = saveTarget;
340 for (
G4int i = 0; i < vecLen; i++)
delete vec[i];
343 for (
G4int i = 0; i <
G4int(savevec.size()); i++) {
363 incidentHasChanged, originalTarget,
364 targetParticle, targetHasChanged,
365 targetNucleus, currentParticle,
367 leadFlag, leadingStrangeParticle);
376 if (finishedTwoClu)
return;
379 incidentHasChanged, originalTarget,
380 targetParticle, targetHasChanged,
381 targetNucleus, currentParticle,
383 false, leadingStrangeParticle);
407 G4bool& incidentHasChanged )
417 incidentHasChanged =
true;
426 incidentHasChanged =
true;
439 for (i = 0; i < vecLen; ++i) {
443 vec[i]->SetDefinitionAndUpdateE(aKaonZL);
445 vec[i]->SetDefinitionAndUpdateE(aKaonZS);
450 if (incidentHasChanged) {
467 if (std::fabs(aE)<.1*
eV) aE=.1*
eV;
471 if (targetParticle.
GetMass() > 0.0)
480 dir = momentum/momentum.
mag();
489 for (i = 0; i < vecLen; ++i) {
490 G4double secKE = vec[i]->GetKineticEnergy();
496 dir = momentum/momentum.
mag();
505 std::pair<G4int, G4double>
511 for (
G4int i = 1; i < 30; i++) {
518 return std::pair<G4int, G4double>(index, fraction);
527 for (i = 0; i <
G4int(sigma.size()); i++) sum += sigma[i];
533 for (i = 0; i <
G4int(sigma.size()); i++) {
534 partialSum += sigma[i];
535 if (fsum < partialSum) {
561 for (
G4int i = 0; i < vecLen; i++) {
562 secDef = vec[i]->GetDefinition();
570 if (chargeSum != Q) {
574 if (baryonSum != B) {
578 if (strangenessSum != S) {
586 for (
G4int i = 0; i < vecLen; i++) {
587 secDef = vec[i]->GetDefinition();
597 0.0, 0.01, 0.013, 0.018, 0.024, 0.032, 0.042, 0.056, 0.075, 0.1,
598 0.13, 0.18, 0.24, 0.32, 0.42, 0.56, 0.75, 1.0, 1.3, 1.8,
599 2.4, 3.2, 4.2, 5.6, 7.5, 10.0, 13.0, 18.0, 24.0, 32.0 };