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

Last change on this file since 1199 was 819, checked in by garnier, 17 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.