| 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