Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // neutron_hp -- source file 26 // neutron_hp -- source file 27 // J.P. Wellisch, Nov-1996 27 // J.P. Wellisch, Nov-1996 28 // A prototype of the low energy neutron trans 28 // A prototype of the low energy neutron transport model. 29 // P. Arce, June-2014 Conversion neutron_hp to 29 // P. Arce, June-2014 Conversion neutron_hp to particle_hp 30 // 30 // 31 #include "G4ParticleHPMadlandNixSpectrum.hh" 31 #include "G4ParticleHPMadlandNixSpectrum.hh" 32 << 33 #include "G4SystemOfUnits.hh" 32 #include "G4SystemOfUnits.hh" 34 33 35 G4double G4ParticleHPMadlandNixSpectrum::Madla << 34 G4double G4ParticleHPMadlandNixSpectrum::Madland(G4double aSecEnergy, G4double tm) 36 { << 35 { 37 G4Pow* Pow = G4Pow::GetInstance(); << 36 G4Pow* Pow=G4Pow::GetInstance(); 38 G4double result; << 37 G4double result; 39 G4double energy = aSecEnergy / eV; << 38 G4double energy = aSecEnergy/eV; 40 G4double EF; << 39 G4double EF; 41 << 40 42 EF = theAvarageKineticPerNucleonForLightFrag << 41 EF = theAvarageKineticPerNucleonForLightFragments/eV; 43 G4double lightU1 = std::sqrt(energy) - std:: << 42 G4double lightU1 = std::sqrt(energy)-std::sqrt(EF); 44 lightU1 *= lightU1 / tm; << 43 lightU1 *= lightU1/tm; 45 G4double lightU2 = std::sqrt(energy) + std:: << 44 G4double lightU2 = std::sqrt(energy)+std::sqrt(EF); 46 lightU2 *= lightU2 / tm; << 45 lightU2 *= lightU2/tm; 47 G4double lightTerm = 0; << 46 G4double lightTerm=0; 48 if (theAvarageKineticPerNucleonForLightFragm << 47 if(theAvarageKineticPerNucleonForLightFragments>1*eV) 49 lightTerm = Pow->powA(lightU2, 1.5) * E1(l << 48 { 50 lightTerm -= Pow->powA(lightU1, 1.5) * E1( << 49 lightTerm = Pow->powA(lightU2, 1.5)*E1(lightU2); 51 lightTerm += Gamma15(lightU2) - Gamma15(li << 50 lightTerm -= Pow->powA(lightU1, 1.5)*E1(lightU1); 52 lightTerm /= 3. * std::sqrt(tm * EF); << 51 lightTerm += Gamma15(lightU2)-Gamma15(lightU1); >> 52 lightTerm /= 3.*std::sqrt(tm*EF); >> 53 } >> 54 >> 55 EF = theAvarageKineticPerNucleonForHeavyFragments/eV; >> 56 G4double heavyU1 = std::sqrt(energy)-std::sqrt(EF); >> 57 heavyU1 *= heavyU1/tm; >> 58 G4double heavyU2 = std::sqrt(energy)+std::sqrt(EF); >> 59 heavyU2 *= heavyU2/tm; >> 60 G4double heavyTerm=0 ; >> 61 if(theAvarageKineticPerNucleonForHeavyFragments> 1*eV) >> 62 { >> 63 heavyTerm = Pow->powA(heavyU2, 1.5)*E1(heavyU2); >> 64 heavyTerm -= Pow->powA(heavyU1, 1.5)*E1(heavyU1); >> 65 heavyTerm += Gamma15(heavyU2)-Gamma15(heavyU1); >> 66 heavyTerm /= 3.*std::sqrt(tm*EF); >> 67 } >> 68 >> 69 result = 0.5*(lightTerm+heavyTerm); >> 70 >> 71 return result; 53 } 72 } 54 73 55 EF = theAvarageKineticPerNucleonForHeavyFrag << 74 G4double G4ParticleHPMadlandNixSpectrum::Sample(G4double anEnergy) 56 G4double heavyU1 = std::sqrt(energy) - std:: << 75 { 57 heavyU1 *= heavyU1 / tm; << 76 G4double tm = theMaxTemp.GetY(anEnergy); 58 G4double heavyU2 = std::sqrt(energy) + std:: << 77 G4double last=0, buff, current = 100*MeV; 59 heavyU2 *= heavyU2 / tm; << 78 G4double precision = 0.001; 60 G4double heavyTerm = 0; << 79 G4double newValue = 0., oldValue=0.; 61 if (theAvarageKineticPerNucleonForHeavyFragm << 80 G4double random = G4UniformRand(); 62 heavyTerm = Pow->powA(heavyU2, 1.5) * E1(h << 81 63 heavyTerm -= Pow->powA(heavyU1, 1.5) * E1( << 82 G4int icounter=0; 64 heavyTerm += Gamma15(heavyU2) - Gamma15(he << 83 G4int icounter_max=1024; 65 heavyTerm /= 3. * std::sqrt(tm * EF); << 84 do >> 85 { >> 86 icounter++; >> 87 if ( icounter > icounter_max ) { >> 88 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl; >> 89 break; >> 90 } >> 91 oldValue = newValue; >> 92 newValue = FissionIntegral(tm, current); >> 93 if(newValue < random) >> 94 { >> 95 buff = current; >> 96 current+=std::abs(current-last)/2.; >> 97 last = buff; >> 98 if(current>190*MeV) throw G4HadronicException(__FILE__, __LINE__, "Madland-Nix Spectrum has not converged in sampling"); >> 99 } >> 100 else >> 101 { >> 102 buff = current; >> 103 current-=std::abs(current-last)/2.; >> 104 last = buff; >> 105 } >> 106 } >> 107 while (std::abs(oldValue-newValue)>precision*newValue); // Loop checking, 11.05.2015, T. Koi >> 108 return current; 66 } 109 } 67 110 68 result = 0.5 * (lightTerm + heavyTerm); << 111 G4double G4ParticleHPMadlandNixSpectrum:: 69 << 112 GIntegral(G4double tm, G4double anEnergy, G4double aMean) 70 return result; << 113 { 71 } << 114 G4Pow* Pow=G4Pow::GetInstance(); 72 << 115 if(aMean<1*eV) return 0; 73 G4double G4ParticleHPMadlandNixSpectrum::Sampl << 116 G4double b = anEnergy/eV; 74 { << 117 G4double sb = std::sqrt(b); 75 G4double tm = theMaxTemp.GetY(anEnergy); << 118 G4double EF = aMean/eV; 76 G4double last = 0, buff, current = 100 * MeV << 119 77 G4double precision = 0.001; << 120 G4double alpha = std::sqrt(tm); 78 G4double newValue = 0., oldValue = 0.; << 121 G4double beta = std::sqrt(EF); 79 G4double random = G4UniformRand(); << 122 G4double A = EF/tm; 80 << 123 G4double B = (sb+beta)*(sb+beta)/tm; 81 G4int icounter = 0; << 124 G4double Ap = A; 82 G4int icounter_max = 1024; << 125 G4double Bp = (sb-beta)*(sb-beta)/tm; 83 do { << 126 84 icounter++; << 127 G4double result; 85 if (icounter > icounter_max) { << 128 G4double alpha2 = alpha*alpha; 86 G4cout << "Loop-counter exceeded the thr << 129 G4double alphabeta = alpha*beta; 87 << __FILE__ << "." << G4endl; << 130 if(b<EF) 88 break; << 131 { 89 } << 132 result = 90 oldValue = newValue; << 133 ( 91 newValue = FissionIntegral(tm, current); << 134 (0.4*alpha2*Pow->powA(B,2.5) - 0.5*alphabeta*B*B)*E1(B) - 92 if (newValue < random) { << 135 (0.4*alpha2*Pow->powA(A,2.5) - 0.5*alphabeta*A*A)*E1(A) 93 buff = current; << 136 ) 94 current += std::abs(current - last) / 2. << 137 - 95 last = buff; << 138 ( 96 if (current > 190 * MeV) << 139 (0.4*alpha2*Pow->powA(Bp,2.5) + 0.5*alphabeta*Bp*Bp)*E1(Bp) - 97 throw G4HadronicException(__FILE__, __ << 140 (0.4*alpha2*Pow->powA(Ap,2.5) + 0.5*alphabeta*Ap*Ap)*E1(Ap) 98 "Madland-Nix << 141 ) >> 142 + >> 143 ( >> 144 (alpha2*B-2*alphabeta*std::sqrt(B))*Gamma15(B) - >> 145 (alpha2*A-2*alphabeta*std::sqrt(A))*Gamma15(A) >> 146 ) >> 147 - >> 148 ( >> 149 (alpha2*Bp-2*alphabeta*std::sqrt(Bp))*Gamma15(Bp) - >> 150 (alpha2*Ap-2*alphabeta*std::sqrt(Ap))*Gamma15(Ap) >> 151 ) >> 152 - 0.6*alpha2*(Gamma25(B) - Gamma25(A) - Gamma25(Bp) + Gamma25(Ap)) >> 153 - 1.5*alphabeta*(G4Exp(-B)*(1+B) - G4Exp(-A)*(1+A) + G4Exp(-Bp)*(1+Bp) + G4Exp(-Ap)*(1+Ap)) ; 99 } 154 } 100 else { << 155 else 101 buff = current; << 156 { 102 current -= std::abs(current - last) / 2. << 157 result = 103 last = buff; << 158 ( >> 159 (0.4*alpha2*Pow->powA(B,2.5) - 0.5*alphabeta*B*B)*E1(B) - >> 160 (0.4*alpha2*Pow->powA(A,2.5) - 0.5*alphabeta*A*A)*E1(A) >> 161 ); >> 162 result -= >> 163 ( >> 164 (0.4*alpha2*Pow->powA(Bp,2.5) + 0.5*alphabeta*Bp*Bp)*E1(Bp) - >> 165 (0.4*alpha2*Pow->powA(Ap,2.5) + 0.5*alphabeta*Ap*Ap)*E1(Ap) >> 166 ); >> 167 result += >> 168 ( >> 169 (alpha2*B-2*alphabeta*std::sqrt(B))*Gamma15(B) - >> 170 (alpha2*A-2*alphabeta*std::sqrt(A))*Gamma15(A) >> 171 ); >> 172 result -= >> 173 ( >> 174 (alpha2*Bp+2*alphabeta*std::sqrt(Bp))*Gamma15(Bp) - >> 175 (alpha2*Ap+2*alphabeta*std::sqrt(Ap))*Gamma15(Ap) >> 176 ); >> 177 result -= 0.6*alpha2*(Gamma25(B) - Gamma25(A) - Gamma25(Bp) + Gamma25(Ap)); >> 178 result -= 1.5*alphabeta*(G4Exp(-B)*(1+B) - G4Exp(-A)*(1+A) + G4Exp(-Bp)*(1+Bp) + G4Exp(-Ap)*(1+Ap) - 2.) ; 104 } 179 } 105 } while (std::abs(oldValue - newValue) << 180 result = result / (3.*std::sqrt(tm*EF)); 106 > precision * newValue); // Loop c << 181 return result; 107 return current; << 108 } << 109 << 110 G4double G4ParticleHPMadlandNixSpectrum::GInte << 111 { << 112 G4Pow* Pow = G4Pow::GetInstance(); << 113 if (aMean < 1 * eV) return 0; << 114 G4double b = anEnergy / eV; << 115 G4double sb = std::sqrt(b); << 116 G4double EF = aMean / eV; << 117 << 118 G4double alpha = std::sqrt(tm); << 119 G4double beta = std::sqrt(EF); << 120 G4double A = EF / tm; << 121 G4double B = (sb + beta) * (sb + beta) / tm; << 122 G4double Ap = A; << 123 G4double Bp = (sb - beta) * (sb - beta) / tm << 124 << 125 G4double result; << 126 G4double alpha2 = alpha * alpha; << 127 G4double alphabeta = alpha * beta; << 128 if (b < EF) { << 129 result = ((0.4 * alpha2 * Pow->powA(B, 2.5 << 130 - (0.4 * alpha2 * Pow->powA(A, 2 << 131 - ((0.4 * alpha2 * Pow->powA(Bp, << 132 - (0.4 * alpha2 * Pow->powA(Ap << 133 + ((alpha2 * B - 2 * alphabeta * << 134 - (alpha2 * A - 2 * alphabeta << 135 - ((alpha2 * Bp - 2 * alphabeta * << 136 - (alpha2 * Ap - 2 * alphabeta << 137 - 0.6 * alpha2 * (Gamma25(B) - Ga << 138 - 1.5 * alphabeta << 139 * (G4Exp(-B) * (1 + B) - G4Ex << 140 + G4Exp(-Ap) * (1 + Ap)); << 141 } << 142 else { << 143 result = ((0.4 * alpha2 * Pow->powA(B, 2.5 << 144 - (0.4 * alpha2 * Pow->powA(A, 2 << 145 result -= ((0.4 * alpha2 * Pow->powA(Bp, 2 << 146 - (0.4 * alpha2 * Pow->powA(Ap, << 147 result += ((alpha2 * B - 2 * alphabeta * s << 148 - (alpha2 * A - 2 * alphabeta * << 149 result -= ((alpha2 * Bp + 2 * alphabeta * << 150 - (alpha2 * Ap + 2 * alphabeta << 151 result -= 0.6 * alpha2 * (Gamma25(B) - Gam << 152 result -= 1.5 * alphabeta << 153 * (G4Exp(-B) * (1 + B) - G4Exp(- << 154 + G4Exp(-Ap) * (1 + Ap) - 2.) << 155 } 182 } 156 result = result / (3. * std::sqrt(tm * EF)); << 157 return result; << 158 } << 159 183