Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/util/src/G4Nucleus.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/util/src/G4Nucleus.cc (Version 11.3.0) and /processes/hadronic/util/src/G4Nucleus.cc (Version 9.4.p2)


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