| 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 // 26 //
27 // 27 //
28 // original by H.P. Wellisch << 28 // original by H.P. Wellisch
29 // modified by J.L. Chuma, TRIUMF, 19-Nov-1996 << 29 // modified by J.L. Chuma, TRIUMF, 19-Nov-1996
30 // last modified: 27-Mar-1997 << 30 // last modified: 27-Mar-1997
31 // J.P.Wellisch: 23-Apr-97: minor simplificati << 31 // J.P.Wellisch: 23-Apr-97: minor simplifications
32 // modified by J.L.Chuma 24-Jul-97 to set the << 32 // modified by J.L.Chuma 24-Jul-97 to set the total momentum in Cinema and
33 // Evaporatio << 33 // EvaporationEffects
34 // modified by J.L.Chuma 21-Oct-97 put std::a << 34 // modified by J.L.Chuma 21-Oct-97 put std::abs() around the totalE^2-mass^2
35 // in calcula << 35 // in calculation of total momentum in
36 // Cinema and << 36 // Cinema and EvaporationEffects
37 // Chr. Volcker, 10-Nov-1997: new methods and << 37 // Chr. Volcker, 10-Nov-1997: new methods and class variables.
38 // HPW added utilities for low energy neutron << 38 // HPW added utilities for low energy neutron transport. (12.04.1998)
39 // M.G. Pia, 2 Oct 1998: modified GetFermiMome << 39 // M.G. Pia, 2 Oct 1998: modified GetFermiMomentum to avoid memory leaks
40 // G.Folger, spring 2010: add integer A/Z int << 40 // G.Folger, spring 2010: add integer A/Z interface
41 // A. Ribon, summer 2015: migrated to G4Exp a << 41 // A. Ribon, 6 August 2015: migrated to G4Exp and G4Log.
42 // A. Ribon, autumn 2021: extended to hypernu <<
43 42
44 #include "G4Nucleus.hh" 43 #include "G4Nucleus.hh"
45 #include "G4NucleiProperties.hh" 44 #include "G4NucleiProperties.hh"
46 #include "G4PhysicalConstants.hh" 45 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh" 46 #include "G4SystemOfUnits.hh"
48 #include "Randomize.hh" 47 #include "Randomize.hh"
49 #include "G4HadronicException.hh" 48 #include "G4HadronicException.hh"
>> 49
50 #include "G4Exp.hh" 50 #include "G4Exp.hh"
51 #include "G4Log.hh" 51 #include "G4Log.hh"
52 #include "G4HyperNucleiProperties.hh" <<
53 #include "G4HadronicParameters.hh" <<
54 <<
55 52
>> 53
56 G4Nucleus::G4Nucleus() 54 G4Nucleus::G4Nucleus()
57 : theA(0), theZ(0), theL(0), aEff(0.0), zEff << 55 : theA(0), theZ(0), aEff(0.0), zEff(0)
58 { 56 {
59 pnBlackTrackEnergy = 0.0; 57 pnBlackTrackEnergy = 0.0;
60 dtaBlackTrackEnergy = 0.0; 58 dtaBlackTrackEnergy = 0.0;
61 pnBlackTrackEnergyfromAnnihilation = 0.0; 59 pnBlackTrackEnergyfromAnnihilation = 0.0;
62 dtaBlackTrackEnergyfromAnnihilation = 0.0; 60 dtaBlackTrackEnergyfromAnnihilation = 0.0;
63 excitationEnergy = 0.0; 61 excitationEnergy = 0.0;
64 momentum = G4ThreeVector(0.,0.,0.); 62 momentum = G4ThreeVector(0.,0.,0.);
65 fermiMomentum = 1.52*hbarc/fermi; 63 fermiMomentum = 1.52*hbarc/fermi;
66 theTemp = 293.16*kelvin; 64 theTemp = 293.16*kelvin;
67 fIsotope = 0; 65 fIsotope = 0;
68 } 66 }
69 67
70 G4Nucleus::G4Nucleus( const G4double A, const << 68 G4Nucleus::G4Nucleus( const G4double A, const G4double Z )
71 { 69 {
72 SetParameters( A, Z, std::max(numberOfLambda << 70 SetParameters( A, Z );
73 pnBlackTrackEnergy = 0.0; 71 pnBlackTrackEnergy = 0.0;
74 dtaBlackTrackEnergy = 0.0; 72 dtaBlackTrackEnergy = 0.0;
75 pnBlackTrackEnergyfromAnnihilation = 0.0; 73 pnBlackTrackEnergyfromAnnihilation = 0.0;
76 dtaBlackTrackEnergyfromAnnihilation = 0.0; 74 dtaBlackTrackEnergyfromAnnihilation = 0.0;
77 excitationEnergy = 0.0; 75 excitationEnergy = 0.0;
78 momentum = G4ThreeVector(0.,0.,0.); 76 momentum = G4ThreeVector(0.,0.,0.);
79 fermiMomentum = 1.52*hbarc/fermi; 77 fermiMomentum = 1.52*hbarc/fermi;
80 theTemp = 293.16*kelvin; 78 theTemp = 293.16*kelvin;
81 fIsotope = 0; 79 fIsotope = 0;
82 } 80 }
83 81
84 G4Nucleus::G4Nucleus( const G4int A, const G4i << 82 G4Nucleus::G4Nucleus( const G4int A, const G4int Z )
85 { 83 {
86 SetParameters( A, Z, std::max(numberOfLambda << 84 SetParameters( A, Z );
87 pnBlackTrackEnergy = 0.0; 85 pnBlackTrackEnergy = 0.0;
88 dtaBlackTrackEnergy = 0.0; 86 dtaBlackTrackEnergy = 0.0;
89 pnBlackTrackEnergyfromAnnihilation = 0.0; 87 pnBlackTrackEnergyfromAnnihilation = 0.0;
90 dtaBlackTrackEnergyfromAnnihilation = 0.0; 88 dtaBlackTrackEnergyfromAnnihilation = 0.0;
91 excitationEnergy = 0.0; 89 excitationEnergy = 0.0;
92 momentum = G4ThreeVector(0.,0.,0.); 90 momentum = G4ThreeVector(0.,0.,0.);
93 fermiMomentum = 1.52*hbarc/fermi; 91 fermiMomentum = 1.52*hbarc/fermi;
94 theTemp = 293.16*kelvin; 92 theTemp = 293.16*kelvin;
95 fIsotope = 0; 93 fIsotope = 0;
96 } 94 }
97 95
98 G4Nucleus::G4Nucleus( const G4Material *aMater 96 G4Nucleus::G4Nucleus( const G4Material *aMaterial )
99 { 97 {
100 ChooseParameters( aMaterial ); 98 ChooseParameters( aMaterial );
101 pnBlackTrackEnergy = 0.0; 99 pnBlackTrackEnergy = 0.0;
102 dtaBlackTrackEnergy = 0.0; 100 dtaBlackTrackEnergy = 0.0;
103 pnBlackTrackEnergyfromAnnihilation = 0.0; 101 pnBlackTrackEnergyfromAnnihilation = 0.0;
104 dtaBlackTrackEnergyfromAnnihilation = 0.0; 102 dtaBlackTrackEnergyfromAnnihilation = 0.0;
105 excitationEnergy = 0.0; 103 excitationEnergy = 0.0;
106 momentum = G4ThreeVector(0.,0.,0.); 104 momentum = G4ThreeVector(0.,0.,0.);
107 fermiMomentum = 1.52*hbarc/fermi; 105 fermiMomentum = 1.52*hbarc/fermi;
108 theTemp = aMaterial->GetTemperature(); 106 theTemp = aMaterial->GetTemperature();
109 fIsotope = 0; 107 fIsotope = 0;
110 } 108 }
111 109
112 G4Nucleus::~G4Nucleus() {} 110 G4Nucleus::~G4Nucleus() {}
113 111
114 << 112 G4ReactionProduct G4Nucleus::
115 //-------------------------------------------- << 113 GetBiasedThermalNucleus(G4double aMass, G4ThreeVector aVelocity, G4double temp) const
116 // SVT (Sampling of the Velocity of the Target <<
117 //-------------------------------------------- <<
118 G4ReactionProduct <<
119 G4Nucleus::GetBiasedThermalNucleus(G4double aM <<
120 { 114 {
121 // If E_neutron <= E_threshold, Then apply t << 115 G4double velMag = aVelocity.mag();
122 // Else consider the target nucleus being wi <<
123 G4double E_threshold = G4HadronicParameters: <<
124 if ( E_threshold == -1. ) { <<
125 E_threshold = 400.0*8.617333262E-11*temp; <<
126 } <<
127 G4double E_neutron = 0.5*aVelocity.mag2()*G4 <<
128 <<
129 G4ReactionProduct result; 116 G4ReactionProduct result;
130 result.SetMass(aMass*G4Neutron::Neutron()->G << 117 G4double value = 0;
131 << 118 G4double random = 1;
132 if ( E_neutron <= E_threshold ) { << 119 G4double norm = 3.*std::sqrt(k_Boltzmann*temp*aMass*G4Neutron::Neutron()->GetPDGMass());
133 << 120 norm /= G4Neutron::Neutron()->GetPDGMass();
134 // Beta = sqrt(m/2kT) << 121 norm *= 5.;
135 G4double beta = std::sqrt(result.GetMass() << 122 norm += velMag;
136 << 123 norm /= velMag;
137 // Neutron speed vn << 124 const G4int maxNumberOfLoops = 1000000;
138 G4double vN_norm = aVelocity.mag(); << 125 G4int loopCounter = -1;
139 G4double vN_norm2 = vN_norm*vN_norm; << 126 while ( (value/norm<random) && ++loopCounter < maxNumberOfLoops ) /* Loop checking, 02.11.2015, A.Ribon */
140 G4double y = beta*vN_norm; << 127 {
141 << 128 result = GetThermalNucleus(aMass, temp);
142 // Normalize neutron velocity << 129 G4ThreeVector targetVelocity = 1./result.GetMass()*result.GetMomentum();
143 aVelocity = (1./vN_norm)*aVelocity; << 130 value = (targetVelocity+aVelocity).mag()/velMag;
144 << 131 random = G4UniformRand();
145 // Sample target speed << 132 }
146 G4double x2; << 133 if ( loopCounter >= maxNumberOfLoops ) {
147 G4double randThreshold; << 134 G4ExceptionDescription ed;
148 G4double vT_norm, vT_norm2, mu; //theta, v << 135 ed << " Failed sampling after maxNumberOfLoops attempts : forced exit! " << G4endl;
149 G4double acceptThreshold; << 136 G4Exception( " G4Nucleus::GetBiasedThermalNucleus ", "HAD_NUCLEUS_001", JustWarning, ed );
150 G4double vRelativeSpeed; << 137 result = GetThermalNucleus(aMass, temp);
151 G4double cdf0 = 2./(2.+std::sqrt(CLHEP::pi <<
152 <<
153 do { <<
154 // Sample the target velocity vT in the <<
155 if ( G4UniformRand() < cdf0 ) { <<
156 // Sample in C45 from https://laws.lan <<
157 x2 = -std::log(G4UniformRand()*G4Unifo <<
158 } else { <<
159 // Sample in C61 from https://laws.lan <<
160 G4double ampl = std::cos(CLHEP::pi/2.0 <<
161 x2 = -std::log(G4UniformRand()) - std: <<
162 } <<
163 <<
164 vT_norm = std::sqrt(x2)/beta; <<
165 vT_norm2 = vT_norm*vT_norm; <<
166 <<
167 // Sample cosine between the incident ne <<
168 mu = 2*G4UniformRand() - 1; <<
169 <<
170 // Define acceptance threshold <<
171 vRelativeSpeed = std::sqrt(vN_norm2 + vT <<
172 acceptThreshold = vRelativeSpeed/(vN_nor <<
173 randThreshold = G4UniformRand(); <<
174 } while ( randThreshold >= acceptThreshold <<
175 <<
176 DoKinematicsOfThermalNucleus(mu, vT_norm, <<
177 <<
178 } else { // target nucleus considered as bei <<
179 <<
180 result.SetMomentum(0., 0., 0.); <<
181 result.SetKineticEnergy(0.); <<
182 <<
183 } 138 }
184 <<
185 return result; 139 return result;
186 } 140 }
187 141
188 <<
189 void <<
190 G4Nucleus::DoKinematicsOfThermalNucleus(const <<
191 G4Reac <<
192 <<
193 // Get target nucleus direction from the neu <<
194 G4double cosTh = mu; <<
195 G4ThreeVector uNorm = aVelocity; <<
196 <<
197 G4double sinTh = std::sqrt(1. - cosTh*cosTh) <<
198 <<
199 // Sample randomly the phi angle between the <<
200 G4double phi = CLHEP::twopi*G4UniformRand(); <<
201 G4double sinPhi = std::sin(phi); <<
202 G4double cosPhi = std::cos(phi); <<
203 <<
204 // Find orthogonal vector to aVelocity - sol <<
205 G4ThreeVector ortho(1., 1., 1.); <<
206 if ( uNorm[0] ) ortho[0] = -(uNorm[1]+ <<
207 else if ( uNorm[1] ) ortho[1] = -(uNorm[0]+ <<
208 else if ( uNorm[2] ) ortho[2] = -(uNorm[0]+ <<
209 <<
210 // Normalize the vector <<
211 ortho = (1/ortho.mag())*ortho; <<
212 <<
213 // Find vector to draw a plan perpendicular <<
214 G4ThreeVector orthoComp( uNorm[1]*ortho[2] - <<
215 uNorm[2]*ortho[0] - <<
216 uNorm[0]*ortho[1] - <<
217 <<
218 // Find the direction of the target velocity <<
219 G4ThreeVector directionTarget( cosTh*uNorm[0 <<
220 cosTh*uNorm[1 <<
221 cosTh*uNorm[2 <<
222 <<
223 // Normalize directionTarget <<
224 directionTarget = ( 1./directionTarget.mag() <<
225 <<
226 // Set momentum <<
227 G4double px = result.GetMass()*vT_norm*direc <<
228 G4double py = result.GetMass()*vT_norm*direc <<
229 G4double pz = result.GetMass()*vT_norm*direc <<
230 result.SetMomentum(px, py, pz); <<
231 <<
232 G4double tMom = std::sqrt(px*px+py*py+pz*pz) <<
233 G4double tEtot = std::sqrt( (tMom+result.Get <<
234 - 2.*tMom*result.GetMass() <<
235 <<
236 if ( tEtot/result.GetMass() - 1. > 0.001 ) { <<
237 // use relativistic energy for higher ener <<
238 result.SetTotalEnergy(tEtot); <<
239 } else { <<
240 // use p**2/2M for lower energies (to pres <<
241 result.SetKineticEnergy(tMom*tMom/(2.*resu <<
242 } <<
243 <<
244 } <<
245 <<
246 <<
247 G4ReactionProduct 142 G4ReactionProduct
248 G4Nucleus::GetThermalNucleus(G4double targetMa 143 G4Nucleus::GetThermalNucleus(G4double targetMass, G4double temp) const
249 { 144 {
250 G4double currentTemp = temp; 145 G4double currentTemp = temp;
251 if (currentTemp < 0) currentTemp = theTemp; 146 if (currentTemp < 0) currentTemp = theTemp;
252 G4ReactionProduct theTarget; 147 G4ReactionProduct theTarget;
253 theTarget.SetMass(targetMass*G4Neutron::Neut 148 theTarget.SetMass(targetMass*G4Neutron::Neutron()->GetPDGMass());
254 G4double px, py, pz; 149 G4double px, py, pz;
255 px = GetThermalPz(theTarget.GetMass(), curre 150 px = GetThermalPz(theTarget.GetMass(), currentTemp);
256 py = GetThermalPz(theTarget.GetMass(), curre 151 py = GetThermalPz(theTarget.GetMass(), currentTemp);
257 pz = GetThermalPz(theTarget.GetMass(), curre 152 pz = GetThermalPz(theTarget.GetMass(), currentTemp);
258 theTarget.SetMomentum(px, py, pz); 153 theTarget.SetMomentum(px, py, pz);
259 G4double tMom = std::sqrt(px*px+py*py+pz*pz) 154 G4double tMom = std::sqrt(px*px+py*py+pz*pz);
260 G4double tEtot = std::sqrt((tMom+theTarget.G 155 G4double tEtot = std::sqrt((tMom+theTarget.GetMass())*
261 (tMom+theTarget.G 156 (tMom+theTarget.GetMass())-
262 2.*tMom*theTarge 157 2.*tMom*theTarget.GetMass());
263 // if(1-tEtot/theTarget.GetMass()>0.001) t 158 // if(1-tEtot/theTarget.GetMass()>0.001) this line incorrect (Bug report 1911)
264 if (tEtot/theTarget.GetMass() - 1. > 0.001) 159 if (tEtot/theTarget.GetMass() - 1. > 0.001) {
265 // use relativistic energy for higher ener 160 // use relativistic energy for higher energies
266 theTarget.SetTotalEnergy(tEtot); 161 theTarget.SetTotalEnergy(tEtot);
267 162
268 } else { 163 } else {
269 // use p**2/2M for lower energies (to pres 164 // use p**2/2M for lower energies (to preserve precision?)
270 theTarget.SetKineticEnergy(tMom*tMom/(2.*t 165 theTarget.SetKineticEnergy(tMom*tMom/(2.*theTarget.GetMass()));
271 } 166 }
272 return theTarget; 167 return theTarget;
273 } 168 }
274 169
275 170
276 void 171 void
277 G4Nucleus::ChooseParameters(const G4Material* 172 G4Nucleus::ChooseParameters(const G4Material* aMaterial)
278 { 173 {
279 G4double random = G4UniformRand(); 174 G4double random = G4UniformRand();
280 G4double sum = aMaterial->GetTotNbOfAtomsPer 175 G4double sum = aMaterial->GetTotNbOfAtomsPerVolume();
281 const G4ElementVector* theElementVector = aM 176 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
282 G4double running(0); 177 G4double running(0);
283 // G4Element* element(0); 178 // G4Element* element(0);
284 const G4Element* element = (*theElementVecto << 179 G4Element* element = (*theElementVector)[aMaterial->GetNumberOfElements()-1];
285 180
286 for (unsigned int i = 0; i < aMaterial->GetN 181 for (unsigned int i = 0; i < aMaterial->GetNumberOfElements(); ++i) {
287 running += aMaterial->GetVecNbOfAtomsPerVo 182 running += aMaterial->GetVecNbOfAtomsPerVolume()[i];
288 if (running > random*sum) { 183 if (running > random*sum) {
289 element = (*theElementVector)[i]; 184 element = (*theElementVector)[i];
290 break; 185 break;
291 } 186 }
292 } 187 }
293 188
294 if (element->GetNumberOfIsotopes() > 0) { 189 if (element->GetNumberOfIsotopes() > 0) {
295 G4double randomAbundance = G4UniformRand() 190 G4double randomAbundance = G4UniformRand();
296 G4double sumAbundance = element->GetRelati 191 G4double sumAbundance = element->GetRelativeAbundanceVector()[0];
297 unsigned int iso=0; 192 unsigned int iso=0;
298 while (iso < element->GetNumberOfIsotopes( 193 while (iso < element->GetNumberOfIsotopes() && /* Loop checking, 02.11.2015, A.Ribon */
299 sumAbundance < randomAbundance) { 194 sumAbundance < randomAbundance) {
300 ++iso; 195 ++iso;
301 sumAbundance += element->GetRelativeAbun 196 sumAbundance += element->GetRelativeAbundanceVector()[iso];
302 } 197 }
303 theA=element->GetIsotope(iso)->GetN(); 198 theA=element->GetIsotope(iso)->GetN();
304 theZ=element->GetIsotope(iso)->GetZ(); 199 theZ=element->GetIsotope(iso)->GetZ();
305 theL=0; <<
306 aEff=theA; 200 aEff=theA;
307 zEff=theZ; 201 zEff=theZ;
308 } else { 202 } else {
309 aEff = element->GetN(); 203 aEff = element->GetN();
310 zEff = element->GetZ(); 204 zEff = element->GetZ();
311 theZ = G4int(zEff + 0.5); 205 theZ = G4int(zEff + 0.5);
312 theA = G4int(aEff + 0.5); << 206 theA = G4int(aEff + 0.5);
313 theL=0; <<
314 } 207 }
315 } 208 }
316 209
317 210
318 void 211 void
319 G4Nucleus::SetParameters( const G4double A, co << 212 G4Nucleus::SetParameters(G4double A, G4double Z)
320 { 213 {
321 theZ = G4lrint(Z); 214 theZ = G4lrint(Z);
322 theA = G4lrint(A); << 215 theA = G4lrint(A);
323 theL = std::max(numberOfLambdas, 0); <<
324 if (theA<1 || theZ<0 || theZ>theA) { 216 if (theA<1 || theZ<0 || theZ>theA) {
325 throw G4HadronicException(__FILE__, __LINE 217 throw G4HadronicException(__FILE__, __LINE__,
326 "G4Nucleus::SetParameters called w 218 "G4Nucleus::SetParameters called with non-physical parameters");
327 } 219 }
328 aEff = A; // atomic weight 220 aEff = A; // atomic weight
329 zEff = Z; // atomic number 221 zEff = Z; // atomic number
330 fIsotope = 0; 222 fIsotope = 0;
331 } 223 }
332 224
333 <<
334 void 225 void
335 G4Nucleus::SetParameters( const G4int A, const << 226 G4Nucleus::SetParameters(G4int A, const G4int Z )
336 { 227 {
337 theZ = Z; 228 theZ = Z;
338 theA = A; << 229 theA = A;
339 theL = std::max(numberOfLambdas, 0); <<
340 if( theA<1 || theZ<0 || theZ>theA ) 230 if( theA<1 || theZ<0 || theZ>theA )
341 { 231 {
342 throw G4HadronicException(__FILE__, __LI 232 throw G4HadronicException(__FILE__, __LINE__,
343 "G4Nucleus::SetParameters called with 233 "G4Nucleus::SetParameters called with non-physical parameters");
344 } 234 }
345 aEff = A; // atomic weight 235 aEff = A; // atomic weight
346 zEff = Z; // atomic number 236 zEff = Z; // atomic number
347 fIsotope = 0; 237 fIsotope = 0;
348 } 238 }
349 239
350 << 240 G4DynamicParticle *
351 G4DynamicParticle * << 241 G4Nucleus::ReturnTargetParticle() const
352 G4Nucleus::ReturnTargetParticle() const << 242 {
353 { << 243 // choose a proton or a neutron as the target particle
354 // choose a proton or a neutron (or a lamba << 244
355 G4DynamicParticle *targetParticle = new G4Dy << 245 G4DynamicParticle *targetParticle = new G4DynamicParticle;
356 const G4double rnd = G4UniformRand(); << 246 if( G4UniformRand() < zEff/aEff )
357 if ( rnd < zEff/aEff ) { << 247 targetParticle->SetDefinition( G4Proton::Proton() );
358 targetParticle->SetDefinition( G4Proton::P << 248 else
359 } else if ( rnd < (zEff + theL*1.0)/aEff ) { << 249 targetParticle->SetDefinition( G4Neutron::Neutron() );
360 targetParticle->SetDefinition( G4Lambda::L << 250 return targetParticle;
361 } else { <<
362 targetParticle->SetDefinition( G4Neutron:: <<
363 } 251 }
364 return targetParticle; <<
365 } <<
366 <<
367 252
368 G4double << 253 G4double
369 G4Nucleus::AtomicMass( const G4double A, const << 254 G4Nucleus::AtomicMass( const G4double A, const G4double Z ) const
370 { << 255 {
371 // Now returns (atomic mass - electron masse << 256 // Now returns (atomic mass - electron masses)
372 if ( numberOfLambdas > 0 ) { <<
373 return G4HyperNucleiProperties::GetNuclear <<
374 } else { <<
375 return G4NucleiProperties::GetNuclearMass( 257 return G4NucleiProperties::GetNuclearMass(A, Z);
376 } 258 }
377 } <<
378 <<
379 259
380 G4double << 260 G4double
381 G4Nucleus::AtomicMass( const G4int A, const G4 << 261 G4Nucleus::AtomicMass( const G4int A, const G4int Z ) const
382 { << 262 {
383 // Now returns (atomic mass - electron masse << 263 // Now returns (atomic mass - electron masses)
384 if ( numberOfLambdas > 0 ) { <<
385 return G4HyperNucleiProperties::GetNuclear <<
386 } else { <<
387 return G4NucleiProperties::GetNuclearMass( 264 return G4NucleiProperties::GetNuclearMass(A, Z);
388 } 265 }
389 } <<
390 266
>> 267 G4double
>> 268 G4Nucleus::GetThermalPz( const G4double mass, const G4double temp ) const
>> 269 {
>> 270 G4double result = G4RandGauss::shoot();
>> 271 result *= std::sqrt(k_Boltzmann*temp*mass); // Das ist impuls (Pz),
>> 272 // nichtrelativistische rechnung
>> 273 // Maxwell verteilung angenommen
>> 274 return result;
>> 275 }
391 276
392 G4double << 277 G4double
393 G4Nucleus::GetThermalPz( const G4double mass, << 278 G4Nucleus::EvaporationEffects( G4double kineticEnergy )
394 { << 279 {
395 G4double result = G4RandGauss::shoot(); << 280 // derived from original FORTRAN code EXNU by H. Fesefeldt (10-Dec-1986)
396 result *= std::sqrt(k_Boltzmann*temp*mass); << 281 //
397 // ni << 282 // Nuclear evaporation as function of atomic number
398 // Ma << 283 // and kinetic energy (MeV) of primary particle
399 return result; << 284 //
400 } << 285 // returns kinetic energy (MeV)
>> 286 //
>> 287 if( aEff < 1.5 )
>> 288 {
>> 289 pnBlackTrackEnergy = dtaBlackTrackEnergy = 0.0;
>> 290 return 0.0;
>> 291 }
>> 292 G4double ek = kineticEnergy/GeV;
>> 293 G4float ekin = std::min( 4.0, std::max( 0.1, ek ) );
>> 294 const G4float atno = std::min( 120., aEff );
>> 295 const G4float gfa = 2.0*((aEff-1.0)/70.)*G4Exp(-(aEff-1.0)/70.);
>> 296 //
>> 297 // 0.35 value at 1 GeV
>> 298 // 0.05 value at 0.1 GeV
>> 299 //
>> 300 G4float cfa = std::max( 0.15, 0.35 + ((0.35-0.05)/2.3)*G4Log(ekin) );
>> 301 G4float exnu = 7.716 * cfa * G4Exp(-cfa)
>> 302 * ((atno-1.0)/120.)*G4Exp(-(atno-1.0)/120.);
>> 303 G4float fpdiv = std::max( 0.5, 1.0-0.25*ekin*ekin );
>> 304 //
>> 305 // pnBlackTrackEnergy is the kinetic energy (in GeV) available for
>> 306 // proton/neutron black track particles
>> 307 // dtaBlackTrackEnergy is the kinetic energy (in GeV) available for
>> 308 // deuteron/triton/alpha black track particles
>> 309 //
>> 310 pnBlackTrackEnergy = exnu*fpdiv;
>> 311 dtaBlackTrackEnergy = exnu*(1.0-fpdiv);
>> 312
>> 313 if( G4int(zEff+0.1) != 82 )
>> 314 {
>> 315 G4double ran1 = -6.0;
>> 316 G4double ran2 = -6.0;
>> 317 for( G4int i=0; i<12; ++i )
>> 318 {
>> 319 ran1 += G4UniformRand();
>> 320 ran2 += G4UniformRand();
>> 321 }
>> 322 pnBlackTrackEnergy *= 1.0 + ran1*gfa;
>> 323 dtaBlackTrackEnergy *= 1.0 + ran2*gfa;
>> 324 }
>> 325 pnBlackTrackEnergy = std::max( 0.0, pnBlackTrackEnergy );
>> 326 dtaBlackTrackEnergy = std::max( 0.0, dtaBlackTrackEnergy );
>> 327 while( pnBlackTrackEnergy+dtaBlackTrackEnergy >= ek ) /* Loop checking, 02.11.2015, A.Ribon */
>> 328 {
>> 329 pnBlackTrackEnergy *= 1.0 - 0.5*G4UniformRand();
>> 330 dtaBlackTrackEnergy *= 1.0 - 0.5*G4UniformRand();
>> 331 }
>> 332 // G4cout << "EvaporationEffects "<<kineticEnergy<<" "
>> 333 // <<pnBlackTrackEnergy+dtaBlackTrackEnergy<<endl;
>> 334 return (pnBlackTrackEnergy+dtaBlackTrackEnergy)*GeV;
>> 335 }
401 336
>> 337 G4double G4Nucleus::AnnihilationEvaporationEffects(G4double kineticEnergy, G4double ekOrg)
>> 338 {
>> 339 // Nuclear evaporation as a function of atomic number and kinetic
>> 340 // energy (MeV) of primary particle. Modified for annihilation effects.
>> 341 //
>> 342 if( aEff < 1.5 || ekOrg < 0.)
>> 343 {
>> 344 pnBlackTrackEnergyfromAnnihilation = 0.0;
>> 345 dtaBlackTrackEnergyfromAnnihilation = 0.0;
>> 346 return 0.0;
>> 347 }
>> 348 G4double ek = kineticEnergy/GeV;
>> 349 G4float ekin = std::min( 4.0, std::max( 0.1, ek ) );
>> 350 const G4float atno = std::min( 120., aEff );
>> 351 const G4float gfa = 2.0*((aEff-1.0)/70.)*G4Exp(-(aEff-1.0)/70.);
>> 352
>> 353 G4float cfa = std::max( 0.15, 0.35 + ((0.35-0.05)/2.3)*G4Log(ekin) );
>> 354 G4float exnu = 7.716 * cfa * G4Exp(-cfa)
>> 355 * ((atno-1.0)/120.)*G4Exp(-(atno-1.0)/120.);
>> 356 G4float fpdiv = std::max( 0.5, 1.0-0.25*ekin*ekin );
402 357
403 G4double << 358 pnBlackTrackEnergyfromAnnihilation = exnu*fpdiv;
404 G4Nucleus::EvaporationEffects( G4double kineti << 359 dtaBlackTrackEnergyfromAnnihilation = exnu*(1.0-fpdiv);
405 { << 360
406 // derived from original FORTRAN code EXNU b <<
407 // <<
408 // Nuclear evaporation as function of atomic <<
409 // and kinetic energy (MeV) of primary parti <<
410 // <<
411 // returns kinetic energy (MeV) <<
412 // <<
413 if( aEff < 1.5 ) <<
414 { <<
415 pnBlackTrackEnergy = dtaBlackTrackEnergy = <<
416 return 0.0; <<
417 } <<
418 G4double ek = kineticEnergy/GeV; <<
419 G4float ekin = std::min( 4.0, std::max( 0.1, <<
420 const G4float atno = std::min( 120., aEff ); <<
421 const G4float gfa = 2.0*((aEff-1.0)/70.)*G4E <<
422 // <<
423 // 0.35 value at 1 GeV <<
424 // 0.05 value at 0.1 GeV <<
425 // <<
426 G4float cfa = std::max( 0.15, 0.35 + ((0.35- <<
427 G4float exnu = 7.716 * cfa * G4Exp(-cfa) <<
428 * ((atno-1.0)/120.)*G4Exp(-(atno-1.0)/120. <<
429 G4float fpdiv = std::max( 0.5, 1.0-0.25*ekin <<
430 // <<
431 // pnBlackTrackEnergy is the kinetic energy <<
432 // proton/neutron black <<
433 // dtaBlackTrackEnergy is the kinetic energy <<
434 // deuteron/triton/alpha <<
435 // <<
436 pnBlackTrackEnergy = exnu*fpdiv; <<
437 dtaBlackTrackEnergy = exnu*(1.0-fpdiv); <<
438 <<
439 if( G4int(zEff+0.1) != 82 ) <<
440 { <<
441 G4double ran1 = -6.0; 361 G4double ran1 = -6.0;
442 G4double ran2 = -6.0; 362 G4double ran2 = -6.0;
443 for( G4int i=0; i<12; ++i ) << 363 for( G4int i=0; i<12; ++i ) {
444 { <<
445 ran1 += G4UniformRand(); 364 ran1 += G4UniformRand();
446 ran2 += G4UniformRand(); 365 ran2 += G4UniformRand();
447 } 366 }
448 pnBlackTrackEnergy *= 1.0 + ran1*gfa; << 367 pnBlackTrackEnergyfromAnnihilation *= 1.0 + ran1*gfa;
449 dtaBlackTrackEnergy *= 1.0 + ran2*gfa; << 368 dtaBlackTrackEnergyfromAnnihilation *= 1.0 + ran2*gfa;
450 } <<
451 pnBlackTrackEnergy = std::max( 0.0, pnBlackT <<
452 dtaBlackTrackEnergy = std::max( 0.0, dtaBlac <<
453 while( pnBlackTrackEnergy+dtaBlackTrackEnerg <<
454 { <<
455 pnBlackTrackEnergy *= 1.0 - 0.5*G4UniformR <<
456 dtaBlackTrackEnergy *= 1.0 - 0.5*G4Uniform <<
457 } <<
458 //G4cout << "EvaporationEffects "<<kineticEn <<
459 // <<pnBlackTrackEnergy+dtaBlackTrackE <<
460 return (pnBlackTrackEnergy+dtaBlackTrackEner <<
461 } <<
462 <<
463 <<
464 G4double <<
465 G4Nucleus::AnnihilationEvaporationEffects(G4do <<
466 { <<
467 // Nuclear evaporation as a function of atom <<
468 // energy (MeV) of primary particle. Modifi <<
469 // <<
470 if( aEff < 1.5 || ekOrg < 0.) <<
471 { <<
472 pnBlackTrackEnergyfromAnnihilation = 0.0; <<
473 dtaBlackTrackEnergyfromAnnihilation = 0.0; <<
474 return 0.0; <<
475 } <<
476 G4double ek = kineticEnergy/GeV; <<
477 G4float ekin = std::min( 4.0, std::max( 0.1, <<
478 const G4float atno = std::min( 120., aEff ); <<
479 const G4float gfa = 2.0*((aEff-1.0)/70.)*G4E <<
480 <<
481 G4float cfa = std::max( 0.15, 0.35 + ((0.35- <<
482 G4float exnu = 7.716 * cfa * G4Exp(-cfa) <<
483 * ((atno-1.0)/120.)*G4Exp(-(atno-1.0)/120. <<
484 G4float fpdiv = std::max( 0.5, 1.0-0.25*ekin <<
485 <<
486 pnBlackTrackEnergyfromAnnihilation = exnu*fp <<
487 dtaBlackTrackEnergyfromAnnihilation = exnu*( <<
488 <<
489 G4double ran1 = -6.0; <<
490 G4double ran2 = -6.0; <<
491 for( G4int i=0; i<12; ++i ) { <<
492 ran1 += G4UniformRand(); <<
493 ran2 += G4UniformRand(); <<
494 } <<
495 pnBlackTrackEnergyfromAnnihilation *= 1.0 + <<
496 dtaBlackTrackEnergyfromAnnihilation *= 1.0 + <<
497 <<
498 pnBlackTrackEnergyfromAnnihilation = std::ma <<
499 dtaBlackTrackEnergyfromAnnihilation = std::m <<
500 G4double blackSum = pnBlackTrackEnergyfromAn <<
501 if (blackSum >= ekOrg/GeV) { <<
502 pnBlackTrackEnergyfromAnnihilation *= ekOr <<
503 dtaBlackTrackEnergyfromAnnihilation *= ekO <<
504 } <<
505 369
506 return (pnBlackTrackEnergyfromAnnihilation+d << 370 pnBlackTrackEnergyfromAnnihilation = std::max( 0.0, pnBlackTrackEnergyfromAnnihilation);
507 } << 371 dtaBlackTrackEnergyfromAnnihilation = std::max( 0.0, dtaBlackTrackEnergyfromAnnihilation);
>> 372 G4double blackSum = pnBlackTrackEnergyfromAnnihilation+dtaBlackTrackEnergyfromAnnihilation;
>> 373 if (blackSum >= ekOrg/GeV) {
>> 374 pnBlackTrackEnergyfromAnnihilation *= ekOrg/GeV/blackSum;
>> 375 dtaBlackTrackEnergyfromAnnihilation *= ekOrg/GeV/blackSum;
>> 376 }
508 377
>> 378 return (pnBlackTrackEnergyfromAnnihilation+dtaBlackTrackEnergyfromAnnihilation)*GeV;
>> 379 }
509 380
510 G4double << 381 G4double
511 G4Nucleus::Cinema( G4double kineticEnergy ) << 382 G4Nucleus::Cinema( G4double kineticEnergy )
512 { << 383 {
513 // derived from original FORTRAN code CINEMA << 384 // derived from original FORTRAN code CINEMA by H. Fesefeldt (14-Oct-1987)
514 // << 385 //
515 // input: kineticEnergy (MeV) << 386 // input: kineticEnergy (MeV)
516 // returns modified kinetic energy (MeV) << 387 // returns modified kinetic energy (MeV)
517 // << 388 //
518 static const G4double expxu = 82.; << 389 static const G4double expxu = 82.; // upper bound for arg. of exp
519 static const G4double expxl = -expxu; << 390 static const G4double expxl = -expxu; // lower bound for arg. of exp
520 << 391
521 G4double ek = kineticEnergy/GeV; << 392 G4double ek = kineticEnergy/GeV;
522 G4double ekLog = G4Log( ek ); << 393 G4double ekLog = G4Log( ek );
523 G4double aLog = G4Log( aEff ); << 394 G4double aLog = G4Log( aEff );
524 G4double em = std::min( 1.0, 0.2390 + 0.0408 << 395 G4double em = std::min( 1.0, 0.2390 + 0.0408*aLog*aLog );
525 G4double temp1 = -ek * std::min( 0.15, 0.001 << 396 G4double temp1 = -ek * std::min( 0.15, 0.0019*aLog*aLog*aLog );
526 G4double temp2 = G4Exp( std::max( expxl, std << 397 G4double temp2 = G4Exp( std::max( expxl, std::min( expxu, -(ekLog-em)*(ekLog-em)*2.0 ) ) );
527 G4double result = 0.0; << 398 G4double result = 0.0;
528 if( std::abs( temp1 ) < 1.0 ) << 399 if( std::abs( temp1 ) < 1.0 )
529 { << 400 {
530 if( temp2 > 1.0e-10 )result = temp1*temp2; << 401 if( temp2 > 1.0e-10 )result = temp1*temp2;
531 } << 402 }
532 else result = temp1*temp2; << 403 else result = temp1*temp2;
533 if( result < -ek )result = -ek; << 404 if( result < -ek )result = -ek;
534 return result*GeV; << 405 return result*GeV;
535 } << 406 }
536 407
>> 408 //
>> 409 // methods for class G4Nucleus ... by Christian Volcker
>> 410 //
537 411
538 G4ThreeVector G4Nucleus::GetFermiMomentum() << 412 G4ThreeVector G4Nucleus::GetFermiMomentum()
539 { << 413 {
540 // chv: .. we assume zero temperature! << 414 // chv: .. we assume zero temperature!
541 << 415
542 // momentum is equally distributed in each p << 416 // momentum is equally distributed in each phasespace volume dpx, dpy, dpz.
543 G4double ranflat1= << 417 G4double ranflat1=
544 G4RandFlat::shoot((G4double)0.,(G4double)f << 418 G4RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);
545 G4double ranflat2= << 419 G4double ranflat2=
546 G4RandFlat::shoot((G4double)0.,(G4double)f << 420 G4RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);
547 G4double ranflat3= << 421 G4double ranflat3=
548 G4RandFlat::shoot((G4double)0.,(G4double)f << 422 G4RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);
549 G4double ranmax = (ranflat1>ranflat2? ranfla << 423 G4double ranmax = (ranflat1>ranflat2? ranflat1: ranflat2);
550 ranmax = (ranmax>ranflat3? ranmax : ranflat3 << 424 ranmax = (ranmax>ranflat3? ranmax : ranflat3);
551 << 425
552 // Isotropic momentum distribution << 426 // Isotropic momentum distribution
553 G4double costheta = 2.*G4UniformRand() - 1.0 << 427 G4double costheta = 2.*G4UniformRand() - 1.0;
554 G4double sintheta = std::sqrt(1.0 - costheta << 428 G4double sintheta = std::sqrt(1.0 - costheta*costheta);
555 G4double phi = 2.0*pi*G4UniformRand(); << 429 G4double phi = 2.0*pi*G4UniformRand();
556 << 430
557 G4double pz=costheta*ranmax; << 431 G4double pz=costheta*ranmax;
558 G4double px=sintheta*std::cos(phi)*ranmax; << 432 G4double px=sintheta*std::cos(phi)*ranmax;
559 G4double py=sintheta*std::sin(phi)*ranmax; << 433 G4double py=sintheta*std::sin(phi)*ranmax;
560 G4ThreeVector p(px,py,pz); << 434 G4ThreeVector p(px,py,pz);
561 return p; << 435 return p;
562 } << 436 }
563 437
564 << 438 G4ReactionProductVector* G4Nucleus::Fragmentate()
565 G4ReactionProductVector* G4Nucleus::Fragmentat << 439 {
566 { << 440 // needs implementation!
567 // needs implementation! << 441 return NULL;
568 return nullptr; << 442 }
569 } <<
570 443
571 << 444 void G4Nucleus::AddMomentum(const G4ThreeVector aMomentum)
572 void G4Nucleus::AddMomentum(const G4ThreeVecto << 445 {
573 { << 446 momentum+=(aMomentum);
574 momentum+=(aMomentum); << 447 }
575 } <<
576 <<
577 448
578 void G4Nucleus::AddExcitationEnergy( G4double << 449 void G4Nucleus::AddExcitationEnergy( G4double anEnergy )
579 { << 450 {
580 excitationEnergy+=anEnergy; << 451 excitationEnergy+=anEnergy;
581 } << 452 }
582 453
583 /* end of file */ 454 /* end of file */
584 455
585 456