Geant4  10.02.p02
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;
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;
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 
82 G4Nucleus::G4Nucleus( const G4int A, const G4int Z )
83 {
84  SetParameters( A, Z );
85  pnBlackTrackEnergy = 0.0;
86  dtaBlackTrackEnergy = 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;
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();
116  G4ReactionProduct result;
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)
159  {
160  theTarget.SetTotalEnergy(tEtot);
161  }
162  else
163  {
164  theTarget.SetKineticEnergy(tMom*tMom/(2.*theTarget.GetMass()));
165  }
166  return theTarget;
167 }
168 
169 
170 void
172 {
173  G4double random = G4UniformRand();
174  G4double sum = aMaterial->GetTotNbOfAtomsPerVolume();
175  const G4ElementVector* theElementVector = aMaterial->GetElementVector();
176  G4double running(0);
177  // G4Element* element(0);
178  G4Element* element = (*theElementVector)[aMaterial->GetNumberOfElements()-1];
179 
180  for (unsigned int i = 0; i < aMaterial->GetNumberOfElements(); ++i) {
181  running += aMaterial->GetVecNbOfAtomsPerVolume()[i];
182  if (running > random*sum) {
183  element = (*theElementVector)[i];
184  break;
185  }
186  }
187 
188  if (element->GetNumberOfIsotopes() > 0) {
189  G4double randomAbundance = G4UniformRand();
190  G4double sumAbundance = element->GetRelativeAbundanceVector()[0];
191  unsigned int iso=0;
192  while (iso < element->GetNumberOfIsotopes() && /* Loop checking, 02.11.2015, A.Ribon */
193  sumAbundance < randomAbundance) {
194  ++iso;
195  sumAbundance += element->GetRelativeAbundanceVector()[iso];
196  }
197  theA=element->GetIsotope(iso)->GetN();
198  theZ=element->GetIsotope(iso)->GetZ();
199  aEff=theA;
200  zEff=theZ;
201  } else {
202  aEff = element->GetN();
203  zEff = element->GetZ();
204  theZ = G4int(zEff + 0.5);
205  theA = G4int(aEff + 0.5);
206  }
207 }
208 
209 
210 void
212 {
213  theZ = G4lrint(Z);
214  theA = G4lrint(A);
215  if (theA<1 || theZ<0 || theZ>theA) {
216  throw G4HadronicException(__FILE__, __LINE__,
217  "G4Nucleus::SetParameters called with non-physical parameters");
218  }
219  aEff = A; // atomic weight
220  zEff = Z; // atomic number
221  fIsotope = 0;
222 }
223 
224 void
226 {
227  theZ = Z;
228  theA = A;
229  if( theA<1 || theZ<0 || theZ>theA )
230  {
231  throw G4HadronicException(__FILE__, __LINE__,
232  "G4Nucleus::SetParameters called with non-physical parameters");
233  }
234  aEff = A; // atomic weight
235  zEff = Z; // atomic number
236  fIsotope = 0;
237 }
238 
241  {
242  // choose a proton or a neutron as the target particle
243 
244  G4DynamicParticle *targetParticle = new G4DynamicParticle;
245  if( G4UniformRand() < zEff/aEff )
246  targetParticle->SetDefinition( G4Proton::Proton() );
247  else
248  targetParticle->SetDefinition( G4Neutron::Neutron() );
249  return targetParticle;
250  }
251 
252  G4double
253  G4Nucleus::AtomicMass( const G4double A, const G4double Z ) const
254  {
255  // Now returns (atomic mass - electron masses)
257  }
258 
259  G4double
260  G4Nucleus::AtomicMass( const G4int A, const G4int Z ) const
261  {
262  // Now returns (atomic mass - electron masses)
264  }
265 
266  G4double
267  G4Nucleus::GetThermalPz( const G4double mass, const G4double temp ) const
268  {
269  G4double result = G4RandGauss::shoot();
270  result *= std::sqrt(k_Boltzmann*temp*mass); // Das ist impuls (Pz),
271  // nichtrelativistische rechnung
272  // Maxwell verteilung angenommen
273  return result;
274  }
275 
276  G4double
278  {
279  // derived from original FORTRAN code EXNU by H. Fesefeldt (10-Dec-1986)
280  //
281  // Nuclear evaporation as function of atomic number
282  // and kinetic energy (MeV) of primary particle
283  //
284  // returns kinetic energy (MeV)
285  //
286  if( aEff < 1.5 )
287  {
289  return 0.0;
290  }
291  G4double ek = kineticEnergy/GeV;
292  G4float ekin = std::min( 4.0, std::max( 0.1, ek ) );
293  const G4float atno = std::min( 120., aEff );
294  const G4float gfa = 2.0*((aEff-1.0)/70.)*G4Exp(-(aEff-1.0)/70.);
295  //
296  // 0.35 value at 1 GeV
297  // 0.05 value at 0.1 GeV
298  //
299  G4float cfa = std::max( 0.15, 0.35 + ((0.35-0.05)/2.3)*G4Log(ekin) );
300  G4float exnu = 7.716 * cfa * G4Exp(-cfa)
301  * ((atno-1.0)/120.)*G4Exp(-(atno-1.0)/120.);
302  G4float fpdiv = std::max( 0.5, 1.0-0.25*ekin*ekin );
303  //
304  // pnBlackTrackEnergy is the kinetic energy (in GeV) available for
305  // proton/neutron black track particles
306  // dtaBlackTrackEnergy is the kinetic energy (in GeV) available for
307  // deuteron/triton/alpha black track particles
308  //
309  pnBlackTrackEnergy = exnu*fpdiv;
310  dtaBlackTrackEnergy = exnu*(1.0-fpdiv);
311 
312  if( G4int(zEff+0.1) != 82 )
313  {
314  G4double ran1 = -6.0;
315  G4double ran2 = -6.0;
316  for( G4int i=0; i<12; ++i )
317  {
318  ran1 += G4UniformRand();
319  ran2 += G4UniformRand();
320  }
321  pnBlackTrackEnergy *= 1.0 + ran1*gfa;
322  dtaBlackTrackEnergy *= 1.0 + ran2*gfa;
323  }
326  while( pnBlackTrackEnergy+dtaBlackTrackEnergy >= ek ) /* Loop checking, 02.11.2015, A.Ribon */
327  {
328  pnBlackTrackEnergy *= 1.0 - 0.5*G4UniformRand();
329  dtaBlackTrackEnergy *= 1.0 - 0.5*G4UniformRand();
330  }
331 // G4cout << "EvaporationEffects "<<kineticEnergy<<" "
332 // <<pnBlackTrackEnergy+dtaBlackTrackEnergy<<endl;
334  }
335 
337  {
338  // Nuclear evaporation as a function of atomic number and kinetic
339  // energy (MeV) of primary particle. Modified for annihilation effects.
340  //
341  if( aEff < 1.5 || ekOrg < 0.)
342  {
345  return 0.0;
346  }
347  G4double ek = kineticEnergy/GeV;
348  G4float ekin = std::min( 4.0, std::max( 0.1, ek ) );
349  const G4float atno = std::min( 120., aEff );
350  const G4float gfa = 2.0*((aEff-1.0)/70.)*G4Exp(-(aEff-1.0)/70.);
351 
352  G4float cfa = std::max( 0.15, 0.35 + ((0.35-0.05)/2.3)*G4Log(ekin) );
353  G4float exnu = 7.716 * cfa * G4Exp(-cfa)
354  * ((atno-1.0)/120.)*G4Exp(-(atno-1.0)/120.);
355  G4float fpdiv = std::max( 0.5, 1.0-0.25*ekin*ekin );
356 
358  dtaBlackTrackEnergyfromAnnihilation = exnu*(1.0-fpdiv);
359 
360  G4double ran1 = -6.0;
361  G4double ran2 = -6.0;
362  for( G4int i=0; i<12; ++i ) {
363  ran1 += G4UniformRand();
364  ran2 += G4UniformRand();
365  }
366  pnBlackTrackEnergyfromAnnihilation *= 1.0 + ran1*gfa;
367  dtaBlackTrackEnergyfromAnnihilation *= 1.0 + ran2*gfa;
368 
372  if (blackSum >= ekOrg/GeV) {
373  pnBlackTrackEnergyfromAnnihilation *= ekOrg/GeV/blackSum;
374  dtaBlackTrackEnergyfromAnnihilation *= ekOrg/GeV/blackSum;
375  }
376 
377  return (pnBlackTrackEnergyfromAnnihilation+dtaBlackTrackEnergyfromAnnihilation)*GeV;
378  }
379 
380  G4double
381  G4Nucleus::Cinema( G4double kineticEnergy )
382  {
383  // derived from original FORTRAN code CINEMA by H. Fesefeldt (14-Oct-1987)
384  //
385  // input: kineticEnergy (MeV)
386  // returns modified kinetic energy (MeV)
387  //
388  static const G4double expxu = 82.; // upper bound for arg. of exp
389  static const G4double expxl = -expxu; // lower bound for arg. of exp
390 
391  G4double ek = kineticEnergy/GeV;
392  G4double ekLog = G4Log( ek );
393  G4double aLog = G4Log( aEff );
394  G4double em = std::min( 1.0, 0.2390 + 0.0408*aLog*aLog );
395  G4double temp1 = -ek * std::min( 0.15, 0.0019*aLog*aLog*aLog );
396  G4double temp2 = G4Exp( std::max( expxl, std::min( expxu, -(ekLog-em)*(ekLog-em)*2.0 ) ) );
397  G4double result = 0.0;
398  if( std::abs( temp1 ) < 1.0 )
399  {
400  if( temp2 > 1.0e-10 )result = temp1*temp2;
401  }
402  else result = temp1*temp2;
403  if( result < -ek )result = -ek;
404  return result*GeV;
405  }
406 
407  //
408  // methods for class G4Nucleus ... by Christian Volcker
409  //
410 
412  {
413  // chv: .. we assume zero temperature!
414 
415  // momentum is equally distributed in each phasespace volume dpx, dpy, dpz.
416  G4double ranflat1=
418  G4double ranflat2=
419  G4RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);
420  G4double ranflat3=
421  G4RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);
422  G4double ranmax = (ranflat1>ranflat2? ranflat1: ranflat2);
423  ranmax = (ranmax>ranflat3? ranmax : ranflat3);
424 
425  // Isotropic momentum distribution
426  G4double costheta = 2.*G4UniformRand() - 1.0;
427  G4double sintheta = std::sqrt(1.0 - costheta*costheta);
428  G4double phi = 2.0*pi*G4UniformRand();
429 
430  G4double pz=costheta*ranmax;
431  G4double px=sintheta*std::cos(phi)*ranmax;
432  G4double py=sintheta*std::sin(phi)*ranmax;
433  G4ThreeVector p(px,py,pz);
434  return p;
435  }
436 
438  {
439  // needs implementation!
440  return NULL;
441  }
442 
443  void G4Nucleus::AddMomentum(const G4ThreeVector aMomentum)
444  {
445  momentum+=(aMomentum);
446  }
447 
449  {
450  excitationEnergy+=anEnergy;
451  }
452 
453  /* end of file */
454 
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:253
G4double EvaporationEffects(G4double kineticEnergy)
Definition: G4Nucleus.cc:277
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4double GetThermalPz(const G4double mass, const G4double temp) const
Definition: G4Nucleus.cc:267
std::vector< G4Element * > G4ElementVector
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
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)
float G4float
Definition: G4Types.hh:77
G4double GetZ() const
Definition: G4Element.hh:131
G4double zEff
Definition: G4Nucleus.hh:196
void ChooseParameters(const G4Material *aMaterial)
Definition: G4Nucleus.cc:171
G4double dtaBlackTrackEnergyfromAnnihilation
Definition: G4Nucleus.hh:207
G4double pnBlackTrackEnergyfromAnnihilation
Definition: G4Nucleus.hh:204
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:240
G4ReactionProductVector * Fragmentate()
Definition: G4Nucleus.cc:437
std::vector< G4ReactionProduct * > G4ReactionProductVector
G4int theA
Definition: G4Nucleus.hh:193
G4double fermiMomentum
Definition: G4Nucleus.hh:225
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)
void SetTotalEnergy(const G4double en)
G4ErrorTarget * theTarget
Definition: errprop.cc:59
void AddMomentum(const G4ThreeVector aMomentum)
Definition: G4Nucleus.cc:443
G4double * GetRelativeAbundanceVector() const
Definition: G4Element.hh:166
static G4Proton * Proton()
Definition: G4Proton.cc:93
static const double GeV
Definition: G4SIunits.hh:214
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4double dtaBlackTrackEnergy
Definition: G4Nucleus.hh:202
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static const double kelvin
Definition: G4SIunits.hh:278
G4double theTemp
Definition: G4Nucleus.hh:226
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:411
G4double AnnihilationEvaporationEffects(G4double kineticEnergy, G4double ekOrg)
Definition: G4Nucleus.cc:336
G4double GetPDGMass() const
static const double pi
Definition: G4SIunits.hh:74
int G4lrint(double ad)
Definition: templates.hh:163
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4int theZ
Definition: G4Nucleus.hh:194
G4double Cinema(G4double kineticEnergy)
Definition: G4Nucleus.cc:381
G4double excitationEnergy
Definition: G4Nucleus.hh:214
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
const G4Isotope * fIsotope
Definition: G4Nucleus.hh:198
G4double GetTemperature() const
Definition: G4Material.hh:182
#define G4endl
Definition: G4ios.hh:61
void AddExcitationEnergy(G4double anEnergy)
Definition: G4Nucleus.cc:448
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
double G4double
Definition: G4Types.hh:76
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
G4double aEff
Definition: G4Nucleus.hh:195
void SetParameters(const G4double A, const G4double Z)
Definition: G4Nucleus.cc:211
G4ThreeVector momentum
Definition: G4Nucleus.hh:217
G4double GetMass() const
static const double fermi
Definition: G4SIunits.hh:102
G4double pnBlackTrackEnergy
Definition: G4Nucleus.hh:200