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

Last change on this file since 1315 was 1315, checked in by garnier, 14 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 14.6 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{
49  pnBlackTrackEnergy = 0.0;
50  dtaBlackTrackEnergy = 0.0;
51  pnBlackTrackEnergyfromAnnihilation = 0.0;
52  dtaBlackTrackEnergyfromAnnihilation = 0.0;
53  excitationEnergy = 0.0;
54  momentum = G4ThreeVector(0.,0.,0.);
55  fermiMomentum = 1.52*hbarc/fermi;
56  theTemp = 293.16*kelvin;
57}
58
59G4Nucleus::G4Nucleus( const G4double A, const G4double Z )
60{
61  SetParameters( A, Z );
62  pnBlackTrackEnergy = 0.0;
63  dtaBlackTrackEnergy = 0.0;
64  pnBlackTrackEnergyfromAnnihilation = 0.0;
65  dtaBlackTrackEnergyfromAnnihilation = 0.0;
66  excitationEnergy = 0.0;
67  momentum = G4ThreeVector(0.,0.,0.);
68  fermiMomentum = 1.52*hbarc/fermi;
69  theTemp = 293.16*kelvin;
70}
71
72G4Nucleus::G4Nucleus( const G4int A, const G4int Z )
73{
74  SetParameters( A, Z );
75  pnBlackTrackEnergy = 0.0;
76  dtaBlackTrackEnergy = 0.0;
77  pnBlackTrackEnergyfromAnnihilation = 0.0;
78  dtaBlackTrackEnergyfromAnnihilation = 0.0;
79  excitationEnergy = 0.0;
80  momentum = G4ThreeVector(0.,0.,0.);
81  fermiMomentum = 1.52*hbarc/fermi;
82  theTemp = 293.16*kelvin;
83}
84
85G4Nucleus::G4Nucleus( const G4Material *aMaterial )
86{
87  ChooseParameters( aMaterial );
88  pnBlackTrackEnergy = 0.0;
89  dtaBlackTrackEnergy = 0.0;
90  pnBlackTrackEnergyfromAnnihilation = 0.0;
91  dtaBlackTrackEnergyfromAnnihilation = 0.0;
92  excitationEnergy = 0.0;
93  momentum = G4ThreeVector(0.,0.,0.);
94  fermiMomentum = 1.52*hbarc/fermi;
95  theTemp = aMaterial->GetTemperature();
96}
97
98G4Nucleus::~G4Nucleus() {}
99
100G4ReactionProduct G4Nucleus::
101GetBiasedThermalNucleus(G4double aMass, G4ThreeVector aVelocity, G4double temp) const
102{
103  G4double velMag = aVelocity.mag();
104  G4ReactionProduct result;
105  G4double value = 0;
106  G4double random = 1;
107  G4double norm = 3.*std::sqrt(k_Boltzmann*temp*aMass*G4Neutron::Neutron()->GetPDGMass());
108  norm /= G4Neutron::Neutron()->GetPDGMass();
109  norm *= 5.;
110  norm += velMag;
111  norm /= velMag;
112  while(value/norm<random)
113  {
114     result = GetThermalNucleus(aMass, temp);
115     G4ThreeVector targetVelocity = 1./result.GetMass()*result.GetMomentum();
116     value = (targetVelocity+aVelocity).mag()/velMag;
117     random = G4UniformRand();
118  }
119  return result;
120}
121
122G4ReactionProduct G4Nucleus::GetThermalNucleus(G4double targetMass, G4double temp) const
123  {
124    G4double currentTemp = temp;
125    if(currentTemp < 0) currentTemp = theTemp;
126    G4ReactionProduct theTarget;   
127    theTarget.SetMass(targetMass*G4Neutron::Neutron()->GetPDGMass());
128    G4double px, py, pz;
129    px = GetThermalPz(theTarget.GetMass(), currentTemp);
130    py = GetThermalPz(theTarget.GetMass(), currentTemp);
131    pz = GetThermalPz(theTarget.GetMass(), currentTemp);
132    theTarget.SetMomentum(px, py, pz);
133    G4double tMom = std::sqrt(px*px+py*py+pz*pz);
134    G4double tEtot = std::sqrt((tMom+theTarget.GetMass())*
135                          (tMom+theTarget.GetMass())-
136                          2.*tMom*theTarget.GetMass());
137    if(1-tEtot/theTarget.GetMass()>0.001)
138    {
139      theTarget.SetTotalEnergy(tEtot);
140    }
141    else
142    {
143      theTarget.SetKineticEnergy(tMom*tMom/(2.*theTarget.GetMass()));
144    }   
145    return theTarget;
146  }
147 
148 void
149  G4Nucleus::ChooseParameters( const G4Material *aMaterial )
150  {
151    G4double random = G4UniformRand();
152    G4double sum = aMaterial->GetTotNbOfAtomsPerVolume();
153    const G4ElementVector *theElementVector = aMaterial->GetElementVector();
154    G4double running(0);
155    G4Element* element(0);
156    for(unsigned int i=0; i<aMaterial->GetNumberOfElements(); ++i )
157    {
158      running += aMaterial->GetVecNbOfAtomsPerVolume()[i];
159      if( running > random*sum ) {
160        element=(*theElementVector)[i];
161        break;
162      }
163    }
164    if ( element->GetNumberOfIsotopes() > 0 ) {
165      G4double randomAbundance = G4UniformRand();
166      G4double sumAbundance = element->GetRelativeAbundanceVector()[0];
167      unsigned int iso=0;
168      while ( iso < element->GetNumberOfIsotopes() &&
169               sumAbundance < randomAbundance ) {
170         ++iso;
171         sumAbundance += element->GetRelativeAbundanceVector()[iso];
172      }
173      theA=element->GetIsotope(iso)->GetN();
174      theZ=element->GetIsotope(iso)->GetZ();
175      aEff=theA;
176      zEff=theZ;
177    } else {   
178      aEff = element->GetN();
179      zEff = element->GetZ();
180      theZ = G4int(zEff + 0.5);
181      theA = G4int(aEff + 0.5);   
182    }
183  }
184 
185 void
186  G4Nucleus::SetParameters( const G4double A, const G4double Z )
187  {
188    theZ = G4int(Z + 0.5);
189    theA = G4int(A + 0.5);   
190    if( theA<1 || theZ<0 || theZ>theA )
191    {
192      throw G4HadronicException(__FILE__, __LINE__,
193                               "G4Nucleus::SetParameters called with non-physical parameters");
194    }
195    aEff = A;  // atomic weight
196    zEff = Z;  // atomic number
197  }
198
199 void
200  G4Nucleus::SetParameters( const G4int A, const G4int Z )
201  {
202    theZ = Z;
203    theA = A;   
204    if( theA<1 || theZ<0 || theZ>theA )
205    {
206      throw G4HadronicException(__FILE__, __LINE__,
207                               "G4Nucleus::SetParameters called with non-physical parameters");
208    }
209    aEff = A;  // atomic weight
210    zEff = Z;  // atomic number
211  }
212
213 G4DynamicParticle *
214  G4Nucleus::ReturnTargetParticle() const
215  {
216    // choose a proton or a neutron as the target particle
217   
218    G4DynamicParticle *targetParticle = new G4DynamicParticle;
219    if( G4UniformRand() < zEff/aEff )
220      targetParticle->SetDefinition( G4Proton::Proton() );
221    else
222      targetParticle->SetDefinition( G4Neutron::Neutron() );
223    return targetParticle;
224  }
225 
226 G4double
227  G4Nucleus::AtomicMass( const G4double A, const G4double Z ) const
228  {
229    // Now returns (atomic mass - electron masses)
230    return G4NucleiProperties::GetNuclearMass(A, Z);
231  }
232 
233 G4double
234  G4Nucleus::AtomicMass( const G4int A, const G4int Z ) const
235  {
236    // Now returns (atomic mass - electron masses)
237    return G4NucleiProperties::GetNuclearMass(A, Z);
238  }
239 
240 G4double
241  G4Nucleus::GetThermalPz( const G4double mass, const G4double temp ) const
242  {
243    G4double result = G4RandGauss::shoot();
244    result *= std::sqrt(k_Boltzmann*temp*mass); // Das ist impuls (Pz),
245                                           // nichtrelativistische rechnung
246                                           // Maxwell verteilung angenommen
247    return result;
248  }
249 
250 G4double
251  G4Nucleus::EvaporationEffects( G4double kineticEnergy )
252  {
253    // derived from original FORTRAN code EXNU by H. Fesefeldt (10-Dec-1986)
254    //
255    // Nuclear evaporation as function of atomic number
256    // and kinetic energy (MeV) of primary particle
257    //
258    // returns kinetic energy (MeV)
259    //
260    if( aEff < 1.5 )
261    {
262      pnBlackTrackEnergy = dtaBlackTrackEnergy = 0.0;
263      return 0.0;
264    }
265    G4double ek = kineticEnergy/GeV;
266    G4float ekin = std::min( 4.0, std::max( 0.1, ek ) );
267    const G4float atno = std::min( 120., aEff ); 
268    const G4float gfa = 2.0*((aEff-1.0)/70.)*std::exp(-(aEff-1.0)/70.);
269    //
270    // 0.35 value at 1 GeV
271    // 0.05 value at 0.1 GeV
272    //
273    G4float cfa = std::max( 0.15, 0.35 + ((0.35-0.05)/2.3)*std::log(ekin) );
274    G4float exnu = 7.716 * cfa * std::exp(-cfa)
275      * ((atno-1.0)/120.)*std::exp(-(atno-1.0)/120.);
276    G4float fpdiv = std::max( 0.5, 1.0-0.25*ekin*ekin );
277    //
278    // pnBlackTrackEnergy  is the kinetic energy (in GeV) available for
279    //                     proton/neutron black track particles
280    // dtaBlackTrackEnergy is the kinetic energy (in GeV) available for
281    //                     deuteron/triton/alpha black track particles
282    //
283    pnBlackTrackEnergy = exnu*fpdiv;
284    dtaBlackTrackEnergy = exnu*(1.0-fpdiv);
285   
286    if( G4int(zEff+0.1) != 82 )
287    { 
288      G4double ran1 = -6.0;
289      G4double ran2 = -6.0;
290      for( G4int i=0; i<12; ++i )
291      {
292        ran1 += G4UniformRand();
293        ran2 += G4UniformRand();
294      }
295      pnBlackTrackEnergy *= 1.0 + ran1*gfa;
296      dtaBlackTrackEnergy *= 1.0 + ran2*gfa;
297    }
298    pnBlackTrackEnergy = std::max( 0.0, pnBlackTrackEnergy );
299    dtaBlackTrackEnergy = std::max( 0.0, dtaBlackTrackEnergy );
300    while( pnBlackTrackEnergy+dtaBlackTrackEnergy >= ek )
301    {
302      pnBlackTrackEnergy *= 1.0 - 0.5*G4UniformRand();
303      dtaBlackTrackEnergy *= 1.0 - 0.5*G4UniformRand();
304    }
305//    G4cout << "EvaporationEffects "<<kineticEnergy<<" "
306//           <<pnBlackTrackEnergy+dtaBlackTrackEnergy<<endl;
307    return (pnBlackTrackEnergy+dtaBlackTrackEnergy)*GeV;
308  }
309 
310 G4double G4Nucleus::AnnihilationEvaporationEffects(G4double kineticEnergy, G4double ekOrg)
311  {
312    // Nuclear evaporation as a function of atomic number and kinetic
313    // energy (MeV) of primary particle.  Modified for annihilation effects.
314    //
315    if( aEff < 1.5 || ekOrg < 0.)
316    {
317      pnBlackTrackEnergyfromAnnihilation = 0.0;
318      dtaBlackTrackEnergyfromAnnihilation = 0.0;
319      return 0.0;
320    }
321    G4double ek = kineticEnergy/GeV;
322    G4float ekin = std::min( 4.0, std::max( 0.1, ek ) );
323    const G4float atno = std::min( 120., aEff ); 
324    const G4float gfa = 2.0*((aEff-1.0)/70.)*std::exp(-(aEff-1.0)/70.);
325
326    G4float cfa = std::max( 0.15, 0.35 + ((0.35-0.05)/2.3)*std::log(ekin) );
327    G4float exnu = 7.716 * cfa * std::exp(-cfa)
328      * ((atno-1.0)/120.)*std::exp(-(atno-1.0)/120.);
329    G4float fpdiv = std::max( 0.5, 1.0-0.25*ekin*ekin );
330
331    pnBlackTrackEnergyfromAnnihilation = exnu*fpdiv;
332    dtaBlackTrackEnergyfromAnnihilation = exnu*(1.0-fpdiv);
333   
334    G4double ran1 = -6.0;
335    G4double ran2 = -6.0;
336    for( G4int i=0; i<12; ++i ) {
337      ran1 += G4UniformRand();
338      ran2 += G4UniformRand();
339    }
340    pnBlackTrackEnergyfromAnnihilation *= 1.0 + ran1*gfa;
341    dtaBlackTrackEnergyfromAnnihilation *= 1.0 + ran2*gfa;
342
343    pnBlackTrackEnergyfromAnnihilation = std::max( 0.0, pnBlackTrackEnergyfromAnnihilation);
344    dtaBlackTrackEnergyfromAnnihilation = std::max( 0.0, dtaBlackTrackEnergyfromAnnihilation);
345    G4double blackSum = pnBlackTrackEnergyfromAnnihilation+dtaBlackTrackEnergyfromAnnihilation;
346    if (blackSum >= ekOrg/GeV) {
347      pnBlackTrackEnergyfromAnnihilation *= ekOrg/GeV/blackSum;
348      dtaBlackTrackEnergyfromAnnihilation *= ekOrg/GeV/blackSum;
349    }
350
351    return (pnBlackTrackEnergyfromAnnihilation+dtaBlackTrackEnergyfromAnnihilation)*GeV;
352  }
353 
354 G4double
355  G4Nucleus::Cinema( G4double kineticEnergy )
356  {
357    // derived from original FORTRAN code CINEMA by H. Fesefeldt (14-Oct-1987)
358    //
359    // input: kineticEnergy (MeV)
360    // returns modified kinetic energy (MeV)
361    //
362    static const G4double expxu =  82.;           // upper bound for arg. of exp
363    static const G4double expxl = -expxu;         // lower bound for arg. of exp
364   
365    G4double ek = kineticEnergy/GeV;
366    G4double ekLog = std::log( ek );
367    G4double aLog = std::log( aEff );
368    G4double em = std::min( 1.0, 0.2390 + 0.0408*aLog*aLog );
369    G4double temp1 = -ek * std::min( 0.15, 0.0019*aLog*aLog*aLog );
370    G4double temp2 = std::exp( std::max( expxl, std::min( expxu, -(ekLog-em)*(ekLog-em)*2.0 ) ) );
371    G4double result = 0.0;
372    if( std::abs( temp1 ) < 1.0 )
373    {
374      if( temp2 > 1.0e-10 )result = temp1*temp2;
375    }
376    else result = temp1*temp2;
377    if( result < -ek )result = -ek;
378    return result*GeV;
379  }
380
381 //
382 // methods for class G4Nucleus  ... by Christian Volcker
383 //
384
385 G4ThreeVector G4Nucleus::GetFermiMomentum()
386  {
387    // chv: .. we assume zero temperature!
388   
389    // momentum is equally distributed in each phasespace volume dpx, dpy, dpz.
390    G4double ranflat1=
391      CLHEP::RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);   
392    G4double ranflat2=
393      CLHEP::RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);   
394    G4double ranflat3=
395      CLHEP::RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);   
396    G4double ranmax = (ranflat1>ranflat2? ranflat1: ranflat2);
397    ranmax = (ranmax>ranflat3? ranmax : ranflat3);
398   
399    // Isotropic momentum distribution
400    G4double costheta = 2.*G4UniformRand() - 1.0;
401    G4double sintheta = std::sqrt(1.0 - costheta*costheta);
402    G4double phi = 2.0*pi*G4UniformRand();
403   
404    G4double pz=costheta*ranmax;
405    G4double px=sintheta*std::cos(phi)*ranmax;
406    G4double py=sintheta*std::sin(phi)*ranmax;
407    G4ThreeVector p(px,py,pz);
408    return p;
409  }
410 
411 G4ReactionProductVector* G4Nucleus::Fragmentate()
412  {
413    // needs implementation!
414    return NULL;
415  }
416 
417 void G4Nucleus::AddMomentum(const G4ThreeVector aMomentum)
418  {
419    momentum+=(aMomentum);
420  }
421 
422 void G4Nucleus::AddExcitationEnergy( G4double anEnergy )
423  {
424    excitationEnergy+=anEnergy;
425  }
426
427 /* end of file */
428
Note: See TracBrowser for help on using the repository browser.