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