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.3.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 << 
 41 // A. Ribon, summer 2015:  migrated to G4Exp a << 
 42 // A. Ribon, autumn 2021:  extended to hypernu << 
 43                                                    40  
 44 #include "G4Nucleus.hh"                            41 #include "G4Nucleus.hh"
 45 #include "G4NucleiProperties.hh"                   42 #include "G4NucleiProperties.hh"
 46 #include "G4PhysicalConstants.hh"              << 
 47 #include "G4SystemOfUnits.hh"                  << 
 48 #include "Randomize.hh"                            43 #include "Randomize.hh"
 49 #include "G4HadronicException.hh"                  44 #include "G4HadronicException.hh"
 50 #include "G4Exp.hh"                            <<  45  
 51 #include "G4Log.hh"                            << 
 52 #include "G4HyperNucleiProperties.hh"          << 
 53 #include "G4HadronicParameters.hh"             << 
 54                                                << 
 55                                                << 
 56 G4Nucleus::G4Nucleus()                             46 G4Nucleus::G4Nucleus()
 57   : theA(0), theZ(0), theL(0), aEff(0.0), zEff << 
 58 {                                                  47 {
 59   pnBlackTrackEnergy = 0.0;                        48   pnBlackTrackEnergy = 0.0;
 60   dtaBlackTrackEnergy = 0.0;                       49   dtaBlackTrackEnergy = 0.0;
 61   pnBlackTrackEnergyfromAnnihilation = 0.0;        50   pnBlackTrackEnergyfromAnnihilation = 0.0;
 62   dtaBlackTrackEnergyfromAnnihilation = 0.0;       51   dtaBlackTrackEnergyfromAnnihilation = 0.0;
 63   excitationEnergy = 0.0;                          52   excitationEnergy = 0.0;
 64   momentum = G4ThreeVector(0.,0.,0.);              53   momentum = G4ThreeVector(0.,0.,0.);
 65   fermiMomentum = 1.52*hbarc/fermi;                54   fermiMomentum = 1.52*hbarc/fermi;
 66   theTemp = 293.16*kelvin;                         55   theTemp = 293.16*kelvin;
 67   fIsotope = 0;                                << 
 68 }                                                  56 }
 69                                                    57 
 70 G4Nucleus::G4Nucleus( const G4double A, const  <<  58 G4Nucleus::G4Nucleus( const G4double A, const G4double Z )
 71 {                                                  59 {
 72   SetParameters( A, Z, std::max(numberOfLambda <<  60   SetParameters( A, Z );
 73   pnBlackTrackEnergy = 0.0;                        61   pnBlackTrackEnergy = 0.0;
 74   dtaBlackTrackEnergy = 0.0;                       62   dtaBlackTrackEnergy = 0.0;
 75   pnBlackTrackEnergyfromAnnihilation = 0.0;        63   pnBlackTrackEnergyfromAnnihilation = 0.0;
 76   dtaBlackTrackEnergyfromAnnihilation = 0.0;       64   dtaBlackTrackEnergyfromAnnihilation = 0.0;
 77   excitationEnergy = 0.0;                          65   excitationEnergy = 0.0;
 78   momentum = G4ThreeVector(0.,0.,0.);              66   momentum = G4ThreeVector(0.,0.,0.);
 79   fermiMomentum = 1.52*hbarc/fermi;                67   fermiMomentum = 1.52*hbarc/fermi;
 80   theTemp = 293.16*kelvin;                         68   theTemp = 293.16*kelvin;
 81   fIsotope = 0;                                << 
 82 }                                              << 
 83                                                << 
 84 G4Nucleus::G4Nucleus( const G4int A, const G4i << 
 85 {                                              << 
 86   SetParameters( A, Z, std::max(numberOfLambda << 
 87   pnBlackTrackEnergy = 0.0;                    << 
 88   dtaBlackTrackEnergy = 0.0;                   << 
 89   pnBlackTrackEnergyfromAnnihilation = 0.0;    << 
 90   dtaBlackTrackEnergyfromAnnihilation = 0.0;   << 
 91   excitationEnergy = 0.0;                      << 
 92   momentum = G4ThreeVector(0.,0.,0.);          << 
 93   fermiMomentum = 1.52*hbarc/fermi;            << 
 94   theTemp = 293.16*kelvin;                     << 
 95   fIsotope = 0;                                << 
 96 }                                                  69 }
 97                                                    70 
 98 G4Nucleus::G4Nucleus( const G4Material *aMater     71 G4Nucleus::G4Nucleus( const G4Material *aMaterial )
 99 {                                                  72 {
100   ChooseParameters( aMaterial );                   73   ChooseParameters( aMaterial );
101   pnBlackTrackEnergy = 0.0;                        74   pnBlackTrackEnergy = 0.0;
102   dtaBlackTrackEnergy = 0.0;                       75   dtaBlackTrackEnergy = 0.0;
103   pnBlackTrackEnergyfromAnnihilation = 0.0;        76   pnBlackTrackEnergyfromAnnihilation = 0.0;
104   dtaBlackTrackEnergyfromAnnihilation = 0.0;       77   dtaBlackTrackEnergyfromAnnihilation = 0.0;
105   excitationEnergy = 0.0;                          78   excitationEnergy = 0.0;
106   momentum = G4ThreeVector(0.,0.,0.);              79   momentum = G4ThreeVector(0.,0.,0.);
107   fermiMomentum = 1.52*hbarc/fermi;                80   fermiMomentum = 1.52*hbarc/fermi;
108   theTemp = aMaterial->GetTemperature();           81   theTemp = aMaterial->GetTemperature();
109   fIsotope = 0;                                << 
110 }                                                  82 }
111                                                    83 
112 G4Nucleus::~G4Nucleus() {}                         84 G4Nucleus::~G4Nucleus() {}
113                                                    85 
114                                                <<  86 G4ReactionProduct G4Nucleus::
115 //-------------------------------------------- <<  87 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 {                                                  88 {
121   // If E_neutron <= E_threshold, Then apply t <<  89   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;                        90   G4ReactionProduct result;
130   result.SetMass(aMass*G4Neutron::Neutron()->G <<  91   G4double value = 0;
131                                                <<  92   G4double random = 1;
132   if ( E_neutron <= E_threshold ) {            <<  93   G4double norm = 3.*std::sqrt(k_Boltzmann*temp*aMass*G4Neutron::Neutron()->GetPDGMass());
133                                                <<  94   norm /= G4Neutron::Neutron()->GetPDGMass();
134     // Beta = sqrt(m/2kT)                      <<  95   norm *= 5.;
135     G4double beta = std::sqrt(result.GetMass() <<  96   norm += velMag;
136                                                <<  97   norm /= velMag;
137     // Neutron speed vn                        <<  98   while(value/norm<random)
138     G4double vN_norm = aVelocity.mag();        <<  99   {
139     G4double vN_norm2 = vN_norm*vN_norm;       << 100      result = GetThermalNucleus(aMass, temp);
140     G4double y = beta*vN_norm;                 << 101      G4ThreeVector targetVelocity = 1./result.GetMass()*result.GetMomentum();
141                                                << 102      value = (targetVelocity+aVelocity).mag()/velMag;
142     // Normalize neutron velocity              << 103      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   }                                               104   }
184                                                << 
185   return result;                                  105   return result;
186 }                                                 106 }
187                                                   107 
188                                                << 108 G4ReactionProduct G4Nucleus::GetThermalNucleus(G4double targetMass, G4double temp) const
189 void                                           << 109   {
190 G4Nucleus::DoKinematicsOfThermalNucleus(const  << 110     G4double currentTemp = temp;
191                                         G4Reac << 111     if(currentTemp < 0) currentTemp = theTemp;
192                                                << 112     G4ReactionProduct theTarget;    
193   // Get target nucleus direction from the neu << 113     theTarget.SetMass(targetMass*G4Neutron::Neutron()->GetPDGMass());
194   G4double cosTh = mu;                         << 114     G4double px, py, pz;
195   G4ThreeVector uNorm = aVelocity;             << 115     px = GetThermalPz(theTarget.GetMass(), currentTemp);
196                                                << 116     py = GetThermalPz(theTarget.GetMass(), currentTemp);
197   G4double sinTh = std::sqrt(1. - cosTh*cosTh) << 117     pz = GetThermalPz(theTarget.GetMass(), currentTemp);
198                                                << 118     theTarget.SetMomentum(px, py, pz);
199   // Sample randomly the phi angle between the << 119     G4double tMom = std::sqrt(px*px+py*py+pz*pz);
200   G4double phi = CLHEP::twopi*G4UniformRand(); << 120     G4double tEtot = std::sqrt((tMom+theTarget.GetMass())*
201   G4double sinPhi = std::sin(phi);             << 121                           (tMom+theTarget.GetMass())-
202   G4double cosPhi = std::cos(phi);             << 122                           2.*tMom*theTarget.GetMass());
203                                                << 123     if(1-tEtot/theTarget.GetMass()>0.001)
204   // Find orthogonal vector to aVelocity - sol << 124     {
205   G4ThreeVector ortho(1., 1., 1.);             << 125       theTarget.SetTotalEnergy(tEtot);
206   if      ( uNorm[0] )  ortho[0] = -(uNorm[1]+ << 126     }
207   else if ( uNorm[1] )  ortho[1] = -(uNorm[0]+ << 127     else
208   else if ( uNorm[2] )  ortho[2] = -(uNorm[0]+ << 128     {
209                                                << 129       theTarget.SetKineticEnergy(tMom*tMom/(2.*theTarget.GetMass()));
210   // Normalize the vector                      << 130     }    
211   ortho = (1/ortho.mag())*ortho;               << 131     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   }                                               132   }
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                                                   133  
276 void                                           << 134  void
277 G4Nucleus::ChooseParameters(const G4Material*  << 135   G4Nucleus::ChooseParameters( const G4Material *aMaterial )
278 {                                              << 136   {
279   G4double random = G4UniformRand();           << 137     G4double random = G4UniformRand();
280   G4double sum = aMaterial->GetTotNbOfAtomsPer << 138     G4double sum = 0;
281   const G4ElementVector* theElementVector = aM << 139     const G4ElementVector *theElementVector = aMaterial->GetElementVector();
282   G4double running(0);                         << 140     unsigned int i;
283   //  G4Element* element(0);                   << 141     for(i=0; i<aMaterial->GetNumberOfElements(); ++i )
284   const G4Element* element = (*theElementVecto << 142     {
285                                                << 143       sum += aMaterial->GetAtomicNumDensityVector()[i];
286   for (unsigned int i = 0; i < aMaterial->GetN << 144     }
287     running += aMaterial->GetVecNbOfAtomsPerVo << 145     G4double running = 0;
288     if (running > random*sum) {                << 146     for(i=0; i<aMaterial->GetNumberOfElements(); ++i )
289       element = (*theElementVector)[i];        << 147     {
290       break;                                   << 148       running += aMaterial->GetAtomicNumDensityVector()[i];
                                                   >> 149       if( running/sum > random ) {
                                                   >> 150         aEff = (*theElementVector)[i]->GetA()*mole/g;
                                                   >> 151         zEff = (*theElementVector)[i]->GetZ();
                                                   >> 152         break;
                                                   >> 153       }
291     }                                             154     }
292   }                                               155   }
293                                                << 156  
294   if (element->GetNumberOfIsotopes() > 0) {    << 157  void
295     G4double randomAbundance = G4UniformRand() << 158   G4Nucleus::SetParameters( const G4double A, const G4double Z )
296     G4double sumAbundance = element->GetRelati << 159   {
297     unsigned int iso=0;                        << 160     G4int myZ = G4int(Z + 0.5);
298     while (iso < element->GetNumberOfIsotopes( << 161     G4int myA = G4int(A + 0.5);   
299            sumAbundance < randomAbundance) {   << 162     if( myA<1 || myZ<0 || myZ>myA )
300       ++iso;                                   << 
301       sumAbundance += element->GetRelativeAbun << 
302     }                                          << 
303     theA=element->GetIsotope(iso)->GetN();     << 
304     theZ=element->GetIsotope(iso)->GetZ();     << 
305     theL=0;                                    << 
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   }                                            << 
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                                                << 
334 void                                           << 
335 G4Nucleus::SetParameters( const G4int A, const << 
336 {                                              << 
337   theZ = Z;                                    << 
338   theA = A;                                    << 
339   theL = std::max(numberOfLambdas, 0);         << 
340   if( theA<1 || theZ<0 || theZ>theA )          << 
341     {                                             163     {
342       throw G4HadronicException(__FILE__, __LI    164       throw G4HadronicException(__FILE__, __LINE__,
343         "G4Nucleus::SetParameters called with  << 165                                "G4Nucleus::SetParameters called with non-physical parameters");
344     }                                             166     }
345   aEff = A;  // atomic weight                  << 167     aEff = A;  // atomic weight
346   zEff = Z;  // atomic number                  << 168     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   }                                               169   }
364   return targetParticle;                       << 
365 }                                              << 
366                                                   170 
367                                                << 171  G4DynamicParticle *
368 G4double                                       << 172   G4Nucleus::ReturnTargetParticle() const
369 G4Nucleus::AtomicMass( const G4double A, const << 173   {
370 {                                              << 174     // choose a proton or a neutron as the target particle
371   // Now returns (atomic mass - electron masse << 175     
372   if ( numberOfLambdas > 0 ) {                 << 176     G4DynamicParticle *targetParticle = new G4DynamicParticle;
373     return G4HyperNucleiProperties::GetNuclear << 177     if( G4UniformRand() < zEff/aEff )
374   } else {                                     << 178       targetParticle->SetDefinition( G4Proton::Proton() );
375     return G4NucleiProperties::GetNuclearMass( << 179     else
                                                   >> 180       targetParticle->SetDefinition( G4Neutron::Neutron() );
                                                   >> 181     return targetParticle;
376   }                                               182   }
377 }                                              << 
378                                                << 
379                                                   183  
380 G4double                                       << 184  G4double
381 G4Nucleus::AtomicMass( const G4int A, const G4 << 185   G4Nucleus::AtomicMass( const G4double A, const G4double Z ) const
382 {                                              << 186   {
383   // Now returns (atomic mass - electron masse << 187     // Now returns (atomic mass - electron masses) 
384   if ( numberOfLambdas > 0 ) {                 << 
385     return G4HyperNucleiProperties::GetNuclear << 
386   } else {                                     << 
387     return G4NucleiProperties::GetNuclearMass(    188     return G4NucleiProperties::GetNuclearMass(A, Z);
388   }                                               189   }
389 }                                              << 
390                                                   190  
                                                   >> 191  G4double
                                                   >> 192   G4Nucleus::GetThermalPz( const G4double mass, const G4double temp ) const
                                                   >> 193   {
                                                   >> 194     G4double result = G4RandGauss::shoot();
                                                   >> 195     result *= std::sqrt(k_Boltzmann*temp*mass); // Das ist impuls (Pz),
                                                   >> 196                                            // nichtrelativistische rechnung
                                                   >> 197                                            // Maxwell verteilung angenommen
                                                   >> 198     return result;
                                                   >> 199   }
391                                                   200  
392 G4double                                       << 201  G4double 
393 G4Nucleus::GetThermalPz( const G4double mass,  << 202   G4Nucleus::EvaporationEffects( G4double kineticEnergy )
394 {                                              << 203   {
395   G4double result = G4RandGauss::shoot();      << 204     // derived from original FORTRAN code EXNU by H. Fesefeldt (10-Dec-1986)
396   result *= std::sqrt(k_Boltzmann*temp*mass);  << 205     //
397                                          // ni << 206     // Nuclear evaporation as function of atomic number
398                                          // Ma << 207     // and kinetic energy (MeV) of primary particle
399   return result;                               << 208     //
400 }                                              << 209     // returns kinetic energy (MeV)
                                                   >> 210     //
                                                   >> 211     if( aEff < 1.5 )
                                                   >> 212     {
                                                   >> 213       pnBlackTrackEnergy = dtaBlackTrackEnergy = 0.0;
                                                   >> 214       return 0.0;
                                                   >> 215     }
                                                   >> 216     G4double ek = kineticEnergy/GeV;
                                                   >> 217     G4float ekin = std::min( 4.0, std::max( 0.1, ek ) );
                                                   >> 218     const G4float atno = std::min( 120., aEff ); 
                                                   >> 219     const G4float gfa = 2.0*((aEff-1.0)/70.)*std::exp(-(aEff-1.0)/70.);
                                                   >> 220     //
                                                   >> 221     // 0.35 value at 1 GeV
                                                   >> 222     // 0.05 value at 0.1 GeV
                                                   >> 223     //
                                                   >> 224     G4float cfa = std::max( 0.15, 0.35 + ((0.35-0.05)/2.3)*std::log(ekin) );
                                                   >> 225     G4float exnu = 7.716 * cfa * std::exp(-cfa)
                                                   >> 226       * ((atno-1.0)/120.)*std::exp(-(atno-1.0)/120.);
                                                   >> 227     G4float fpdiv = std::max( 0.5, 1.0-0.25*ekin*ekin );
                                                   >> 228     //
                                                   >> 229     // pnBlackTrackEnergy  is the kinetic energy (in GeV) available for
                                                   >> 230     //                     proton/neutron black track particles
                                                   >> 231     // dtaBlackTrackEnergy is the kinetic energy (in GeV) available for
                                                   >> 232     //                     deuteron/triton/alpha black track particles
                                                   >> 233     //
                                                   >> 234     pnBlackTrackEnergy = exnu*fpdiv;
                                                   >> 235     dtaBlackTrackEnergy = exnu*(1.0-fpdiv);
                                                   >> 236     
                                                   >> 237     if( G4int(zEff+0.1) != 82 )
                                                   >> 238     { 
                                                   >> 239       G4double ran1 = -6.0;
                                                   >> 240       G4double ran2 = -6.0;
                                                   >> 241       for( G4int i=0; i<12; ++i )
                                                   >> 242       {
                                                   >> 243         ran1 += G4UniformRand();
                                                   >> 244         ran2 += G4UniformRand();
                                                   >> 245       }
                                                   >> 246       pnBlackTrackEnergy *= 1.0 + ran1*gfa;
                                                   >> 247       dtaBlackTrackEnergy *= 1.0 + ran2*gfa;
                                                   >> 248     }
                                                   >> 249     pnBlackTrackEnergy = std::max( 0.0, pnBlackTrackEnergy );
                                                   >> 250     dtaBlackTrackEnergy = std::max( 0.0, dtaBlackTrackEnergy );
                                                   >> 251     while( pnBlackTrackEnergy+dtaBlackTrackEnergy >= ek )
                                                   >> 252     {
                                                   >> 253       pnBlackTrackEnergy *= 1.0 - 0.5*G4UniformRand();
                                                   >> 254       dtaBlackTrackEnergy *= 1.0 - 0.5*G4UniformRand();
                                                   >> 255     }
                                                   >> 256 //    G4cout << "EvaporationEffects "<<kineticEnergy<<" "
                                                   >> 257 //           <<pnBlackTrackEnergy+dtaBlackTrackEnergy<<endl;
                                                   >> 258     return (pnBlackTrackEnergy+dtaBlackTrackEnergy)*GeV;
                                                   >> 259   }
401                                                   260  
                                                   >> 261  G4double G4Nucleus::AnnihilationEvaporationEffects(G4double kineticEnergy, G4double ekOrg)
                                                   >> 262   {
                                                   >> 263     // Nuclear evaporation as a function of atomic number and kinetic 
                                                   >> 264     // energy (MeV) of primary particle.  Modified for annihilation effects. 
                                                   >> 265     //
                                                   >> 266     if( aEff < 1.5 || ekOrg < 0.)
                                                   >> 267     {
                                                   >> 268       pnBlackTrackEnergyfromAnnihilation = 0.0;
                                                   >> 269       dtaBlackTrackEnergyfromAnnihilation = 0.0;
                                                   >> 270       return 0.0;
                                                   >> 271     }
                                                   >> 272     G4double ek = kineticEnergy/GeV;
                                                   >> 273     G4float ekin = std::min( 4.0, std::max( 0.1, ek ) );
                                                   >> 274     const G4float atno = std::min( 120., aEff ); 
                                                   >> 275     const G4float gfa = 2.0*((aEff-1.0)/70.)*std::exp(-(aEff-1.0)/70.);
                                                   >> 276 
                                                   >> 277     G4float cfa = std::max( 0.15, 0.35 + ((0.35-0.05)/2.3)*std::log(ekin) );
                                                   >> 278     G4float exnu = 7.716 * cfa * std::exp(-cfa)
                                                   >> 279       * ((atno-1.0)/120.)*std::exp(-(atno-1.0)/120.);
                                                   >> 280     G4float fpdiv = std::max( 0.5, 1.0-0.25*ekin*ekin );
402                                                   281 
403 G4double                                       << 282     pnBlackTrackEnergyfromAnnihilation = exnu*fpdiv;
404 G4Nucleus::EvaporationEffects( G4double kineti << 283     dtaBlackTrackEnergyfromAnnihilation = exnu*(1.0-fpdiv);
405 {                                              << 284     
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;                         285     G4double ran1 = -6.0;
442     G4double ran2 = -6.0;                         286     G4double ran2 = -6.0;
443     for( G4int i=0; i<12; ++i )                << 287     for( G4int i=0; i<12; ++i ) {
444     {                                          << 
445       ran1 += G4UniformRand();                    288       ran1 += G4UniformRand();
446       ran2 += G4UniformRand();                    289       ran2 += G4UniformRand();
447     }                                             290     }
448     pnBlackTrackEnergy *= 1.0 + ran1*gfa;      << 291     pnBlackTrackEnergyfromAnnihilation *= 1.0 + ran1*gfa;
449     dtaBlackTrackEnergy *= 1.0 + ran2*gfa;     << 292     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                                                << 
463                                                << 
464 G4double                                       << 
465 G4Nucleus::AnnihilationEvaporationEffects(G4do << 
466 {                                              << 
467   // Nuclear evaporation as a function of atom << 
468   // energy (MeV) of primary particle.  Modifi << 
469   //                                           << 
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                                                   293 
506   return (pnBlackTrackEnergyfromAnnihilation+d << 294     pnBlackTrackEnergyfromAnnihilation = std::max( 0.0, pnBlackTrackEnergyfromAnnihilation);
507 }                                              << 295     dtaBlackTrackEnergyfromAnnihilation = std::max( 0.0, dtaBlackTrackEnergyfromAnnihilation);
                                                   >> 296     G4double blackSum = pnBlackTrackEnergyfromAnnihilation+dtaBlackTrackEnergyfromAnnihilation;
                                                   >> 297     if (blackSum >= ekOrg/GeV) {
                                                   >> 298       pnBlackTrackEnergyfromAnnihilation *= ekOrg/GeV/blackSum;
                                                   >> 299       dtaBlackTrackEnergyfromAnnihilation *= ekOrg/GeV/blackSum;
                                                   >> 300     }
508                                                   301 
                                                   >> 302     return (pnBlackTrackEnergyfromAnnihilation+dtaBlackTrackEnergyfromAnnihilation)*GeV;
                                                   >> 303   }
509                                                   304  
510 G4double                                       << 305  G4double 
511 G4Nucleus::Cinema( G4double kineticEnergy )    << 306   G4Nucleus::Cinema( G4double kineticEnergy )
512 {                                              << 307   {
513   // derived from original FORTRAN code CINEMA << 308     // derived from original FORTRAN code CINEMA by H. Fesefeldt (14-Oct-1987)
514   //                                           << 309     //
515   // input: kineticEnergy (MeV)                << 310     // input: kineticEnergy (MeV)
516   // returns modified kinetic energy (MeV)     << 311     // returns modified kinetic energy (MeV)
517   //                                           << 312     //
518   static const G4double expxu =  82.;          << 313     static const G4double expxu =  82.;           // upper bound for arg. of exp
519   static const G4double expxl = -expxu;        << 314     static const G4double expxl = -expxu;         // lower bound for arg. of exp
520                                                << 315     
521   G4double ek = kineticEnergy/GeV;             << 316     G4double ek = kineticEnergy/GeV;
522   G4double ekLog = G4Log( ek );                << 317     G4double ekLog = std::log( ek );
523   G4double aLog = G4Log( aEff );               << 318     G4double aLog = std::log( aEff );
524   G4double em = std::min( 1.0, 0.2390 + 0.0408 << 319     G4double em = std::min( 1.0, 0.2390 + 0.0408*aLog*aLog );
525   G4double temp1 = -ek * std::min( 0.15, 0.001 << 320     G4double temp1 = -ek * std::min( 0.15, 0.0019*aLog*aLog*aLog );
526   G4double temp2 = G4Exp( std::max( expxl, std << 321     G4double temp2 = std::exp( std::max( expxl, std::min( expxu, -(ekLog-em)*(ekLog-em)*2.0 ) ) );
527   G4double result = 0.0;                       << 322     G4double result = 0.0;
528   if( std::abs( temp1 ) < 1.0 )                << 323     if( std::abs( temp1 ) < 1.0 )
529   {                                            << 324     {
530     if( temp2 > 1.0e-10 )result = temp1*temp2; << 325       if( temp2 > 1.0e-10 )result = temp1*temp2;
531   }                                            << 326     }
532   else result = temp1*temp2;                   << 327     else result = temp1*temp2;
533   if( result < -ek )result = -ek;              << 328     if( result < -ek )result = -ek;
534   return result*GeV;                           << 329     return result*GeV;
535 }                                              << 330   }
536                                                   331 
                                                   >> 332  //
                                                   >> 333  // methods for class G4Nucleus  ... by Christian Volcker
                                                   >> 334  //
537                                                   335 
538 G4ThreeVector G4Nucleus::GetFermiMomentum()    << 336  G4ThreeVector G4Nucleus::GetFermiMomentum()
539 {                                              << 337   {
540   // chv: .. we assume zero temperature!       << 338     // chv: .. we assume zero temperature!
541                                                << 339     
542   // momentum is equally distributed in each p << 340     // momentum is equally distributed in each phasespace volume dpx, dpy, dpz.
543   G4double ranflat1=                           << 341     G4double ranflat1=
544     G4RandFlat::shoot((G4double)0.,(G4double)f << 342       CLHEP::RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);   
545   G4double ranflat2=                           << 343     G4double ranflat2=
546     G4RandFlat::shoot((G4double)0.,(G4double)f << 344       CLHEP::RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);   
547   G4double ranflat3=                           << 345     G4double ranflat3=
548     G4RandFlat::shoot((G4double)0.,(G4double)f << 346       CLHEP::RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);   
549   G4double ranmax = (ranflat1>ranflat2? ranfla << 347     G4double ranmax = (ranflat1>ranflat2? ranflat1: ranflat2);
550   ranmax = (ranmax>ranflat3? ranmax : ranflat3 << 348     ranmax = (ranmax>ranflat3? ranmax : ranflat3);
551                                                << 349     
552   // Isotropic momentum distribution           << 350     // Isotropic momentum distribution
553   G4double costheta = 2.*G4UniformRand() - 1.0 << 351     G4double costheta = 2.*G4UniformRand() - 1.0;
554   G4double sintheta = std::sqrt(1.0 - costheta << 352     G4double sintheta = std::sqrt(1.0 - costheta*costheta);
555   G4double phi = 2.0*pi*G4UniformRand();       << 353     G4double phi = 2.0*pi*G4UniformRand();
556                                                << 354     
557   G4double pz=costheta*ranmax;                 << 355     G4double pz=costheta*ranmax;
558   G4double px=sintheta*std::cos(phi)*ranmax;   << 356     G4double px=sintheta*std::cos(phi)*ranmax;
559   G4double py=sintheta*std::sin(phi)*ranmax;   << 357     G4double py=sintheta*std::sin(phi)*ranmax;
560   G4ThreeVector p(px,py,pz);                   << 358     G4ThreeVector p(px,py,pz);
561   return p;                                    << 359     return p;
562 }                                              << 360   }
563                                                   361  
564                                                << 362  G4ReactionProductVector* G4Nucleus::Fragmentate()
565 G4ReactionProductVector* G4Nucleus::Fragmentat << 363   {
566 {                                              << 364     // needs implementation!
567   // needs implementation!                     << 365     return NULL;
568   return nullptr;                              << 366   }
569 }                                              << 
570                                                   367  
571                                                << 368  void G4Nucleus::AddMomentum(const G4ThreeVector aMomentum)
572 void G4Nucleus::AddMomentum(const G4ThreeVecto << 369   {
573 {                                              << 370     momentum+=(aMomentum);
574   momentum+=(aMomentum);                       << 371   }
575 }                                              << 
576                                                << 
577                                                   372  
578 void G4Nucleus::AddExcitationEnergy( G4double  << 373  void G4Nucleus::AddExcitationEnergy( G4double anEnergy )
579 {                                              << 374   {
580   excitationEnergy+=anEnergy;                  << 375     excitationEnergy+=anEnergy;
581 }                                              << 376   }
582                                                   377 
583  /* end of file */                                378  /* end of file */
584                                                   379 
585                                                   380