00001 // 00002 // ******************************************************************** 00003 // * License and Disclaimer * 00004 // * * 00005 // * The Geant4 software is copyright of the Copyright Holders of * 00006 // * the Geant4 Collaboration. It is provided under the terms and * 00007 // * conditions of the Geant4 Software License, included in the file * 00008 // * LICENSE and available at http://cern.ch/geant4/license . These * 00009 // * include a list of copyright holders. * 00010 // * * 00011 // * Neither the authors of this software system, nor their employing * 00012 // * institutes,nor the agencies providing financial support for this * 00013 // * work make any representation or warranty, express or implied, * 00014 // * regarding this software system or assume any liability for its * 00015 // * use. Please see the license in the file LICENSE and URL above * 00016 // * for the full disclaimer and the limitation of liability. * 00017 // * * 00018 // * This code implementation is the result of the scientific and * 00019 // * technical work of the GEANT4 collaboration. * 00020 // * By using, copying, modifying or distributing the software (or * 00021 // * any work based on the software) you agree to acknowledge its * 00022 // * use in resulting scientific publications, and indicate your * 00023 // * acceptance of all terms of the Geant4 Software license. * 00024 // ******************************************************************** 00025 // 00026 // 00027 // Implementation of the HETC88 code into Geant4. 00028 // Evaporation and De-excitation parts 00029 // T. Lampen, Helsinki Institute of Physics, May-2000 00030 // 00031 // 20120608 M. Kelsey -- Change vars 's','m','m2' to avoid name collisions 00032 00033 #include "globals.hh" 00034 #include "G4ios.hh" 00035 #include "Randomize.hh" 00036 #include "G4SystemOfUnits.hh" 00037 #include "G4Neutron.hh" 00038 #include "G4Proton.hh" 00039 #include "G4Deuteron.hh" 00040 #include "G4Triton.hh" 00041 #include "G4Alpha.hh" 00042 #include "G4ParticleTable.hh" 00043 #include "G4Nucleus.hh" 00044 #include "G4BENeutronChannel.hh" 00045 00046 00047 G4BENeutronChannel::G4BENeutronChannel() 00048 { 00049 name = "neutron"; 00050 particleA = 1; 00051 particleZ = 0; 00052 verboseLevel = 0; 00053 rho = 0; 00054 } 00055 00056 00057 G4BENeutronChannel::~G4BENeutronChannel() 00058 { 00059 } 00060 00061 00062 void G4BENeutronChannel::calculateProbability() 00063 { 00064 const G4int residualZ = nucleusZ - particleZ; 00065 const G4int residualA = nucleusA - particleA; 00066 00067 if ( nucleusA < 2.0 * particleA || 00068 nucleusZ < 2.0 * particleZ || 00069 residualA <= residualZ || 00070 excitationEnergy - getThresh() - correction < 0 ) 00071 { 00072 if ( verboseLevel >= 6 ) 00073 G4cout << "G4BENeutronChannel : calculateProbability = 0 " << G4endl; 00074 emissionProbability = 0; 00075 return; 00076 } 00077 00078 // In HETC88 s-s0 was used in std::exp( s ), in which s0 was either 50 or 00079 // max(s_i), where i goes over all channels. 00080 00081 const G4double levelParam = getLevelDensityParameter(); 00082 00083 const G4double slevel = 2 * std::sqrt( levelParam * ( excitationEnergy - getThresh() - correction ) ); 00084 const G4double eye0 = std::exp( slevel ) * ( slevel - 1 ) / ( 2 * levelParam ); 00085 const G4double eye1 = ( slevel*slevel - 3*slevel +3 ) * std::exp( slevel ) / ( 4 * levelParam*levelParam ) ; 00086 00087 emissionProbability = std::pow( G4double(residualA), 0.666666 ) * alpha() * ( eye1 + beta() * eye0 ); 00088 00089 if ( verboseLevel >= 6 ) 00090 G4cout << "G4BENeutronChannel : calculateProbability " << G4endl 00091 << " res A = " << residualA << G4endl 00092 << " res Z = " << residualZ << G4endl 00093 << " alpha = " << alpha() << G4endl 00094 << " beta = " << beta() << G4endl 00095 << " E = " << excitationEnergy << G4endl 00096 << " correction = " << correction << G4endl 00097 << " eye1 = " << eye1 << G4endl 00098 << " eye0 = " << eye0 << G4endl 00099 << " levelParam = " << levelParam << G4endl 00100 << " thresh = " << getThresh() << G4endl 00101 << " s = " << s << G4endl 00102 << " probability = " << emissionProbability << G4endl; 00103 00104 return; 00105 } 00106 00107 00108 G4double G4BENeutronChannel::sampleKineticEnergy() 00109 { 00111 // A random number is sampled from the density function 00112 // P(x) = x * std::exp ( 2 std::sqrt ( a ( xMax - x ) ) ) [not normalized], 00113 // x belongs to [ 0, xMax ] 00114 // with the 'Hit or Miss' -method 00115 // Kinetic energy is this energy scaled properly 00116 00117 G4double levelParam; 00118 levelParam = getLevelDensityParameter(); 00119 00120 const G4double xMax = excitationEnergy - getThresh() - correction + beta(); // maximum number 00121 const G4double xProb = ( - 1 + std::sqrt ( 1 + 4 * levelParam * xMax ) ) / ( 2 * levelParam ); // most probable value 00122 const G4double maxProb = xProb * std::exp ( 2 * std::sqrt ( levelParam * ( xMax - xProb ) ) ); // maximum value of P(x) 00123 00124 // Sample x according to density function P(x) with rejection method 00125 G4double r1; 00126 G4double r2; 00127 G4int koe=0; 00128 do 00129 { 00130 r1 = beta() + G4UniformRand() * ( xMax - beta() ); 00131 r2 = G4UniformRand() * maxProb; 00132 koe++; 00133 } 00134 while ( r1 * std::exp ( 2 * std::sqrt ( levelParam * ( xMax - r1 ) ) ) < r2 ); 00135 00136 // G4cout << koe << G4endl; 00137 G4double kineticEnergy = r1 - beta(); 00138 00139 if ( verboseLevel >= 10 ) 00140 G4cout << " G4BENeutronChannel : sampleKineticEnergy() " << G4endl 00141 << " kinetic n e = " << kineticEnergy << G4endl 00142 << " levelParam = " << levelParam << G4endl 00143 << " thresh= " << getThresh() << G4endl 00144 << " beta= " << beta() << G4endl; 00145 00146 return kineticEnergy; 00147 } 00148 00149 00150 G4DynamicParticle * G4BENeutronChannel::emit() 00151 { 00152 G4double u; 00153 G4double v; 00154 G4double w; 00155 G4DynamicParticle * pParticle = new G4DynamicParticle; 00156 00157 pParticle -> SetDefinition( G4Neutron::Neutron() ); 00158 pParticle -> SetKineticEnergy( sampleKineticEnergy() ); 00159 isotropicCosines( u, v, w ); 00160 pParticle -> SetMomentumDirection( u , v , w ); 00161 00162 return pParticle; 00163 } 00164 00165 00166 G4double G4BENeutronChannel::alpha() 00167 { 00168 const G4double residualA = nucleusA - particleA; 00169 return 0.76 + 1.93 * std::pow( residualA, -0.33333 ); 00170 } 00171 00172 00173 G4double G4BENeutronChannel::beta() 00174 { 00175 G4double residualA = nucleusA - particleA; 00176 return ( 1.66 * std::pow ( residualA, -0.66666 ) - 0.05 )/alpha()*MeV; 00177 } 00178