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 // $Id$ 00027 // 00028 //--------------------------------------------------------------------- 00029 // 00030 // Geant4 class G4E1Probability 00031 // 00032 // by V. Lara (May 2003) 00033 // 00034 // Modifications: 00035 // 18.05.2010 V.Ivanchenko trying to speedup the most slow method 00036 // by usage of G4Pow, integer A and introduction of const members 00037 // 17.11.2010 V.Ivanchenko perform general cleanup and simplification 00038 // of integration method; low-limit of integration is defined 00039 // by gamma energy or is zero (was always zero before) 00040 // 00041 00042 #include "G4E1Probability.hh" 00043 #include "Randomize.hh" 00044 #include "G4Pow.hh" 00045 #include "G4SystemOfUnits.hh" 00046 00047 // Constructors and operators 00048 // 00049 00050 G4E1Probability::G4E1Probability():G4VEmissionProbability() 00051 { 00052 G4double x = CLHEP::pi*CLHEP::hbarc; 00053 normC = 1.0 / (x*x); 00054 theLevelDensityParameter = 0.125/MeV; 00055 fG4pow = G4Pow::GetInstance(); 00056 } 00057 00058 G4E1Probability::~G4E1Probability() 00059 {} 00060 00061 // Calculate the emission probability 00062 // 00063 00064 G4double G4E1Probability::EmissionProbDensity(const G4Fragment& frag, 00065 G4double gammaE) 00066 { 00067 00068 // Calculate the probability density here 00069 00070 // From nuclear fragment properties and the excitation energy, calculate 00071 // the probability density for photon evaporation from U to U - gammaE 00072 // (U = nucleus excitation energy, gammaE = total evaporated photon 00073 // energy). Fragment = nuclear fragment BEFORE de-excitation 00074 00075 G4double theProb = 0.0; 00076 00077 G4int Afrag = frag.GetA_asInt(); 00078 G4double Uexcite = frag.GetExcitationEnergy(); 00079 G4double U = std::max(0.0,Uexcite-gammaE); 00080 00081 if(gammaE < 0.0) { return theProb; } 00082 00083 // Need a level density parameter. 00084 // For now, just use the constant approximation (not reliable near magic 00085 // nuclei). 00086 00087 G4double aLevelDensityParam = Afrag*theLevelDensityParameter; 00088 00089 // G4double levelDensBef = std::exp(2*std::sqrt(aLevelDensityParam*Uexcite)); 00090 // G4double levelDensAft = std::exp(2*std::sqrt(aLevelDensityParam*(Uexcite-gammaE))); 00091 // VI reduce number of calls to exp 00092 G4double levelDens = 00093 std::exp(2*(std::sqrt(aLevelDensityParam*U)-std::sqrt(aLevelDensityParam*Uexcite))); 00094 // Now form the probability density 00095 00096 // Define constants for the photoabsorption cross-section (the reverse 00097 // process of our de-excitation) 00098 00099 G4double sigma0 = 2.5 * Afrag * millibarn; // millibarns 00100 00101 G4double Egdp = (40.3 / fG4pow->powZ(Afrag,0.2) )*MeV; 00102 G4double GammaR = 0.30 * Egdp; 00103 00104 // CD 00105 //cout<<" PROB TESTS "<<G4endl; 00106 //cout<<" hbarc = "<<hbarc<<G4endl; 00107 //cout<<" pi = "<<pi<<G4endl; 00108 //cout<<" Uexcite, gammaE = "<<Uexcite<<" "<<gammaE<<G4endl; 00109 //cout<<" Uexcite, gammaE = "<<Uexcite*MeV<<" "<<gammaE*MeV<<G4endl; 00110 //cout<<" lev density param = "<<aLevelDensityParam<<G4endl; 00111 //cout<<" level densities = "<<levelDensBef<<" "<<levelDensAft<<G4endl; 00112 //cout<<" sigma0 = "<<sigma0<<G4endl; 00113 //cout<<" Egdp, GammaR = "<<Egdp<<" "<<GammaR<<G4endl; 00114 //cout<<" normC = "<<normC<<G4endl; 00115 00116 // VI implementation 18.05.2010 00117 G4double gammaE2 = gammaE*gammaE; 00118 G4double gammaR2 = gammaE2*GammaR*GammaR; 00119 G4double egdp2 = gammaE2 - Egdp*Egdp; 00120 G4double sigmaAbs = sigma0*gammaR2/(egdp2*egdp2 + gammaR2); 00121 theProb = normC * sigmaAbs * gammaE2 * levelDens; 00122 00123 // old implementation 00124 // G4double numerator = sigma0 * gammaE*gammaE * GammaR*GammaR; 00125 // G4double denominator = (gammaE*gammaE - Egdp*Egdp)* 00126 // (gammaE*gammaE - Egdp*Egdp) + GammaR*GammaR*gammaE*gammaE; 00127 00128 //G4double sigmaAbs = numerator/denominator; 00129 //theProb = normC * sigmaAbs * gammaE2 * levelDensAft/levelDensBef; 00130 00131 // CD 00132 //cout<<" sigmaAbs = "<<sigmaAbs<<G4endl; 00133 //cout<<" Probability = "<<theProb<<G4endl; 00134 00135 return theProb; 00136 00137 } 00138 00139 G4double G4E1Probability::EmissionProbability(const G4Fragment& frag, 00140 G4double gammaE) 00141 { 00142 // From nuclear fragment properties and the excitation energy, calculate 00143 // the probability for photon evaporation down to last ground level. 00144 // fragment = nuclear fragment BEFORE de-excitation 00145 00146 G4double upperLim = gammaE; 00147 G4double lowerLim = 0.0; 00148 00149 //G4cout << "G4E1Probability::EmissionProbability: Emin= " << lowerLim 00150 // << " Emax= " << upperLim << G4endl; 00151 if( upperLim - lowerLim <= CLHEP::keV ) { return 0.0; } 00152 00153 // Need to integrate EmissionProbDensity from lowerLim to upperLim 00154 // and multiply by factor 3 (?!) 00155 00156 G4double integ = 3.0 * EmissionIntegration(frag,lowerLim,upperLim); 00157 00158 return integ; 00159 00160 } 00161 00162 G4double G4E1Probability::EmissionIntegration(const G4Fragment& frag, 00163 G4double lowLim, G4double upLim) 00164 00165 { 00166 // Simple integration 00167 // VI replace by direct integration over 100 point 00168 00169 const G4int numIters = 100; 00170 G4double Step = (upLim-lowLim)/G4double(numIters); 00171 00172 G4double res = 0.0; 00173 G4double x = lowLim - 0.5*Step; 00174 00175 for(G4int i = 0; i < numIters; ++i) { 00176 x += Step; 00177 res += EmissionProbDensity(frag, x); 00178 } 00179 00180 if(res > 0.0) { res /= G4double(numIters); } 00181 else { res = 0.0; } 00182 00183 return res; 00184 00185 } 00186 00187