Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LowEnergyGammaConversion.cc
Go to the documentation of this file.
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 // --------------------------------------------------------------------
28 // $Id$
29 // GEANT4 tag $Name: geant4-09-01-ref-00 $
30 //
31 //
32 // --------------------------------------------------------------
33 //
34 // Author: A. Forti
35 // Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
36 //
37 // History:
38 // --------
39 // 02/03/1999 A. Forti 1st implementation
40 // 14.03.2000 Veronique Lefebure;
41 // Change initialisation of lowestEnergyLimit from 1.22 to 1.022.
42 // Note that the hard coded value 1.022 should be used instead of
43 // 2*electron_mass_c2 in order to agree with the value of the data bank EPDL97
44 // 24.04.01 V.Ivanchenko remove RogueWave
45 // 27.07.01 F.Longo correct bug in energy distribution
46 // 21.01.03 V.Ivanchenko Cut per region
47 // 25.03.03 F.Longo fix in angular distribution of e+/e-
48 // 24.04.03 V.Ivanchenko - Cut per region mfpt
49 //
50 // --------------------------------------------------------------
51 
53 
54 #include "Randomize.hh"
55 #include "G4PhysicalConstants.hh"
56 #include "G4SystemOfUnits.hh"
57 #include "G4ParticleDefinition.hh"
58 #include "G4Track.hh"
59 #include "G4Step.hh"
60 #include "G4ForceCondition.hh"
61 #include "G4Gamma.hh"
62 #include "G4Electron.hh"
63 #include "G4DynamicParticle.hh"
64 #include "G4VParticleChange.hh"
65 #include "G4ThreeVector.hh"
66 #include "G4Positron.hh"
67 #include "G4IonisParamElm.hh"
68 #include "G4Material.hh"
71 #include "G4RDVEMDataSet.hh"
72 #include "G4RDVDataSetAlgorithm.hh"
74 #include "G4RDVRangeTest.hh"
75 #include "G4RDRangeTest.hh"
76 #include "G4MaterialCutsCouple.hh"
77 
79  : G4VDiscreteProcess(processName),
80  lowEnergyLimit(1.022000*MeV),
81  highEnergyLimit(100*GeV),
82  intrinsicLowEnergyLimit(1.022000*MeV),
83  intrinsicHighEnergyLimit(100*GeV),
84  smallEnergy(2.*MeV)
85 
86 {
87  if (lowEnergyLimit < intrinsicLowEnergyLimit ||
88  highEnergyLimit > intrinsicHighEnergyLimit)
89  {
90  G4Exception("G4LowEnergyGammaConversion::G4LowEnergyGammaConversion()",
91  "OutOfRange", FatalException,
92  "Energy limit outside intrinsic process validity range!");
93  }
94 
95  // The following pointer is owned by G4DataHandler
96 
97  crossSectionHandler = new G4RDCrossSectionHandler();
98  crossSectionHandler->Initialise(0,1.0220*MeV,100.*GeV,400);
99  meanFreePathTable = 0;
100  rangeTest = new G4RDRangeTest;
101 
102  if (verboseLevel > 0)
103  {
104  G4cout << GetProcessName() << " is created " << G4endl
105  << "Energy range: "
106  << lowEnergyLimit / MeV << " MeV - "
107  << highEnergyLimit / GeV << " GeV"
108  << G4endl;
109  }
110 }
111 
113 {
114  delete meanFreePathTable;
115  delete crossSectionHandler;
116  delete rangeTest;
117 }
118 
120 {
121 
122  crossSectionHandler->Clear();
123  G4String crossSectionFile = "pair/pp-cs-";
124  crossSectionHandler->LoadData(crossSectionFile);
125 
126  delete meanFreePathTable;
127  meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
128 }
129 
131  const G4Step& aStep)
132 {
133 // The energies of the e+ e- secondaries are sampled using the Bethe - Heitler
134 // cross sections with Coulomb correction. A modified version of the random
135 // number techniques of Butcher & Messel is used (Nuc Phys 20(1960),15).
136 
137 // Note 1 : Effects due to the breakdown of the Born approximation at low
138 // energy are ignored.
139 // Note 2 : The differential cross section implicitly takes account of
140 // pair creation in both nuclear and atomic electron fields. However triplet
141 // prodution is not generated.
142 
143  aParticleChange.Initialize(aTrack);
144 
145  const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
146 
147  const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
148  G4double photonEnergy = incidentPhoton->GetKineticEnergy();
149  G4ParticleMomentum photonDirection = incidentPhoton->GetMomentumDirection();
150 
151  G4double epsilon ;
152  G4double epsilon0 = electron_mass_c2 / photonEnergy ;
153 
154  // Do it fast if photon energy < 2. MeV
155  if (photonEnergy < smallEnergy )
156  {
157  epsilon = epsilon0 + (0.5 - epsilon0) * G4UniformRand();
158  }
159  else
160  {
161  // Select randomly one element in the current material
162  const G4Element* element = crossSectionHandler->SelectRandomElement(couple,photonEnergy);
163 
164  if (element == 0)
165  {
166  G4cout << "G4LowEnergyGammaConversion::PostStepDoIt - element = 0" << G4endl;
167  }
168  G4IonisParamElm* ionisation = element->GetIonisation();
169  if (ionisation == 0)
170  {
171  G4cout << "G4LowEnergyGammaConversion::PostStepDoIt - ionisation = 0" << G4endl;
172  }
173 
174  // Extract Coulomb factor for this Element
175  G4double fZ = 8. * (ionisation->GetlogZ3());
176  if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb());
177 
178  // Limits of the screening variable
179  G4double screenFactor = 136. * epsilon0 / (element->GetIonisation()->GetZ3()) ;
180  G4double screenMax = std::exp ((42.24 - fZ)/8.368) - 0.952 ;
181  G4double screenMin = std::min(4.*screenFactor,screenMax) ;
182 
183  // Limits of the energy sampling
184  G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ;
185  G4double epsilonMin = std::max(epsilon0,epsilon1);
186  G4double epsilonRange = 0.5 - epsilonMin ;
187 
188  // Sample the energy rate of the created electron (or positron)
189  G4double screen;
190  G4double gReject ;
191 
192  G4double f10 = ScreenFunction1(screenMin) - fZ;
193  G4double f20 = ScreenFunction2(screenMin) - fZ;
194  G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
195  G4double normF2 = std::max(1.5 * f20,0.);
196 
197  do {
198  if (normF1 / (normF1 + normF2) > G4UniformRand() )
199  {
200  epsilon = 0.5 - epsilonRange * std::pow(G4UniformRand(), 0.3333) ;
201  screen = screenFactor / (epsilon * (1. - epsilon));
202  gReject = (ScreenFunction1(screen) - fZ) / f10 ;
203  }
204  else
205  {
206  epsilon = epsilonMin + epsilonRange * G4UniformRand();
207  screen = screenFactor / (epsilon * (1 - epsilon));
208  gReject = (ScreenFunction2(screen) - fZ) / f20 ;
209  }
210  } while ( gReject < G4UniformRand() );
211 
212  } // End of epsilon sampling
213 
214  // Fix charges randomly
215 
216  G4double electronTotEnergy;
217  G4double positronTotEnergy;
218 
220  {
221  electronTotEnergy = (1. - epsilon) * photonEnergy;
222  positronTotEnergy = epsilon * photonEnergy;
223  }
224  else
225  {
226  positronTotEnergy = (1. - epsilon) * photonEnergy;
227  electronTotEnergy = epsilon * photonEnergy;
228  }
229 
230  // Scattered electron (positron) angles. ( Z - axis along the parent photon)
231  // Universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211),
232  // derived from Tsai distribution (Rev. Mod. Phys. 49, 421 (1977)
233 
234  G4double u;
235  const G4double a1 = 0.625;
236  G4double a2 = 3. * a1;
237  // G4double d = 27. ;
238 
239  // if (9. / (9. + d) > G4UniformRand())
240  if (0.25 > G4UniformRand())
241  {
242  u = - std::log(G4UniformRand() * G4UniformRand()) / a1 ;
243  }
244  else
245  {
246  u = - std::log(G4UniformRand() * G4UniformRand()) / a2 ;
247  }
248 
249  G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
250  G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
251  G4double phi = twopi * G4UniformRand();
252 
253  G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
254  G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
255 
256 
257  // Kinematics of the created pair:
258  // the electron and positron are assumed to have a symetric angular
259  // distribution with respect to the Z axis along the parent photon
260 
261  G4double localEnergyDeposit = 0. ;
262 
264  G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
265 
266  // Generate the electron only if with large enough range w.r.t. cuts and safety
267 
268  G4double safety = aStep.GetPostStepPoint()->GetSafety();
269 
270  if (rangeTest->Escape(G4Electron::Electron(),couple,electronKineEnergy,safety))
271  {
272  G4ThreeVector electronDirection (dxEle, dyEle, dzEle);
273  electronDirection.rotateUz(photonDirection);
274 
276  electronDirection,
277  electronKineEnergy);
278  aParticleChange.AddSecondary(particle1) ;
279  }
280  else
281  {
282  localEnergyDeposit += electronKineEnergy ;
283  }
284 
285  // The e+ is always created (even with kinetic energy = 0) for further annihilation
286  G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
287 
288  // Is the local energy deposit correct, if the positron is always created?
289  if (! (rangeTest->Escape(G4Positron::Positron(),couple,positronKineEnergy,safety)))
290  {
291  localEnergyDeposit += positronKineEnergy ;
292  positronKineEnergy = 0. ;
293  }
294 
295  G4ThreeVector positronDirection (dxPos, dyPos, dzPos);
296  positronDirection.rotateUz(photonDirection);
297 
298  // Create G4DynamicParticle object for the particle2
300  positronDirection, positronKineEnergy);
301  aParticleChange.AddSecondary(particle2) ;
302 
303  aParticleChange.ProposeLocalEnergyDeposit(localEnergyDeposit) ;
304 
305  // Kill the incident photon
309 
310  // Reset NbOfInteractionLengthLeft and return aParticleChange
311  return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
312 }
313 
315 {
316  return ( &particle == G4Gamma::Gamma() );
317 }
318 
320  G4double, // previousStepSize
322 {
324  G4double energy = photon->GetKineticEnergy();
325  const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
326  size_t materialIndex = couple->GetIndex();
327 
328  G4double meanFreePath;
329  if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
330  else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
331  else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
332  return meanFreePath;
333 }
334 
335 G4double G4LowEnergyGammaConversion::ScreenFunction1(G4double screenVariable)
336 {
337  // Compute the value of the screening function 3*phi1 - phi2
338 
339  G4double value;
340 
341  if (screenVariable > 1.)
342  value = 42.24 - 8.368 * std::log(screenVariable + 0.952);
343  else
344  value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
345 
346  return value;
347 }
348 
349 G4double G4LowEnergyGammaConversion::ScreenFunction2(G4double screenVariable)
350 {
351  // Compute the value of the screening function 1.5*phi1 - 0.5*phi2
352 
353  G4double value;
354 
355  if (screenVariable > 1.)
356  value = 42.24 - 8.368 * std::log(screenVariable + 0.952);
357  else
358  value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
359 
360  return value;
361 }