source: trunk/source/processes/hadronic/util/src/G4Nucleus.cc @ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

geant4 tag 9.4
File size: 14.7 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26//
27//
28 // original by H.P. Wellisch
29 // modified by J.L. Chuma, TRIUMF, 19-Nov-1996
30 // last modified: 27-Mar-1997
31 // J.P.Wellisch: 23-Apr-97: minor simplifications
32 // modified by J.L.Chuma 24-Jul-97  to set the total momentum in Cinema and
33 //                                  EvaporationEffects
34 // modified by J.L.Chuma 21-Oct-97  put std::abs() around the totalE^2-mass^2
35 //                                  in calculation of total momentum in
36 //                                  Cinema and EvaporationEffects
37 // Chr. Volcker, 10-Nov-1997: new methods and class variables.
38 // HPW added utilities for low energy neutron transport. (12.04.1998)
39 // M.G. Pia, 2 Oct 1998: modified GetFermiMomentum to avoid memory leaks
40 // G.Folger, spring 2010:  add integer A/Z interface
41 
42#include "G4Nucleus.hh"
43#include "G4NucleiProperties.hh"
44#include "Randomize.hh"
45#include "G4HadronicException.hh"
46 
47G4Nucleus::G4Nucleus()
48  : theA(0), theZ(0), aEff(0.0), zEff(0)
49{
50  pnBlackTrackEnergy = 0.0;
51  dtaBlackTrackEnergy = 0.0;
52  pnBlackTrackEnergyfromAnnihilation = 0.0;
53  dtaBlackTrackEnergyfromAnnihilation = 0.0;
54  excitationEnergy = 0.0;
55  momentum = G4ThreeVector(0.,0.,0.);
56  fermiMomentum = 1.52*hbarc/fermi;
57  theTemp = 293.16*kelvin;
58}
59
60G4Nucleus::G4Nucleus( const G4double A, const G4double Z )
61{
62  SetParameters( A, Z );
63  pnBlackTrackEnergy = 0.0;
64  dtaBlackTrackEnergy = 0.0;
65  pnBlackTrackEnergyfromAnnihilation = 0.0;
66  dtaBlackTrackEnergyfromAnnihilation = 0.0;
67  excitationEnergy = 0.0;
68  momentum = G4ThreeVector(0.,0.,0.);
69  fermiMomentum = 1.52*hbarc/fermi;
70  theTemp = 293.16*kelvin;
71}
72
73G4Nucleus::G4Nucleus( const G4int A, const G4int Z )
74{
75  SetParameters( A, Z );
76  pnBlackTrackEnergy = 0.0;
77  dtaBlackTrackEnergy = 0.0;
78  pnBlackTrackEnergyfromAnnihilation = 0.0;
79  dtaBlackTrackEnergyfromAnnihilation = 0.0;
80  excitationEnergy = 0.0;
81  momentum = G4ThreeVector(0.,0.,0.);
82  fermiMomentum = 1.52*hbarc/fermi;
83  theTemp = 293.16*kelvin;
84}
85
86G4Nucleus::G4Nucleus( const G4Material *aMaterial )
87{
88  ChooseParameters( aMaterial );
89  pnBlackTrackEnergy = 0.0;
90  dtaBlackTrackEnergy = 0.0;
91  pnBlackTrackEnergyfromAnnihilation = 0.0;
92  dtaBlackTrackEnergyfromAnnihilation = 0.0;
93  excitationEnergy = 0.0;
94  momentum = G4ThreeVector(0.,0.,0.);
95  fermiMomentum = 1.52*hbarc/fermi;
96  theTemp = aMaterial->GetTemperature();
97}
98
99G4Nucleus::~G4Nucleus() {}
100
101G4ReactionProduct G4Nucleus::
102GetBiasedThermalNucleus(G4double aMass, G4ThreeVector aVelocity, G4double temp) const
103{
104  G4double velMag = aVelocity.mag();
105  G4ReactionProduct result;
106  G4double value = 0;
107  G4double random = 1;
108  G4double norm = 3.*std::sqrt(k_Boltzmann*temp*aMass*G4Neutron::Neutron()->GetPDGMass());
109  norm /= G4Neutron::Neutron()->GetPDGMass();
110  norm *= 5.;
111  norm += velMag;
112  norm /= velMag;
113  while(value/norm<random)
114  {
115     result = GetThermalNucleus(aMass, temp);
116     G4ThreeVector targetVelocity = 1./result.GetMass()*result.GetMomentum();
117     value = (targetVelocity+aVelocity).mag()/velMag;
118     random = G4UniformRand();
119  }
120  return result;
121}
122
123G4ReactionProduct G4Nucleus::GetThermalNucleus(G4double targetMass, G4double temp) const
124  {
125    G4double currentTemp = temp;
126    if(currentTemp < 0) currentTemp = theTemp;
127    G4ReactionProduct theTarget;   
128    theTarget.SetMass(targetMass*G4Neutron::Neutron()->GetPDGMass());
129    G4double px, py, pz;
130    px = GetThermalPz(theTarget.GetMass(), currentTemp);
131    py = GetThermalPz(theTarget.GetMass(), currentTemp);
132    pz = GetThermalPz(theTarget.GetMass(), currentTemp);
133    theTarget.SetMomentum(px, py, pz);
134    G4double tMom = std::sqrt(px*px+py*py+pz*pz);
135    G4double tEtot = std::sqrt((tMom+theTarget.GetMass())*
136                          (tMom+theTarget.GetMass())-
137                          2.*tMom*theTarget.GetMass());
138    if(1-tEtot/theTarget.GetMass()>0.001)
139    {
140      theTarget.SetTotalEnergy(tEtot);
141    }
142    else
143    {
144      theTarget.SetKineticEnergy(tMom*tMom/(2.*theTarget.GetMass()));
145    }   
146    return theTarget;
147  }
148 
149 void
150  G4Nucleus::ChooseParameters( const G4Material *aMaterial )
151  {
152    G4double random = G4UniformRand();
153    G4double sum = aMaterial->GetTotNbOfAtomsPerVolume();
154    const G4ElementVector *theElementVector = aMaterial->GetElementVector();
155    G4double running(0);
156    G4Element* element(0);
157    for(unsigned int i=0; i<aMaterial->GetNumberOfElements(); ++i )
158    {
159      running += aMaterial->GetVecNbOfAtomsPerVolume()[i];
160      if( running > random*sum ) {
161        element=(*theElementVector)[i];
162        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);   
183    }
184  }
185 
186 void
187  G4Nucleus::SetParameters( const G4double A, const G4double Z )
188  {
189    theZ = G4int(Z + 0.5);
190    theA = G4int(A + 0.5);   
191    if( theA<1 || theZ<0 || theZ>theA )
192    {
193      throw G4HadronicException(__FILE__, __LINE__,
194                               "G4Nucleus::SetParameters called with non-physical parameters");
195    }
196    aEff = A;  // atomic weight
197    zEff = Z;  // atomic number
198  }
199
200 void
201  G4Nucleus::SetParameters( const G4int A, const G4int Z )
202  {
203    theZ = Z;
204    theA = A;   
205    if( theA<1 || theZ<0 || theZ>theA )
206    {
207      throw G4HadronicException(__FILE__, __LINE__,
208                               "G4Nucleus::SetParameters called with non-physical parameters");
209    }
210    aEff = A;  // atomic weight
211    zEff = Z;  // atomic number
212  }
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  }
226 
227 G4double
228  G4Nucleus::AtomicMass( const G4double A, const G4double Z ) const
229  {
230    // Now returns (atomic mass - electron masses)
231    return G4NucleiProperties::GetNuclearMass(A, Z);
232  }
233 
234 G4double
235  G4Nucleus::AtomicMass( const G4int A, const G4int Z ) const
236  {
237    // Now returns (atomic mass - electron masses)
238    return G4NucleiProperties::GetNuclearMass(A, Z);
239  }
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  }
250 
251 G4double
252  G4Nucleus::EvaporationEffects( G4double kineticEnergy )
253  {
254    // derived from original FORTRAN code EXNU by H. Fesefeldt (10-Dec-1986)
255    //
256    // Nuclear evaporation as function of atomic number
257    // and kinetic energy (MeV) of primary particle
258    //
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  }
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 );
331
332    pnBlackTrackEnergyfromAnnihilation = exnu*fpdiv;
333    dtaBlackTrackEnergyfromAnnihilation = exnu*(1.0-fpdiv);
334   
335    G4double ran1 = -6.0;
336    G4double ran2 = -6.0;
337    for( G4int i=0; i<12; ++i ) {
338      ran1 += G4UniformRand();
339      ran2 += G4UniformRand();
340    }
341    pnBlackTrackEnergyfromAnnihilation *= 1.0 + ran1*gfa;
342    dtaBlackTrackEnergyfromAnnihilation *= 1.0 + ran2*gfa;
343
344    pnBlackTrackEnergyfromAnnihilation = std::max( 0.0, pnBlackTrackEnergyfromAnnihilation);
345    dtaBlackTrackEnergyfromAnnihilation = std::max( 0.0, dtaBlackTrackEnergyfromAnnihilation);
346    G4double blackSum = pnBlackTrackEnergyfromAnnihilation+dtaBlackTrackEnergyfromAnnihilation;
347    if (blackSum >= ekOrg/GeV) {
348      pnBlackTrackEnergyfromAnnihilation *= ekOrg/GeV/blackSum;
349      dtaBlackTrackEnergyfromAnnihilation *= ekOrg/GeV/blackSum;
350    }
351
352    return (pnBlackTrackEnergyfromAnnihilation+dtaBlackTrackEnergyfromAnnihilation)*GeV;
353  }
354 
355 G4double
356  G4Nucleus::Cinema( G4double kineticEnergy )
357  {
358    // derived from original FORTRAN code CINEMA by H. Fesefeldt (14-Oct-1987)
359    //
360    // input: kineticEnergy (MeV)
361    // returns modified kinetic energy (MeV)
362    //
363    static const G4double expxu =  82.;           // upper bound for arg. of exp
364    static const G4double expxl = -expxu;         // lower bound for arg. of exp
365   
366    G4double ek = kineticEnergy/GeV;
367    G4double ekLog = std::log( ek );
368    G4double aLog = std::log( aEff );
369    G4double em = std::min( 1.0, 0.2390 + 0.0408*aLog*aLog );
370    G4double temp1 = -ek * std::min( 0.15, 0.0019*aLog*aLog*aLog );
371    G4double temp2 = std::exp( std::max( expxl, std::min( expxu, -(ekLog-em)*(ekLog-em)*2.0 ) ) );
372    G4double result = 0.0;
373    if( std::abs( temp1 ) < 1.0 )
374    {
375      if( temp2 > 1.0e-10 )result = temp1*temp2;
376    }
377    else result = temp1*temp2;
378    if( result < -ek )result = -ek;
379    return result*GeV;
380  }
381
382 //
383 // methods for class G4Nucleus  ... by Christian Volcker
384 //
385
386 G4ThreeVector G4Nucleus::GetFermiMomentum()
387  {
388    // chv: .. we assume zero temperature!
389   
390    // momentum is equally distributed in each phasespace volume dpx, dpy, dpz.
391    G4double ranflat1=
392      CLHEP::RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);   
393    G4double ranflat2=
394      CLHEP::RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);   
395    G4double ranflat3=
396      CLHEP::RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);   
397    G4double ranmax = (ranflat1>ranflat2? ranflat1: ranflat2);
398    ranmax = (ranmax>ranflat3? ranmax : ranflat3);
399   
400    // Isotropic momentum distribution
401    G4double costheta = 2.*G4UniformRand() - 1.0;
402    G4double sintheta = std::sqrt(1.0 - costheta*costheta);
403    G4double phi = 2.0*pi*G4UniformRand();
404   
405    G4double pz=costheta*ranmax;
406    G4double px=sintheta*std::cos(phi)*ranmax;
407    G4double py=sintheta*std::sin(phi)*ranmax;
408    G4ThreeVector p(px,py,pz);
409    return p;
410  }
411 
412 G4ReactionProductVector* G4Nucleus::Fragmentate()
413  {
414    // needs implementation!
415    return NULL;
416  }
417 
418 void G4Nucleus::AddMomentum(const G4ThreeVector aMomentum)
419  {
420    momentum+=(aMomentum);
421  }
422 
423 void G4Nucleus::AddExcitationEnergy( G4double anEnergy )
424  {
425    excitationEnergy+=anEnergy;
426  }
427
428 /* end of file */
429
Note: See TracBrowser for help on using the repository browser.