Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Nucleus.cc
Go to the documentation of this file.
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 
42 #include "G4Nucleus.hh"
43 #include "G4NucleiProperties.hh"
44 #include "G4PhysicalConstants.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "Randomize.hh"
47 #include "G4HadronicException.hh"
48 
50  : theA(0), theZ(0), aEff(0.0), zEff(0)
51 {
52  pnBlackTrackEnergy = 0.0;
53  dtaBlackTrackEnergy = 0.0;
54  pnBlackTrackEnergyfromAnnihilation = 0.0;
55  dtaBlackTrackEnergyfromAnnihilation = 0.0;
56  excitationEnergy = 0.0;
57  momentum = G4ThreeVector(0.,0.,0.);
58  fermiMomentum = 1.52*hbarc/fermi;
59  theTemp = 293.16*kelvin;
60  fIsotope = 0;
61 }
62 
64 {
65  SetParameters( A, Z );
66  pnBlackTrackEnergy = 0.0;
67  dtaBlackTrackEnergy = 0.0;
68  pnBlackTrackEnergyfromAnnihilation = 0.0;
69  dtaBlackTrackEnergyfromAnnihilation = 0.0;
70  excitationEnergy = 0.0;
71  momentum = G4ThreeVector(0.,0.,0.);
72  fermiMomentum = 1.52*hbarc/fermi;
73  theTemp = 293.16*kelvin;
74  fIsotope = 0;
75 }
76 
77 G4Nucleus::G4Nucleus( const G4int A, const G4int Z )
78 {
79  SetParameters( A, Z );
80  pnBlackTrackEnergy = 0.0;
81  dtaBlackTrackEnergy = 0.0;
82  pnBlackTrackEnergyfromAnnihilation = 0.0;
83  dtaBlackTrackEnergyfromAnnihilation = 0.0;
84  excitationEnergy = 0.0;
85  momentum = G4ThreeVector(0.,0.,0.);
86  fermiMomentum = 1.52*hbarc/fermi;
87  theTemp = 293.16*kelvin;
88  fIsotope = 0;
89 }
90 
91 G4Nucleus::G4Nucleus( const G4Material *aMaterial )
92 {
93  ChooseParameters( aMaterial );
94  pnBlackTrackEnergy = 0.0;
95  dtaBlackTrackEnergy = 0.0;
96  pnBlackTrackEnergyfromAnnihilation = 0.0;
97  dtaBlackTrackEnergyfromAnnihilation = 0.0;
98  excitationEnergy = 0.0;
99  momentum = G4ThreeVector(0.,0.,0.);
100  fermiMomentum = 1.52*hbarc/fermi;
101  theTemp = aMaterial->GetTemperature();
102  fIsotope = 0;
103 }
104 
106 
109 {
110  G4double velMag = aVelocity.mag();
111  G4ReactionProduct result;
112  G4double value = 0;
113  G4double random = 1;
114  G4double norm = 3.*std::sqrt(k_Boltzmann*temp*aMass*G4Neutron::Neutron()->GetPDGMass());
115  norm /= G4Neutron::Neutron()->GetPDGMass();
116  norm *= 5.;
117  norm += velMag;
118  norm /= velMag;
119  while(value/norm<random)
120  {
121  result = GetThermalNucleus(aMass, temp);
122  G4ThreeVector targetVelocity = 1./result.GetMass()*result.GetMomentum();
123  value = (targetVelocity+aVelocity).mag()/velMag;
124  random = G4UniformRand();
125  }
126  return result;
127 }
128 
131 {
132  G4double currentTemp = temp;
133  if(currentTemp < 0) currentTemp = theTemp;
135  theTarget.SetMass(targetMass*G4Neutron::Neutron()->GetPDGMass());
136  G4double px, py, pz;
137  px = GetThermalPz(theTarget.GetMass(), currentTemp);
138  py = GetThermalPz(theTarget.GetMass(), currentTemp);
139  pz = GetThermalPz(theTarget.GetMass(), currentTemp);
140  theTarget.SetMomentum(px, py, pz);
141  G4double tMom = std::sqrt(px*px+py*py+pz*pz);
142  G4double tEtot = std::sqrt((tMom+theTarget.GetMass())*
143  (tMom+theTarget.GetMass())-
144  2.*tMom*theTarget.GetMass());
145  if(1-tEtot/theTarget.GetMass()>0.001)
146  {
147  theTarget.SetTotalEnergy(tEtot);
148  }
149  else
150  {
151  theTarget.SetKineticEnergy(tMom*tMom/(2.*theTarget.GetMass()));
152  }
153  return theTarget;
154 }
155 
156 
157 void
159 {
160  G4double random = G4UniformRand();
161  G4double sum = aMaterial->GetTotNbOfAtomsPerVolume();
162  const G4ElementVector* theElementVector = aMaterial->GetElementVector();
163  G4double running(0);
164  // G4Element* element(0);
165  G4Element* element = (*theElementVector)[aMaterial->GetNumberOfElements()-1];
166 
167  for (unsigned int i = 0; i < aMaterial->GetNumberOfElements(); ++i) {
168  running += aMaterial->GetVecNbOfAtomsPerVolume()[i];
169  if (running > random*sum) {
170  element = (*theElementVector)[i];
171  break;
172  }
173  }
174 
175  if (element->GetNumberOfIsotopes() > 0) {
176  G4double randomAbundance = G4UniformRand();
177  G4double sumAbundance = element->GetRelativeAbundanceVector()[0];
178  unsigned int iso=0;
179  while (iso < element->GetNumberOfIsotopes() &&
180  sumAbundance < randomAbundance) {
181  ++iso;
182  sumAbundance += element->GetRelativeAbundanceVector()[iso];
183  }
184  theA=element->GetIsotope(iso)->GetN();
185  theZ=element->GetIsotope(iso)->GetZ();
186  aEff=theA;
187  zEff=theZ;
188  } else {
189  aEff = element->GetN();
190  zEff = element->GetZ();
191  theZ = G4int(zEff + 0.5);
192  theA = G4int(aEff + 0.5);
193  }
194 }
195 
196 
197 void
199 {
200  theZ = G4lrint(Z);
201  theA = G4lrint(A);
202  if (theA<1 || theZ<0 || theZ>theA) {
203  throw G4HadronicException(__FILE__, __LINE__,
204  "G4Nucleus::SetParameters called with non-physical parameters");
205  }
206  aEff = A; // atomic weight
207  zEff = Z; // atomic number
208  fIsotope = 0;
209 }
210 
211 void
213 {
214  theZ = Z;
215  theA = A;
216  if( theA<1 || theZ<0 || theZ>theA )
217  {
218  throw G4HadronicException(__FILE__, __LINE__,
219  "G4Nucleus::SetParameters called with non-physical parameters");
220  }
221  aEff = A; // atomic weight
222  zEff = Z; // atomic number
223  fIsotope = 0;
224 }
225 
228  {
229  // choose a proton or a neutron as the target particle
230 
231  G4DynamicParticle *targetParticle = new G4DynamicParticle;
232  if( G4UniformRand() < zEff/aEff )
233  targetParticle->SetDefinition( G4Proton::Proton() );
234  else
235  targetParticle->SetDefinition( G4Neutron::Neutron() );
236  return targetParticle;
237  }
238 
239  G4double
240  G4Nucleus::AtomicMass( const G4double A, const G4double Z ) const
241  {
242  // Now returns (atomic mass - electron masses)
244  }
245 
246  G4double
247  G4Nucleus::AtomicMass( const G4int A, const G4int Z ) const
248  {
249  // Now returns (atomic mass - electron masses)
251  }
252 
253  G4double
254  G4Nucleus::GetThermalPz( const G4double mass, const G4double temp ) const
255  {
256  G4double result = G4RandGauss::shoot();
257  result *= std::sqrt(k_Boltzmann*temp*mass); // Das ist impuls (Pz),
258  // nichtrelativistische rechnung
259  // Maxwell verteilung angenommen
260  return result;
261  }
262 
263  G4double
265  {
266  // derived from original FORTRAN code EXNU by H. Fesefeldt (10-Dec-1986)
267  //
268  // Nuclear evaporation as function of atomic number
269  // and kinetic energy (MeV) of primary particle
270  //
271  // returns kinetic energy (MeV)
272  //
273  if( aEff < 1.5 )
274  {
275  pnBlackTrackEnergy = dtaBlackTrackEnergy = 0.0;
276  return 0.0;
277  }
278  G4double ek = kineticEnergy/GeV;
279  G4float ekin = std::min( 4.0, std::max( 0.1, ek ) );
280  const G4float atno = std::min( 120., aEff );
281  const G4float gfa = 2.0*((aEff-1.0)/70.)*std::exp(-(aEff-1.0)/70.);
282  //
283  // 0.35 value at 1 GeV
284  // 0.05 value at 0.1 GeV
285  //
286  G4float cfa = std::max( 0.15, 0.35 + ((0.35-0.05)/2.3)*std::log(ekin) );
287  G4float exnu = 7.716 * cfa * std::exp(-cfa)
288  * ((atno-1.0)/120.)*std::exp(-(atno-1.0)/120.);
289  G4float fpdiv = std::max( 0.5, 1.0-0.25*ekin*ekin );
290  //
291  // pnBlackTrackEnergy is the kinetic energy (in GeV) available for
292  // proton/neutron black track particles
293  // dtaBlackTrackEnergy is the kinetic energy (in GeV) available for
294  // deuteron/triton/alpha black track particles
295  //
296  pnBlackTrackEnergy = exnu*fpdiv;
297  dtaBlackTrackEnergy = exnu*(1.0-fpdiv);
298 
299  if( G4int(zEff+0.1) != 82 )
300  {
301  G4double ran1 = -6.0;
302  G4double ran2 = -6.0;
303  for( G4int i=0; i<12; ++i )
304  {
305  ran1 += G4UniformRand();
306  ran2 += G4UniformRand();
307  }
308  pnBlackTrackEnergy *= 1.0 + ran1*gfa;
309  dtaBlackTrackEnergy *= 1.0 + ran2*gfa;
310  }
311  pnBlackTrackEnergy = std::max( 0.0, pnBlackTrackEnergy );
312  dtaBlackTrackEnergy = std::max( 0.0, dtaBlackTrackEnergy );
313  while( pnBlackTrackEnergy+dtaBlackTrackEnergy >= ek )
314  {
315  pnBlackTrackEnergy *= 1.0 - 0.5*G4UniformRand();
316  dtaBlackTrackEnergy *= 1.0 - 0.5*G4UniformRand();
317  }
318 // G4cout << "EvaporationEffects "<<kineticEnergy<<" "
319 // <<pnBlackTrackEnergy+dtaBlackTrackEnergy<<endl;
320  return (pnBlackTrackEnergy+dtaBlackTrackEnergy)*GeV;
321  }
322 
324  {
325  // Nuclear evaporation as a function of atomic number and kinetic
326  // energy (MeV) of primary particle. Modified for annihilation effects.
327  //
328  if( aEff < 1.5 || ekOrg < 0.)
329  {
330  pnBlackTrackEnergyfromAnnihilation = 0.0;
331  dtaBlackTrackEnergyfromAnnihilation = 0.0;
332  return 0.0;
333  }
334  G4double ek = kineticEnergy/GeV;
335  G4float ekin = std::min( 4.0, std::max( 0.1, ek ) );
336  const G4float atno = std::min( 120., aEff );
337  const G4float gfa = 2.0*((aEff-1.0)/70.)*std::exp(-(aEff-1.0)/70.);
338 
339  G4float cfa = std::max( 0.15, 0.35 + ((0.35-0.05)/2.3)*std::log(ekin) );
340  G4float exnu = 7.716 * cfa * std::exp(-cfa)
341  * ((atno-1.0)/120.)*std::exp(-(atno-1.0)/120.);
342  G4float fpdiv = std::max( 0.5, 1.0-0.25*ekin*ekin );
343 
344  pnBlackTrackEnergyfromAnnihilation = exnu*fpdiv;
345  dtaBlackTrackEnergyfromAnnihilation = exnu*(1.0-fpdiv);
346 
347  G4double ran1 = -6.0;
348  G4double ran2 = -6.0;
349  for( G4int i=0; i<12; ++i ) {
350  ran1 += G4UniformRand();
351  ran2 += G4UniformRand();
352  }
353  pnBlackTrackEnergyfromAnnihilation *= 1.0 + ran1*gfa;
354  dtaBlackTrackEnergyfromAnnihilation *= 1.0 + ran2*gfa;
355 
356  pnBlackTrackEnergyfromAnnihilation = std::max( 0.0, pnBlackTrackEnergyfromAnnihilation);
357  dtaBlackTrackEnergyfromAnnihilation = std::max( 0.0, dtaBlackTrackEnergyfromAnnihilation);
358  G4double blackSum = pnBlackTrackEnergyfromAnnihilation+dtaBlackTrackEnergyfromAnnihilation;
359  if (blackSum >= ekOrg/GeV) {
360  pnBlackTrackEnergyfromAnnihilation *= ekOrg/GeV/blackSum;
361  dtaBlackTrackEnergyfromAnnihilation *= ekOrg/GeV/blackSum;
362  }
363 
364  return (pnBlackTrackEnergyfromAnnihilation+dtaBlackTrackEnergyfromAnnihilation)*GeV;
365  }
366 
367  G4double
368  G4Nucleus::Cinema( G4double kineticEnergy )
369  {
370  // derived from original FORTRAN code CINEMA by H. Fesefeldt (14-Oct-1987)
371  //
372  // input: kineticEnergy (MeV)
373  // returns modified kinetic energy (MeV)
374  //
375  static const G4double expxu = 82.; // upper bound for arg. of exp
376  static const G4double expxl = -expxu; // lower bound for arg. of exp
377 
378  G4double ek = kineticEnergy/GeV;
379  G4double ekLog = std::log( ek );
380  G4double aLog = std::log( aEff );
381  G4double em = std::min( 1.0, 0.2390 + 0.0408*aLog*aLog );
382  G4double temp1 = -ek * std::min( 0.15, 0.0019*aLog*aLog*aLog );
383  G4double temp2 = std::exp( std::max( expxl, std::min( expxu, -(ekLog-em)*(ekLog-em)*2.0 ) ) );
384  G4double result = 0.0;
385  if( std::abs( temp1 ) < 1.0 )
386  {
387  if( temp2 > 1.0e-10 )result = temp1*temp2;
388  }
389  else result = temp1*temp2;
390  if( result < -ek )result = -ek;
391  return result*GeV;
392  }
393 
394  //
395  // methods for class G4Nucleus ... by Christian Volcker
396  //
397 
399  {
400  // chv: .. we assume zero temperature!
401 
402  // momentum is equally distributed in each phasespace volume dpx, dpy, dpz.
403  G4double ranflat1=
404  G4RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);
405  G4double ranflat2=
406  G4RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);
407  G4double ranflat3=
408  G4RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);
409  G4double ranmax = (ranflat1>ranflat2? ranflat1: ranflat2);
410  ranmax = (ranmax>ranflat3? ranmax : ranflat3);
411 
412  // Isotropic momentum distribution
413  G4double costheta = 2.*G4UniformRand() - 1.0;
414  G4double sintheta = std::sqrt(1.0 - costheta*costheta);
415  G4double phi = 2.0*pi*G4UniformRand();
416 
417  G4double pz=costheta*ranmax;
418  G4double px=sintheta*std::cos(phi)*ranmax;
419  G4double py=sintheta*std::sin(phi)*ranmax;
420  G4ThreeVector p(px,py,pz);
421  return p;
422  }
423 
425  {
426  // needs implementation!
427  return NULL;
428  }
429 
430  void G4Nucleus::AddMomentum(const G4ThreeVector aMomentum)
431  {
432  momentum+=(aMomentum);
433  }
434 
436  {
437  excitationEnergy+=anEnergy;
438  }
439 
440  /* end of file */
441 
size_t GetNumberOfIsotopes() const
Definition: G4Element.hh:158
ThreeVector shoot(const G4int Ap, const G4int Af)
G4double AtomicMass(const G4double A, const G4double Z) const
Definition: G4Nucleus.cc:240
G4double EvaporationEffects(G4double kineticEnergy)
Definition: G4Nucleus.cc:264
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4double GetThermalPz(const G4double mass, const G4double temp) const
Definition: G4Nucleus.cc:254
std::vector< G4Element * > G4ElementVector
CLHEP::Hep3Vector G4ThreeVector
G4double GetN() const
Definition: G4Element.hh:134
void SetKineticEnergy(const G4double en)
void SetMomentum(const G4double x, const G4double y, const G4double z)
const char * p
Definition: xmltok.h:285
float G4float
Definition: G4Types.hh:77
G4double GetZ() const
Definition: G4Element.hh:131
void ChooseParameters(const G4Material *aMaterial)
Definition: G4Nucleus.cc:158
G4ReactionProduct GetThermalNucleus(G4double aMass, G4double temp=-1) const
Definition: G4Nucleus.cc:130
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:188
int G4int
Definition: G4Types.hh:78
G4DynamicParticle * ReturnTargetParticle() const
Definition: G4Nucleus.cc:227
G4ReactionProductVector * Fragmentate()
Definition: G4Nucleus.cc:424
std::vector< G4ReactionProduct * > G4ReactionProductVector
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:204
G4ReactionProduct GetBiasedThermalNucleus(G4double aMass, G4ThreeVector aVelocity, G4double temp=-1) const
Definition: G4Nucleus.cc:108
void SetMass(const G4double mas)
G4int GetN() const
Definition: G4Isotope.hh:94
#define G4UniformRand()
Definition: Randomize.hh:87
float k_Boltzmann
Definition: hepunit.py:299
void SetTotalEnergy(const G4double en)
G4ErrorTarget * theTarget
Definition: errprop.cc:59
void AddMomentum(const G4ThreeVector aMomentum)
Definition: G4Nucleus.cc:430
G4double * GetRelativeAbundanceVector() const
Definition: G4Element.hh:166
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4int GetZ() const
Definition: G4Isotope.hh:91
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:207
G4ThreeVector GetFermiMomentum()
Definition: G4Nucleus.cc:398
G4double AnnihilationEvaporationEffects(G4double kineticEnergy, G4double ekOrg)
Definition: G4Nucleus.cc:323
G4double GetPDGMass() const
int G4lrint(double ad)
Definition: templates.hh:163
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double Cinema(G4double kineticEnergy)
Definition: G4Nucleus.cc:368
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4ThreeVector GetMomentum() const
const G4Isotope * GetIsotope(G4int iso) const
Definition: G4Element.hh:169
G4double GetTemperature() const
Definition: G4Material.hh:180
const XML_Char int const XML_Char * value
void AddExcitationEnergy(G4double anEnergy)
Definition: G4Nucleus.cc:435
size_t GetNumberOfElements() const
Definition: G4Material.hh:184
double G4double
Definition: G4Types.hh:76
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
double mag() const
void SetParameters(const G4double A, const G4double Z)
Definition: G4Nucleus.cc:198
G4double GetMass() const