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

Last change on this file since 1350 was 1347, checked in by garnier, 15 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.