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

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