Geant4  10.03.p03
 All Classes 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  // A. Ribon, 6 August 2015: migrated to G4Exp and G4Log.
42 
43 #include "G4Nucleus.hh"
44 #include "G4NucleiProperties.hh"
45 #include "G4PhysicalConstants.hh"
46 #include "G4SystemOfUnits.hh"
47 #include "Randomize.hh"
48 #include "G4HadronicException.hh"
49 
50 #include "G4Exp.hh"
51 #include "G4Log.hh"
52 
53 
55  : theA(0), theZ(0), aEff(0.0), zEff(0)
56 {
57  pnBlackTrackEnergy = 0.0;
58  dtaBlackTrackEnergy = 0.0;
59  pnBlackTrackEnergyfromAnnihilation = 0.0;
60  dtaBlackTrackEnergyfromAnnihilation = 0.0;
61  excitationEnergy = 0.0;
62  momentum = G4ThreeVector(0.,0.,0.);
63  fermiMomentum = 1.52*hbarc/fermi;
64  theTemp = 293.16*kelvin;
65  fIsotope = 0;
66 }
67 
69 {
70  SetParameters( A, Z );
71  pnBlackTrackEnergy = 0.0;
72  dtaBlackTrackEnergy = 0.0;
73  pnBlackTrackEnergyfromAnnihilation = 0.0;
74  dtaBlackTrackEnergyfromAnnihilation = 0.0;
75  excitationEnergy = 0.0;
76  momentum = G4ThreeVector(0.,0.,0.);
77  fermiMomentum = 1.52*hbarc/fermi;
78  theTemp = 293.16*kelvin;
79  fIsotope = 0;
80 }
81 
83 {
84  SetParameters( A, Z );
85  pnBlackTrackEnergy = 0.0;
86  dtaBlackTrackEnergy = 0.0;
87  pnBlackTrackEnergyfromAnnihilation = 0.0;
88  dtaBlackTrackEnergyfromAnnihilation = 0.0;
89  excitationEnergy = 0.0;
90  momentum = G4ThreeVector(0.,0.,0.);
91  fermiMomentum = 1.52*hbarc/fermi;
92  theTemp = 293.16*kelvin;
93  fIsotope = 0;
94 }
95 
96 G4Nucleus::G4Nucleus( const G4Material *aMaterial )
97 {
98  ChooseParameters( aMaterial );
99  pnBlackTrackEnergy = 0.0;
100  dtaBlackTrackEnergy = 0.0;
101  pnBlackTrackEnergyfromAnnihilation = 0.0;
102  dtaBlackTrackEnergyfromAnnihilation = 0.0;
103  excitationEnergy = 0.0;
104  momentum = G4ThreeVector(0.,0.,0.);
105  fermiMomentum = 1.52*hbarc/fermi;
106  theTemp = aMaterial->GetTemperature();
107  fIsotope = 0;
108 }
109 
111 
114 {
115  G4double velMag = aVelocity.mag();
117  G4double value = 0;
118  G4double random = 1;
119  G4double norm = 3.*std::sqrt(k_Boltzmann*temp*aMass*G4Neutron::Neutron()->GetPDGMass());
120  norm /= G4Neutron::Neutron()->GetPDGMass();
121  norm *= 5.;
122  norm += velMag;
123  norm /= velMag;
124  const G4int maxNumberOfLoops = 1000000;
125  G4int loopCounter = -1;
126  while ( (value/norm<random) && ++loopCounter < maxNumberOfLoops ) /* Loop checking, 02.11.2015, A.Ribon */
127  {
128  result = GetThermalNucleus(aMass, temp);
129  G4ThreeVector targetVelocity = 1./result.GetMass()*result.GetMomentum();
130  value = (targetVelocity+aVelocity).mag()/velMag;
131  random = G4UniformRand();
132  }
133  if ( loopCounter >= maxNumberOfLoops ) {
135  ed << " Failed sampling after maxNumberOfLoops attempts : forced exit! " << G4endl;
136  G4Exception( " G4Nucleus::GetBiasedThermalNucleus ", "HAD_NUCLEUS_001", JustWarning, ed );
137  result = GetThermalNucleus(aMass, temp);
138  }
139  return result;
140 }
141 
144 {
145  G4double currentTemp = temp;
146  if (currentTemp < 0) currentTemp = theTemp;
148  theTarget.SetMass(targetMass*G4Neutron::Neutron()->GetPDGMass());
149  G4double px, py, pz;
150  px = GetThermalPz(theTarget.GetMass(), currentTemp);
151  py = GetThermalPz(theTarget.GetMass(), currentTemp);
152  pz = GetThermalPz(theTarget.GetMass(), currentTemp);
153  theTarget.SetMomentum(px, py, pz);
154  G4double tMom = std::sqrt(px*px+py*py+pz*pz);
155  G4double tEtot = std::sqrt((tMom+theTarget.GetMass())*
156  (tMom+theTarget.GetMass())-
157  2.*tMom*theTarget.GetMass());
158  // if(1-tEtot/theTarget.GetMass()>0.001) this line incorrect (Bug report 1911)
159  if (tEtot/theTarget.GetMass() - 1. > 0.001) {
160  // use relativistic energy for higher energies
161  theTarget.SetTotalEnergy(tEtot);
162 
163  } else {
164  // use p**2/2M for lower energies (to preserve precision?)
165  theTarget.SetKineticEnergy(tMom*tMom/(2.*theTarget.GetMass()));
166  }
167  return theTarget;
168 }
169 
170 
171 void
173 {
174  G4double random = G4UniformRand();
175  G4double sum = aMaterial->GetTotNbOfAtomsPerVolume();
176  const G4ElementVector* theElementVector = aMaterial->GetElementVector();
177  G4double running(0);
178  // G4Element* element(0);
179  G4Element* element = (*theElementVector)[aMaterial->GetNumberOfElements()-1];
180 
181  for (unsigned int i = 0; i < aMaterial->GetNumberOfElements(); ++i) {
182  running += aMaterial->GetVecNbOfAtomsPerVolume()[i];
183  if (running > random*sum) {
184  element = (*theElementVector)[i];
185  break;
186  }
187  }
188 
189  if (element->GetNumberOfIsotopes() > 0) {
190  G4double randomAbundance = G4UniformRand();
191  G4double sumAbundance = element->GetRelativeAbundanceVector()[0];
192  unsigned int iso=0;
193  while (iso < element->GetNumberOfIsotopes() && /* Loop checking, 02.11.2015, A.Ribon */
194  sumAbundance < randomAbundance) {
195  ++iso;
196  sumAbundance += element->GetRelativeAbundanceVector()[iso];
197  }
198  theA=element->GetIsotope(iso)->GetN();
199  theZ=element->GetIsotope(iso)->GetZ();
200  aEff=theA;
201  zEff=theZ;
202  } else {
203  aEff = element->GetN();
204  zEff = element->GetZ();
205  theZ = G4int(zEff + 0.5);
206  theA = G4int(aEff + 0.5);
207  }
208 }
209 
210 
211 void
213 {
214  theZ = G4lrint(Z);
215  theA = G4lrint(A);
216  if (theA<1 || theZ<0 || theZ>theA) {
217  throw G4HadronicException(__FILE__, __LINE__,
218  "G4Nucleus::SetParameters called with non-physical parameters");
219  }
220  aEff = A; // atomic weight
221  zEff = Z; // atomic number
222  fIsotope = 0;
223 }
224 
225 void
227 {
228  theZ = Z;
229  theA = A;
230  if( theA<1 || theZ<0 || theZ>theA )
231  {
232  throw G4HadronicException(__FILE__, __LINE__,
233  "G4Nucleus::SetParameters called with non-physical parameters");
234  }
235  aEff = A; // atomic weight
236  zEff = Z; // atomic number
237  fIsotope = 0;
238 }
239 
242  {
243  // choose a proton or a neutron as the target particle
244 
245  G4DynamicParticle *targetParticle = new G4DynamicParticle;
246  if( G4UniformRand() < zEff/aEff )
247  targetParticle->SetDefinition( G4Proton::Proton() );
248  else
249  targetParticle->SetDefinition( G4Neutron::Neutron() );
250  return targetParticle;
251  }
252 
253  G4double
254  G4Nucleus::AtomicMass( const G4double A, const G4double Z ) const
255  {
256  // Now returns (atomic mass - electron masses)
258  }
259 
260  G4double
261  G4Nucleus::AtomicMass( const G4int A, const G4int Z ) const
262  {
263  // Now returns (atomic mass - electron masses)
265  }
266 
267  G4double
268  G4Nucleus::GetThermalPz( const G4double mass, const G4double temp ) const
269  {
271  result *= std::sqrt(k_Boltzmann*temp*mass); // Das ist impuls (Pz),
272  // nichtrelativistische rechnung
273  // Maxwell verteilung angenommen
274  return result;
275  }
276 
277  G4double
279  {
280  // derived from original FORTRAN code EXNU by H. Fesefeldt (10-Dec-1986)
281  //
282  // Nuclear evaporation as function of atomic number
283  // and kinetic energy (MeV) of primary particle
284  //
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  }
336 
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 );
357 
358  pnBlackTrackEnergyfromAnnihilation = exnu*fpdiv;
359  dtaBlackTrackEnergyfromAnnihilation = exnu*(1.0-fpdiv);
360 
361  G4double ran1 = -6.0;
362  G4double ran2 = -6.0;
363  for( G4int i=0; i<12; ++i ) {
364  ran1 += G4UniformRand();
365  ran2 += G4UniformRand();
366  }
367  pnBlackTrackEnergyfromAnnihilation *= 1.0 + ran1*gfa;
368  dtaBlackTrackEnergyfromAnnihilation *= 1.0 + ran2*gfa;
369 
370  pnBlackTrackEnergyfromAnnihilation = std::max( 0.0, pnBlackTrackEnergyfromAnnihilation);
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  }
377 
378  return (pnBlackTrackEnergyfromAnnihilation+dtaBlackTrackEnergyfromAnnihilation)*GeV;
379  }
380 
381  G4double
382  G4Nucleus::Cinema( G4double kineticEnergy )
383  {
384  // derived from original FORTRAN code CINEMA by H. Fesefeldt (14-Oct-1987)
385  //
386  // input: kineticEnergy (MeV)
387  // returns modified kinetic energy (MeV)
388  //
389  static const G4double expxu = 82.; // upper bound for arg. of exp
390  static const G4double expxl = -expxu; // lower bound for arg. of exp
391 
392  G4double ek = kineticEnergy/GeV;
393  G4double ekLog = G4Log( ek );
394  G4double aLog = G4Log( aEff );
395  G4double em = std::min( 1.0, 0.2390 + 0.0408*aLog*aLog );
396  G4double temp1 = -ek * std::min( 0.15, 0.0019*aLog*aLog*aLog );
397  G4double temp2 = G4Exp( std::max( expxl, std::min( expxu, -(ekLog-em)*(ekLog-em)*2.0 ) ) );
398  G4double result = 0.0;
399  if( std::abs( temp1 ) < 1.0 )
400  {
401  if( temp2 > 1.0e-10 )result = temp1*temp2;
402  }
403  else result = temp1*temp2;
404  if( result < -ek )result = -ek;
405  return result*GeV;
406  }
407 
408  //
409  // methods for class G4Nucleus ... by Christian Volcker
410  //
411 
413  {
414  // chv: .. we assume zero temperature!
415 
416  // momentum is equally distributed in each phasespace volume dpx, dpy, dpz.
417  G4double ranflat1=
418  G4RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);
419  G4double ranflat2=
420  G4RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);
421  G4double ranflat3=
422  G4RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);
423  G4double ranmax = (ranflat1>ranflat2? ranflat1: ranflat2);
424  ranmax = (ranmax>ranflat3? ranmax : ranflat3);
425 
426  // Isotropic momentum distribution
427  G4double costheta = 2.*G4UniformRand() - 1.0;
428  G4double sintheta = std::sqrt(1.0 - costheta*costheta);
429  G4double phi = 2.0*pi*G4UniformRand();
430 
431  G4double pz=costheta*ranmax;
432  G4double px=sintheta*std::cos(phi)*ranmax;
433  G4double py=sintheta*std::sin(phi)*ranmax;
434  G4ThreeVector p(px,py,pz);
435  return p;
436  }
437 
439  {
440  // needs implementation!
441  return NULL;
442  }
443 
444  void G4Nucleus::AddMomentum(const G4ThreeVector aMomentum)
445  {
446  momentum+=(aMomentum);
447  }
448 
450  {
451  excitationEnergy+=anEnergy;
452  }
453 
454  /* end of file */
455 
G4double G4ParticleHPJENDLHEData::G4double result
size_t GetNumberOfIsotopes() const
Definition: G4Element.hh:159
ThreeVector shoot(const G4int Ap, const G4int Af)
G4double AtomicMass(const G4double A, const G4double Z) const
Definition: G4Nucleus.cc:254
G4double EvaporationEffects(G4double kineticEnergy)
Definition: G4Nucleus.cc:278
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4double GetThermalPz(const G4double mass, const G4double temp) const
Definition: G4Nucleus.cc:268
std::vector< G4Element * > G4ElementVector
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
CLHEP::Hep3Vector G4ThreeVector
G4double GetN() const
Definition: G4Element.hh:135
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:172
G4ReactionProduct GetThermalNucleus(G4double aMass, G4double temp=-1) const
Definition: G4Nucleus.cc:143
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
G4DynamicParticle * ReturnTargetParticle() const
Definition: G4Nucleus.cc:241
G4ReactionProductVector * Fragmentate()
Definition: G4Nucleus.cc:438
std::vector< G4ReactionProduct * > G4ReactionProductVector
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:206
G4ReactionProduct GetBiasedThermalNucleus(G4double aMass, G4ThreeVector aVelocity, G4double temp=-1) const
Definition: G4Nucleus.cc:113
void SetMass(const G4double mas)
G4int GetN() const
Definition: G4Isotope.hh:94
#define G4UniformRand()
Definition: Randomize.hh:97
double A(double temperature)
float k_Boltzmann
Definition: hepunit.py:299
const XML_Char int const XML_Char * value
Definition: expat.h:331
void SetTotalEnergy(const G4double en)
G4ErrorTarget * theTarget
Definition: errprop.cc:59
G4double ek
void AddMomentum(const G4ThreeVector aMomentum)
Definition: G4Nucleus.cc:444
G4double * GetRelativeAbundanceVector() const
Definition: G4Element.hh:167
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
static constexpr double kelvin
Definition: G4SIunits.hh:281
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4int GetZ() const
Definition: G4Isotope.hh:91
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:209
G4ThreeVector GetFermiMomentum()
Definition: G4Nucleus.cc:412
G4double AnnihilationEvaporationEffects(G4double kineticEnergy, G4double ekOrg)
Definition: G4Nucleus.cc:337
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:382
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static constexpr double GeV
Definition: G4SIunits.hh:217
G4ThreeVector GetMomentum() const
const G4Isotope * GetIsotope(G4int iso) const
Definition: G4Element.hh:170
G4double GetTemperature() const
Definition: G4Material.hh:182
#define G4endl
Definition: G4ios.hh:61
static constexpr double pi
Definition: G4SIunits.hh:75
void AddExcitationEnergy(G4double anEnergy)
Definition: G4Nucleus.cc:449
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
double G4double
Definition: G4Types.hh:76
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
static constexpr double fermi
Definition: G4SIunits.hh:103
double mag() const
void SetParameters(const G4double A, const G4double Z)
Definition: G4Nucleus.cc:212
G4double GetMass() const