00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include "G4BEChargedChannel.hh"
00033 #include "G4SystemOfUnits.hh"
00034
00035 G4BEChargedChannel::G4BEChargedChannel()
00036 {
00037 verboseLevel = 0;
00038 }
00039
00040
00041 G4BEChargedChannel::~G4BEChargedChannel()
00042 {
00043 }
00044
00045
00046 void G4BEChargedChannel::calculateProbability()
00047 {
00048 G4int residualZ = nucleusZ - particleZ;
00049 G4int residualA = nucleusA - particleA;
00050
00051
00052
00053
00054 if ( nucleusA < 2.0 * particleA ||
00055 nucleusZ < 2.0 * particleZ ||
00056 residualA <= residualZ ||
00057 excitationEnergy - getThresh() - correction < 0 )
00058 {
00059 if ( verboseLevel >= 6 )
00060 G4cout << "G4BEChargedChannel : calculateProbability for " << getName() << " = 0 " << G4endl;
00061 emissionProbability = 0;
00062 return;
00063 }
00064
00065
00066
00067
00068 G4double levelParam = getLevelDensityParameter();
00069 G4double slevel = 2 * std::sqrt( levelParam * ( excitationEnergy - getThresh() - correction ) );
00070 G4double constant = A / 2 * ( 2 * spin + 1 ) * ( 1 + coulombFactor() );
00071 G4double eye1 = ( slevel*slevel - 3 * slevel + 3 ) / ( 4 * levelParam*levelParam ) * std::exp( slevel );
00072
00073 emissionProbability = constant * std::pow( G4double(residualA), 0.6666666 ) * eye1;
00074
00075 if ( verboseLevel >= 6 )
00076 G4cout << "G4BEChargedChannel : calculateProbability for " << getName() << G4endl
00077 << " res A = " << residualA << G4endl
00078 << " res Z = " << residualZ << G4endl
00079 << " c*(c_i+1) = "<< constant << G4endl
00080 << " qmfactor = "<< qmFactor() << G4endl
00081 << " coulombfactor = "<< coulombFactor() << G4endl
00082 << " E = " << excitationEnergy << G4endl
00083 << " correction = " << correction << G4endl
00084 << " eye1 = " << eye1 << G4endl
00085 << " levelParam = " << levelParam << G4endl
00086 << " thresh = " << getThresh() << G4endl
00087 << " s = " << s << G4endl
00088 << " probability = " << emissionProbability << G4endl;
00089
00090 return;
00091 }
00092
00093
00094 G4double G4BEChargedChannel::sampleKineticEnergy()
00095 {
00096 G4double levelParam;
00097 levelParam = getLevelDensityParameter();
00098
00099 const G4double xMax = excitationEnergy - getThresh() - correction;
00100 const G4double xProb = ( - 1 + std::sqrt ( 1 + 4 * levelParam * xMax ) ) / ( 2 * levelParam );
00101 const G4double maxProb = xProb * std::exp ( 2 * std::sqrt ( levelParam * ( xMax - xProb ) ) );
00102
00103
00104 G4double r1;
00105 G4double r2;
00106 G4int koe=0;
00107 do
00108 {
00109 r1 = G4UniformRand() * xMax;
00110 r2 = G4UniformRand() * maxProb;
00111 koe++;
00112 }
00113 while ( r1 * std::exp ( 2 * std::sqrt ( levelParam * ( xMax - r1 ) ) ) < r2 );
00114
00115
00116 G4double kineticEnergy = r1 + getCoulomb();
00117
00118 if ( verboseLevel >= 10 )
00119 G4cout << " G4BENeutronChannel : sampleKineticEnergy() " << G4endl
00120 << " kinetic n e = " << kineticEnergy << G4endl
00121 << " levelParam = " << levelParam << G4endl
00122 << " thresh= " << getThresh() << G4endl;
00123
00124 return kineticEnergy;
00125 }
00126
00127
00128 G4double G4BEChargedChannel::coulombFactorForProton()
00129 {
00130
00131
00132
00133 G4double t[7] = { 0.08 , 0 , -0.06 , -0.1 , -0.1 , -0.1 , -0.1 };
00134 G4int Z = nucleusZ - particleZ;
00135
00136 if ( Z >= 70.0 ) return t[6];
00137 if ( Z <= 10.0 ) return t[0];
00138
00139
00140 G4int n = G4int( 0.1 * Z + 1.0 );
00141 G4float x = ( 10 * n - Z ) * 0.1;
00142 G4double ret_val = x * t[n - 2] + ( 1.0 - x ) * t[n-1];
00143
00144 return ret_val;
00145 }
00146
00147
00148 G4double G4BEChargedChannel::qmFactorForProton()
00149 {
00150
00151
00152
00153 G4double t[7] = { 0.36, 0.51, 0.60, 0.66, 0.68, 0.69, 0.69 };
00154 G4int Z = nucleusZ - particleZ;
00155
00156 if ( Z >= 70.0 ) return t[6];
00157 if ( Z <= 10.0 ) return t[0];
00158
00159
00160 G4int n = G4int( 0.1 * Z + 1.0 );
00161 G4float x = ( 10 * n - Z ) * 0.1;
00162 return x * t[n - 2] + ( 1.0 - x ) * t[n-1];
00163 }
00164
00165
00166 G4double G4BEChargedChannel::qmFactorForAlpha()
00167 {
00168
00169
00170
00171 G4double t[7] = { 0.77, 0.81, 0.85, 0.89, 0.93, 0.97, 1.00 };
00172 G4int Z = nucleusZ - particleZ;
00173
00174 if ( Z >= 70.0 ) return t[6];
00175 if ( Z <= 10.0 ) return t[0];
00176
00177
00178 G4int n = G4int( 0.1 * Z + 1.0 );
00179 G4float x = ( 10 * n - Z ) * 0.1;
00180 return x * t[n - 2] + ( 1.0 - x ) * t[n-1];
00181 }