| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 //
28 // original by H.P. Wellisch
29 // modified by J.L. Chuma, TRIUMF, 19-Nov-1996
30 // last modified: 27-Mar-1997
31 // J.P.Wellisch: 23-Apr-97: minor simplifications
32 // modified by J.L.Chuma 24-Jul-97 to set the total momentum in Cinema and
33 // EvaporationEffects
34 // modified by J.L.Chuma 21-Oct-97 put std::abs() around the totalE^2-mass^2
35 // in calculation of total momentum in
36 // Cinema and EvaporationEffects
37 // Chr. Volcker, 10-Nov-1997: new methods and class variables.
38 // HPW added utilities for low energy neutron transport. (12.04.1998)
39 // M.G. Pia, 2 Oct 1998: modified GetFermiMomentum to avoid memory leaks
40 // G.Folger, spring 2010: add integer A/Z interface
41 // A. Ribon, summer 2015: migrated to G4Exp and G4Log
42 // A. Ribon, autumn 2021: extended to hypernuclei
43
44 #include "G4Nucleus.hh"
45 #include "G4NucleiProperties.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "Randomize.hh"
49 #include "G4HadronicException.hh"
50 #include "G4Exp.hh"
51 #include "G4Log.hh"
52 #include "G4HyperNucleiProperties.hh"
53 #include "G4HadronicParameters.hh"
54
55
56 G4Nucleus::G4Nucleus()
57 : theA(0), theZ(0), theL(0), aEff(0.0), zEff(0)
58 {
59 pnBlackTrackEnergy = 0.0;
60 dtaBlackTrackEnergy = 0.0;
61 pnBlackTrackEnergyfromAnnihilation = 0.0;
62 dtaBlackTrackEnergyfromAnnihilation = 0.0;
63 excitationEnergy = 0.0;
64 momentum = G4ThreeVector(0.,0.,0.);
65 fermiMomentum = 1.52*hbarc/fermi;
66 theTemp = 293.16*kelvin;
67 fIsotope = 0;
68 }
69
70 G4Nucleus::G4Nucleus( const G4double A, const G4double Z, const G4int numberOfLambdas )
71 {
72 SetParameters( A, Z, std::max(numberOfLambdas, 0) );
73 pnBlackTrackEnergy = 0.0;
74 dtaBlackTrackEnergy = 0.0;
75 pnBlackTrackEnergyfromAnnihilation = 0.0;
76 dtaBlackTrackEnergyfromAnnihilation = 0.0;
77 excitationEnergy = 0.0;
78 momentum = G4ThreeVector(0.,0.,0.);
79 fermiMomentum = 1.52*hbarc/fermi;
80 theTemp = 293.16*kelvin;
81 fIsotope = 0;
82 }
83
84 G4Nucleus::G4Nucleus( const G4int A, const G4int Z, const G4int numberOfLambdas )
85 {
86 SetParameters( A, Z, std::max(numberOfLambdas, 0) );
87 pnBlackTrackEnergy = 0.0;
88 dtaBlackTrackEnergy = 0.0;
89 pnBlackTrackEnergyfromAnnihilation = 0.0;
90 dtaBlackTrackEnergyfromAnnihilation = 0.0;
91 excitationEnergy = 0.0;
92 momentum = G4ThreeVector(0.,0.,0.);
93 fermiMomentum = 1.52*hbarc/fermi;
94 theTemp = 293.16*kelvin;
95 fIsotope = 0;
96 }
97
98 G4Nucleus::G4Nucleus( const G4Material *aMaterial )
99 {
100 ChooseParameters( aMaterial );
101 pnBlackTrackEnergy = 0.0;
102 dtaBlackTrackEnergy = 0.0;
103 pnBlackTrackEnergyfromAnnihilation = 0.0;
104 dtaBlackTrackEnergyfromAnnihilation = 0.0;
105 excitationEnergy = 0.0;
106 momentum = G4ThreeVector(0.,0.,0.);
107 fermiMomentum = 1.52*hbarc/fermi;
108 theTemp = aMaterial->GetTemperature();
109 fIsotope = 0;
110 }
111
112 G4Nucleus::~G4Nucleus() {}
113
114
115 //-------------------------------------------------------------------------------------------------
116 // SVT (Sampling of the Velocity of the Target nucleus) method, L. Thulliez (CEA-Saclay) 2021/05/04
117 //-------------------------------------------------------------------------------------------------
118 G4ReactionProduct
119 G4Nucleus::GetBiasedThermalNucleus(G4double aMass, G4ThreeVector aVelocity, G4double temp) const
120 {
121 // If E_neutron <= E_threshold, Then apply the Sampling ot the Velocity of the Target (SVT) method;
122 // Else consider the target nucleus being without motion.
123 G4double E_threshold = G4HadronicParameters::Instance()->GetNeutronKineticEnergyThresholdForSVT();
124 if ( E_threshold == -1. ) {
125 E_threshold = 400.0*8.617333262E-11*temp;
126 }
127 G4double E_neutron = 0.5*aVelocity.mag2()*G4Neutron::Neutron()->GetPDGMass(); // E=0.5*m*v2
128
129 G4ReactionProduct result;
130 result.SetMass(aMass*G4Neutron::Neutron()->GetPDGMass());
131
132 if ( E_neutron <= E_threshold ) {
133
134 // Beta = sqrt(m/2kT)
135 G4double beta = std::sqrt(result.GetMass()/(2.*8.617333262E-11*temp)); // kT E-5[eV] mass E-11[MeV] => beta in [m/s]-1
136
137 // Neutron speed vn
138 G4double vN_norm = aVelocity.mag();
139 G4double vN_norm2 = vN_norm*vN_norm;
140 G4double y = beta*vN_norm;
141
142 // Normalize neutron velocity
143 aVelocity = (1./vN_norm)*aVelocity;
144
145 // Sample target speed
146 G4double x2;
147 G4double randThreshold;
148 G4double vT_norm, vT_norm2, mu; //theta, val1, val2,
149 G4double acceptThreshold;
150 G4double vRelativeSpeed;
151 G4double cdf0 = 2./(2.+std::sqrt(CLHEP::pi)*y);
152
153 do {
154 // Sample the target velocity vT in the laboratory frame
155 if ( G4UniformRand() < cdf0 ) {
156 // Sample in C45 from https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-9721.pdf
157 x2 = -std::log(G4UniformRand()*G4UniformRand());
158 } else {
159 // Sample in C61 from https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-9721.pdf
160 G4double ampl = std::cos(CLHEP::pi/2.0 * G4UniformRand());
161 x2 = -std::log(G4UniformRand()) - std::log(G4UniformRand())*ampl*ampl;
162 }
163
164 vT_norm = std::sqrt(x2)/beta;
165 vT_norm2 = vT_norm*vT_norm;
166
167 // Sample cosine between the incident neutron and the target in the laboratory frame
168 mu = 2*G4UniformRand() - 1;
169
170 // Define acceptance threshold
171 vRelativeSpeed = std::sqrt(vN_norm2 + vT_norm2 - 2*vN_norm*vT_norm*mu);
172 acceptThreshold = vRelativeSpeed/(vN_norm + vT_norm);
173 randThreshold = G4UniformRand();
174 } while ( randThreshold >= acceptThreshold );
175
176 DoKinematicsOfThermalNucleus(mu, vT_norm, aVelocity, result);
177
178 } else { // target nucleus considered as being without motion
179
180 result.SetMomentum(0., 0., 0.);
181 result.SetKineticEnergy(0.);
182
183 }
184
185 return result;
186 }
187
188
189 void
190 G4Nucleus::DoKinematicsOfThermalNucleus(const G4double mu, const G4double vT_norm, const G4ThreeVector& aVelocity,
191 G4ReactionProduct& result) const {
192
193 // Get target nucleus direction from the neutron direction and the relative angle between target nucleus and neutron (mu)
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 neutron veloicty and the target velocity
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 - solve equation xx' + yy' + zz' = 0
205 G4ThreeVector ortho(1., 1., 1.);
206 if ( uNorm[0] ) ortho[0] = -(uNorm[1]+uNorm[2])/uNorm[0];
207 else if ( uNorm[1] ) ortho[1] = -(uNorm[0]+uNorm[2])/uNorm[1];
208 else if ( uNorm[2] ) ortho[2] = -(uNorm[0]+uNorm[1])/uNorm[2];
209
210 // Normalize the vector
211 ortho = (1/ortho.mag())*ortho;
212
213 // Find vector to draw a plan perpendicular to uNorm (i.e neutron velocity) with vectors ortho & orthoComp
214 G4ThreeVector orthoComp( uNorm[1]*ortho[2] - ortho[1]*uNorm[2],
215 uNorm[2]*ortho[0] - ortho[2]*uNorm[0],
216 uNorm[0]*ortho[1] - ortho[0]*uNorm[1] );
217
218 // Find the direction of the target velocity in the laboratory frame
219 G4ThreeVector directionTarget( cosTh*uNorm[0] + sinTh*(cosPhi*orthoComp[0] + sinPhi*ortho[0]),
220 cosTh*uNorm[1] + sinTh*(cosPhi*orthoComp[1] + sinPhi*ortho[1]),
221 cosTh*uNorm[2] + sinTh*(cosPhi*orthoComp[2] + sinPhi*ortho[2]) );
222
223 // Normalize directionTarget
224 directionTarget = ( 1./directionTarget.mag() )*directionTarget;
225
226 // Set momentum
227 G4double px = result.GetMass()*vT_norm*directionTarget[0];
228 G4double py = result.GetMass()*vT_norm*directionTarget[1];
229 G4double pz = result.GetMass()*vT_norm*directionTarget[2];
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.GetMass())*(tMom+result.GetMass())
234 - 2.*tMom*result.GetMass() );
235
236 if ( tEtot/result.GetMass() - 1. > 0.001 ) {
237 // use relativistic energy for higher energies
238 result.SetTotalEnergy(tEtot);
239 } else {
240 // use p**2/2M for lower energies (to preserve precision?)
241 result.SetKineticEnergy(tMom*tMom/(2.*result.GetMass()));
242 }
243
244 }
245
246
247 G4ReactionProduct
248 G4Nucleus::GetThermalNucleus(G4double targetMass, G4double temp) const
249 {
250 G4double currentTemp = temp;
251 if (currentTemp < 0) currentTemp = theTemp;
252 G4ReactionProduct theTarget;
253 theTarget.SetMass(targetMass*G4Neutron::Neutron()->GetPDGMass());
254 G4double px, py, pz;
255 px = GetThermalPz(theTarget.GetMass(), currentTemp);
256 py = GetThermalPz(theTarget.GetMass(), currentTemp);
257 pz = GetThermalPz(theTarget.GetMass(), currentTemp);
258 theTarget.SetMomentum(px, py, pz);
259 G4double tMom = std::sqrt(px*px+py*py+pz*pz);
260 G4double tEtot = std::sqrt((tMom+theTarget.GetMass())*
261 (tMom+theTarget.GetMass())-
262 2.*tMom*theTarget.GetMass());
263 // if(1-tEtot/theTarget.GetMass()>0.001) this line incorrect (Bug report 1911)
264 if (tEtot/theTarget.GetMass() - 1. > 0.001) {
265 // use relativistic energy for higher energies
266 theTarget.SetTotalEnergy(tEtot);
267
268 } else {
269 // use p**2/2M for lower energies (to preserve precision?)
270 theTarget.SetKineticEnergy(tMom*tMom/(2.*theTarget.GetMass()));
271 }
272 return theTarget;
273 }
274
275
276 void
277 G4Nucleus::ChooseParameters(const G4Material* aMaterial)
278 {
279 G4double random = G4UniformRand();
280 G4double sum = aMaterial->GetTotNbOfAtomsPerVolume();
281 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
282 G4double running(0);
283 // G4Element* element(0);
284 const G4Element* element = (*theElementVector)[aMaterial->GetNumberOfElements()-1];
285
286 for (unsigned int i = 0; i < aMaterial->GetNumberOfElements(); ++i) {
287 running += aMaterial->GetVecNbOfAtomsPerVolume()[i];
288 if (running > random*sum) {
289 element = (*theElementVector)[i];
290 break;
291 }
292 }
293
294 if (element->GetNumberOfIsotopes() > 0) {
295 G4double randomAbundance = G4UniformRand();
296 G4double sumAbundance = element->GetRelativeAbundanceVector()[0];
297 unsigned int iso=0;
298 while (iso < element->GetNumberOfIsotopes() && /* Loop checking, 02.11.2015, A.Ribon */
299 sumAbundance < randomAbundance) {
300 ++iso;
301 sumAbundance += element->GetRelativeAbundanceVector()[iso];
302 }
303 theA=element->GetIsotope(iso)->GetN();
304 theZ=element->GetIsotope(iso)->GetZ();
305 theL=0;
306 aEff=theA;
307 zEff=theZ;
308 } else {
309 aEff = element->GetN();
310 zEff = element->GetZ();
311 theZ = G4int(zEff + 0.5);
312 theA = G4int(aEff + 0.5);
313 theL=0;
314 }
315 }
316
317
318 void
319 G4Nucleus::SetParameters( const G4double A, const G4double Z, const G4int numberOfLambdas )
320 {
321 theZ = G4lrint(Z);
322 theA = G4lrint(A);
323 theL = std::max(numberOfLambdas, 0);
324 if (theA<1 || theZ<0 || theZ>theA) {
325 throw G4HadronicException(__FILE__, __LINE__,
326 "G4Nucleus::SetParameters called with non-physical parameters");
327 }
328 aEff = A; // atomic weight
329 zEff = Z; // atomic number
330 fIsotope = 0;
331 }
332
333
334 void
335 G4Nucleus::SetParameters( const G4int A, const G4int Z, const G4int numberOfLambdas )
336 {
337 theZ = Z;
338 theA = A;
339 theL = std::max(numberOfLambdas, 0);
340 if( theA<1 || theZ<0 || theZ>theA )
341 {
342 throw G4HadronicException(__FILE__, __LINE__,
343 "G4Nucleus::SetParameters called with non-physical parameters");
344 }
345 aEff = A; // atomic weight
346 zEff = Z; // atomic number
347 fIsotope = 0;
348 }
349
350
351 G4DynamicParticle *
352 G4Nucleus::ReturnTargetParticle() const
353 {
354 // choose a proton or a neutron (or a lamba if a hypernucleus) as the target particle
355 G4DynamicParticle *targetParticle = new G4DynamicParticle;
356 const G4double rnd = G4UniformRand();
357 if ( rnd < zEff/aEff ) {
358 targetParticle->SetDefinition( G4Proton::Proton() );
359 } else if ( rnd < (zEff + theL*1.0)/aEff ) {
360 targetParticle->SetDefinition( G4Lambda::Lambda() );
361 } else {
362 targetParticle->SetDefinition( G4Neutron::Neutron() );
363 }
364 return targetParticle;
365 }
366
367
368 G4double
369 G4Nucleus::AtomicMass( const G4double A, const G4double Z, const G4int numberOfLambdas ) const
370 {
371 // Now returns (atomic mass - electron masses)
372 if ( numberOfLambdas > 0 ) {
373 return G4HyperNucleiProperties::GetNuclearMass(G4int(A), G4int(Z), numberOfLambdas);
374 } else {
375 return G4NucleiProperties::GetNuclearMass(A, Z);
376 }
377 }
378
379
380 G4double
381 G4Nucleus::AtomicMass( const G4int A, const G4int Z, const G4int numberOfLambdas ) const
382 {
383 // Now returns (atomic mass - electron masses)
384 if ( numberOfLambdas > 0 ) {
385 return G4HyperNucleiProperties::GetNuclearMass(A, Z, numberOfLambdas);
386 } else {
387 return G4NucleiProperties::GetNuclearMass(A, Z);
388 }
389 }
390
391
392 G4double
393 G4Nucleus::GetThermalPz( const G4double mass, const G4double temp ) const
394 {
395 G4double result = G4RandGauss::shoot();
396 result *= std::sqrt(k_Boltzmann*temp*mass); // Das ist impuls (Pz),
397 // nichtrelativistische rechnung
398 // Maxwell verteilung angenommen
399 return result;
400 }
401
402
403 G4double
404 G4Nucleus::EvaporationEffects( G4double kineticEnergy )
405 {
406 // derived from original FORTRAN code EXNU by H. Fesefeldt (10-Dec-1986)
407 //
408 // Nuclear evaporation as function of atomic number
409 // and kinetic energy (MeV) of primary particle
410 //
411 // returns kinetic energy (MeV)
412 //
413 if( aEff < 1.5 )
414 {
415 pnBlackTrackEnergy = dtaBlackTrackEnergy = 0.0;
416 return 0.0;
417 }
418 G4double ek = kineticEnergy/GeV;
419 G4float ekin = std::min( 4.0, std::max( 0.1, ek ) );
420 const G4float atno = std::min( 120., aEff );
421 const G4float gfa = 2.0*((aEff-1.0)/70.)*G4Exp(-(aEff-1.0)/70.);
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-0.05)/2.3)*G4Log(ekin) );
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*ekin );
430 //
431 // pnBlackTrackEnergy is the kinetic energy (in GeV) available for
432 // proton/neutron black track particles
433 // dtaBlackTrackEnergy is the kinetic energy (in GeV) available for
434 // deuteron/triton/alpha black track particles
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;
442 G4double ran2 = -6.0;
443 for( G4int i=0; i<12; ++i )
444 {
445 ran1 += G4UniformRand();
446 ran2 += G4UniformRand();
447 }
448 pnBlackTrackEnergy *= 1.0 + ran1*gfa;
449 dtaBlackTrackEnergy *= 1.0 + ran2*gfa;
450 }
451 pnBlackTrackEnergy = std::max( 0.0, pnBlackTrackEnergy );
452 dtaBlackTrackEnergy = std::max( 0.0, dtaBlackTrackEnergy );
453 while( pnBlackTrackEnergy+dtaBlackTrackEnergy >= ek ) /* Loop checking, 02.11.2015, A.Ribon */
454 {
455 pnBlackTrackEnergy *= 1.0 - 0.5*G4UniformRand();
456 dtaBlackTrackEnergy *= 1.0 - 0.5*G4UniformRand();
457 }
458 //G4cout << "EvaporationEffects "<<kineticEnergy<<" "
459 // <<pnBlackTrackEnergy+dtaBlackTrackEnergy<< G4endl;
460 return (pnBlackTrackEnergy+dtaBlackTrackEnergy)*GeV;
461 }
462
463
464 G4double
465 G4Nucleus::AnnihilationEvaporationEffects(G4double kineticEnergy, G4double ekOrg)
466 {
467 // Nuclear evaporation as a function of atomic number and kinetic
468 // energy (MeV) of primary particle. Modified for annihilation effects.
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, ek ) );
478 const G4float atno = std::min( 120., aEff );
479 const G4float gfa = 2.0*((aEff-1.0)/70.)*G4Exp(-(aEff-1.0)/70.);
480
481 G4float cfa = std::max( 0.15, 0.35 + ((0.35-0.05)/2.3)*G4Log(ekin) );
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*ekin );
485
486 pnBlackTrackEnergyfromAnnihilation = exnu*fpdiv;
487 dtaBlackTrackEnergyfromAnnihilation = exnu*(1.0-fpdiv);
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 + ran1*gfa;
496 dtaBlackTrackEnergyfromAnnihilation *= 1.0 + ran2*gfa;
497
498 pnBlackTrackEnergyfromAnnihilation = std::max( 0.0, pnBlackTrackEnergyfromAnnihilation);
499 dtaBlackTrackEnergyfromAnnihilation = std::max( 0.0, dtaBlackTrackEnergyfromAnnihilation);
500 G4double blackSum = pnBlackTrackEnergyfromAnnihilation+dtaBlackTrackEnergyfromAnnihilation;
501 if (blackSum >= ekOrg/GeV) {
502 pnBlackTrackEnergyfromAnnihilation *= ekOrg/GeV/blackSum;
503 dtaBlackTrackEnergyfromAnnihilation *= ekOrg/GeV/blackSum;
504 }
505
506 return (pnBlackTrackEnergyfromAnnihilation+dtaBlackTrackEnergyfromAnnihilation)*GeV;
507 }
508
509
510 G4double
511 G4Nucleus::Cinema( G4double kineticEnergy )
512 {
513 // derived from original FORTRAN code CINEMA by H. Fesefeldt (14-Oct-1987)
514 //
515 // input: kineticEnergy (MeV)
516 // returns modified kinetic energy (MeV)
517 //
518 static const G4double expxu = 82.; // upper bound for arg. of exp
519 static const G4double expxl = -expxu; // lower bound for arg. of exp
520
521 G4double ek = kineticEnergy/GeV;
522 G4double ekLog = G4Log( ek );
523 G4double aLog = G4Log( aEff );
524 G4double em = std::min( 1.0, 0.2390 + 0.0408*aLog*aLog );
525 G4double temp1 = -ek * std::min( 0.15, 0.0019*aLog*aLog*aLog );
526 G4double temp2 = G4Exp( std::max( expxl, std::min( expxu, -(ekLog-em)*(ekLog-em)*2.0 ) ) );
527 G4double result = 0.0;
528 if( std::abs( temp1 ) < 1.0 )
529 {
530 if( temp2 > 1.0e-10 )result = temp1*temp2;
531 }
532 else result = temp1*temp2;
533 if( result < -ek )result = -ek;
534 return result*GeV;
535 }
536
537
538 G4ThreeVector G4Nucleus::GetFermiMomentum()
539 {
540 // chv: .. we assume zero temperature!
541
542 // momentum is equally distributed in each phasespace volume dpx, dpy, dpz.
543 G4double ranflat1=
544 G4RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);
545 G4double ranflat2=
546 G4RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);
547 G4double ranflat3=
548 G4RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);
549 G4double ranmax = (ranflat1>ranflat2? ranflat1: ranflat2);
550 ranmax = (ranmax>ranflat3? ranmax : ranflat3);
551
552 // Isotropic momentum distribution
553 G4double costheta = 2.*G4UniformRand() - 1.0;
554 G4double sintheta = std::sqrt(1.0 - costheta*costheta);
555 G4double phi = 2.0*pi*G4UniformRand();
556
557 G4double pz=costheta*ranmax;
558 G4double px=sintheta*std::cos(phi)*ranmax;
559 G4double py=sintheta*std::sin(phi)*ranmax;
560 G4ThreeVector p(px,py,pz);
561 return p;
562 }
563
564
565 G4ReactionProductVector* G4Nucleus::Fragmentate()
566 {
567 // needs implementation!
568 return nullptr;
569 }
570
571
572 void G4Nucleus::AddMomentum(const G4ThreeVector aMomentum)
573 {
574 momentum+=(aMomentum);
575 }
576
577
578 void G4Nucleus::AddExcitationEnergy( G4double anEnergy )
579 {
580 excitationEnergy+=anEnergy;
581 }
582
583 /* end of file */
584
585