lars/geant4
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
fbd4999cf7355ec0b49f3aa4c3a443c94f2c0529
geant4/source/processes/hadronic/models/cascade/pre_compound/src/G4PreEquilibrium.cc
T
Gabriele Cosmo fbd4999cf7 Import Geant4 5.1.0 source tree
2016-06-09 10:15:15 +02:00

3611 lines
123 KiB
C++
Raw Blame History

//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4PreEquilibrium.cc,v 1.9 2001/07/11 10:03:55 gunter Exp $
// GEANT4 tag $Name: geant4-05-01 $
//
// Hadronic Process: Pre-equilibrium HETC
// Joseph L. Chuma, TRIUMF, 24-Mar-2000
#include "G4PreEquilibrium.hh"
#include "Randomize.hh"
G4double oneThird = 1./3.;
G4double twoThirds = 2./3.;
G4double fourThirds = 4./3.;
// S(Z)+P(Z) from Tab. 1 from A.G.W. Cameron, Canad. J. Phys., 35(1957)1021
// or Delta M(Z) from Tab. 97 of book [1]
G4double G4PreEquilibrium::t1y[130] =
{ 20.80, 15.80, 21.00, 16.80, 19.80,
16.50, 18.80, 16.50, 18.50, 17.20,
18.26, 15.05, 16.01, 12.04, 13.27,
11.09, 12.17, 10.26, 11.04, 8.41,
9.79, 7.36, 8.15, 5.63, 5.88,
3.17, 3.32, 0.82, 1.83, 0.97,
2.33, 1.27, 2.92, 1.61, 2.91,
1.35, 2.40, 0.89, 1.74, 0.36,
0.95, -0.65, -0.04, -1.73, -0.96,
-2.87, -2.05, -4.05, -3.40, -5.72,
-3.75, -4.13, -2.42, -2.85, -1.01,
-1.33, 0.54, -0.02, 1.74, 0.75,
2.24, 1.00, 1.98, 0.79, 1.54,
0.39, 1.08, 0.00, 0.78, -0.35,
0.58, -0.55, 0.59, -0.61, 0.59,
-0.35, 0.32, -0.96, -0.52, -2.08,
-2.46, -3.64, -1.55, -0.96, 0.97,
0.88, 2.37, 1.75, 2.72, 1.90,
2.55, 1.46, 1.93, 0.86, 1.17,
0.08, 0.39, -0.76, -0.39, -1.51,
-1.17, -2.36, -1.95, -3.06, -2.62,
-3.55, -2.95, -3.75, -3.07, -3.79,
-3.06, -3.77, -3.05, -3.78, -3.12,
-3.90, -3.35, -4.24, -3.86, -4.92,
-5.06, -6.77, -7.41, -9.18, -10.16,
-11.12, -9.76, -9.23, -7.96, -7.65 };
// S(N)+P(N) from Tab. 1 from A.G.W. Cameron, Canad. J. Phys., 35(1957)1021
// or Delta M(N) from Tab. 97 of book [1]
G4double G4PreEquilibrium::t2xy[] =
{ -8.40, -12.90, -8.00, 11.90, -9.20,
-12.50, -10.80, -13.60, -11.20, -12.20,
-12.81, -15.40, -13.07, -15.80, -13.81,
-14.98, -12.63, -13.76, -11.37, -12.38,
-9.23, -9.65, -7.64, -9.17, -8.05,
-9.72, -8.87, -10.76, -8.64, -8.89,
-6.60, -7.13, -4.77, -5.33, -3.06,
-3.79, -1.72, -2.79, -0.93, -2.19,
-0.52, -1.90, -0.45, -2.20, -1.22,
-3.07, -2.42, -4.37, -3.94, -6.08,
-4.49, -4.50, -3.14, -2.93, -1.04,
-1.36, 0.69, 0.21, 2.11, 1.33,
3.29, 2.46, 4.30, 3.32, 4.79,
3.62, 4.97, 3.64, 4.63, 3.07,
4.06, 2.49, 3.30, 1.46, 2.06,
0.51, 0.74, -1.18, -1.26, -3.54,
-3.97, -5.26, -4.18, -3.71, -2.10,
-1.70, -0.08, -0.18, 0.94, 0.27,
1.13, 0.08, 0.91, -0.31, 0.49,
-0.78, 0.08, -1.15, -0.23, -1.41,
-0.42, -1.55, -0.55, -1.66, -0.66,
-1.73, -0.75, -1.74, -0.78, -1.69,
-0.78, -1.60, -0.75, -1.46, -0.67,
-1.26, -0.51, -1.04, -0.53, -1.84,
-2.42, -4.52, -4.76, -6.33, -6.76,
-7.81, -5.80, -5.37, -3.63, -3.35,
-1.75, -1.88, -0.61, -0.90, 0.09,
-0.32, 0.55, -0.13, 0.70, -0.06,
0.49, -0.20, 0.40, -0.22, 0.36,
-0.09, 0.58, 0.12, 0.75, 0.15,
0.70, 0.17, 1.11, 0.89, 1.85,
1.62, 2.54, 2.29, 3.20, 2.91,
3.84, 3.53, 4.48, 4.15, 5.12,
4.78, 5.75, 5.39, 6.31, 5.91,
6.87, 6.33, 7.13, 6.61, 7.30,
6.31, 6.27, 4.83, 4.49, 2.85,
2.32, 0.58, -0.11, -0.98, 0.81,
1.77, 3.37, 4.13, 5.60, 6.15,
7.29, 7.35, 7.95, 7.67, 8.16,
7.83, 8.31, 8.01, 8.53, 8.27
};
G4double G4PreEquilibrium::cemgeo( const G4int nuclideTypeNumber, const G4int mediumNumber )
{
// this routine was called by input
bertcem( nuclideTypeNumber, mediumNumber );
return GetA_FCOMON(nuclideTypeNumber,mediumNumber);
}
G4double G4PreEquilibrium::erupcem( DPvector &neutrons, DPvector &protons, DPvector &deuterons,
DPvector &tritons, DPvector &he3s, DPvector &he4s )
{
precof();
G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
G4ParticleDefinition *aProton = G4Proton::Proton();
G4ParticleDefinition *aDeuteron = G4Deuteron::Deuteron();
G4ParticleDefinition *aTriton = G4Triton::Triton();
G4ParticleDefinition *aHe3 = G4He3::He3();
//G4ParticleDefinition *aHe4 = G4He4::He4();
G4double protonMass = aProton->GetPDGMass()/MeV;
G4double neutronMass = aNeutron->GetPDGMass()/MeV;
G4double deuteronMass = aDeuteron->GetPDGMass()/MeV;
G4double tritonMass = aTriton->GetPDGMass()/MeV;
G4double he3Mass = aHe3->GetPDGMass()/MeV;
//G4double he4Mass = aHe4->GetPDGMass()/MeV;
G4int numberOfProtons = G4int(theNucleus->GetZ()+0.5);
G4int numberOfNucleons = G4int(theNucleus->GetN()+0.5);
G4int numberOfNeutrons = numberOfNucleons - numberOfProtons;
G4double nuclearMass = protonMass*numberOfProtons + neutronMass*numberOfNeutrons;
G4double P = theNucleus->GetMomentum().mag()/MeV;
G4double totalEnergy = sqrt( P*P + nuclearMass*nuclearMass );
G4double kineticEnergyAfterEvaporation = totalEnergy - nuclearMass;
for( G4int k = 0; k < KTOT; ++k )
{
Dvector sptK( *(spt[k]) );
G4double sint = sptK[0];
G4double cost = sptK[1];
G4double sinf = sin( parz[k].angle2 );
G4double cosf = cos( parz[k].angle2 );
G4double alc = sint*cosf;
G4double bec = sint*sinf;
G4double gac = cost;
G4double energy = parz[k].energy*GeV/MeV;
G4DynamicParticle *dp = new G4DynamicParticle();
G4ThreeVector m;
switch ( parz[k].particleType )
{
case 1: // neutrons
dp->SetDefinition( aNeutron );
dp->SetKineticEnergy( energy-neutronMass );
m.setX( alc ); m.setY( bec ); m.setZ( gac );
dp->SetMomentum( m );
neutrons.push_back( dp );
break;
case 2: // protons
dp->SetDefinition( aProton );
dp->SetKineticEnergy( energy-protonMass );
m.setX( alc ); m.setY( bec ); m.setZ( gac );
dp->SetMomentum( m );
protons.push_back( dp );
break;
case 3: // deuterons
dp->SetDefinition( aDeuteron );
dp->SetKineticEnergy( energy-deuteronMass );
m.setX( alc ); m.setY( bec ); m.setZ( gac );
dp->SetMomentum( m );
deuterons.push_back( dp );
break;
case 4: // tritons
dp->SetDefinition( aTriton );
dp->SetKineticEnergy( energy-tritonMass );
m.setX( alc ); m.setY( bec ); m.setZ( gac );
dp->SetMomentum( m );
tritons.push_back( dp );
break;
case 5: // he-3
dp->SetDefinition( aHe3 );
dp->SetKineticEnergy( energy-he3Mass );
m.setX( alc ); m.setY( bec ); m.setZ( gac );
dp->SetMomentum( m );
he3s.push_back( dp );
break;
//case 6: // he-4
//dp->SetDefinition( aHe4 );
//dp->SetKineticEnergy( energy-he4Mass );
//m.setX( alc ); m.setY( bec ); m.setZ( gac );
//dp->SetMomentum( m );
//dp->SetMomentum( alc, bec, gac );
//he4s.push_back( dp );
//break;
}
}
return kineticEnergyAfterEvaporation;
}
void G4PreEquilibrium::precof()
{
const G4double RM = 1.5;
G4double pevapj[6], gj[7], ami[6];
G4int i, j;
SetR0( RM );
G4double wt = 1.0;
G4double protonMass = G4Proton::Proton()->GetPDGMass()/GeV;
G4double neutronMass = G4Neutron::Neutron()->GetPDGMass()/GeV;
G4int numberOfProtons = G4int(theNucleus->GetZ()+0.5);
G4int numberOfNucleons = G4int(theNucleus->GetN()+0.5);
G4int numberOfNeutrons = numberOfNucleons - numberOfProtons;
G4double nuclearMass = protonMass*numberOfProtons + neutronMass*numberOfNeutrons;
G4double P = theNucleus->GetMomentum().mag()/GeV;
G4double totalEnergy = sqrt( P*P + nuclearMass*nuclearMass );
G4double kineticEnergy = (totalEnergy-nuclearMass)*GeV/MeV;
G4ThreeVector momentum = theNucleus->GetMomentum() * (1/totalEnergy*GeV);
G4ThreeVector uglmom( GetLXYZ() );
KTOT = 0;
G4double pz = GetPZ0();
G4int n = GetN0();
G4double p = GetP0();
G4double h = GetH0();
theNucleus->AddExcitationEnergy( -kineticEnergy );
//
if( theNucleus->GetEnergyDeposit() < 3 ||
theNucleus->GetN() <= 4 ||
theNucleus->GetZ() <= 2 )return;
//
G4double dl = massDefect( theNucleus->GetN(), theNucleus->GetZ() );
//
// auxiliary code for nuclear data extraction
//
const G4double z1[5] = { 10.00, 20.00, 30.00, 50.00, 70.00 }; // Z:Tab.96
const G4double a1[5] = { 0.42, 0.58, 0.68, 0.77, 0.80 }; // Kp -//-
const G4double c1[5] = { 0.50, 0.28, 0.20, 0.15, 0.10 }; // Cp -//-
const G4double a2[5] = { 0.68, 0.82, 0.91, 0.97, 0.98 }; // Ka -//-
const G4double c2[5] = { 0.10, 0.10, 0.10, 0.08, 0.06 }; // Ca -//-
//
// particle b: n p d t He-3 He-4
const G4double aj[6] = { 1.0, 1.0, 2.0, 3.0, 3.0, 4.0 }; // A(b)
const G4double zj[6] = { 0.0, 1.0, 1.0, 1.0, 2.0, 2.0 }; // Z(b)
const G4double dlm[6] = { 8.368, 7.569, 13.835, 15.835, 15.817, 3.607 }; // q(b)
//
G4double cc[6], vk[6];
alj.reserve(6);
gb.reserve(6);
afj.reserve(7);
zfj.reserve(6);
rj.reserve(7);
vj.reserve(6);
bj.reserve(6);
r0j.reserve(6);
for( G4int k = 0; k < 100; ++k )
{
cc[0] = 0;
cc[1] = quadraticInterpolation( theNucleus->GetZ(), z1, c1, 4 );
cc[5] = quadraticInterpolation( theNucleus->GetZ(), z1, c2, 4 );
cc[2] = cc[1]/2;
cc[3] = cc[1]/3;
cc[4] = cc[5]*fourThirds;
//
vk[0] = 0;
vk[1] = quadraticInterpolation( theNucleus->GetZ(), z1, a1, 4 );
vk[5] = quadraticInterpolation( theNucleus->GetZ(), z1, a2, 4 );
vk[2] = vk[1] + 0.06;
vk[3] = vk[1] + 0.12;
vk[4] = vk[5] - 0.06;
//
G4double um2 = uglmom.mag2();
G4double un = theNucleus->GetN() - theNucleus->GetZ();
G4double ue = theNucleus->GetEnergyDeposit()-
12*((1-theNucleus->GetZ()+2*(G4int(theNucleus->GetZ()+0.5)/2))+
(1-un+2*(G4int(un)/2)))/sqrt(theNucleus->GetN());
if( ue <= 0.1 )return;
G4double am = fam( theNucleus->GetN(), theNucleus->GetZ(), ue );
afj.erase( afj.begin(), afj.end() );
zfj.erase( zfj.begin(), zfj.end() );
rj.erase( rj.begin(), rj.end() );
bj.erase( bj.begin(), bj.end() );
vj.erase( vj.begin(), vj.end() );
for( i = 0; i < 6; ++i )
{
G4double naj = theNucleus->GetN() - aj[i];
G4double zzj = theNucleus->GetZ() - zj[i];
afj.push_back( naj );
zfj.push_back( zzj );
G4double pairj = (1-zzj+2*(G4int(zzj)/2)) + (1-(naj-zzj)+2*(G4int(naj-zzj)/2));
pevapj[i] = 12 * pairj / sqrt(naj);
G4double uej = theNucleus->GetEnergyDeposit() - pevapj[i];
//
if( uej <= 1 )rj.push_back( 0.0 );
else
{
ami[i] = fam( naj, zzj, uej );
//
// calculation of coulomb energy
//
G4double tmp = vk[i]*1.44/RM*zj[i]*zzj/
(pow(aj[i],oneThird)+pow(naj,oneThird))*
(1-theNucleus->GetEnergyDeposit()/(81*theNucleus->GetN()*am));
tmp < 0 ? vj.push_back( 0.0 ) : vj.push_back( tmp );
//
bj.push_back( massDefect(naj,zzj) - (dl-dlm[i]) );
rj.push_back( uej-bj[i]-vj[i] );
}
}
afj.push_back( theNucleus->GetN() );
//
r0j.erase( r0j.begin(), r0j.end() );
r0j.push_back( 0.76+2.2/pow(afj[1],oneThird) );
SetBN( (2.12/pow(afj[1],twoThirds)-0.05)/r0j[0] );
for( i = 1; i < 6; ++i )r0j.push_back( 1+cc[i] );
G4double ep1, ep2, ep3;
G4bool ind;
G4ThreeVector p12;
G4int lm;
LABEL:
if( n < sqrt(1.19*am*theNucleus->GetN()*theNucleus->GetEnergyDeposit()+0.5) )
{
// pre-equilibrium emission
//
while ( p < 1 )
{
++p;
++h;
n += 2;
}
exn = p + h;
alj.erase( alj.begin(), alj.end() );
alj.push_back( p*(exn-1) );
alj.push_back( alj[0] );
alj.push_back( alj[1]*(p-1)*(exn-2)/2 );
alj.push_back( alj[2]*(p-2)*(exn-3)/6 );
alj.push_back( alj[3] );
alj.push_back( alj[4]*(p-3)*(exn-4)/12 );
//
G4double nucleusN = theNucleus->GetN();
gb.erase( gb.begin(), gb.end() );
gb.push_back( 1. );
gb.push_back( 1. );
gb.push_back( 16/nucleusN );
gb.push_back( 243/nucleusN/nucleusN );
gb.push_back( gb[3] );
gb.push_back( 4096/nucleusN/nucleusN/nucleusN );
for( i = 0; i < 6; ++i )
{
if( p <= aj[i] - 0.01 )gj[i] = 0;
else
{
if( p+h <= aj[i] + 0.01 )gj[i] = 0;
else
{
if( pz <= zj[i] - 0.01 )gj[i] = 0;
else
{
if( rj[i] )gj[i] = 0;
else
{
SetAC( 0.585*ami[i] );
gj[i] = gamagu(i);
}
}
}
}
}
G4double g = 0;
for( i = 0; i < 6; ++i )g += gj[i];
if( n <= 0 )
{
// write(16,20)
// 20 format(20x,20hnumber of exitons = 0)
return;
}
SetAC( 0.595*am );
G4double c1, c2, c3;
transitionRates( p, h, c1, c2, c3 );
G4double c = c1+c2+c3;
G4double b1 = G4UniformRand();
G4double p2 = 250*pow(G4UniformRand(),oneThird);
G4double ct2 = 1-2*G4UniformRand();
G4double st2 = sqrt(1-ct2*ct2);
G4double fi2 = 2*pi*G4UniformRand();
p12.setX( p2*st2*cos(fi2)+theNucleus->GetMomentum().x()*GeV/p );
p12.setY( p2*st2*sin(fi2)+theNucleus->GetMomentum().y()*GeV/p );
p12.setZ( p2*ct2+theNucleus->GetMomentum().z()*GeV/p );
if( b1 > g/(c+g) )
{
G4double ran = G4UniformRand();
if( ran <= c1/c )
{
++p; ++h; ++(++n);
if( G4UniformRand() <= theNucleus->GetZ()/theNucleus->GetN() )++pz;
}
else
{
if( ran <= (c2+c1)/c )
{
--(--n); --p; --h;
if( pz > 0 && G4UniformRand() <= theNucleus->GetZ()/theNucleus->GetN() )--pz;
}
}
goto LABEL;
}
for( j = 1; j < 6; ++j )gj[j] += gj[j-1];
G4double b = g*G4UniformRand();
for( j = 0; j < 6; ++j )
{
if( b <= gj[j] )
{
lm = j;
break;
}
}
ep1 = tkinm1( lm, p, h );
ep2 = zj[lm];
ep3 = 940*aj[lm];
ind = true;
}
else
{
// equilibrium emission
//
rj.push_back( 0.0 );
G4double per = arfaf( ami );
for( i = 0; i < 7; ++i )
{
rj[i] <= 0 ? gj[i] = 0 : gj[i] = gameqf( i, cc[i], per, ami[i], RM );
}
G4double gt = 0;
for( i = 0; i < 7; ++i )gt += gj[i];
if( gt <= 0 )return;
wt *= 1-gj[7]/gt;
SetWF( 1-wt );
for( j = 1; j < 7; ++j )gj[j] += gj[j-1];
G4double b = G4UniformRand()*gt;
for( j = 1; j <= 7; ++j )
{
if( b < gj[j] )
{
lm = j;
break;
}
}
if( lm == 7 )return;
ep1 = tkin( lm, am );
ep2 = zj[lm];
ep3 = 940*aj[lm];
ind = false;
}
++KTOT;
p -= aj[lm];
n -= G4int(aj[lm]);
pz -= ep2;
parz[k].charge = ep2;
parz[k].particleType = lm;
if( !ind )parz[k].what = 1000.0;
else parz[k].what = 100.0;
theNucleus->SetParameters( afj[lm], zfj[lm] );
theNucleus->AddExcitationEnergy(
-theNucleus->GetEnergyDeposit()+rj[lm]-ep1+vj[lm]+pevapj[lm] );
G4double tl, cot, st, cf, sf;
G4ThreeVector pl;
Dvector *vtmp = new Dvector( 5 );
vtmp->insert( vtmp->begin()+3, ep2 );
vtmp->insert( vtmp->begin()+4, ep3/GeV );
if( ind && (momentum.mag2()/GeV > 0.0000001) )
{
G4double cx = sqrt(G4UniformRand());
G4double sx = sqrt(1-cx*cx);
G4double fi = 2*pi*G4UniformRand();
G4double pm = sqrt(ep1*(ep1+2*ep3));
pl.setX( pm*sx*cos(fi) );
pl.setY( pm*sx*sin(fi) );
pl.setZ( pm*cx );
G4ThreeVector ps = rotation( p12, momentum, pl );
tl = cinema( ps, momentum, pl, cot, st, cf, sf, ep3 );
vtmp->insert( vtmp->begin(), st );
vtmp->insert( vtmp->begin()+1, cot );
vtmp->insert( vtmp->begin()+2, tl/GeV );
}
else
{
G4double angl[4];
//
// choose isotropic distributed angle for particle emitted
//
angl[0] = 1-2*G4UniformRand();
angl[3] = sqrt(1-angl[0]*angl[0]);
G4double f = 2*pi*G4UniformRand();
angl[1] = sin(f);
angl[2] = cos(f);
//
G4double pm = sqrt(ep1*(ep1+2*ep3));
G4ThreeVector ps( pm*angl[3]*angl[2], pm*angl[3]*angl[1], pm*angl[0] );
if( momentum.mag2() <= 0.0001 )
{
vtmp->insert( vtmp->begin(), angl[3] );
vtmp->insert( vtmp->begin()+1, angl[0] );
vtmp->insert( vtmp->begin()+2, ep1/GeV );
cot = angl[0];
st = angl[2];
cf = angl[1];
sf = angl[0];
pl = ps;
}
else
{
tl = cinema( ps, momentum, pl, cot, st, cf, sf, ep3 );
vtmp->insert( vtmp->begin(), st );
vtmp->insert( vtmp->begin()+1, cot );
vtmp->insert( vtmp->begin()+2, tl/GeV );
}
}
spt.push_back( vtmp );
theNucleus->GetMomentum().setX( theNucleus->GetMomentum().x()-pl.x()*GeV );
theNucleus->GetMomentum().setY( theNucleus->GetMomentum().y()-pl.y()*GeV );
theNucleus->GetMomentum().setZ( theNucleus->GetMomentum().z()-pl.z()*GeV );
numberOfProtons = G4int(theNucleus->GetZ()+0.5);
numberOfNucleons = G4int(theNucleus->GetN()+0.5);
numberOfNeutrons = numberOfNucleons - numberOfProtons;
nuclearMass = protonMass*numberOfProtons + neutronMass*numberOfNeutrons;
P = theNucleus->GetMomentum().mag()/GeV;
totalEnergy = sqrt( P*P + nuclearMass*nuclearMass );
momentum = theNucleus->GetMomentum() * (1/totalEnergy/GeV);
parz[k].angle1 = atan2(st,cot);
G4double fi1 = atan2(sf,cf);
if( fi1 < 0 )fi1 = pi-fi1;
parz[k].angle2 = fi1;
parz[k].energy = (*vtmp)[2];
G4double almax = 0.219327*RM*(pow(theNucleus->GetN(),oneThird)+pow(aj[lm],oneThird))*
sqrt(theNucleus->GetN()*aj[lm]*(ep1-vj[lm])/(theNucleus->GetN()+aj[lm]));
G4double alp = almax * sqrt(G4UniformRand());
uglmom.setX( uglmom.x() - alp*st*cf );
uglmom.setY( uglmom.y() - alp*st*sf );
uglmom.setZ( uglmom.z() - alp*cot );
}
// write(16,49)u,a,z
// format('massiv spt exceeded after evaporation',40x,'u=',f10.5,' a=',f5.1,' z=',f4.1)
return;
}
void G4PreEquilibrium::mashnk( G4int &icem, const G4double anucc, const G4double znucc,
const G4double ecno, const G4double tipno,
const G4double ehipi, const G4double ehin,
const G4double ehicut )
{
// this routine is called by cascade
//
// ehicem is the value of max energy for which siginb is calculated for ifircl > 0
// siginb is inelastic cross-section for either bert or bertcem
//
// icem is trigger for use of bertcem of s.mashik cem code.
// if imash == 0, bertcem is not used and icem = 0
// if imash == 1, ehicem is max energy for which icem = 1 and bertcem is called
//
icem = 0;
G4double ehicem; // ehicem is in common/callcm/ in cascad.f
G4int mat, icm; // mat is in COMON.F, icm is in common/callcm/ in cascad.f
G4double ehic[7]; // ehic is in common/callcm/ in cascad.f
G4int imash = 0, ihecc = 0;
if( imash == 0 ) // imash is in common/callcm/ in cascad.f
{
if( ihecc == 0 ) // ihecc in in common/hecc1/ in cascad.f
{
icem = 0;
ehicem = ehicut;
icm = icem;
ehic[mat-1] = ehicem;
return;
}
else if( ihecc == 1 )
{
icem = 0;
ehicem = 10000;
icm = icem;
ehic[mat-1] = ehicem;
return;
}
}
if( anucc < 14 ) ehicem = 0;
else if( anucc < 17 )ehicem = 1000;
else
{
// icem = 0 for oxygen and nitrogen if ecno > ehicem
// ehicem can be set by trial and error for this case
//
// next statement is tentative
//
tipno <= 1 ? ehicem = ehin : ehicem = ehipi;
//
// icem = 0 for energies > ehicem, which can be set at will
// may want to set it lower if ihecc == 1
//
ecno <= ehicem ? icem = 1 : icem = 0; // if icem == 1, bertcem is called
}
icm = icem;
ehic[mat-1] = ehicem;
return;
}
G4double G4PreEquilibrium::bf( const G4double a, const G4double z, const G4double e )
{
const G4double xx1[51] =
{ 0.00, 0.02, 0.04, 0.06, 0.08, 0.10, 0.12, 0.14, 0.16, 0.18,
0.20, 0.22, 0.24, 0.26, 0.28, 0.30, 0.32, 0.34, 0.36, 0.38,
0.40, 0.42, 0.44, 0.46, 0.48, 0.50, 0.52, 0.54, 0.56, 0.58,
0.60, 0.62, 0.64, 0.66, 0.68, 0.70, 0.72, 0.74, 0.76, 0.78,
0.80, 0.82, 0.84, 0.86, 0.88, 0.90, 0.92, 0.94, 0.96, 0.98, 1.0 };
//
// b(x) from Krappe, Nix & Sierk 1979 for a/r0 = 0.5508
//
const G4double yy2[51] =
{ 0.1495, 0.1475, 0.1450, 0.1425, 0.1400, 0.1375, 0.1350, 0.1328, 0.1298, 0.1270,
0.1245, 0.1215, 0.1185, 0.1150, 0.1118, 0.1083, 0.1040, 0.0997, 0.0949, 0.0900,
0.0849, 0.0795, 0.0739, 0.0685, 0.0630, 0.0578, 0.0520, 0.0460, 0.0408, 0.0342,
0.0289, 0.0236, 0.0188, 0.0146, 0.0111, 0.0082, 0.0059, 0.0043, 0.0029, 0.0020,
0.0017, 0.0012, 0.0010, 0.0007, 0.0004, 0.0002, 0.0001, 0.0000, 0.0000, 0.0000, 0.0 };
//
const G4double r0m = 1.2;
//
// Yukawa-plus-exponential macroscopic model of Krappe, Nix & Sierk (1979)
//
const G4double gamma = 3;
const G4double a2 = 21.7;
const G4double a3 = 0.7322;
G4double a3rt = pow(a,oneThird);
G4double sufnuc = a2*(1-gamma*(1-2*z/a)*(1-2*z/a))*a3rt*a3rt;
G4double x = 0.5*a3*z*z/a3rt/sufnuc;
if ( x >= 1 )return 0.0;
//
// Cameron (Can.J.Phys.35(1957)1021) shell and pairing corr. for g.s. mass
// calculate saddle-point shell and/or pairing corrections
//
return sufnuc*quadraticInterpolation( x, xx1, yy2, 50 ) -
GetT1Y(G4int(z)-1) - GetT2XY(G4int(a-z)-1);
}
G4double G4PreEquilibrium::quadraticInterpolation( const G4double u, const G4double *e,
const G4double *f, const G4int n )
{
G4double x1, x2, x3, y1, y2, y3;
if( u <= e[0] )
{
x1 = e[0];
x2 = e[1];
x3 = e[2];
y1 = f[0];
y2 = f[1];
y3 = f[2];
}
else
{
if( u >= e[n-1] )
{
if( u > e[n] )return f[n];
x1 = e[n-2];
x2 = e[n-1];
x3 = e[n];
y1 = f[n-2];
y2 = f[n-1];
y3 = f[n];
}
else
{
for( G4int j = 0; j < n; ++j )
{
if( u < e[j] )
{
x1 = e[j-1];
x2 = e[j];
x3 = e[j+1];
y1 = f[j-1];
y2 = f[j];
y3 = f[j+1];
break;
}
}
}
}
return y1*(u-x2)*(u-x3)/((x1-x2)*(x1-x3)) + y2*(u-x1)*(u-x3)/((x2-x1)*(x2-x3)) +
y3*(u-x1)*(u-x2)/((x3-x1)*(x3-x2));
}
G4double G4PreEquilibrium::massDefect( const G4double x, const G4double y ) const
{
// calculation of mass defect
//
G4double c = pow(x,oneThird);
G4double d = pow(x,twoThirds);
G4double e = pow(x,fourThirds);
G4double f = pow(y,fourThirds);
G4double a = 1 - 0.62025/d;
G4double b = (x-2*y)/x;
G4double es = (25.8357-44.2355*b*b)*a*a*d;
G4double ec = 0.779*y*(y-1)/c*( 1 - 1.5849/d + 1.2273/x + 1.5772/e );
G4double ealfa = -0.4323*f/c*( 1 + 0.49597/x - 0.57811/c - 0.14518/d );
G4double edob = 8.367*x + 31.4506*x*b*b - 0.783*y - 17.0354*x;
G4int i = G4int(x);
G4int j = G4int(y);
return es + ec + ealfa + edob + GetT1Y(j-1) + GetT2XY(i-j-1);
}
G4double G4PreEquilibrium::gamagu( const G4int j )
{
// particle b: n p d t He-3 He-4
const G4double aj[6] = { 1.0, 1.0, 2.0, 3.0, 3.0, 4.0 }; // A(b)
const G4double w[8] = { 0.1012285363, 0.2223810345, 0.3137066459, 0.3626837834,
0.3626837834, 0.3137066459, 0.2223810345, 0.1012285363 };
const G4double fiks[8] = { 0.9602898565, 0.7966664774, 0.5255324099, 0.1834346425,
-0.1834346425, -0.5255324099, -0.7966664774, -0.9602898565 };
//
G4double a = 0;
if( j > 0 )a = vj[j];
G4double b = theNucleus->GetEnergyDeposit() - bj[j];
G4double y = 0;
for( G4int k = 0; k < 8; ++k )
{
G4double e = 0.5*( (b-a)*fiks[k] + b + a );
G4double result;
if( j > 1 )
result = gb[j]*r0j[j]*0.104/
(GetR0()*pow(afj[j],oneThird)*sqrt(aj[j]*theNucleus->GetEnergyDeposit()))*
alj[j]*((e-vj[j])/theNucleus->GetEnergyDeposit())*
pow(abs((e+bj[j])/theNucleus->GetEnergyDeposit()),aj[j]-1.5)*
pow(abs(1-(e+bj[j])/theNucleus->GetEnergyDeposit()),exn-aj[j]-1);
else if( j == 1 )
result = 0.000234*GetR0()*GetR0()*pow(afj[1],twoThirds)*r0j[1]*
alj[1]/(GetAC()*theNucleus->GetEnergyDeposit()*afj[1])*
pow(1-(e+bj[1])/theNucleus->GetEnergyDeposit(),exn-2)*(e-vj[1]);
else // j must be 0
result = 0.000234*GetR0()*GetR0()*pow(afj[0],twoThirds)*r0j[0]*
alj[0]/(GetAC()*theNucleus->GetEnergyDeposit()*afj[0])*
pow(1-(e+bj[0])/theNucleus->GetEnergyDeposit(),exn-2)*(e+GetBN());
y += w[k]*(b-a)*result/2;
}
return y;
}
G4double G4PreEquilibrium::tkin( const G4int i, const G4double am )
{
//
// kinetic energy for particles in equilibrium decay
//
G4double rk, frk;
G4double rb = 4*am*afj[i]*rj[i];
do
{
G4double b1 = G4UniformRand();
rk = 1 + 1/sqrt(rb)*log(b1+(1-b1)*exp(-sqrt(rb)));
if( i==0 )
{
G4double q1 = 1 + (2.12/pow(afj[0],twoThirds)-0.05)/
(0.76+2.2/pow(afj[0],oneThird))/rj[0];
frk = 3*sqrt(3.)/2/(q1*sqrt(q1))*(q1*rk-rk*rk*rk);
}
else
{
frk = 3*sqrt(3.)/2*(rk-rk*rk*rk);
}
} while( G4UniformRand() > frk );
return rj[i]*(1-rk*rk) + vj[i];
}
void G4PreEquilibrium::transitionRates( const G4double p, const G4double h,
G4double &c1, G4double &c2, G4double &c3 )
{
// calculation of transition rates (was trncem)
//
G4double est = 1.6*45+theNucleus->GetEnergyDeposit()/(p+h);
G4double b = sqrt(2*est/940);
G4double sf = (10.63/(b*b)-29.93/b+42.9+34.10/(b*b)-82.20/b+82.2)/2;
G4double t;
if( 45/est <= 0.5 )
t = 1-7*45/est/5;
else
t = 1-7*45/est/5+0.4*45/est*pow(2-1/45*est,2.5);
c1 = 0.00332*sf*t*sqrt(est)/pow(1.2+1/(4.7*b),3);
c2 = (c1*p*h*(p+h+1)*(p+h-2))/(GetAC()*theNucleus->GetN()*theNucleus->GetEnergyDeposit()*
GetAC()*theNucleus->GetN()*theNucleus->GetEnergyDeposit());
c3 = c1*(p+h+1)*(p*(p+1)+4*p*h+h*(h-1))/((p+h)*GetAC()*theNucleus->GetN()*
theNucleus->GetEnergyDeposit());
if( c2 < 0 )c2 = 0;
return;
}
G4double G4PreEquilibrium::tkinm1( const G4int j, const G4double p, const G4double h )
{
// kinetic energy for particles in pre-equilibrium decay
//
// particle b: n p d t He-3 He-4
const G4double aj[6] = { 1.0, 1.0, 2.0, 3.0, 3.0, 4.0 }; // A(b)
//
G4double e, e1, dj, t3;
j == 0 ?
dj = (2.12/pow(afj[0],twoThirds)-0.05)/(0.76+2.2/pow(afj[0],oneThird)) :
dj = -vj[j];
G4double t = p + h - aj[j] - 1;
G4double r2 = rj[j];
G4double r1 = r2 + vj[j];
if( j < 2 )
{
if( t <= -0.01 )return r1;
if( t <= 0.1 )
{
G4double result;
j <= 1 ?
result = -dj+sqrt(dj*dj+(G4UniformRand()*(r2*r2+2*dj*r2))) :
result = sqrt(G4UniformRand())*r2+vj[j];
return result;
}
e1 = (r1-dj*t)/(t+1);
do
{
e = vj[j] + G4UniformRand()*r2;
t3 = (e+dj)/(e1+dj)*pow(abs((r1-e)/(r1-e1)),t);
}
while ( G4UniformRand() > t3 );
return e;
}
if( t <= -0.1 )return r1;
if( t <= 0.1 )
{
do
{
e = vj[j] + G4UniformRand()*r2;
t3 = pow( abs((theNucleus->GetEnergyDeposit()-r1+e)/
(theNucleus->GetEnergyDeposit()-r1+r1)), aj[j]-1.5 ) *
(e+dj)/(r1+dj);
}
while ( G4UniformRand() > t3 );
return e;
}
G4double es = theNucleus->GetEnergyDeposit()*(aj[j]-0.5) +
(theNucleus->GetEnergyDeposit()-r2)*(p+h-2.5);
e1 = (es+sqrt(es*es-(theNucleus->GetEnergyDeposit()-r2)*(aj[j]-1.5)*
(p+h-1.5)*4.*theNucleus->GetEnergyDeposit()))/
((p+h-1.5)*2)-theNucleus->GetEnergyDeposit()+r1;
do
{
e = vj[j] + G4UniformRand()*r2;
if( e1+0.001 > r1 )break;
t3 = pow( abs((theNucleus->GetEnergyDeposit()-r1+e)/
(theNucleus->GetEnergyDeposit()-r1+e1)), aj[j]-1.5 ) *
(e+dj)/(e1+dj) * pow( abs((r1-e)/(r1-e1)), t );
}
while ( G4UniformRand() > t3 );
return e;
}
void G4PreEquilibrium::lpoly( const G4double x, const G4int n, G4double *pl )
{
// This subroutine calculates the ordinary Legendre polynomials of
// order 0 to n-1 of argument x and stores them in the vector pl.
// They are calculated by recursion relation from the first two polynomials.
//
// written by A. J. Sierk LANL T-9 February,1984
//
pl[0] = 1;
pl[1] = x;
for( G4int i = 2; i < n; ++i )pl[i] = ((2*i-3)*x*pl[i-1]-(i-2)*pl[i-2])/(i-1);
return;
}
G4double G4PreEquilibrium::barfit( const G4double a, const G4double z, const G4int il )
{
// This subroutine returns the barrier height, bfis.
// Arguments: z, the atomic number
// a, the atomic mass number
// il, the angular momentum in units of h-bar,
// where h-bar is Plancks constant divided by 2*pi
//
// The fission barrier for il = 0 is calculated from a 7th order
// fit in two variables to 638 calculated fission barriers for z values
// from 20 to 110. These 638 barriers are fit with an rms deviation of
// 0.10 MeV by this 49-parameter function.
// If barfit is called with (z,a) values outside the range of the fit,
// the barrier height is set to 0
//
// For il values not equal to zero, the values of il at which the barrier
// is 80% and 20% of the il = 0 value are respectively fit to 20-parameter
// functions of z and a, over a more restricted range of a values, than is
// the case for il = 0. The value of il where the barrier disappears, lmax,
// for 61 nuclei, is fit to a 35-parameter function of z and a, with the
// same range of z and a values as il-80 and il-20.
// Once again, if a (z,a) pair is outside of the range of validity of the
// fit, the barrier value is set to 0. These three values
// (bfis(il=0), il-80, and l-20) and the constraints of bfis = 0 and
// d(bfis)/dl = 0 at il = lmax and il = 0 lead to a fifth-order fit to
// bfis(il) for il > l-20. the first three constraints lead to a third-order
// fit for the region il < il-20.
//
// The ground-state energies are calculated from a 175-parameter
// fit in z, a, and il to 329 ground-state energies for 36 different
// z and a values.
// (the range of z and a is the same as for l-80, l-20, and l-max)
//
// The calculated barriers from which the fits were made were calculated
// in 1983-1985 by A. J. Sierk of Los Alamos National Laboratory, group t-9,
// using Yukawa-plus-exponential double folded nuclear energy, exact Couloub
// diffuseness corrections, and diffuse-matter moments of inertia. The
// parameters of the model are those derived by Moller and Nix in 1979:
// r-0 = 1.16 fm, as = 21.13 MeV, kappa-s = 2.3 a = 0.68 fm.
// The diffuseness of the matter and charge distributions used
// corresponds to a surface diffuseness parameter (defined by Myers)
// of 0.99 fm. The calculated barriers for il = 0 are
// accurate to a little less than 0.1 MeV; the output from this
// routine is a little less accurate. Worst errors may be as large
// as 0.5 MeV; characteristic uncertainty is in the range of 0.1-0.2
// MeV. The values of egs are generally approximated to within
// about 0.1-0.2 MeV; the largest deviation is about 0.5 MeV,
// near il-i for light nuclei.
//
// The rms deviation of lmax from the 61 input values is 0.31
// h-bar. the approximate value is nearly always within
// 0.5 h-bar of the calculated one.
//
// Below is a table of test values to check implementation of the program
// z, a, l egnd st fiss bar moments of inertia lmax
//
// 28, 58, 0 0.00 33.14 0.816 3.603 3.608 46.1
// ,25 21.36 19.50 0.778 3.662 3.662 46.1
// ,40 49.66 2.97 0.724 3.648 3.650 46.1
// ,46.1 59.14 0.00 0.746 3.160 3.160 46.1
// 65,153, 0 0.00 28.88 0.621 3.698 3.698 82.3
// ,50 19.00 16.16 0.615 3.639 3.639 82.3
// ,80 45.24 0.26 0.616 2.765 2.788 82.3
// ,82.3 47.04 0.00 0.682 2.231 2.276 82.3
// 93,229, 0 0.00 3.76 0.715 1.747 1.747 68.1
// ,45 8.21 1.26 0.765 1.578 1.578 68.1
// ,68.1 17.96 0.00 1.053 1.053 1.236 68.1
//
// written by A. J. Sierk, LANL t-9
// version 1.0 February, 1984
// version 1.1 January, 1985 improved coefficients in egs and lmax
// version 1.2 September, 1985 improved lmax, egs coefficients
// version 1.21 June, 1986 minor changes made
//
// Copyright, 1986, the regents of the University of California.
// This software was produced under a U. s. government contract
// (w-7405-eng-36) by the Los Alamos National Laboratory, which is
// operated by the University of California for the U. S. Department
// of Energy. The U. S. government is licensed to use, reproduce,
// and distribute this software. Permission is granted to the public
// to copy and use this software without charge, provided that this
// notice and any statement of authorship are reproduced on all
// copies. Neither the government nor the University makes any
// warranty, expressed or implied, or assumes any liability
// or responsibility for the use of this software.
//
G4double pa[7], pz[7], pl[10];
//
const G4double emncof[4][5] =
{ { -901.100, -1408.18, 2770.00, -706.695, 889.867 },
{ 13535.5, -20384.7, 10938.4, -4862.97, -618.603 },
{ -3263.67, 1624.47, 1368.56, 1317.31, 153.372 },
{ 7488.63, -12158.1, 5502.81, -1336.30, 0.0505367 } };
const G4double elmcof[4][5] =
{ { 1845.42, -5640.02, 5667.30, -3151.50, 954.160 },
{ -2245.77, 8561.33, -9673.48, 5817.44, -1869.97 },
{ 2797.72, -8730.73, 9197.06, -4919.00, 1372.83 },
{ -30.1866, 1411.61, -2859.19, 2130.16, -649.072 } };
const G4double emxcof[5][7] = {
{-4106527.32, 10006494.7,-10953375.1, 7847972.52,-3785749.26, 1122379.45,-177561.170},
{10876333.0,-26375824.5, 28547240.0,-20110746.7, 9483736.41,-2734385.28, 413247.256},
{-8765309.03, 21425051.3,-23579959.5, 17016134.7,-8237381.90, 2424479.57,-365427.239},
{6302589.54,-15299900.4, 16564020.0,-11669577.6, 5473691.53,-1549863.42, 215409.246},
{-1455398.91, 3649618.35,-4212674.23, 3243125.55,-1679279.04, 523795.062,-76657.6599}};
const G4double elzcof[7][7] = {
{511819.909,-1303031.86, 1901198.70,-1206282.42, 568208.488, 54834.6483,-24588.3052},
{-1132694.53, 2977645.90,-4543263.26, 3004648.70,-1449892.74,-102026.610, 62795.9815},
{1375433.04,-3658089.88, 5477989.99,-3781092.83, 1841317.65, 15366.9695,-69681.7834},
{-856559.835, 2488722.66,-4073491.28, 3128358.99,-1623940.90, 119797.378, 42573.7058},
{328723.311,-1098921.75, 2039972.69,-1771857.18, 996051.545,-153305.699,-11298.2954},
{41585.0238, 72965.3408,-493776.346, 601254.680,-401308.292, 96596.8391,-349960.27},
{-182751.044, 391386.300,-303639.248, 115782.417,-424992.80,-611772.47, 36682.647}};
G4double bfis = 0;
if( G4int(z) < 19 || G4int(z) > 111 )
{
// write(16,1000)
// 1000 format(/10x,'* * * * barfit called with z less than 19 or ',
// ' greater than 111. bfis is set to 0.0. * * * *')
return bfis;
}
if( G4int(z) > 102 && il > 0 )
{
// write(16,1010)
// 1010 format(/10x,'* * * * barfit called with z greater than 102',
// ' and l not equal to zero. bfis is set to 0.0. * * * *')
return bfis;
}
if( a < 1.2*z+0.01*z*z || a > 5.8*z-0.024*z*z )
{
// write(16,1020)a
// 1020 format(/10x,'* * * * barfit called with a = ',i3,', outside ',
// 'the allowed values for z = ',i3,' * * * *')
return bfis;
}
lpoly( 0.01*z, 7, pz );
lpoly( 0.0025*a, 7, pa );
G4int i, j;
for( i = 0; i < 7; ++i )
{
for( j = 0; j < 7; ++j )bfis += elzcof[j][i]*pz[j]*pa[i];
}
if( (a < 1.4*z+0.009*z*z-5 || a > 20+3*z+10) && il > 0 )
{
// write(16,1030)a,il
// 1030 format(/10x,'* * * * barfit called with a = ',i3,', outside',
// ' the allowed values for z = ',i3/26x,'for nonzero l = ',i3,
// ' * * * *')
bfis = 0;
return bfis;
}
G4double el80 = 0;
G4double el20 = 0;
G4double elmax = 0;
for( i = 0; i < 4; ++i )
{
for( j = 0; j < 5; ++j )
{
el80 += elmcof[j][i]*pz[j]*pa[i];
el20 += emncof[j][i]*pz[j]*pa[i];
}
}
for( i = 0; i < 5; ++i )
{
for( j = 0; j < 7; ++j )elmax += emxcof[j][i]*pz[j]*pa[i];
}
if( il < 1 )return bfis;
G4double x = el20/elmax;
G4double y = el80/elmax;
if( il <= el20 )
{
G4double q = 0.2/(el20*el20*el80*el80*(el20-el80));
bfis *=
1+q*(4*el80*el80*el80-el20*el20*el20)*il*il-q*(4*el80*el80-el20*el20)*il*il*il;
}
else
{
G4double aj = (-20*pow(x,5)+25*pow(x,4)-4)*(y-1)*(y-1)*y*y;
G4double ak = (-20*pow(y,5)+25*pow(y,4)-1)*(x-1)*(x-1)*x*x;
G4double q = 0.2/((y-x)*((1-x)*(1-y)*x*y)*((1-x)*(1-y)*x*y));
G4double z = il/elmax;
bfis *= 4*pow(z,5)-5*pow(z,4)+1+
(z-1)*(q*(aj*y-ak*x)*(2*z+1)-q*(aj*(2*y+1)-ak*(2*x+1))*z)*z*z*(z-1);
}
if( bfis <= 0 || il > elmax )bfis = 0;
return bfis;
}
G4bool G4PreEquilibrium::cascem( G4ThreeVector amnucl )
{
const G4int N = 7;
const G4double EPS = 0.007;
const G4double VPI = 0.025;
const G4double CM0 = 0.94;
const G4int ME0 = 1;
const G4int MQ0 = 1;
G4int mv = 0;
G4ThreeVector v;
G4double tin1, sigp, sign, sigabs;
SetINDI( false );
SetING( 0 );
exitons.protons = 0;
exitons.neutrons = 0;
exitons.hols = 0;
//
// calculation of entering point of particle in nucleus
//
SetT1( 0 );
SetT2( 0 );
G4double temp1 = G4UniformRand();
G4double temp2 = 2*pi*G4UniformRand();
//
G4ThreeVector tmp3vec( rbig[N] * sqrt(temp1) * cos(temp2),
rbig[N] * sqrt(temp1) * sin(temp2),
-rbig[N] * sqrt(1-temp1) );
MASTRUCT p = { tmp3vec, { 0.0, 1.0, 0.0, 1.0, GetT0(), CM0 } };
//
IPSTRUCT ip = { ME0, 0, 0, MQ0, N };
//
// the following appears to be a bug in the original cem.f
// T1 is set to 0 just above
//
if( GetT0() == 0 )p.array[2] = GetT1()*pow(GetT2()/GetT1(),G4UniformRand());
const G4double p0 = sqrt(p.array[2]*(p.array[2]+2*p.array[5])); // not the data member p0
//
G4ThreeVector am0( p0*(p.momentum.z()*p.array[0]*p.array[2]-p.momentum.y()*p.array[1]),
p0*(p.momentum.x()*p.array[1]-p.momentum.z()*p.array[0]*p.array[3]),
p0*(p.momentum.y()*p.array[0]*p.array[3]-p.momentum.x()*p.array[0]*p.array[2]) );
//
const G4double t3 = p.array[4];
p.array[4] += poten( N, ip );
theNucleus->AddExcitationEnergy( -theNucleus->GetEnergyDeposit() + t3 + MQ0*EPS );
G4double atwght = aNucl+MQ0;
G4double charge = zNucl+ME0;
amnucl = am0;
const G4double obr = af[6]/2;
LABEL:
G4double cutof1 = ip.nuclearZone*( poten(ip.ip4,ip) + ip.proton*obr ) + 0.001;
if( p.array[4] <= cutof1 )
{
exitons.protons += ip.proton;
exitons.neutrons += 1-ip.proton;
cascem1( amnucl, p, ip, mv );
goto LABEL;
}
G4int ipe[5];
G4double pe[9];
if( pointe( p, ip, pe, ipe, v, u, tin1, sigp, sign, sigabs, t3 ) )
{
G4int np;
if( ip.nuclearZone != 0 )
{
G4double cutof3 = cutof1-ip.nuclearZone*obr+0.05;
if( p.array[4] <= cutof3 )
{
G4double r = rsm[N]*p.momentum.mag();
if( rsm[N] >= r+0.001 )
{
G4double wm = 0;
for( G4int i = 0; i < 50; ++i )
wm += wim( 1, p, ip ) + wim( 0, p, ip );
G4double w3 = wopt( p.array[4]-cutof3+0.05, r, ip.proton );
if( abs((wm/100-w3)/w3) >= 0.3 )
{
typint( p, ip, pe, ipe, v, u, tin1, sigp, sign, sigabs, mv, np );
if( np <= 0 )
{
exitons.protons += ip.proton;
exitons.neutrons += 1-ip.proton;
cascem1( amnucl, p, ip, mv );
}
else
{
pauliPrinciple( p, ip, v, mv, np, 0 );
exitons.protons += ip.proton;
exitons.neutrons += 1-ip.proton;
cascem1( amnucl, p, ip, mv );
}
goto LABEL;
}
}
}
}
typint( p, ip, pe, ipe, v, u, tin1, sigp, sign, sigabs, mv, np );
if( np <= 0 )
{
exitons.protons += ip.proton;
exitons.neutrons += 1-ip.proton;
cascem1( amnucl, p, ip, mv );
}
else
pauliPrinciple( p, ip, v, mv, np, 1 );
goto LABEL;
}
ip.proton <= 0 ? temp1 = 0 : temp1 = ip.proton;
G4double cutof2 = poten( ip.ip4, ip ) + temp1*obr + 0.001;
if( p.array[4] <= cutof2 )
{
exitons.protons += ip.proton;
exitons.neutrons += 1-ip.proton;
cascem1( amnucl, p, ip, mv );
goto LABEL;
}
theNucleus->AddExcitationEnergy(
-(p.array[4]+ip.nuclearZone*EPS+(1-ip.nuclearZone)*(0.14+VPI)*(1-ip.ip1)) );
atwght -= ip.nuclearZone;
charge -= ip.proton;
temp1 = sqrt(p.array[4]*(p.array[4]+2*p.array[5]));
amnucl.setX( amnucl.x()-temp1*(p.momentum.z()*p.array[0]*p.array[2]-p.momentum.y()*p.array[1]) );
amnucl.setY( amnucl.y()-temp1*(p.momentum.x()*p.array[1]-p.momentum.z()*p.array[0]*p.array[3]) );
amnucl.setZ(
amnucl.z()-temp1*(p.momentum.y()*p.array[0]*p.array[3]-p.momentum.x()*p.array[0]*p.array[2]) );
if( theNucleus->GetEnergyDeposit() <= 0.0001 &&
atwght == aNucl && charge == zNucl )return true;
if( theNucleus->GetEnergyDeposit() < 0 )return true;
if( KTOT > 100 )
{
//write(16,27)
//27 format (40x,32hmassiv spt exceeded after cascad)
return false;
}
Dvector vtmp( 5 );
vtmp[0] = p.array[0];
vtmp[1] = p.array[1];
vtmp[2] = p.array[4];
vtmp[3] = ip.proton;
vtmp[4] = p.array[5];
spt[KTOT]->swap( vtmp );
PARZS *pztmp = new PARZS;
pztmp->charge = ip.proton;
pztmp->what = GetING();
if( p.array[5] <= 0.5 )
{
if( ip.proton < 0 )
pztmp->particleType = 7;
else if( ip.proton == 0 )
pztmp->particleType = 8;
else
pztmp->particleType = 9;
}
else
{
if( ip.proton <= 0.099 )pztmp->particleType = 1; // neutron
else pztmp->particleType = 2; // proton
if( GetINDI() )pztmp->what = -pztmp->what;
}
pztmp->angle1 = atan2( p.array[0], p.array[1] );
temp1 = atan2( p.array[2], p.array[3] );
if( temp1 < 0 )temp1 = pi - temp1;
pztmp->angle2 = temp1;
pztmp->energy = p.array[4];
parz.push_back( *pztmp );
++KTOT;
cascem1( amnucl, p, ip, mv );
goto LABEL;
}
void G4PreEquilibrium::cascem1( G4ThreeVector a, MASTRUCT &p, IPSTRUCT &ip, G4int &m )
{
const G4int N = 7;
if( m <= 0 )
{
G4double temp = 5.07*rsm[N];
a *= temp;
}
else
{
G4int i;
p.momentum = pmemo[m].momentum;
for( i = 0; i < 6; ++i )p.array[i] = pmemo[m].array[i];
ip.proton = imemo[m].proton;
ip.ip1 = imemo[m].ip1;
ip.ip2 = imemo[m].ip2;
ip.nuclearZone = imemo[m].nuclearZone;
ip.ip4 = imemo[m].ip4;
SetING( ngen[m--] );
}
return;
}
G4bool G4PreEquilibrium::pointe( MASTRUCT &p, IPSTRUCT &ip,
G4double *pe, G4int *ipe,
G4ThreeVector v, G4double u,
G4double tin1, G4double sigp,
G4double sign, G4double sigabs, const G4double t3 )
{
// determine interaction point inside nucleus
//
const G4int N = 7;
G4bool nout;
const G4double sk = 3;
G4double temp1 = log(G4UniformRand());
do
{
G4double s = geometricalParticlePath( p, ip );
do
{
partnerSelection( p, ip, pe, ipe );
//
// calculation of t, v, u
//
G4double pin = sqrt( p.array[4]*(p.array[4]+2*p.array[5]) );
G4double pn = sqrt( pe[7]*(pe[7]+2*pe[8]) );
G4double denom = p.array[4]+p.array[5]+pe[7]+pe[8];
v.setX( (pin*p.array[0]*p.array[3]+pn*pe[3]*pe[6])/denom );
v.setY( (pin*p.array[0]*p.array[2]+pn*pe[3]*pe[5])/denom );
v.setZ( (pin*p.array[1]+pn*pe[4])/denom );
u = sqrt( 1-v.mag2() )*denom;
tin1 = (u*u-(p.array[5]+pe[8])*(p.array[5]+pe[8]))/(2*pe[8]);
//
G4int i, ms, mq, ksin, ksip, me;
slqek( i, ms, mq, ksin, me, ip.ip1, ip.ip2, ip.nuclearZone,
ip.proton, ipe[1], ipe[2], ipe[3], 0 );
sign = sigmat( i, ms, mq, ksin, 0, tin1 );
slqek( i, ms, mq, ksip, me, ip.ip1, ip.ip2, ip.nuclearZone,
ip.proton, ipe[1], ipe[2], ipe[3], 1 );
sigp = sigmat( i, ms, mq, ksip, 0, tin1 );
sigabs = sigmat( i, ms, mq, ksip, 3, p.array[4] );
G4double plambi = 10/((rhon[ip.ip4]*sign+rhop[ip.ip4]*sigp+
rhop[ip.ip4]*sigabs)*rsm[N-1]);
if( -plambi*temp1-6*rhon[0]/rhon[N-1] > 0 )
{
p.array[4] = t3;
return false;
}
G4double temp2 = sk;
G4double deltsi = plambi/temp2;
G4double temp5 = min( s, deltsi );
G4double piks = temp1+temp5/plambi;
if( piks > 0 )
{
G4ThreeVector tmp3vec( plambi*temp1*p.array[0]*p.array[3],
plambi*temp1*p.array[0]*p.array[2],
plambi*temp1*p.array[1] );
p.momentum -= tmp3vec;
pe[0] = p.momentum.x();
pe[1] = p.momentum.y();
pe[2] = p.momentum.z();
return nout;
}
temp1 = piks;
s -= temp5;
G4ThreeVector tmp3vec( temp5*p.array[0]*p.array[3],
temp5*p.array[0]*p.array[2],
temp5*p.array[1] );
p.momentum += tmp3vec;
} while ( s != 0 );
nout = refrac( N, p, ip );
} while ( ip.ip4 < N+1 );
G4double temp3 = geometricalParticlePath( p, ip );
G4ThreeVector tmp3vec( temp3*p.array[0]*p.array[3],
temp3*p.array[0]*p.array[2],
temp3*p.array[1] );
p.momentum += tmp3vec;
nout = refrac( N, p, ip );
return false;
}
G4double G4PreEquilibrium::geometricalParticlePath( const MASTRUCT p, const IPSTRUCT ip )
{
G4double rin = p.momentum.mag();
G4double c = (p.momentum.x()*p.array[0]*p.array[3] +
p.momentum.y()*p.array[0]*p.array[2] +
p.momentum.z()*p.array[1])/rin;
G4double g, r;
if( c < 0 )
{
g = -1;
ip.ip4 < 2 ? r = 0 : r = rbig[ip.ip4-2];
}
else
{
g = 1;
r = rbig[ip.ip4-1];
}
G4double d = r*r - rin*rin*(1-c*c);
if( d < 0 )
{
g = 1;
r = rbig[ip.ip4-1];
d = r*r - rin*rin*(1-c*c);
}
return g*sqrt(d) - rin*c;
}
G4bool G4PreEquilibrium::refrac( const G4int N, MASTRUCT &p, IPSTRUCT &ip ) const
{
// calculation of energy and direction change of
// particle while entering another zone of the nucleus
//
G4int j;
G4bool result;
G4double tmp = p.momentum.x()*p.array[0]*p.array[3] +
p.momentum.y()*p.array[0]*p.array[2] +
p.momentum.z()*p.array[1];
tmp < 0 ? j = ip.ip4-1 : j = ip.ip4+1;
p.array[4] += poten( j, ip ) - poten( ip.ip4, ip );
ip.ip4 = j;
j >= N+2 ? result = false : result = true;
return result;
}
G4double G4PreEquilibrium::wopt( const G4double e1, const G4double r, const G4int i )
{
const G4double a = aNucl;
const G4double z = zNucl;
G4double x = (a-2*z)/a;
G4double e = e1*GeV;
G4double rm, am, wv, wsf; // this rm is not the data member rm
if( i <= 0 )
{
if( e <= 25 )
{
rm = 1.26*pow(a,oneThird);
am = 0.58;
wv = 0.22*e - 1.56;
wsf = 13 - 0.25*e - 12*x;
}
else
{
rm = 1.21*pow(a,oneThird);
am = 0.6448;
wv = 0.459 + 0.111*e;
wsf = 4.28 - 0.0414*e;
}
}
else
{
if( e <= 25 )
{
rm = 1.32*pow(a,oneThird);
am = 0.51 + 0.7*x;
wv = 0.22*e - 2.7;
wsf = 11.8 - 0.25*e + 12*x;
}
else
{
rm = 1.37*pow(a,oneThird);
am = 0.74 - 0.008*e + x;
wv = 1.2 + 0.09*e;
wsf = 4.2 - 0.05*e + 15.5*x;
}
}
if( wv < 0 )wv = 0;
if( wsf < 0 )wsf = 0;
G4double f = 1/(1+exp((r-rm)/am));
return f*(wv+4*wsf*(1-f));
}
void G4PreEquilibrium::typint( MASTRUCT &p, IPSTRUCT &ip,
G4double *pe, G4int *ipe, G4ThreeVector v,
G4double u, G4double tin1, G4double sigp,
G4double sign, G4double sigabs, G4int mv, G4int &np )
{
// block of determining of interaction type and calculation
// of secondary particles characteristics
//
G4bool nin = true;
while ( nin )
{
G4int mtemp = ip.ip4-1;
G4double betabs = rhop[mtemp]*sigabs/
(rhop[mtemp]*sigp+rhon[mtemp]*sign+rhop[mtemp]*sigabs);
if( G4UniformRand() <= betabs )
{
absorption( p, ip, pe, mv, np, v );
return;
}
G4int i, ms, mq, ksi, me;
slqek( i, ms, mq, ksi, me, ip.ip1, ip.ip2, ip.nuclearZone,
ip.proton, ipe[1], ipe[2], ipe[3], ipe[0] );
G4double betael = (sigmat( i, ms, mq, ksi, 1, tin1 )+
sigmat( i, ms, mq, ksi, 2, tin1 ))/
sigmat( i, ms, mq, ksi, 0, tin1 );
if (G4UniformRand() <= betael )
{
np = elasticScattering( v, u, tin1, p, ip, ipe, mv, i, ms, mq, ksi, me );
return;
}
nin = inelasticScattering( p, ip, ipe, i, ms, mq, ksi, me, v, u, tin1, mv, np );
}
return;
}
G4int G4PreEquilibrium::elasticScattering( const G4ThreeVector v, const G4double u,
const G4double tin, const MASTRUCT p,
const IPSTRUCT ip, const G4int *ipatne,
const G4int mv, const G4int i,
const G4int ms, const G4int mq,
const G4int ksi, const G4int me )
{
// calculate particle characteristics in elastic and charge exchange scattering
//
G4double cmi, cmn;
if( ip.ip1 != 0 )
{
cmi = 0.14;
cmn = 0.94;
}
else
{
cmi = p.array[5];
cmn = 0.94;
}
G4double tmp1 = sigmat( i, ms, mq, ksi, 2, tin );
G4double tmp2 = sigmat( i, ms, mq, ksi, 1, tin );
G4int ie, ne;
G4double ctsti;
if( G4UniformRand() >= tmp1/(tmp1+tmp2) )
{
ie = ip.proton;
ne = ipatne[0];
ctsti = cosel( i, mq, ksi, tin, p.array[5] );
}
else
{
if( ip.proton != 0 )
{
ie = 0;
ne = me;
}
else
{
ne = 1 - ipatne[0];
ie = me - ne;
}
ctsti = cosex( i, tin, p.array[5] );
}
G4ThreeVector pist, pnst;
momentaCalc( p, v, u, pist, pnst, ctsti, cmi, cmn );
//
MASTRUCT p1 = { pist, { 0.0, 0.0, 0.0, 0.0, 0.0, cmi } };
pmemo[mv+2] = p1;
imemo[mv+2].proton = ie;
imemo[mv+2].ip1 = 0;
imemo[mv+2].ip2 = ip.ip2;
imemo[mv+2].nuclearZone = ip.nuclearZone;
MASTRUCT p2 = { pnst, { 0.0, 0.0, 0.0, 0.0, 0.0, cmn } };
pmemo[mv] = p2;
imemo[mv].proton = ne;
imemo[mv].ip1 = 0;
imemo[mv].ip2 = 0;
imemo[mv].nuclearZone = 1;
return 2;
}
G4bool G4PreEquilibrium::inelasticScattering( MASTRUCT &p, IPSTRUCT &ip, G4int *ipatne,
G4int l, G4int ms, G4int mq, G4int ksi,
G4int me, G4ThreeVector v, G4double u,
G4double tin1, G4int mv, G4int np )
{
if( ip.ip1 > 0 )
{
G4double betais = crossSectionInterp(tin1,25)/crossSectionInterp(tin1,26);
//
// there is a bug in the original cem.f routine
// isocem uses TINI as 3rd parameter
// it should be TIN1
//
if( G4UniformRand() <= betais )
isocem( u, v, tin1, p, ipatne, mv );
else
statisticalModel( u, v, p, ipatne, mv );
return true;
}
G4double betath;
if( tin1 >= 4 )
betath = 0;
else
betath = sigmat(l,ms,mq,ksi,7,tin1)/
(sigmat(l,ms,mq,ksi,0,tin1)-sigmat(l,ms,mq,ksi,1,tin1)-sigmat(l,ms,mq,ksi,2,tin1));
G4int ith;
G4double th;
if( G4UniformRand() < betath )
{
ith = 1;
th = 1;
}
else
{
ith = 0;
th = 2;
}
if( u > p.array[5]+0.14*th+0.96 )
{
G4int ik = 0;
G4int kp = 1;
while ( kp != 0 )
{
G4bool lp = vmnsp( p, ip, u, mv, np, ith, mq, tin1 );
if( np == 0 )return false;
if( lp )return true;
direction( v, u, tin1, mq, mv, np, p, kp, ith );
if( ++ik >= 50 )return true;
}
chinel( ip, l, ms, mq, ksi, np, mv, tin1, me, ipatne );
}
return true;
}
G4double G4PreEquilibrium::wim( const G4int i, const MASTRUCT p, const IPSTRUCT ip )
{
const G4double EPS = 0.007;
const G4int N = 7;
G4double r, r1, a1, fac, tf, rho;
r = rsm[N]*p.momentum.mag();
if( i <= 0 )
{
r1 = 1.26*pow(aNucl,oneThird);
a1 = 0.58;
fac = 1/(1+exp((r-r1)/a1));
tf = tfn[ip.ip4-1];
rho = rhon[0]*fac;
}
else
{
r1 = 1.32*pow(aNucl,oneThird);
a1 = 0.51+0.7*(aNucl-2*zNucl)/aNucl;
fac = 1/(1+exp((r-r1)/a1));
tf = tfp[ip.ip4-1];
rho = rhop[0]*fac;
}
G4double tp = tf*pow(G4UniformRand(),twoThirds);
SetT0( p.array[5] );
G4double z = abs(GetT0()-tf-EPS);
G4double y = sqrt(z*(z+2*p.array[5]));
G4double ct = 1.0 - 2*G4UniformRand();
G4double fi = 2*pi*G4UniformRand();
G4double pp = sqrt(tp*(tp+1.88));
G4double et = GetT0()+p.array[5]+tp+0.94;
G4double vx = pp*sqrt(1-ct*ct)*cos(fi)/et;
G4double vy = pp*sqrt(1-ct*ct)*sin(fi)/et;
G4double vz = (pp*ct+sqrt(GetT0()*(GetT0()+2*p.array[5])))/et;
G4double t = (et*et*(1-vx*vx-vy*vy-vz*vz)-(p.array[5]+0.94)*(p.array[5]+0.94))/1.88;
G4double b = sqrt(t*(t+2*p.array[5]))/(t+p.array[5]);
G4double s;
ip.proton == i ? s = 10.63/(b*b)-29.92/b+42.9 : s = 34.10/(b*b)-82.20/b+82.2;
G4double bs = tf/GetT0();
G4double dz = 1 - 1.4*bs;
if( bs > 0.5 )dz += 0.4*bs*pow(2-1/bs,2.5);
return 10*y/(z+p.array[5])*dz*rho*s;
}
void G4PreEquilibrium::pauliPrinciple( MASTRUCT &p, IPSTRUCT &ip, G4ThreeVector v,
G4int mv, G4int np, const G4int ipa )
{
G4int mtemp, i, j = ip.ip4, cntr = 1;
L10:
if( np == 2 )
{
if( cntr == 2 )
{
++cntr;
goto L10;
}
else mtemp = mv+cntr-1;
}
else mtemp = mv+cntr-1;
//
G4ThreeVector pstar( pmemo[mtemp].momentum );
//
G4double ct, st, cfi, sfi, t;
G4ThreeVector ptmp( p.momentum );
t = cinema( pstar, v, ptmp, ct, st, cfi, sfi, pmemo[mtemp].array[5] );
MASTRUCT tmp = { ptmp, { st, ct, sfi, cfi, t, 0.0 } };
pmemo.insert( pmemo.begin()+mtemp, tmp );
imemo[mtemp].ip4 = ip.ip4;
ngen.insert( ngen.begin()+mtemp, GetING()+1 );
if( imemo[mtemp].nuclearZone == 1 &&
imemo[mtemp].ip2 == 0 &&
pmemo[mtemp].array[4]-tfp[j]*imemo[mtemp].proton-tfn[j]*(1-imemo[mtemp].proton) <= 0 )
return;
if( cntr < np )
{
++cntr;
goto L10;
}
if( ipa == 0 )return;
if( p.array[5] == 0.140 && np == 2 &&
pmemo[mv].array[5] == 0.94 && pmemo[mv+2].array[5] == 0.94 )SetINDI( true );
for( i = 0; i < 6; ++i )p.array[i] = pmemo[mv+2].array[i];
ip.proton = imemo[mv+2].proton;
ip.ip1 = imemo[mv+2].ip2;
ip.ip2 = imemo[mv+2].ip2;
ip.nuclearZone = imemo[mv+2].nuclearZone;
SetING( ngen[mv+2] );
++(exitons.hols);
if( np > 2 )
{
G4int ntemp = mv+np;
for( i = 0; i < 6; ++i )pmemo[mv+2].array[i] = pmemo[ntemp].array[i];
imemo[mv+2].proton = imemo[ntemp].proton;
imemo[mv+2].ip1 = imemo[ntemp].ip1;
imemo[mv+2].ip2 = imemo[ntemp].ip2;
imemo[mv+2].nuclearZone = imemo[ntemp].nuclearZone;
}
mv += np-1;
return;
}
void G4PreEquilibrium::absorption( const MASTRUCT p, const IPSTRUCT ip,
G4double *pe, G4int mv, G4int &np, G4ThreeVector v )
{
// calculation of particle characteristics in absorption
//
G4double p1[9];
G4int ip1[5];
//
partnerSelection( p, ip, p1, ip1 );
G4double pn1 = sqrt( p1[7]*(p1[7]+1.88) );
G4double pn = sqrt( pe[7]*(pe[7]+1.88) );
G4ThreeVector paf( pn*pe[3]*pe[6]+pn1*p1[3]*p1[6],
pn*pe[3]*pe[5]+pn1*p1[3]*p1[5],
pn*pe[4]+pn1*p1[4] );
G4double pafm = paf.mag();
G4double ctn1, stn1, sfn1, cfn1;
if( paf.z()*paf.z() >= pafm*pafm )
{
ctn1 = 1;
stn1 = 0;
sfn1 = 0;
cfn1 = 1;
}
else
{
ctn1 = paf.z()/pafm;
stn1 = sqrt(1-paf.z()/pafm*paf.z()/pafm);
sfn1 = paf.y()/(pafm*stn1);
cfn1 = paf.x()/(pafm*stn1);
}
G4double taf = sqrt(pafm*pafm+1.88*1.88)-1.88;
//
// calculation of t, v, u
//
G4double pin = sqrt( p.array[4]*(p.array[4]+2*p.array[5]) );
pn = sqrt( taf*(taf+2*1.88) );
G4double denom = p.array[4]+p.array[5]+taf+1.88;
v.setX( (pin*p.array[0]*p.array[3]+pn*stn1*cfn1)/denom );
v.setY( (pin*p.array[0]*p.array[2]+pn*stn1*sfn1)/denom );
v.setZ( (pin*p.array[1]+pn*ctn1)/denom );
G4double u = sqrt( 1-v.mag2() )*denom;
G4double tin1 = (u*u-(p.array[5]+1.88)*(p.array[5]+1.88))/(2*1.88);
//
G4int ne, ie;
chargeInAbsorption( ip.ip1, ip.proton, ne, ie );
G4double ctst;
if( ip.ip1 == 0 )
ctst = 1-2*G4UniformRand();
else
tin1 <= 0.455 ? ctst = costa( 18, tin1 ) : ctst = 1-2*G4UniformRand();
G4ThreeVector pist, pnst;
momentaCalc( p, v, u, pist, pnst, ctst, 0.94, 0.94 );
if( mv <= 97 )
{
MASTRUCT pm1 = { pist, { 0.0, 0.0, 0.0, 0.0, 0.0, p1[8] } };
pmemo.insert( pmemo.begin()+mv+2, pm1 );
imemo[mv+2].proton = ie;
imemo[mv+2].ip1 = 0;
imemo[mv+2].ip2 = 0;
imemo[mv+2].nuclearZone = 1;
MASTRUCT pm2 = { (-1)*pist, { 0.0, 0.0, 0.0, 0.0, 0.0, pe[8] } };
pmemo.insert( pmemo.begin()+mv, pm2 );
imemo[mv+2].proton = ne;
imemo[mv+2].ip1 = 0;
imemo[mv+2].ip2 = 0;
imemo[mv+2].nuclearZone = 1;
np = 2;
return;
}
np = 0;
//write(16,19)
//19 format(45x,29h memory is exceeded in cascad)
return;
}
void G4PreEquilibrium::chargeInAbsorption( const G4int m, const G4int ine,
G4int &ne1, G4int &ne2 )
{
// determine charge in absorption
//
const G4double a = aNucl;
const G4double z = zNucl;
if( ine < 0 )
{
if( m != 0 )
{
ne1 = 0;
ne2 = 0;
}
else
{
G4double temp2 = (z*(z-1))/(2*z*(a-z)+z*(z-1));
if( G4UniformRand() <= temp2 )
{
ne1 = 0;
ne2 = 1;
}
else
{
ne1 = 0;
ne2 = 0;
}
}
}
else if( ine == 0 )
{
if( m != 0 )
{
ne1 = 0;
ne2 = 1;
}
else
{
G4double temp3 = (2*z*(a-z))/(a*(a-1));
if( G4UniformRand() <= temp3 )
{
ne1 = 0;
ne2 = 1;
}
else
{
G4double temp4 = (z*(z-1))/(a*(a-1)) + temp3;
if( G4UniformRand() <= temp4 )
{
ne1 = 1;
ne2 = 1;
}
else
{
ne1 = 0;
ne2 = 0;
}
}
}
}
else
{
if( m != 0 )
{
ne1 = 1;
ne2 = 1;
}
else
{
G4double temp1 = (z*(a-z))/(z*(a-z)+(a-z)*(a-z-1)/2);
if( G4UniformRand() <= temp1 )
{
ne1 = 1;
ne2 = 1;
}
else
{
ne1 = 0;
ne2 = 1;
}
}
}
return;
}
void G4PreEquilibrium::momentaCalc( const MASTRUCT p, const G4ThreeVector v,
const G4double u, G4ThreeVector pist,
G4ThreeVector pnst, const G4double cti,
const G4double cmi, const G4double cmn )
{
// calculate momenta of secondary particles in centre of mass system
// for absorption and elastic scattering
//
G4double angle = 2*pi*G4UniformRand();
G4double pim = sqrt((u*u+cmi*cmi-cmn*cmn)/(2*u)*(u*u+cmi*cmi-cmn*cmn)/(2*u)-cmi*cmi);
G4double v2 = v.mag2();
G4double temp3 = sqrt(p.array[4]*(p.array[4]+2*p.array[5]));
G4double temp1 = (temp3*p.array[0]*p.array[3]*v.x()+temp3*p.array[1]*v.y()+
temp3*p.array[0]*p.array[2]*v.z())*(1/sqrt(1-v2)-1)/v2;
G4double temp2 = (p.array[4]+p.array[5])/sqrt(1-v2);
G4ThreeVector pins( temp3*p.array[0]*p.array[3]+v.x()*(temp1-temp2),
temp3*p.array[0]*p.array[2]+v.y()*(temp1-temp2),
temp3*p.array[1]+v.z()*(temp1-temp2) );
G4ThreeVector pii( pim*sqrt(1-cti*cti)*cos(angle),
pim*sqrt(1-cti*cti)*sin(angle),
pim*cti );
pist = rotation( pins, v, pii );
pnst = pist * (-1);
return;
}
void G4PreEquilibrium::partnerSelection( const MASTRUCT p, const IPSTRUCT ip,
G4double *pe, G4int *ipe )
{
G4double t = 0;
if( G4UniformRand() >= (aNucl-zNucl)/aNucl )
{
while( t <= 0 )
t = tfp[ip.ip4-1]*pow(G4UniformRand(),twoThirds);
ipe[0] = 1;
}
else
{
while( t <= 0 )
t = tfn[ip.ip4-1]*pow(G4UniformRand(),twoThirds);
ipe[0] = 0;
}
pe[4] = 1-2*G4UniformRand();
pe[3] = sqrt(1-pe[4]*pe[4]);
G4double phin = 2*pi*G4UniformRand();
pe[6] = cos(phin);
pe[5] = sin(phin);
pe[0] = p.momentum.x();
pe[1] = p.momentum.y();
pe[2] = p.momentum.z();
pe[8] = 0.940;
pe[7] = t;
ipe[4] = ip.ip4;
ipe[1] = 0;
ipe[2] = 0;
ipe[3] = 1;
return;
}
void G4PreEquilibrium::chinel( IPSTRUCT &ipatin, const G4int l, const G4int ms,
const G4int mq, const G4int ksi, const G4int np,
const G4int mv, const G4double tin1, const G4int me,
G4int *ipatne )
{
// determine secondary particles charges in inelastic scattering
//
if( np == 3 )
{
G4double spi0 = sigmat( l, ms, mq, ksi, 4, tin1 );
G4double sth = sigmat( l, ms, mq, ksi, 7, tin1 );
G4double ran = G4UniformRand();
if( ran >= spi0/sth )
{
if( ran >= (spi0+sigmat(l,ms,mq,ksi,5,tin1))/sth )
{
imemo[mv].proton = ipatne[1] - (ipatin.ip4-1)*ipatin.ip1;
imemo[mv+2].proton = (ipatin.ip4-1)*ipatin.ip1*ipatin.ip1 - ipatin.ip4*ipatin.ip1 + 1;
imemo[mv+1].proton = me - imemo[mv].proton - imemo[mv+2].proton;
}
else
{
imemo[mv].proton = 1 - ipatne[1];
imemo[mv+2].proton = ipatin.ip1;
imemo[mv+1].proton = me - imemo[mv].proton - imemo[mv+2].proton;
}
}
else
{
imemo[mv].proton = ipatne[1];
imemo[mv+1].proton = 0;
imemo[mv+2].proton = ipatin.ip1;
}
return;
}
G4int migq;
do
{
G4UniformRand() <= 0.5 ? imemo[mv].proton = 1 : imemo[mv].proton = 0;
if( mq >= 1 )
{
G4UniformRand() <= 0.5 ? imemo[mv+2].proton = 1 : imemo[mv+2].proton = 0;
}
G4int lambda = 0;
do
{
++lambda;
if( mq <= 1 || lambda != 2 )
{
G4double ran = G4UniformRand();
if( ran >= oneThird )
{
ran < twoThirds ? imemo[mv+lambda].proton = 0 : imemo[mv+lambda].proton = -1;
}
else imemo[mv+lambda].proton = 1;
}
}
while ( lambda < np );
migq = 0;
for( G4int i = 0; i < np; ++i )migq += imemo[mv+i].proton;
} while( me != migq );
return;
}
void G4PreEquilibrium::direction( G4ThreeVector v, const G4double u, const G4double tin1,
const G4int mq, const G4int mv, const G4int np,
const MASTRUCT partin, G4int &kp, const G4int ith )
{
// determine direction of secondary particle motion
//
G4int nd = 0;
G4int m1, m2;
kp = 0;
if( mq > 1 )
{
if( G4UniformRand() >= 0.5 )
{
if( G4UniformRand() >= 0.5 )
{
m1 = 0;
m2 = 1;
}
else
{
m1 = 1;
m2 = 2;
}
}
else
{
m1 = 0;
m2 = 2;
}
}
else
{
if( G4UniformRand() >= 0.5 )
{
m1 = 0;
if( G4UniformRand() >= 0.5 )m2 = 1;
else m2 = 2;
}
else
{
m1 = 1;
m2 = 2;
}
}
G4ThreeVector pl, pakv;
G4int lambda, m1temp, m2temp;
LABEL:
pakv.setX( 0.0 );
pakv.setY( 0.0 );
pakv.setZ( 0.0 );
m1temp = mv + m1;
m2temp = mv + m2;
lambda = 0;
do
{
if( lambda != m1 && lambda != m2 )
{
G4int ja = coefficientTypeA( ith, mq, lambda );
G4double ctl = costa( ja, tin1 );
G4double fl = 2*pi*G4UniformRand();
G4double stl = sqrt(1-ctl*ctl);
G4double temp1 = cos(fl);
G4double temp2 = sin(fl);
G4double temp3 = pmemo[mv+lambda].array[0];
G4ThreeVector tmp3vec( temp3*stl*temp1, temp3*stl*temp2, temp3*ctl );
MASTRUCT tmp = { tmp3vec, { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } };
pmemo.insert( pmemo.begin()+mv+lambda, tmp );
pakv += pmemo[mv+lambda].momentum;
}
} while( lambda++ < np );
G4double pakvm = pakv.mag();
if( np != 3 && pmemo[m1temp].momentum.z() >= pakvm+pmemo[m2temp].array[0] &&
pmemo[m1temp].momentum.z() <= abs(pakvm-pmemo[m2temp].array[0]) )
{
if( ++nd >= 100 )
{
kp = 2;
return;
}
goto LABEL;
}
G4double v2 = v.mag2();
G4double temp4 = sqrt(partin.array[4]*(partin.array[4]+2*partin.array[5]));
G4double spv = temp4*partin.array[0]*partin.array[3]*v.x() + temp4*partin.array[1]*v.z() +
temp4*partin.array[0]*partin.array[2]*v.y();
G4double temp5 = spv*(1/sqrt(1-v2)-1)/v2;
G4double temp6 = (partin.array[4]+partin.array[5])/sqrt(1-v2);
G4ThreeVector pin( temp4*partin.array[0]*partin.array[3] + v.x()*temp5 - v.x()*temp6,
temp4*partin.array[0]*partin.array[2] + v.y()*temp5 - v.y()*temp6,
temp4*partin.array[1] + v.z()*temp5 - v.z()*temp6 );
//
lambda = 1;
if( lambda == m1 && lambda < np )++lambda;
G4int ltemp = mv + lambda;
G4ThreeVector plst;
if( lambda == m1 )
{
G4ThreeVector pakst;
pakst = rotation( pin, v, pakv );
G4double ctm1 = (pmemo[m2temp].array[0]*pmemo[m2temp].array[0]-
pmemo[m1temp].array[0]*pmemo[m1temp].array[0]-pakvm*pakvm)/
(2*pakvm*pmemo[m1temp].array[0]);
G4double ctm2 = (pmemo[m1temp].array[0]*pmemo[m1temp].array[0]-
pmemo[m2temp].array[0]*pmemo[m2temp].array[0]-pakvm*pakvm)/
(2*pakvm*pmemo[m2temp].array[0]);
G4double fm1 = 2*pi*G4UniformRand();
G4double fm2 = pi+fm1;
G4double stm1 = sqrt(1-ctm1*ctm1);
G4double stm2 = sqrt(1-ctm2*ctm2);
G4double cfm1 = cos(fm1);
G4double sfm1 = sin(fm1);
G4double cfm2 = cos(fm2);
G4double sfm2 = sin(fm2);
pl.setX( pmemo[m1temp].array[0]*stm1*cfm1 );
pl.setY( pmemo[m1temp].array[0]*stm1*sfm1 );
pl.setZ( pmemo[m1temp].array[1]*ctm1 );
plst = rotation( pakst, v, pl );
pmemo[m1temp].momentum = plst;
pl.setX( pmemo[m2temp].array[0]*stm2*cfm2 );
pl.setY( pmemo[m2temp].array[0]*stm2*sfm2 );
pl.setZ( pmemo[m2temp].array[0]*ctm2 );
plst = rotation( pakst, v, pl );
pmemo[m2temp].momentum = plst;
}
else if( lambda != m2 )
{
pl = pmemo[ltemp].momentum;
plst = rotation( pin, v, pl );
pmemo[ltemp].momentum = plst;
}
return;
}
G4int G4PreEquilibrium::coefficientTypeA( const G4int ith, const G4int mq,
const G4int lamb )
{
// determine type of coefficients a(n,k)
//
G4int result;
if( ith != 0 )
{
if( mq <= 1 )
{
lamb > 1 ? result = 25 : result = 24;
}
else
{
if( lamb > 1 )
{
lamb == 3 ? result = 16 : result = 21;
}
else result = 20;
}
}
else
{
if( mq <= 1 )
{
lamb <= 1 ? result = 26 : result = 27;
}
else
{
if( lamb > 1 )
{
lamb == 3 ? result = 22 : result = 23;
}
else result = 22;
}
}
return result;
}
G4bool G4PreEquilibrium::vmnsp( const MASTRUCT partin, const IPSTRUCT ipatin, const G4double u,
const G4int mv, G4int &np, const G4int ith,
const G4int mq, const G4double tin1 )
{
// calculate secondary particle number and
// determine absolute values of momenta in inelastic interaction
//
G4int counter = 0;
LABEL1:
G4double u1 = u;
G4int lambda = 0;
LABEL2:
if( mv+lambda >= 100 )
{
np = 0;
//write(16,39)
//39 format(45x,29h memory is exceeded in cascad)
return false;
}
if( lambda == 0 )
{
G4ThreeVector tmp3vec( 0.0, 0.0, 0.0 );
MASTRUCT matmp = { tmp3vec, { 0.0, 0.0, 0.0, 0.0, 0.0, 0.94 } };
pmemo.insert( pmemo.begin()+mv+lambda, matmp );
IPSTRUCT iptmp = { 0, 0, 0, 1, 0 };
imemo.insert( imemo.begin()+mv+lambda, iptmp );
}
else if( lambda == 2 )
{
G4ThreeVector tmp3vec( 0.0, 0.0, 0.0 );
MASTRUCT matmp = { tmp3vec, { 0.0, 0.0, 0.0, 0.0, 0.0, partin.array[5] } };
pmemo.insert( pmemo.begin()+mv+lambda, matmp );
IPSTRUCT iptmp = { 0, 0, ipatin.ip2, ipatin.nuclearZone, 0 };
imemo.insert( imemo.begin()+mv+lambda, iptmp );
}
else
{
G4ThreeVector tmp3vec( 0.0, 0.0, 0.0 );
MASTRUCT matmp = { tmp3vec, { 0.0, 0.0, 0.0, 0.0, 0.0, 0.14 } };
pmemo.insert( pmemo.begin()+mv+lambda, matmp );
IPSTRUCT iptmp = { 0, 0, 0, 0, 0 };
imemo.insert( imemo.begin()+mv+lambda, iptmp );
}
G4int jb = coefficientTypeB( ith, mq, lambda );
pmemo[mv+lambda].array[0] = secondaryParticleMomentum(jb,tin1);
G4double el = sqrt(pmemo[mv+lambda].array[0]*pmemo[mv+lambda].array[0]+
pmemo[mv+lambda].array[5]*pmemo[mv+lambda].array[5]);
G4double deltu = u1-el;
if( lambda == 0 )
{
if( deltu <= partin.array[5] )
{
if( ++counter == 100 )return true;
goto LABEL1;
}
if( ith == 0 )
{
u1 = deltu;
++lambda;
goto LABEL2;
}
pmemo[mv+2].array[0] = sqrt(deltu*deltu-partin.array[5]*partin.array[5]);
pmemo[mv+2].array[5] = partin.array[5];
imemo[mv+2].ip1 = 0;
imemo[mv+2].ip2 = ipatin.ip2;
imemo[mv+2].nuclearZone = ipatin.nuclearZone;
if( pmemo[mv].array[0] > pmemo[mv+1].array[0]+pmemo[mv+2].array[0] ||
pmemo[mv].array[0] <= abs(pmemo[mv+1].array[0]-pmemo[mv+2].array[0]) )
{
if( ++counter == 100 )return true;
goto LABEL1;
}
np = 3;
return true;
}
if( deltu > 0.14 )
{
u1 = deltu;
++lambda;
goto LABEL2;
}
if( lambda <= 1 || lambda == 3 )
{
if( ++counter == 100 )return true;
goto LABEL1;
}
el += deltu;
pmemo[mv+lambda].array[0] = sqrt(el*el-pmemo[mv+lambda].array[5]*pmemo[mv+lambda].array[5]);
np = lambda;
G4double c = pmemo[mv].array[0];
G4int i;
for( i = 0; i < np; ++i )
{
if( c < pmemo[mv+i].array[0] )c = pmemo[mv+i].array[0];
}
G4double pmax = c;
G4double sigma = 0;
for( i = 0; i < np; ++i )sigma += pmemo[mv+i].array[0];
if( 2*pmax >= sigma )
{
if( ++counter == 100 )return true;
goto LABEL1;
}
return false;
}
G4double G4PreEquilibrium::secondaryParticleMomentum( const G4int j, const G4double t )
{
// calculation of secondary particle momentum
//
const G4double bnkj[8][4][4] =
{ { { 0.50278, 3.1442, -7.8172, 8.1667 },
{ 0.93482, -10.590, 29.227, -34.550 },
{-0.096685, 4.7335, -14.298, 17.685 },
{-0.025041, -0.62478, 2.0282, -2.5895 } },
{ { 1.1965, -0.82889, 1.0426, -1.9090 },
{ 0.28703, -4.9065, 16.264, -19.904 },
{-0.24492, 2.9191, -9.5776, 11.938 },
{ 0.037297, -0.42200, 1.3883, -1.7476 } },
{ { 1.3508, -4.3139, 12.291, -15.288 },
{-0.20086, 1.3641, -3.4030, 3.8559 },
{ 0.012583, -0.083492, 0.18600, -0.20043 },
{-0.00023628, 0.0013514, -0.0024324, 0.0021906 } },
{ { 1.2419, -4.3633, 13.743, -18.592 },
{-0.24404, 1.3158, -3.5691, 4.3867 },
{ 0.015693, -0.082579, 0.21427, -0.25846 },
{-0.00029386, 0.0014060, -0.0033835, 0.0038664 } },
{ { 0.63054, -3.7333, 13.464, -18.594 },
{ 2.1801, 1.5163, -16.380, 27.944 },
{-1.2886, -2.4570, 15.129, -23.295 },
{ 0.20915, 0.52279, -2.8687, 4.2688 } },
{ { 0.93363, -1.8181, 5.5157, -8.5216 },
{ 1.7811, -8.2927, 20.607, -20.827 },
{-1.5264, 6.8433, -16.067, 16.845 },
{ 0.27128, -1.1944, 2.7495, -2.9045 } },
{ { 1.9439, -4.6268, 9.7879, -9.6074 },
{-0.34640, 1.1093, -1.9313, 1.7064 },
{ 0.027054, -0.11638, 0.26969, -0.31853 },
{-0.00066092, 0.0050728, -0.014995, 0.019605 } },
{ { 1.8693, -5.5678, 14.795, -16.903 },
{-0.49965, 1.7874, -4.1330, 3.8393 },
{ 0.046194, -0.18536, 0.45315, -0.46273 },
{-0.0013341, 0.0057710, -0.014554, 0.015554 } } };
const G4double ckj[8][3] = {
0.14509, 0.46520, -0.033005, 0.15376,
0.27436, -0.014604, 0.62959, 0.17866,
-0.0026216, 0.83810, 0.0086137, 0.0032946,
0.092852, 0.53886, -0.054493, 0.13032,
0.40709, -0.028782, 0.14909, 0.38502,
-0.012775, 0.18024, 0.33022, -0.0094491 };
G4double bnk[4][4];
G4int n, k;
for( k = 0; k < 4; ++k )
{
for( n = 0; n < 4; ++n )bnk[n][k] = bnkj[n][k][j];
}
G4double s1 = 0;
G4double r1 = G4UniformRand();
G4double s2 = 0;
G4double s3 = 0;
for( n = 0; n < 4; ++n )
{
for( k = 0; k < 4; ++k )s1 += bnk[n][k]*pow(t,k)*pow(r1,n);
}
for( n = 0; n < 4; ++n )
{
for( k = 0; k < 4; ++k )s2 += bnk[n][k]*pow(t,k);
}
for( k = 0; k < 3; ++k )s3 += ckj[k][j]*pow(t,k);
return s3*sqrt(r1)*(s1+(1-s2)*r1*r1*r1*r1);
}
G4int G4PreEquilibrium::coefficientTypeB( const G4int i, const G4int mq,
const G4int lamb )
{
// determine type of coefficients b(n,k)
//
G4int result;
if( i != 0 )
{
if( mq <= 1 )
{
lamb > 1 ? result = 5 : result = 4;
}
else
{
if( lamb > 1 )
{
lamb == 3 ? result = 0 : result = 1;
}
else result = 0;
}
}
else
{
if( mq <= 1 )
{
lamb <= 1 ? result = 6 : result = 7;
}
else
{
if( lamb > 1 )
{
lamb == 3 ? result = 2 : result = 3;
}
else result = 2;
}
}
return result;
}
void G4PreEquilibrium::statisticalModel( const G4double u, const G4ThreeVector v,
const MASTRUCT partin, const G4int *ipatne,
const G4int mv )
{
// determine secondary particle characteristics for
// gamma-n interaction with statistical model
//
G4double tpim = (u*u+0.0196-1.08*1.08)/(2*u) - 0.14;
LABEL:
G4double t1 = G4UniformRand()*tpim;
G4double t2 = G4UniformRand()*tpim;
G4double e1 = t1 + 0.14;
G4double e2 = t2 + 0.14;
//
if( G4UniformRand() >= 27*e1*e2*(u-e1-e2)/(u*u*u) )goto LABEL;
//
G4double t3 = u - e1 - e2 - 0.94;
if( t3 <= 0 )goto LABEL;
//
G4double p1 = sqrt(t1*(t1+0.28));
G4double p2 = sqrt(t2*(t2+0.28));
G4double p3 = sqrt(t3*(t3+1.88));
//
if( (p1+p2-p3)*(p1-p2+p3)*(p2+p3-p1) <= 0 )goto LABEL;
//
G4double ct3 = 1-2*G4UniformRand();
G4double fi3 = 2*pi*G4UniformRand();
G4double temp2 = sqrt(1-ct3*ct3);
//
G4ThreeVector pv3( p3*temp2*cos(fi3),
p3*temp2*sin(fi3),
p3*ct3 );
G4double temp3 = sqrt(partin.array[4]*(partin.array[4]+2*partin.array[5]));
//
G4ThreeVector pin( temp3*partin.array[0]*partin.array[3],
temp3*partin.array[0]*partin.array[2],
temp3*partin.array[1] );
//
G4ThreeVector pinst = cms( pin, v, partin.array[4]+partin.array[5] );
//
G4ThreeVector ps3 = rotation( pinst, v, pv3 );
G4double ct1 = (p2*p2-p1*p1-p3*p3)/(2*p3*p1);
G4double ct2 = (p1*p1-p2*p2-p3*p3)/(2*p3*p2);
G4double fi1 = 2*pi*G4UniformRand();
G4double fi2 = pi + fi1;
G4double st1 = sqrt(1-ct1*ct1);
G4double st2 = sqrt(1-ct2*ct2);
//
G4ThreeVector pv1( p1*st1*cos(fi1),
p1*st1*sin(fi1),
p1*ct1 );
G4ThreeVector ps1 = rotation( ps3, v, pv1 );
//
G4ThreeVector pv2( p2*st2*cos(fi2),
p2*st2*sin(fi2),
p2*ct2 );
//
G4ThreeVector ps2 = rotation( ps3, v, pv2 );
pmemo[mv].momentum = ps1;
pmemo[mv].array[5] = 0.14;
pmemo[mv+1].momentum = ps2;
pmemo[mv+1].array[5] = 0.14;
pmemo[mv+2].momentum = ps3;
pmemo[mv+2].array[5] = 0.94;
if( ipatne[1] <= 0 )
{
G4double temp4 = G4UniformRand();
if( temp4 > oneThird)
{
if( temp4 <= twoThirds )
{
imemo[mv].proton = 0;
imemo[mv].ip1 = 0;
imemo[mv].ip2 = 0;
imemo[mv].nuclearZone = 0;
imemo[mv+1].proton = -1;
imemo[mv+1].ip1 = 0;
imemo[mv+1].ip2 = 0;
imemo[mv+1].nuclearZone = 0;
imemo[mv+2].proton = 1;
imemo[mv+2].ip1 = 0;
imemo[mv+2].ip2 = 0;
imemo[mv+2].nuclearZone = 1;
}
else
{
imemo[mv].proton = -1;
imemo[mv].ip1 = 0;
imemo[mv].ip2 = 0;
imemo[mv].nuclearZone = 0;
imemo[mv+1].proton = 1;
imemo[mv+1].ip1 = 0;
imemo[mv+1].ip2 = 0;
imemo[mv+1].nuclearZone = 0;
imemo[mv+2].proton = 0;
imemo[mv+2].ip1 = 0;
imemo[mv+2].ip2 = 0;
imemo[mv+2].nuclearZone = 1;
}
}
else
{
imemo[mv].proton = 0;
imemo[mv].ip1 = 0;
imemo[mv].ip2 = 0;
imemo[mv].nuclearZone = 0;
imemo[mv+1].proton = 0;
imemo[mv+1].ip1 = 0;
imemo[mv+1].ip2 = 0;
imemo[mv+1].nuclearZone = 0;
imemo[mv+2].proton = 0;
imemo[mv+2].ip1 = 0;
imemo[mv+2].ip2 = 0;
imemo[mv+2].nuclearZone = 1;
}
}
else
{
G4double temp4 = G4UniformRand();
if( temp4 > oneThird )
{
if( temp4 <= twoThirds )
{
imemo[mv].proton = 0;
imemo[mv].ip1 = 0;
imemo[mv].ip2 = 0;
imemo[mv].nuclearZone = 0;
imemo[mv+1].proton = 1;
imemo[mv+1].ip1 = 0;
imemo[mv+1].ip2 = 0;
imemo[mv+1].nuclearZone = 0;
imemo[mv+2].proton = 0;
imemo[mv+2].ip1 = 0;
imemo[mv+2].ip2 = 0;
imemo[mv+2].nuclearZone = 1;
}
else
{
imemo[mv].proton = 1;
imemo[mv].ip1 = 0;
imemo[mv].ip2 = 0;
imemo[mv].nuclearZone = 0;
imemo[mv+1].proton = -1;
imemo[mv+1].ip1 = 0;
imemo[mv+1].ip2 = 0;
imemo[mv+1].nuclearZone = 0;
imemo[mv+2].proton = 1;
imemo[mv+2].ip1 = 0;
imemo[mv+2].ip2 = 0;
imemo[mv+2].nuclearZone = 1;
}
}
else
{
imemo[mv].proton = 0;
imemo[mv].ip1 = 0;
imemo[mv].ip2 = 0;
imemo[mv].nuclearZone = 0;
imemo[mv+1].proton = 0;
imemo[mv+1].ip1 = 0;
imemo[mv+1].ip2 = 0;
imemo[mv+1].nuclearZone = 0;
imemo[mv+2].proton = 1;
imemo[mv+2].ip1 = 0;
imemo[mv+2].ip2 = 0;
imemo[mv+2].nuclearZone = 1;
}
}
return;
}
void G4PreEquilibrium::isocem( const G4double u, const G4ThreeVector v,
const G4double tin1, const MASTRUCT partin,
const G4int *ipatne, const G4int mv )
{
// interaction with (3/2,3/2) isobar production
//
//
G4double a1 = (u*u+0.0196-1.08*1.08)/(2*u);
G4double a2 = sqrt(a1*a1-0.0196);
G4double a3 = u-a1;
//
G4double epim, bms, pim, p;
do
{
bms = G4UniformRand()*(u-1.22)+1.08;
epim = (u*u+0.0196-bms*bms)/(2*u);
pim = sqrt(epim*epim-0.0196);
p = ((pim*epim*(u-epim))/u)*
crossSectionInterp((bms*bms-1.08*1.08)/1.88,8)/(200*a1*a2*a3/u);
}
while ( G4UniformRand() >= p );
//
G4double ctpi;
tin1 < 1 ? ctpi = costa( 27, tin1 ) : ctpi = costa( 28, tin1 );
G4double fipi = 2*pi*G4UniformRand();
G4double epit = (bms*bms+0.0196-0.94*0.94)/(2*bms);
G4double temp1 = sqrt(1-ctpi*ctpi);
G4ThreeVector ppim( pim*temp1*cos(fipi),
pim*temp1*sin(fipi),
pim*ctpi );
G4ThreeVector vt = ppim * (1/(epim-u));
G4double ctilpi = 1-2*G4UniformRand();
G4double ftilpi = 2*pi*G4UniformRand();
temp1 = sqrt(1-ctilpi*ctilpi);
G4double pmt = sqrt(epit*epit-0.0196);
G4ThreeVector ppit( pmt*temp1*cos(ftilpi),
pmt*temp1*sin(ftilpi),
pmt*ctilpi );
vt *= -1.0;
G4ThreeVector ppi = cms( ppit, vt, epit );
G4ThreeVector ppt = ppit * (-1.0);
G4ThreeVector pp = cms( ppt, vt, bms-epit );
G4double temp2 = sqrt(partin.array[4]*(partin.array[4]+2*partin.array[5]));
G4ThreeVector pin( temp2*partin.array[0]*partin.array[3],
temp2*partin.array[0]*partin.array[2],
temp2*partin.array[1] );
G4ThreeVector pinst = cms( pin, v, partin.array[4]+partin.array[5] );
G4ThreeVector ppimst;
ppimst = rotation( pinst, v, ppim );
pmemo[mv].momentum = ppimst;
pmemo[mv].array[5] = 0.14;
G4ThreeVector ppist;
ppist = rotation( pinst, v, ppi );
pmemo[mv+1].momentum = ppist;
pmemo[mv+1].array[5] = 0.14;
G4ThreeVector ppst;
ppst = rotation( pinst, v, pp );
pmemo[mv+2].momentum = ppst;
pmemo[mv+2].array[5] = 0.94;
if( ipatne[1] > 0 )
{
imemo[mv].proton = -1;
imemo[mv].ip1 = 0;
imemo[mv].ip2 = 0;
imemo[mv].nuclearZone = 0;
imemo[mv+1].proton = 1;
imemo[mv+1].ip1 = 0;
imemo[mv+1].ip2 = 0;
imemo[mv+1].nuclearZone = 0;
imemo[mv+2].proton = 1;
imemo[mv+2].ip1 = 0;
imemo[mv+2].ip2 = 0;
imemo[mv+2].nuclearZone = 1;
}
else
{
imemo[mv].proton = 1;
imemo[mv].ip1 = 0;
imemo[mv].ip2 = 0;
imemo[mv].nuclearZone = 0;
imemo[mv+1].proton = -1;
imemo[mv+1].ip1 = 0;
imemo[mv+1].ip2 = 0;
imemo[mv+1].nuclearZone = 0;
imemo[mv+2].proton = 0;
imemo[mv+2].ip1 = 0;
imemo[mv+2].ip2 = 0;
imemo[mv+2].nuclearZone = 1;
}
return;
}
G4ThreeVector G4PreEquilibrium::cms( const G4ThreeVector p, const G4ThreeVector v,
const G4double tcm ) const
{
// momentum calculation in system which has velocity v relative to given one
//
G4double temp1 = sqrt( 1 - v.mag2() );
G4double temp2 = p.dot(v)/v.mag2()*(1/temp1-1);
G4ThreeVector result( p.x() + v.x()*(temp2-tcm/temp1),
p.y() + v.y()*(temp2-tcm/temp1),
p.z() + v.z()*(temp2-tcm/temp1) );
return result;
}
void G4PreEquilibrium::bertcem( const G4int nuclideTypeNumber, const G4int mediumNumber )
{
// this routine was called by cascade and by cemgeo
//
aNucl = GetA_FCOMON(nuclideTypeNumber,mediumNumber);
zNucl = GetZ_FCOMON(nuclideTypeNumber,mediumNumber);
SetT0( 3.4995 );
//
// neutron density = proton density = nucleon density
//
const G4double a[10] = { 0.95, 0.8, 0.5, 0.2, 0.1, 0.05, 0.01, 0.0, 0.0, 0.0 };
const G4double r0n = 1.07;
const G4double rmax = 10;
const G4double bnr = 0.545;
//
G4double rn = r0n*pow(aNucl,oneThird);
G4int i;
//
const G4int N = 7;
for( i = 0; i < N; ++i )rsm.push_back( rn+bnr*(log((1-a[i])/a[i])) );
rsm.push_back( rmax );
//
G4double hin[10], fi[10];
fi[0] = fintfis( rsm[0], 0.0, bnr, rn );
hin[0] = fints2( rsm[0], 0.0, bnr, rn );
G4double sumhin = hin[0];
for( i = 1; i < N; ++i )
{
fi[i] = fintfis( rsm[i], rsm[i-1], bnr, rn );
hin[i] = fints2( rsm[i], rsm[i-1], bnr, rn );
sumhin += hin[i];
}
//
G4double rhon0 = aNucl/(4*pi*sumhin);
G4double fi0 = (4*pi*zNucl*0.00144*rhon0)/aNucl;
rhon.push_back( 3*hin[0]*rhon0/rsm[0]/rsm[0]/rsm[0] );
af.push_back( fi[0]*fi0*3/rsm[0]/rsm[0]/rsm[0] );
for( i = 1; i < N; ++i )
{
rhon.push_back( 3*hin[i]*rhon0/(rsm[i]*rsm[i]*rsm[i]-rsm[i-1]*rsm[i-1]*rsm[i-1]) );
tfp.push_back( 0.1985*pow( rhon[i]*zNucl/aNucl, twoThirds ) );
tfn.push_back( 0.1985*pow( rhon[i]*(aNucl-zNucl)/aNucl, twoThirds ) );
af.push_back( 3*fi0*fi[i]/(rsm[i]*rsm[i]*rsm[i]-rsm[i-1]*rsm[i-1]*rsm[i-1]) );
rbig.push_back( rsm[i]/rsm[N] );
}
rbig.push_back( rmax/rsm[N-1] );
af.push_back( zNucl*0.00144/rmax );
for( i = 0; i < N; ++i )
{
rhop.push_back( rhon[i]*zNucl/aNucl );
rhon[i] *= (aNucl-zNucl)/aNucl;
}
//
af.insert( af.begin(), 10, 0.0 );
//
G4ThreeVector amnucl;
if( !cascem( amnucl ) )
{
SetLXYZ( amnucl );
SetN0( exitons.protons + exitons.neutrons + exitons.hols );
SetH0( exitons.hols );
SetPZ0( exitons.protons );
SetP0( exitons.protons + exitons.neutrons );
}
return;
}
G4double G4PreEquilibrium::fintfis( const G4double flim1, const G4double flim2,
const G4double bs, const G4double rs )
{
const G4int N = 7;
const G4double w[8] = { 0.1012285363, 0.2223810345, 0.3137066459, 0.3626837834,
0.3626837834, 0.3137066459, 0.2223810345, 0.1012285363 };
const G4double fiks[8] = { 0.9602898565, 0.7966664774, 0.5255324099, 0.1834346425,
-0.1834346425, -0.5255324099, -0.7966664774, -0.9602898565 };
//
G4double temp = 0;
for( G4int i = 0; i < 8; ++i )
{
G4double fy = ((flim1-flim2)*fiks[i]+flim1+flim2)/2;
G4double tmp1 = 0.0;
G4double tmp2 = 0.0;
G4double tmp3 = 0.0;
for( G4int j = 0; j < 8; ++j )
{
fy = ((flim1-flim2)*fiks[j]+flim1+flim2)/2;
tmp1 += w[j]*(fy*fy)/(1+exp((fy-rs)/bs));
tmp2 += w[j]*rsm[N]/(1+exp((rsm[N]-rs)/bs));
tmp3 += w[j]*fy/(1+exp((fy-rs)/bs));
}
G4double f1 = tmp1*(flim1-flim2)/2;
G4double f2 = tmp2*(flim1-flim2)/2;
G4double f3 = tmp3*(flim1-flim2)/2;
temp += w[i] * (f1/fy + f2 - f3)*fy*fy;
}
return temp*(flim1-flim2)/2;
}
G4double G4PreEquilibrium::fints2( G4double flim1, G4double flim2,
G4double bs, G4double rs )
{
const G4double w[8] = { 0.1012285363, 0.2223810345, 0.3137066459, 0.3626837834,
0.3626837834, 0.3137066459, 0.2223810345, 0.1012285363 };
const G4double fiks[8] = { 0.9602898565, 0.7966664774, 0.5255324099, 0.1834346425,
-0.1834346425, -0.5255324099, -0.7966664774, -0.9602898565 };
G4double temp = 0;
for( G4int i = 0; i < 8; ++i )
{
G4double fy = ((flim1-flim2)*fiks[i]+flim1+flim2)/2;
temp += w[i] * (fy*fy)/(1+exp((fy-rs)/bs));
}
return temp*(flim1-flim2)/2;
}
G4double G4PreEquilibrium::arfaf( G4double *ami )
{
// renormalisation of emission probabilities
//
const G4double AMF = 0.0;
//
G4double sfix = 30;
G4double smx = 0;
for( G4int k = 0; k < 7; ++k )
{
G4double q8;
if( rj[k] <= 0 )q8 = 0;
else
{
G4double tmp;
k == 6 ? tmp = AMF : tmp = ami[k];
q8 = 2*sqrt(tmp*afj[k]*rj[k]);
}
smx = max( smx, q8 );
}
G4double result;
smx <= sfix ? result = 0.0 : result = smx-sfix;
return result;
}
G4double G4PreEquilibrium::fam( const G4double a, const G4double z, const G4double e )
{
//
// third set of Iljinov, Mebel et al. parameters: without collective effects
// Truran, Cameron & Hils shell corrections, without a dependence on ga
//
const G4double al = 0.072;
const G4double be = 0.257;
const G4double ga = 0.059;
// shell is an inline function defined in the header
return (al+be/pow(a,oneThird))*(1+(1-exp(-ga*e))*shell(a,z)/e);
}
G4double G4PreEquilibrium::gameqf( const G4int j, const G4double cc, const G4double per,
const G4double am, const G4double radncl )
{
//
// a possible bug in the original cem.f
// idel set to 0 in inpcem, and nowhere else
// idel checked in precof, bypassing setting amf
//
const G4double AMF = 0.0;
//
// n p d t He-3 He-4
const G4int gam[6] = { 1, 1, 3, 3, 3, 2 }; // gammab/2
G4double result;
G4double q1, q2;
if( j == 6 )
{
q1 = 2*sqrt( AMF * afj[6] * rj[6] );
q2 = 33.4/(pi*radncl*radncl);
result = q2/(AMF*afj[6])*((q1-1)*exp(q1-per)+exp(-per));
}
else
{
G4double alfa, beta, q3, q4, q5;
if( j == 0 )
{
alfa = 0.76 + 2.2/pow(afj[0],oneThird);
beta = (2.12/pow(afj[0],twoThirds)-0.05)/alfa;
}
else
{
alfa = 1 + cc;
beta = 0.0;
}
q1 = am*afj[j];
q2 = q1*rj[j];
q3 = gam[j]*pow(afj[j],twoThirds)*alfa/(q1*q1);
q4 = (2*beta*q1-3)/2 + q2;
q5 = (2*beta*q1-3)*(sqrt(q2)-0.5) + 2*q2;
result = q3*(q4*exp(-per)+q5*exp(2*sqrt(q2)-per));
}
return result;
}
G4double G4PreEquilibrium::poten( const G4int i, const IPSTRUCT ip ) const
{
// calculation of particle potential in i'th nuclear zone
//
const G4int N = 7;
const G4double EPS = 0.007;
const G4double VPI = 0.025;
G4double result;
if( i >= N+2 )
result = 0.0;
else if( i == N+1 )
result = ip.proton*af[N];
else if( ip.ip2 != 0 )
result = ip.proton*af[i-1];
else if( ip.nuclearZone != 0 )
result = ip.proton*af[i-1] + tfp[i-1]*ip.proton + (1-ip.proton)*tfn[i-1] + ip.nuclearZone*EPS;
else
result = ip.proton*af[i-1] - VPI*(ip.ip1-1);
return result;
}
G4double G4PreEquilibrium::cinema( const G4ThreeVector pstar, const G4ThreeVector v,
G4ThreeVector &p, G4double &ct, G4double &st,
G4double &cfi, G4double &sfi, const G4double cm ) const
{
// kinematic code
//
G4double spv = pstar.dot(v);
G4double v2 = v.mag2();
G4double temp1 = sqrt(1-v2);
G4double temp2 = spv*(1/temp1-1)/v2;
G4double pmstar = pstar.mag();
G4double tstar = sqrt(pmstar*pmstar+cm*cm) - cm;
p.setX( pstar.x() + v.x()*(temp2 + (tstar+cm)/temp1) );
p.setY( pstar.y() + v.y()*(temp2 + (tstar+cm)/temp1) );
p.setZ( pstar.z() + v.z()*(temp2 + (tstar+cm)/temp1) );
G4double pm = p.mag();
if( p.z()*p.z() >= pm*pm )
{
ct = 1;
st = 0;
cfi = 1;
sfi = 0;
}
else
{
ct = p.z()/pm;
st = sqrt(1-ct*ct);
cfi = p.x()/(pm*st);
sfi = p.y()/(pm*st);
}
return sqrt(pm*pm+cm*cm)-cm;
}
G4ThreeVector G4PreEquilibrium::rotation( const G4ThreeVector a,
const G4ThreeVector b,
const G4ThreeVector ps ) const
{
G4double sp = a.dot(b);
//
G4double amod = a.mag();
G4double alpha1 = sp/amod;
G4double alpha2 = sqrt( b.mag2() - alpha1*alpha1 );
//
G4ThreeVector an = a.cross(b);
//
G4ThreeVector pr;
pr.setX( ps.x()*b.x()/alpha2 + (ps.z()-alpha1*ps.x()/alpha2)*a.x()/amod + (ps.y()*an.x()) );
pr.setY( ps.x()*b.y()/alpha2 + (ps.z()-alpha1*ps.x()/alpha2)*a.y()/amod + (ps.y()*an.y()) );
pr.setZ( ps.x()*b.z()/alpha2 + (ps.z()-alpha1*ps.x()/alpha2)*a.z()/amod + (ps.y()*an.z()) );
pr *= 1/alpha2/amod;
return pr;
}
G4double G4PreEquilibrium::crossSectionInterp( const G4double x, const G4int lq )
{
const G4double sigma[28][30] =
{ {17613.000, 330.000, 154.000, 96.000, 70.000,
51.000, 38.400, 30.000, 23.600, 22.400,
22.200, 22.600, 23.400, 24.700, 29.500,
40.500, 48.500, 47.400, 47.000, 46.700,
46.000, 45.000, 43.000, 41.200, 40.800,
40.500, 39.800, 39.000, 39.000, 38.500},
{17613.000, 330.000, 154.000, 96.000, 70.000,
51.000, 38.400, 30.000, 23.600, 22.400,
22.200, 22.600, 22.800, 22.900, 25.000,
25.000, 25.000, 25.000, 25.000, 22.000,
19.500, 17.500, 15.000, 13.700, 12.000,
11.000, 9.800, 8.800, 8.500, 6.500},
{20357.000, 950.000, 480.000, 300.000, 200.000,
160.000, 108.000, 74.000, 50.000, 41.000,
36.500, 34.000, 32.500, 32.000, 34.200,
36.100, 37.800, 38.400, 39.000, 40.000,
40.000, 40.500, 41.100, 42.300, 42.300,
42.000, 41.100, 39.600, 39.500, 39.200},
{20357.000, 950.000, 480.000, 300.000, 200.000,
160.000, 108.000, 74.000, 50.000, 41.000,
36.500, 34.000, 32.500, 31.200, 30.800,
25.100, 19.100, 18.000, 16.500, 16.000,
15.200, 14.000, 12.100, 10.800, 10.700,
11.000, 9.600, 8.000, 7.000, 6.200},
{ 6.000, 6.000, 6.500, 7.000, 8.500,
10.500, 16.500, 25.300, 57.500, 68.500,
64.500, 52.000, 40.500, 25.700, 29.000,
32.100, 45.600, 38.000, 44.300, 54.000,
58.000, 45.700, 35.300, 35.000, 34.300,
34.000, 32.400, 30.400, 26.500, 25.000},
{ 2.000, 2.000, 2.250, 2.300, 2.400,
2.500, 4.500, 7.700, 20.000, 25.200,
23.900, 21.000, 16.000, 10.500, 11.200,
14.000, 20.300, 16.000, 19.300, 26.000,
26.500, 18.700, 11.700, 10.500, 9.600,
9.500, 6.800, 5.900, 5.000, 4.000},
{ 4.000, 4.000, 4.250, 4.700, 6.100,
8.000, 12.000, 17.600, 37.500, 43.300,
40.600, 31.000, 24.500, 13.100, 11.000,
9.500, 8.300, 5.000, 5.500, 6.800,
7.000, 3.500, 2.000, 1.900, 1.800,
1.600, 0.220, 0.150, 0.048, 0.009},
{ 1.900, 2.300, 3.500, 5.500, 9.000,
14.000, 28.000, 60.000, 163.000, 195.000,
185.000, 145.000, 113.000, 45.000, 25.200,
21.600, 15.600, 15.200, 19.500, 22.800,
24.500, 27.600, 36.700, 41.000, 39.000,
32.300, 28.900, 27.700, 24.900, 23.500},
{ 1.900, 2.300, 3.500, 5.500, 9.000,
14.000, 28.000, 60.000, 163.000, 195.000,
184.900, 144.800, 112.800, 44.400, 23.200,
18.600, 10.800, 7.700, 9.000, 10.200,
11.300, 13.500, 16.900, 19.000, 16.900,
12.600, 5.700, 5.600, 4.900, 4.000},
{ 10.000, 12.000, 14.000, 16.000, 17.000,
20.000, 30.000, 38.000, 42.000, 38.500,
32.000, 24.000, 18.000, 1.000, 0.000,
0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 0.000, 0.000, 0.000, 0.000},
{ 0.000, 0.000, 0.000, 0.050, 0.100,
0.300, 0.600, 1.200, 2.200, 3.200,
3.400, 3.600, 3.700, 3.800, 3.900,
3.900, 3.900, 4.000, 4.000, 4.000,
4.000, 4.000, 3.900, 3.800, 3.500,
3.200, 2.900, 2.700, 2.400, 2.100},
{ 0.000, 0.000, 0.000, 0.400, 0.800,
1.400, 2.300, 4.400, 8.000, 10.800,
15.000, 16.100, 16.600, 17.000, 17.300,
17.500, 17.600, 17.700, 17.800, 17.900,
17.800, 17.500, 16.800, 16.000, 13.700,
12.600, 11.600, 10.800, 9.200, 8.000},
{ 0.000, 0.000, 0.000, 0.300, 0.600,
1.000, 1.600, 2.700, 5.000, 6.300,
6.900, 7.200, 7.400, 7.500, 7.600,
7.700, 7.800, 8.000, 8.100, 8.100,
8.100, 8.000, 7.900, 7.700, 7.000,
6.600, 6.200, 5.800, 4.600, 3.500},
{ 0.000, 0.000, 0.000, 0.000, 0.100,
0.200, 0.500, 0.800, 1.800, 2.100,
2.400, 2.600, 2.800, 3.000, 3.200,
3.300, 3.400, 3.650, 3.800, 3.850,
3.900, 4.000, 4.000, 3.900, 3.500,
3.200, 2.700, 2.100, 1.000, 0.600},
{ 0.000, 0.000, 0.000, 0.000, 0.200,
0.600, 1.200, 2.100, 3.300, 4.100,
4.900, 6.000, 7.700, 8.900, 10.000,
10.200, 10.100, 9.800, 9.800, 10.100,
10.300, 8.500, 6.900, 5.500, 3.900,
3.300, 3.000, 2.700, 2.200, 2.100},
{ 0.000, 0.000, 0.000, 0.000, 0.100,
0.300, 0.500, 0.700, 0.900, 1.100,
1.300, 1.500, 1.800, 2.200, 2.600,
2.600, 2.400, 2.400, 3.000, 3.400,
3.500, 3.100, 2.700, 2.400, 2.100,
1.900, 1.800, 1.700, 1.500, 1.300},
{ 0.000, 0.000, 0.000, 0.200, 0.700,
1.400, 3.000, 4.900, 5.100, 5.000,
4.400, 3.600, 3.100, 2.700, 2.400,
2.200, 2.000, 1.700, 1.600, 1.400,
1.300, 1.100, 0.900, 0.800, 0.600,
0.600, 0.600, 0.500, 0.400, 0.300},
{ 0.000, 0.000, 0.000, 0.200, 0.500,
0.800, 1.300, 2.100, 3.400, 4.700,
6.400, 6.800, 6.200, 5.200, 7.100,
8.000, 8.900, 9.800, 10.300, 10.500,
10.600, 9.900, 7.700, 5.500, 3.400,
2.800, 2.700, 2.600, 2.200, 1.900},
{ 0.000, 0.100, 0.500, 1.200, 3.000,
4.200, 5.000, 5.200, 6.100, 7.900,
8.900, 9.600, 9.900, 10.000, 10.100,
10.100, 10.000, 9.800, 9.500, 8.800,
8.000, 6.800, 5.900, 5.200, 4.200,
2.900, 2.700, 2.500, 2.200, 1.900},
{ 0.000, 0.011, 0.026, 0.060, 0.120,
0.180, 0.270, 0.280, 0.240, 0.186,
0.130, 0.076, 0.052, 0.040, 0.031,
0.032, 0.040, 0.044, 0.039, 0.033,
0.028, 0.025, 0.022, 0.019, 0.015,
0.012, 0.007, 0.004, 0.002, 0.000},
{ 0.000, 0.070, 0.107, 0.143, 0.183,
0.217, 0.240, 0.214, 0.177, 0.141,
0.117, 0.089, 0.081, 0.080, 0.083,
0.093, 0.103, 0.095, 0.061, 0.050,
0.049, 0.052, 0.056, 0.050, 0.025,
0.017, 0.008, 0.004, 0.002, 0.000},
{ 5.000, 3.000, 1.125, 0.900, 0.700,
0.580, 0.500, 0.415, 0.350, 0.280,
0.260, 0.270, 0.280, 0.295, 0.310,
0.298, 0.250, 0.200, 0.145, 0.115,
0.095, 0.065, 0.040, 0.028, 0.018,
0.010, 0.005, 0.000, 0.000, 0.000},
{ 0.000, 0.011, 0.031, 0.055, 0.067,
0.078, 0.084, 0.088, 0.088, 0.087,
0.085, 0.083, 0.080, 0.078, 0.074,
0.072, 0.069, 0.066, 0.064, 0.060,
0.057, 0.052, 0.038, 0.031, 0.027,
0.024, 0.023, 0.022, 0.022, 0.022},
{ 0.000, 0.002, 0.004, 0.006, 0.007,
0.008, 0.008, 0.009, 0.011, 0.012,
0.013, 0.014, 0.015, 0.016, 0.016,
0.016, 0.016, 0.016, 0.016, 0.015,
0.015, 0.014, 0.012, 0.010, 0.009,
0.008, 0.007, 0.007, 0.007, 0.007},
{ 0.000, 0.002, 0.004, 0.006, 0.007,
0.008, 0.008, 0.009, 0.011, 0.012,
0.013, 0.014, 0.015, 0.016, 0.016,
0.016, 0.016, 0.016, 0.016, 0.015,
0.015, 0.014, 0.012, 0.010, 0.009,
0.008, 0.007, 0.007, 0.007, 0.007},
{ 0.000, 0.010, 0.027, 0.047, 0.057,
0.065, 0.069, 0.071, 0.069, 0.065,
0.061, 0.057, 0.052, 0.048, 0.043,
0.039, 0.035, 0.030, 0.027, 0.023,
0.019, 0.015, 0.008, 0.004, 0.003,
0.002, 0.002, 0.001, 0.001, 0.001},
{ 0.000, 0.015, 0.039, 0.067, 0.081,
0.094, 0.100, 0.106, 0.110, 0.111,
0.111, 0.111, 0.110, 0.110, 0.106,
0.104, 0.101, 0.098, 0.096, 0.090,
0.087, 0.080, 0.062, 0.051, 0.045,
0.040, 0.037, 0.036, 0.036, 0.036},
{ 0.000, 0.081, 0.133, 0.203, 0.303,
0.397, 0.510, 0.494, 0.417, 0.327,
0.247, 0.180, 0.172, 0.187, 0.208,
0.231, 0.253, 0.250, 0.211, 0.194,
0.187, 0.187, 0.184, 0.173, 0.138,
0.119, 0.095, 0.070, 0.049, 0.037} };
const G4double argus[6][30] =
{ {0.000, 0.010, 0.020, 0.030, 0.040,
0.050, 0.070, 0.100, 0.150, 0.200,
0.250, 0.300, 0.350, 0.400, 0.500,
0.650, 0.850, 0.950, 1.100, 1.300,
1.500, 2.000, 3.000, 4.000, 5.000,
7.000, 10.000, 16.000, 22.000, 30.000},
{0.000, 0.010, 0.020, 0.030, 0.040,
0.050, 0.075, 0.100, 0.150, 0.175,
0.200, 0.225, 0.250, 0.350, 0.450,
0.500, 0.600, 0.700, 0.800, 0.850,
0.900, 1.000, 1.200, 1.300, 1.400,
1.600, 3.000, 4.000, 10.000, 20.000},
{0.200, 0.250, 0.300, 0.350, 0.400,
0.450, 0.500, 0.550, 0.600, 0.650,
0.700, 0.750, 0.800, 0.850, 0.900,
0.950, 1.000, 1.100, 1.200, 1.300,
1.400, 1.600, 1.800, 2.000, 2.400,
2.600, 2.800, 3.000, 3.500, 4.000},
{0.150, 0.175, 0.200, 0.225, 0.250,
0.275, 0.300, 0.325, 0.350, 0.375,
0.400, 0.450, 0.500, 0.550, 0.600,
0.650, 0.700, 0.750, 0.800, 0.850,
0.900, 0.950, 1.000, 1.050, 1.150,
1.250, 1.500, 2.000, 3.000, 4.000},
{0.010, 0.020, 0.040, 0.050, 0.060,
0.070, 0.080, 0.090, 0.100, 0.120,
0.150, 0.170, 0.190, 0.210, 0.250,
0.270, 0.300, 0.320, 0.350, 0.370,
0.400, 0.440, 0.500, 0.600, 0.700,
0.800, 0.900, 1.000, 2.000, 3.000},
{0.400, 0.450, 0.500, 0.550, 0.575,
0.600, 0.625, 0.650, 0.700, 0.750,
0.800, 0.850, 0.900, 0.950, 1.000,
1.050, 1.100, 1.150, 1.200, 1.250,
1.300, 1.500, 2.000, 2.500, 3.000,
3.500, 4.000, 4.500, 5.000, 5.500} };
G4int i;
if( lq < 3 )i = 0;
else
{
if( lq < 9 )i = 1;
else
{
if( lq < 18 )i = 2;
else
{
if( lq < 20 )i = 3;
else
{
if( lq < 21 )i = 4;
else
{
if( lq < 26 )i = 5;
else i = 3;
}
}
}
}
}
G4int lpha = 0;
while( x > argus[lpha][i] )++lpha;
if( x == argus[lpha][i] )return sigma[lpha][lq];
if( lpha == 0 )return 0.0;
G4double phi1, psi1, phi2, psi2, phi3, psi3;
if( lpha >= 28 )
{
phi1 = sigma[27][lq];
psi1 = argus[27][i];
phi2 = sigma[28][lq];
psi2 = argus[28][i];
phi3 = sigma[29][lq];
psi3 = argus[29][i];
}
else
{
phi1 = sigma[lpha-1][lq];
psi1 = argus[lpha-1][i];
phi2 = sigma[lpha][lq];
psi2 = argus[lpha][i];
phi3 = sigma[lpha+1][lq];
psi3 = argus[lpha+1][i];
}
G4double delta = (psi2-psi3)*psi1*psi1 + (psi3-psi1)*psi2*psi2 + (psi1-psi2)*psi3*psi3;
G4double deltaa = phi1*(psi2-psi3) + phi2*(psi3-psi1) + phi3*(psi1-psi2);
G4double deltab = (phi2-phi3)*psi1*psi1 + (phi3-phi1)*psi2*psi2 + (phi1-phi2)*psi3*psi3;
G4double deltac = (psi2*phi3-psi3*phi2)*psi1*psi1 + (psi3*phi1-psi1*phi3)*psi2*psi2 +
(psi1*phi2-psi2*phi1)*psi3*psi3;
G4double a = deltaa/delta;
G4double b = deltab/delta;
G4double c = deltac/delta;
return a*x*x + b*x + c;
}
G4double G4PreEquilibrium::sigmat( const G4int i, const G4int ms, const G4int mq,
const G4int ksi, const G4int iks, const G4double t )
{
// choose cross section type and calculate cross section value for given energy
//
const G4int icst[28] = {
210, 211, 220, 221, 120, 121, 122,
110, 111, 123, 214, 215, 224, 225,
114, 115, 126, 124, 125, 10111, 10112,
10113, 10115, 10114, 10116, 10118, 10117, 10110 };
const G4int nsicst[21] = {
112, 113, 116, 117, 127, 130, 131,
132, 133, 134, 135, 136, 137, 212,
213, 216, 217, 222, 223, 226, 227 };
G4int ics = 10000*i + 1000*ms + 100*mq + 10*ksi + iks;
G4int js = 0;
do
{
if( ics == icst[js] )return crossSectionInterp( t, js );
} while ( js++ < 28 );
G4int nsjs = 0;
while( ics != nsicst[nsjs] )++nsjs;
G4double result = 0;
switch (nsjs)
{
case 2:
result = crossSectionInterp( t, 9 );
break;
case 4:
result = crossSectionInterp( t, 14 ) + crossSectionInterp( t, 15 );
break;
case 5:
result = crossSectionInterp( t, 17 ) + crossSectionInterp( t, 18 ) +
crossSectionInterp( t, 16 );
break;
case 6:
result = (crossSectionInterp( t, 7 ) + crossSectionInterp( t, 4 ))/2;
break;
case 7:
result = (crossSectionInterp( t, 8 ) + crossSectionInterp( t, 5 ) -
crossSectionInterp( t, 6 ))/2;
break;
case 8:
result = crossSectionInterp( t, 6 );
break;
case 9:
result = crossSectionInterp( t, 9 );
break;
case 10:
result = (crossSectionInterp( t, 14 ) + crossSectionInterp( t, 17 ))/2;
break;
case 11:
result = (crossSectionInterp( t, 15 ) + crossSectionInterp( t, 18 ))/2;
break;
case 12:
result = crossSectionInterp( t, 16 )/2;
break;
case 13:
result = (crossSectionInterp( t, 14 ) + crossSectionInterp( t, 15 ) +
crossSectionInterp( t, 17 ) + crossSectionInterp( t, 18 ) +
crossSectionInterp( t, 16 ))/2;
break;
case 17:
result = crossSectionInterp( t, 10 ) + crossSectionInterp( t, 11 );
break;
case 20:
result = crossSectionInterp( t, 12 );
break;
case 21:
result = 2*crossSectionInterp( t, 12 ) + crossSectionInterp( t, 13 );
break;
}
return result;
}
void G4PreEquilibrium::slqek( G4int &i, G4int &ms, G4int &mq, G4int &ksi, G4int &me,
const G4int lin, const G4int msin, const G4int mqin,
const G4int mein, const G4int ln, const G4int msn,
const G4int mqn, const G4int men )
{
// form cross-sections type
//
ms = msin + msn;
i = lin + ln; // this ln is not the data member ln
mq = mqin + mqn;
me = mein + men;
if( ms != 0 )return; // strange particle
if( i > 0 )
ksi = 1;
else if( mq > 1 )
me != 1 ? ksi = 1 : ksi = 2;
else if( me == 2 )
ksi = 1;
else if( me == -1 )
ksi = 1;
else if( me != 0 )
mein == 1 ? ksi = 2 : ksi = 3;
else
mein == -1 ? ksi = 2 : ksi = 3;
return;
}
G4double G4PreEquilibrium::costa( const G4int j, const G4double t )
{
// cosinus calculation
//
const G4double ankj[29][4][4] =
{ { { 2.7404, -9.6998, 10.400, 2.3882 },
{-7.5137, 44.096, -74.379, 46.038 },
{ 7.5479, -39.274, 64.835, -41.609 },
{-1.8369, 8.6911, -13.060, 7.1880 } },
{ {-30.853, 106.24, -129.39, 54.339 },
{ 19.465, -68.102, 96.358, -56.827 },
{-3.4831, 12.341, -18.592, 12.024 },
{ 0.18941, -0.67880, 1.0665, -0.72910 } },
{ { 0.10258, -1.0542, 11.389, -16.638 },
{-0.49607, 11.800, -90.857, 164.76 },
{ 1.5437, -33.769, 251.92, -450.71 },
{-1.2021, 25.336, -186.58, 332.54 } },
{ { 0.15789, 2.9671, -5.5251, 6.8925 },
{-7.0218, -205.34, 569.51, -898.58 },
{ 134.96, 4872.2, -14674., 23924. },
{-821.16, -32586., 100980., -165530. } },
{ { 0.31531, -7.4981, 43.295, -76.360 },
{-6.5373, 193.07, -1018.1, 1742.6 },
{ 46.864, -1303.0, 6729.1, -11075. },
{-95.192, 2637.3, -12857., 20294. } },
{ {-17.953, 109.72, -239.54, 228.26 },
{ 91.968, -519.63, 1126.6, -1074.0 },
{-132.70, 741.12, -1600.0, 1524.9 },
{ 58.598, -318.74, 677.51, -640.11 } },
{ { 0.42169, 147.05, -653.35, 915.07 },
{-3.5198, -260.19, 1225.0, -1748.1 },
{ 3.6373, 155.92, -752.01, 1079.6 },
{-0.78041, -30.563, 147.95, -212.50 } },
{ {-0.38288, 3.7587, -6.5144, 6.7740 },
{ 103.81, -272.82, 477.59, -512.22 },
{-1788.2, 4305.2, -7931.4, 9347.1 },
{ 7147.5, -3339.5, -4139.2, -4436.4 } },
{ { 0.24991, 32.028, -118.82, 150.99 },
{-2.6994, -460.45, 1895.9, -2519.0 },
{ 16.268, 2138.4, -9126.2, 12431. },
{-29.654, -3182.3, 13944., -19342. } },
{ { 3.9025, -91.126, 323.73, -400.48 },
{-20.619, 491.70, -1715.5, 2114.3 },
{ 33.004, -766.84, 2700.3, -3352.5 },
{-16.367, 373.94, -1320.2, 1642.3 } },
{ { 19.402, -224.46, 747.33, -935.70 },
{-44.180, 471.94, -1485.6, 1805.5 },
{ 31.567, -301.76, 907.63, -1077.3 },
{-6.8648, 60.476, -175.20, 203.81 } },
{ { 0.40693, -4.1404, 14.044, -17.265 },
{-3.6799, 59.610, -162.69, 188.73 },
{ 14.556, -175.50, 458.39, -533.90 },
{-12.621, 149.64, -381.18, 451.41 } },
{ {-0.47554, 2.2641, -12.528, 24.647 },
{ 5.1620, -9.9236, 55.623, -104.62 },
{-8.1117, 19.315, -84.255, 139.08 },
{ 3.5187, -9.1783, 34.950, -51.243 } },
{ { 0.48173, 5.7726, -13.745, 27.125 },
{-4.4804, -38.582, 111.59, -243.05 },
{ 16.306, 110.46, -330.45, 722.70 },
{-15.968, -80.140, 246.16, -607.53 } },
{ {-5.1646, -6.0776, 78.989, -107.05 },
{ 21.871, 56.915, -401.59, 512.15 },
{-27.993, -94.670, 569.28, -696.21 },
{ 11.587, 45.998, -245.6, 284.52 } },
{ {-53.067, 576.12, -1543.8, 164550. },
{ 147.50, -1638., 4592.3, -4994.9 },
{-134.36, 1578.0, -4446.3, 4902.2 },
{ 40.253, -488.60, 1400.1, -1560.6 } },
{ { 0.14988, 2.8753, -5.3078, 6.2233 },
{-5.9558, -162.03, 430.79, -625.48 },
{ 128.75, 3140.2, -7918.9, 10983. },
{-851.61, -18780., 44607., -58790. } },
{ { 0.53689, -13.216, 81.011, -142.85 },
{-10.550, 296.29, -1695.7, 2893.5 },
{ 69.621, -1924.5, 10620., -17468. },
{-138.65, 3928.1, -20293., 32058. } },
{ { 0.65288, 0.38977, 0.84078, 0.18893 },
{-4.3964, 34.309, -73.692, 84.308 },
{ 14.889, -143.80, 312.27, -350.14 },
{-15.658, 171.60, -372.12, 412.99 } },
{ { 0.085591, 5.0390, -13.782, 14.661 },
{ 0.054284, -9.2324, 36.397, -42.962 },
{-0.051111, 4.6003, -20.534, 27.731 },
{ 0.0074514, -0.62529, 2.9159, -4.1101 } },
{ { 0.071622, 3.0960, -11.125, 18.130 },
{ 0.092581, -3.2186, 20.273, -33.245 },
{-0.051531, 0.89886, -7.5084, 13.188 },
{ 0.0058258, -0.0017288, 0.70224, -1.4854 } },
{ { 0.082300, 0.15854, 3.7716, -4.0562 },
{ 0.010802, -0.33688, 1.1727, -0.67476 },
{-0.0021798, 0.052166, -0.25816, 0.32048 },
{ 0.000065764, -0.0014711, 0.0078209, -0.010580 } },
{ { 0.11138, 0.60396, 3.0174, -4.4190 },
{-0.017709, 0.23015, -1.8187, 3.4518 },
{ 0.0020977, -0.025458, 0.21626, -0.40692 },
{-0.000054799, 0.00059111, -0.0055552, 0.010647 } },
{ { 0.17288, 7.1080, -17.961, 16.403 },
{-0.14504, -13.032, 41.781, -40.799 },
{ 0.045390, 8.3515, -30.260, 32.882 },
{-0.0047961, -1.4095, 5.3505, -6.0946 } },
{ { 0.037596, 1.4331, -3.1350, 6.4864 },
{ 0.23827, 1.8253, 1.7648, -16.735 },
{-0.15410, -1.5201, -1.5692, 17.185 },
{ 0.025037, 0.30588, 0.32520, -3.5277 } },
{ { 0.12489, 1.3573, 0.82338, -1.4595 },
{-0.051577, -0.35778, -1.1690, 1.8078 },
{ 0.0074864, 0.032888, 0.23744, -0.39802 },
{-0.00029880, -0.00075117, -0.011402, 0.019505 } },
{ { 0.18470, 1.9269, -3.2979, 3.6843 },
{-0.073932, 0.27213, 1.0600, -2.3354 },
{ 0.018907, -0.056473, -0.16487, 0.38426 },
{-0.00092984, 0.0025506, 0.0073052, -0.017220 } },
{ {-1.0306, 32.849, -75.052, 60.255 },
{ 7.9586, -125.72, 256.04, -165.47 },
{-14.797, 165.90, -279.91, 113.33 },
{ 8.2309, -67.871, 85.762, 5.9727 } },
{ {-237.22, 968.90, -1621.9, 1363.7 },
{ 658.00, -2694.1, 4548.0, -3846.0 },
{-606.53, 2498.3, -4249.8, 3613.6 },
{ 186.04, -769.33, 1316.6, -1124.2 } } };
G4double ank[4][4];
for( G4int k = 0; k < 4; ++k )
{
for( G4int n = 0; n < 4; ++n )ank[n][k] = ankj[n][k][j];
}
G4double s1 = 0;
G4double s2 = 0;
G4double r = G4UniformRand();
for( G4int n = 0; n < 4; ++n )
{
for( G4int k = 0; k < 4; ++k )
{
s1 += ank[n][k]*pow(t,k)*pow(r,n);
s2 += ank[n][k]*pow(t,k);
}
}
G4double cta = 2*sqrt(r)*(s1+(1-s2)*pow(r,4)) - 1;
if( cta < -1 )cta = -1;
else if( cta > 1 )cta = 1;
return cta;
}
G4double G4PreEquilibrium::cosex( const G4int i, const G4double t, const G4double cm )
{
// cosinus calculation for charge exchange scattering
//
G4double result;
if( i != 0 )
{
if( t <= 0.51 )result = costa( 13, t );
else
t <= 1 ? result = costa( 14, t ) : result = costa( 15, t );
}
else
{
if( t <= 0.08 )result = costa( 16, t );
else
{
if( t <= 0.3 )result = costa( 17, t );
else
{
if( t <= 1 )result = costa( 9, t );
else
{
if( t <= 2.4 )result = costa( 10, t );
else
{
G4double tmax = sqrt(t*(t+2*cm));
result = 1+(2*log(1+G4UniformRand()*(exp(-7.5*tmax)-1)))/(7.5*tmax);
}
}
}
}
}
return result;
}
G4double G4PreEquilibrium::cosel( const G4int i, const G4int mq, const G4int ksi,
const G4double t, const G4double cm )
{
// cosinus calculation for elastic scattering
//
G4double result;
if( i != 0 )
{
t <= 0.45 ? result = costa( 11, t ) : result = costa( 12, t );
}
else
{
if( mq >= 2 )
{
if( ksi >= 2 )
{
if( t <= 0.97 )result = costa( 2, t );
else
{
if( t <= 2.8 )result = (1+costa( 0, t ))/2;
else
{
if( t <= 10 )result = (3-costa( 1, t ))/4;
else
{
G4double tmax = sqrt(t*(t+2*cm));
result = 1+(2*log(1+G4UniformRand()*(exp(-8.7*tmax)-1)))/(8.7*tmax);
}
}
}
}
else
{
if( t <= 0.46 )result = 1-2*G4UniformRand();
else
{
if( t <= 2.8 )result = (1+costa( 0, t ))/2;
else
{
if( t <= 10 )result = (3-costa( 1, t ))/4;
else
{
G4double tmax = sqrt(t*(t+2*cm));
result = 1+(2*log(1+G4UniformRand()*(exp(-8.7*tmax)-1)))/(8.7*tmax);
}
}
}
}
}
else
{
if( ksi < 2 )
{
if( t <= 0.08 )result = costa( 3, t );
else
{
if( t <= 0.3 )result = costa( 4, t );
else
{
if( t <= 1 )result = costa( 5, t );
else
{
if( t <= 2.4 )result = costa( 6, t );
else
{
G4double tmax = sqrt(t*(t+2*cm));
result = 1+(2*log(1+G4UniformRand()*(exp(-7.5*tmax)-1)))/(7.5*tmax);
}
}
}
}
}
else if( ksi == 2 )
{
if( t <= 0.08 )result = costa( 7, t );
else
{
if( t <= 0.3 )result = costa( 8, t );
else
{
if( t <= 1 )result = costa( 9, t );
else
{
if( t <= 2.4 )result = costa( 10, t );
else
{
G4double tmax = sqrt(t*(t+2*cm));
result = 1+(2*log(1+G4UniformRand()*(exp(-7.5*tmax)-1)))/(7.5*tmax);
}
}
}
}
}
else if( ksi > 2 )
{
if( G4UniformRand() <= 0.5 )
{
if( t <= 0.08 )result = costa( 3, t );
else
{
if( t <= 0.3 )result = costa( 4, t );
else
{
if( t <= 1 )result = costa( 5, t );
else
{
if( t <= 2.4 )result = costa( 6, t );
else
{
G4double tmax = sqrt(t*(t+2*cm));
result = 1+(2*log(1+G4UniformRand()*(exp(-7.5*tmax)-1)))/(7.5*tmax);
}
}
}
}
}
else
{
if( t <= 0.08 )result = costa( 7, t );
else
{
if( t <= 0.3 )result = costa( 8, t );
else
{
if( t <= 1 )result = costa( 9, t );
else
{
if( t <= 2.4 )result = costa( 10, t );
else
{
G4double tmax = sqrt(t*(t+2*cm));
result = 1+(2*log(1+G4UniformRand()*(exp(-7.5*tmax)-1)))/(7.5*tmax);
}
}
}
}
}
}
}
}
return result;
}
// end of file
Reference in New Issue View Git Blame Copy Permalink