ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Nucleus.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4Nucleus.cc
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 
54 G4Nucleus::G4Nucleus()
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 
68 G4Nucleus::G4Nucleus( const G4double A, const G4double Z )
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 
82 G4Nucleus::G4Nucleus( const G4int A, const G4int Z )
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 
110 G4Nucleus::~G4Nucleus() {}
111 
112 G4ReactionProduct G4Nucleus::
113 GetBiasedThermalNucleus(G4double aMass, G4ThreeVector aVelocity, G4double temp) const
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 ) {
134  G4ExceptionDescription ed;
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 
142 G4ReactionProduct
143 G4Nucleus::GetThermalNucleus(G4double targetMass, G4double temp) const
144 {
145  G4double currentTemp = temp;
146  if (currentTemp < 0) currentTemp = theTemp;
147  G4ReactionProduct theTarget;
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
172 G4Nucleus::ChooseParameters(const G4Material* aMaterial)
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
212 G4Nucleus::SetParameters(G4double A, G4double Z)
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
226 G4Nucleus::SetParameters(G4int A, const G4int Z )
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 
240  G4DynamicParticle *
241  G4Nucleus::ReturnTargetParticle() const
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)
257  return G4NucleiProperties::GetNuclearMass(A, Z);
258  }
259 
260  G4double
261  G4Nucleus::AtomicMass( const G4int A, const G4int Z ) const
262  {
263  // Now returns (atomic mass - electron masses)
264  return G4NucleiProperties::GetNuclearMass(A, Z);
265  }
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  }
276 
277  G4double
278  G4Nucleus::EvaporationEffects( G4double kineticEnergy )
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 
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 );
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 
412  G4ThreeVector G4Nucleus::GetFermiMomentum()
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 
438  G4ReactionProductVector* G4Nucleus::Fragmentate()
439  {
440  // needs implementation!
441  return NULL;
442  }
443 
444  void G4Nucleus::AddMomentum(const G4ThreeVector aMomentum)
445  {
446  momentum+=(aMomentum);
447  }
448 
449  void G4Nucleus::AddExcitationEnergy( G4double anEnergy )
450  {
451  excitationEnergy+=anEnergy;
452  }
453 
454  /* end of file */
455