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

Last change on this file since 1183 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 13.2 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 
41#include "G4Nucleus.hh"
42#include "G4NucleiProperties.hh"
43#include "Randomize.hh"
44#include "G4HadronicException.hh"
45 
46G4Nucleus::G4Nucleus()
47{
48  pnBlackTrackEnergy = 0.0;
49  dtaBlackTrackEnergy = 0.0;
50  pnBlackTrackEnergyfromAnnihilation = 0.0;
51  dtaBlackTrackEnergyfromAnnihilation = 0.0;
52  excitationEnergy = 0.0;
53  momentum = G4ThreeVector(0.,0.,0.);
54  fermiMomentum = 1.52*hbarc/fermi;
55  theTemp = 293.16*kelvin;
56}
57
58G4Nucleus::G4Nucleus( const G4double A, const G4double Z )
59{
60  SetParameters( A, Z );
61  pnBlackTrackEnergy = 0.0;
62  dtaBlackTrackEnergy = 0.0;
63  pnBlackTrackEnergyfromAnnihilation = 0.0;
64  dtaBlackTrackEnergyfromAnnihilation = 0.0;
65  excitationEnergy = 0.0;
66  momentum = G4ThreeVector(0.,0.,0.);
67  fermiMomentum = 1.52*hbarc/fermi;
68  theTemp = 293.16*kelvin;
69}
70
71G4Nucleus::G4Nucleus( const G4Material *aMaterial )
72{
73  ChooseParameters( aMaterial );
74  pnBlackTrackEnergy = 0.0;
75  dtaBlackTrackEnergy = 0.0;
76  pnBlackTrackEnergyfromAnnihilation = 0.0;
77  dtaBlackTrackEnergyfromAnnihilation = 0.0;
78  excitationEnergy = 0.0;
79  momentum = G4ThreeVector(0.,0.,0.);
80  fermiMomentum = 1.52*hbarc/fermi;
81  theTemp = aMaterial->GetTemperature();
82}
83
84G4Nucleus::~G4Nucleus() {}
85
86G4ReactionProduct G4Nucleus::
87GetBiasedThermalNucleus(G4double aMass, G4ThreeVector aVelocity, G4double temp) const
88{
89  G4double velMag = aVelocity.mag();
90  G4ReactionProduct result;
91  G4double value = 0;
92  G4double random = 1;
93  G4double norm = 3.*std::sqrt(k_Boltzmann*temp*aMass*G4Neutron::Neutron()->GetPDGMass());
94  norm /= G4Neutron::Neutron()->GetPDGMass();
95  norm *= 5.;
96  norm += velMag;
97  norm /= velMag;
98  while(value/norm<random)
99  {
100     result = GetThermalNucleus(aMass, temp);
101     G4ThreeVector targetVelocity = 1./result.GetMass()*result.GetMomentum();
102     value = (targetVelocity+aVelocity).mag()/velMag;
103     random = G4UniformRand();
104  }
105  return result;
106}
107
108G4ReactionProduct G4Nucleus::GetThermalNucleus(G4double targetMass, G4double temp) const
109  {
110    G4double currentTemp = temp;
111    if(currentTemp < 0) currentTemp = theTemp;
112    G4ReactionProduct theTarget;   
113    theTarget.SetMass(targetMass*G4Neutron::Neutron()->GetPDGMass());
114    G4double px, py, pz;
115    px = GetThermalPz(theTarget.GetMass(), currentTemp);
116    py = GetThermalPz(theTarget.GetMass(), currentTemp);
117    pz = GetThermalPz(theTarget.GetMass(), currentTemp);
118    theTarget.SetMomentum(px, py, pz);
119    G4double tMom = std::sqrt(px*px+py*py+pz*pz);
120    G4double tEtot = std::sqrt((tMom+theTarget.GetMass())*
121                          (tMom+theTarget.GetMass())-
122                          2.*tMom*theTarget.GetMass());
123    if(1-tEtot/theTarget.GetMass()>0.001)
124    {
125      theTarget.SetTotalEnergy(tEtot);
126    }
127    else
128    {
129      theTarget.SetKineticEnergy(tMom*tMom/(2.*theTarget.GetMass()));
130    }   
131    return theTarget;
132  }
133 
134 void
135  G4Nucleus::ChooseParameters( const G4Material *aMaterial )
136  {
137    G4double random = G4UniformRand();
138    G4double sum = 0;
139    const G4ElementVector *theElementVector = aMaterial->GetElementVector();
140    unsigned int i;
141    for(i=0; i<aMaterial->GetNumberOfElements(); ++i )
142    {
143      sum += aMaterial->GetAtomicNumDensityVector()[i];
144    }
145    G4double running = 0;
146    for(i=0; i<aMaterial->GetNumberOfElements(); ++i )
147    {
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      }
154    }
155  }
156 
157 void
158  G4Nucleus::SetParameters( const G4double A, const G4double Z )
159  {
160    G4int myZ = G4int(Z + 0.5);
161    G4int myA = G4int(A + 0.5);   
162    if( myA<1 || myZ<0 || myZ>myA )
163    {
164      throw G4HadronicException(__FILE__, __LINE__,
165                               "G4Nucleus::SetParameters called with non-physical parameters");
166    }
167    aEff = A;  // atomic weight
168    zEff = Z;  // atomic number
169  }
170
171 G4DynamicParticle *
172  G4Nucleus::ReturnTargetParticle() const
173  {
174    // choose a proton or a neutron as the target particle
175   
176    G4DynamicParticle *targetParticle = new G4DynamicParticle;
177    if( G4UniformRand() < zEff/aEff )
178      targetParticle->SetDefinition( G4Proton::Proton() );
179    else
180      targetParticle->SetDefinition( G4Neutron::Neutron() );
181    return targetParticle;
182  }
183 
184 G4double
185  G4Nucleus::AtomicMass( const G4double A, const G4double Z ) const
186  {
187    // Now returns (atomic mass - electron masses)
188    return G4NucleiProperties::GetNuclearMass(A, Z);
189  }
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  }
200 
201 G4double
202  G4Nucleus::EvaporationEffects( G4double kineticEnergy )
203  {
204    // derived from original FORTRAN code EXNU by H. Fesefeldt (10-Dec-1986)
205    //
206    // Nuclear evaporation as function of atomic number
207    // and kinetic energy (MeV) of primary particle
208    //
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  }
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 );
281
282    pnBlackTrackEnergyfromAnnihilation = exnu*fpdiv;
283    dtaBlackTrackEnergyfromAnnihilation = exnu*(1.0-fpdiv);
284   
285    G4double ran1 = -6.0;
286    G4double ran2 = -6.0;
287    for( G4int i=0; i<12; ++i ) {
288      ran1 += G4UniformRand();
289      ran2 += G4UniformRand();
290    }
291    pnBlackTrackEnergyfromAnnihilation *= 1.0 + ran1*gfa;
292    dtaBlackTrackEnergyfromAnnihilation *= 1.0 + ran2*gfa;
293
294    pnBlackTrackEnergyfromAnnihilation = std::max( 0.0, pnBlackTrackEnergyfromAnnihilation);
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    }
301
302    return (pnBlackTrackEnergyfromAnnihilation+dtaBlackTrackEnergyfromAnnihilation)*GeV;
303  }
304 
305 G4double
306  G4Nucleus::Cinema( G4double kineticEnergy )
307  {
308    // derived from original FORTRAN code CINEMA by H. Fesefeldt (14-Oct-1987)
309    //
310    // input: kineticEnergy (MeV)
311    // returns modified kinetic energy (MeV)
312    //
313    static const G4double expxu =  82.;           // upper bound for arg. of exp
314    static const G4double expxl = -expxu;         // lower bound for arg. of exp
315   
316    G4double ek = kineticEnergy/GeV;
317    G4double ekLog = std::log( ek );
318    G4double aLog = std::log( aEff );
319    G4double em = std::min( 1.0, 0.2390 + 0.0408*aLog*aLog );
320    G4double temp1 = -ek * std::min( 0.15, 0.0019*aLog*aLog*aLog );
321    G4double temp2 = std::exp( std::max( expxl, std::min( expxu, -(ekLog-em)*(ekLog-em)*2.0 ) ) );
322    G4double result = 0.0;
323    if( std::abs( temp1 ) < 1.0 )
324    {
325      if( temp2 > 1.0e-10 )result = temp1*temp2;
326    }
327    else result = temp1*temp2;
328    if( result < -ek )result = -ek;
329    return result*GeV;
330  }
331
332 //
333 // methods for class G4Nucleus  ... by Christian Volcker
334 //
335
336 G4ThreeVector G4Nucleus::GetFermiMomentum()
337  {
338    // chv: .. we assume zero temperature!
339   
340    // momentum is equally distributed in each phasespace volume dpx, dpy, dpz.
341    G4double ranflat1=
342      CLHEP::RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);   
343    G4double ranflat2=
344      CLHEP::RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);   
345    G4double ranflat3=
346      CLHEP::RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);   
347    G4double ranmax = (ranflat1>ranflat2? ranflat1: ranflat2);
348    ranmax = (ranmax>ranflat3? ranmax : ranflat3);
349   
350    // Isotropic momentum distribution
351    G4double costheta = 2.*G4UniformRand() - 1.0;
352    G4double sintheta = std::sqrt(1.0 - costheta*costheta);
353    G4double phi = 2.0*pi*G4UniformRand();
354   
355    G4double pz=costheta*ranmax;
356    G4double px=sintheta*std::cos(phi)*ranmax;
357    G4double py=sintheta*std::sin(phi)*ranmax;
358    G4ThreeVector p(px,py,pz);
359    return p;
360  }
361 
362 G4ReactionProductVector* G4Nucleus::Fragmentate()
363  {
364    // needs implementation!
365    return NULL;
366  }
367 
368 void G4Nucleus::AddMomentum(const G4ThreeVector aMomentum)
369  {
370    momentum+=(aMomentum);
371  }
372 
373 void G4Nucleus::AddExcitationEnergy( G4double anEnergy )
374  {
375    excitationEnergy+=anEnergy;
376  }
377
378 /* end of file */
379
Note: See TracBrowser for help on using the repository browser.