source: trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeGammaConversion.cc@ 1199

Last change on this file since 1199 was 1055, checked in by garnier, 17 years ago

maj sur la beta de geant 4.9.3
File size: 13.5 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//
29//
30// --------------------------------------------------------------
31//
32// Author: L.Pandola
33// History:
34// --------
35// 02 Dec 2002 L.Pandola 1st implementation
36// 12 Feb 2003 MG Pia Migration to "cuts per region"
37// 10 Mar 2003 V.Ivanchenko Remove CutPerMaterial warning
38// 13 Mar 2003 L.Pandola Code "cleaned"
39// 25 Mar 2003 L.Pandola Changed the name of the database file to read
40// 24 Apr 2003 V.Ivanchenko Cut per region mfpt
41// 17 Mar 2004 L.Pandola Removed unnecessary calls to std::pow(a,b)
42// --------------------------------------------------------------
43
44#include "G4PenelopeGammaConversion.hh"
45
46#include "Randomize.hh"
47#include "G4ParticleDefinition.hh"
48#include "G4Track.hh"
49#include "G4Step.hh"
50#include "G4ForceCondition.hh"
51#include "G4Gamma.hh"
52#include "G4Electron.hh"
53#include "G4DynamicParticle.hh"
54#include "G4VParticleChange.hh"
55#include "G4ThreeVector.hh"
56#include "G4Positron.hh"
57#include "G4IonisParamElm.hh"
58#include "G4Material.hh"
59#include "G4VCrossSectionHandler.hh"
60#include "G4CrossSectionHandler.hh"
61#include "G4VEMDataSet.hh"
62#include "G4VDataSetAlgorithm.hh"
63#include "G4LogLogInterpolation.hh"
64#include "G4VRangeTest.hh"
65#include "G4RangeTest.hh"
66#include "G4MaterialCutsCouple.hh"
67
68
69G4PenelopeGammaConversion::G4PenelopeGammaConversion(const G4String& processName)
70 : G4VDiscreteProcess(processName),
71 lowEnergyLimit(1.022000*MeV),
72 highEnergyLimit(100*GeV),
73 intrinsicLowEnergyLimit(1.022000*MeV),
74 intrinsicHighEnergyLimit(100*GeV),
75 smallEnergy(1.1*MeV)
76
77{
78 if (lowEnergyLimit < intrinsicLowEnergyLimit ||
79 highEnergyLimit > intrinsicHighEnergyLimit)
80 {
81 G4Exception("G4PenelopeGammaConversion::G4PenelopeGammaConversion - energy limit outside intrinsic process validity range");
82 }
83
84 // The following pointer is owned by G4DataHandler
85 crossSectionHandler = new G4CrossSectionHandler();
86 // Log log interpolation (default)
87 crossSectionHandler->Initialise(0,1.0220*MeV,100.*GeV,400);
88 meanFreePathTable = 0;
89 rangeTest = new G4RangeTest;
90
91 if (verboseLevel > 0)
92 {
93 G4cout << GetProcessName() << " is created " << G4endl
94 << "Energy range: "
95 << lowEnergyLimit / MeV << " MeV - "
96 << highEnergyLimit / GeV << " GeV"
97 << G4endl;
98 }
99
100 G4cout << G4endl;
101 G4cout << "*******************************************************************************" << G4endl;
102 G4cout << "*******************************************************************************" << G4endl;
103 G4cout << " The class G4PenelopeGammaConversion is NOT SUPPORTED ANYMORE. " << G4endl;
104 G4cout << " It will be REMOVED with the next major release of Geant4. " << G4endl;
105 G4cout << " Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
106 G4cout << "*******************************************************************************" << G4endl;
107 G4cout << "*******************************************************************************" << G4endl;
108 G4cout << G4endl;
109}
110
111G4PenelopeGammaConversion::~G4PenelopeGammaConversion()
112{
113 delete meanFreePathTable;
114 delete crossSectionHandler;
115 delete rangeTest;
116}
117
118void G4PenelopeGammaConversion::BuildPhysicsTable(const G4ParticleDefinition& )
119{
120
121 crossSectionHandler->Clear();
122 G4String crossSectionFile = "penelope/pp-cs-pen-";
123 crossSectionHandler->LoadData(crossSectionFile);
124 delete meanFreePathTable;
125 meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
126}
127
128G4VParticleChange* G4PenelopeGammaConversion::PostStepDoIt(const G4Track& aTrack,
129 const G4Step& aStep)
130{
131 aParticleChange.Initialize(aTrack);
132
133 const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
134
135 const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
136 G4double photonEnergy = incidentPhoton->GetKineticEnergy();
137 G4ParticleMomentum photonDirection = incidentPhoton->GetMomentumDirection();
138
139 G4double eps ;
140 G4double eki = electron_mass_c2 / photonEnergy ;
141
142 // Do it fast if photon energy < 1.1 MeV
143 if (photonEnergy < smallEnergy )
144 {
145 eps = eki + (1-2*eki) * G4UniformRand();
146 }
147 else
148 {
149 // Select randomly one element in the current material
150 const G4Element* element = crossSectionHandler->SelectRandomElement(couple,photonEnergy);
151
152 if (element == 0)
153 {
154 G4cout << "G4PenelopeGammaConversion::PostStepDoIt - element = 0" << G4endl;
155 }
156 G4IonisParamElm* ionisation = element->GetIonisation();
157 if (ionisation == 0)
158 {
159 G4cout << "G4PenelopeGammaConversion::PostStepDoIt - ionisation = 0" << G4endl;
160 }
161
162 //Low energy and Coulomb corrections
163 G4double Z=ionisation->GetZ();
164 G4double ZAlpha = Z*fine_structure_const;
165 G4double ScreenRadius = GetScreeningRadius(Z);
166 G4double funct1=0,g0=0;
167 G4double g1min=0,g2min=0;
168 funct1 = 4.0*std::log(ScreenRadius);
169 g0 = funct1-4*CoulombCorrection(ZAlpha)+LowEnergyCorrection(ZAlpha,eki);
170 G4double bmin = 2*eki*ScreenRadius;
171 g1min=g0+ScreenFunction(bmin,1);
172 g2min=g0+ScreenFunction(bmin,2);
173 G4double xr,a1,p1;
174 xr=0.5-eki;
175 a1=(2.0/3.0)*g1min*xr*xr;
176 p1=a1/(a1+g2min);
177
178 //Random sampling of eps
179 G4double rand1,rand2,rand3,b;
180 G4double g1;
181
182 do{
183 rand1 = G4UniformRand();
184 if (rand1 < p1) {
185 rand2 = 2.0*G4UniformRand()-1.0;
186 if (rand2 < 0) {
187 eps = 0.5 - xr*std::pow(std::abs(rand2),(1./3.));
188 }
189 else
190 {
191 eps = 0.5 + xr*std::pow(rand2,(1./3.));
192 }
193 b = (eki*ScreenRadius)/(2*eps*(1.0-eps));
194 g1 = g0+ScreenFunction(b,1);
195 if (g1 < 0) g1=0;
196 rand3 = G4UniformRand()*g1min;
197 }
198 else
199 {
200 eps = eki+2.0*xr*G4UniformRand();
201 b = (eki*ScreenRadius)/(2*eps*(1.0-eps));
202 g1 = g0+ScreenFunction(b,2);
203 if (g1 < 0) g1=0;
204 rand3 = G4UniformRand()*g2min;
205 }
206 } while (rand3>g1);
207 } //End of eps sampling
208
209 G4double electronTotEnergy;
210 G4double positronTotEnergy;
211
212 electronTotEnergy = eps*photonEnergy;
213 positronTotEnergy = (1.0-eps)*photonEnergy;
214
215 // Scattered electron (positron) angles. ( Z - axis along the parent photon)
216
217 //electron kinematics
218 G4double costheta_el,costheta_po;
219 G4double phi_el,phi_po;
220 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
221 costheta_el = G4UniformRand()*2.0-1.0;
222 G4double kk = std::sqrt(electronKineEnergy*(electronKineEnergy+2.*electron_mass_c2));
223 costheta_el = (costheta_el*electronTotEnergy+kk)/(electronTotEnergy+costheta_el*kk);
224 phi_el = twopi * G4UniformRand() ;
225 G4double dirX_el = std::sqrt(1.-costheta_el*costheta_el) * std::cos(phi_el);
226 G4double dirY_el = std::sqrt(1.-costheta_el*costheta_el) * std::sin(phi_el);
227 G4double dirZ_el = costheta_el;
228
229 //positron kinematics
230 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
231 costheta_po = G4UniformRand()*2.0-1.0;
232 kk = std::sqrt(positronKineEnergy*(positronKineEnergy+2.*electron_mass_c2));
233 costheta_po = (costheta_po*positronTotEnergy+kk)/(positronTotEnergy+costheta_po*kk);
234 phi_po = twopi * G4UniformRand() ;
235 G4double dirX_po = std::sqrt(1.-costheta_po*costheta_po) * std::cos(phi_po);
236 G4double dirY_po = std::sqrt(1.-costheta_po*costheta_po) * std::sin(phi_po);
237 G4double dirZ_po = costheta_po;
238
239// Kinematics of the created pair:
240// the electron and positron are assumed to have a symetric angular
241// distribution with respect to the Z axis along the parent photon
242
243 G4double localEnergyDeposit = 0. ;
244
245 aParticleChange.SetNumberOfSecondaries(2) ;
246
247
248 // Generate the electron only if with large enough range w.r.t. cuts and safety
249
250 G4double safety = aStep.GetPostStepPoint()->GetSafety();
251
252 if (rangeTest->Escape(G4Electron::Electron(),couple,electronKineEnergy,safety))
253 {
254 G4ThreeVector electronDirection ( dirX_el, dirY_el, dirZ_el);
255 electronDirection.rotateUz(photonDirection);
256 G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
257 electronDirection,
258 electronKineEnergy);
259 aParticleChange.AddSecondary(particle1) ;
260 }
261 else
262 {
263 localEnergyDeposit += electronKineEnergy ;
264 }
265
266
267 if (! (rangeTest->Escape(G4Positron::Positron(),couple,positronKineEnergy,safety)))
268 {
269 localEnergyDeposit += positronKineEnergy ;
270 positronKineEnergy = 0. ;
271 }
272 G4ThreeVector positronDirection(dirX_po,dirY_po,dirZ_po);
273 positronDirection.rotateUz(photonDirection);
274
275 // Create G4DynamicParticle object for the particle2
276 G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
277 positronDirection, positronKineEnergy);
278 aParticleChange.AddSecondary(particle2) ;
279
280 aParticleChange.ProposeLocalEnergyDeposit(localEnergyDeposit) ;
281
282 // Kill the incident photon
283 aParticleChange.ProposeMomentumDirection(0.,0.,0.) ;
284 aParticleChange.ProposeEnergy(0.) ;
285 aParticleChange.ProposeTrackStatus(fStopAndKill) ;
286
287 // Reset NbOfInteractionLengthLeft and return aParticleChange
288 return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
289}
290
291G4bool G4PenelopeGammaConversion::IsApplicable(const G4ParticleDefinition& particle)
292{
293 return ( &particle == G4Gamma::Gamma() );
294}
295
296G4double G4PenelopeGammaConversion::GetMeanFreePath(const G4Track& track,
297 G4double, // previousStepSize
298 G4ForceCondition*)
299{
300 const G4DynamicParticle* photon = track.GetDynamicParticle();
301 G4double energy = photon->GetKineticEnergy();
302 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
303 size_t materialIndex = couple->GetIndex();
304
305 G4double meanFreePath;
306 if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
307 else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
308 else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
309 return meanFreePath;
310}
311
312G4double G4PenelopeGammaConversion::ScreenFunction(G4double b,G4int icase)
313{
314 G4double bsquare=b*b;
315 G4double a0,f1,f2,g1,g2;
316 f1=2.0-2*std::log(1+bsquare);
317 f2=f1-(2.0/3.0);
318 if (b < 1.0e-10)
319 {
320 f1=f1-twopi*b;
321 }
322 else
323 {
324 a0 = 4*b*std::atan(1.0/b);
325 f1 = f1 - a0;
326 f2 = f2+2*bsquare*(4.0-a0-3*std::log((1+bsquare)/bsquare));
327 }
328 g1=0.5*(3*f1-f2);
329 g2=0.25*(3*f1+f2);
330 if (icase==1) {
331 return g1;
332 }
333 else
334 {
335 return g2;
336 }
337}
338
339G4double G4PenelopeGammaConversion::CoulombCorrection(G4double a)
340{
341 G4double fc=0;
342 G4double b[7] = {0.202059,-0.03693,0.00835,-0.00201,0.00049,-0.00012,0.00003};
343 G4double aSquared = a*a;
344 G4double aFourth = aSquared*aSquared;
345 G4double aEighth = aFourth*aFourth;
346
347 fc = ((1.0/(1.0+a*a))+b[0]+b[1]*aSquared+b[2]*aFourth+b[3]*(aSquared*aFourth)+
348 b[4]*aEighth+b[5]*(aEighth*aSquared)+b[6]*(aEighth*aFourth));
349 fc=aSquared*fc;
350 return fc;
351}
352
353G4double G4PenelopeGammaConversion::LowEnergyCorrection(G4double a,G4double eki)
354{
355 G4double f0=0,t=0;
356 G4double b[12] = {-1.744,-12.10,11.18,8.523,73.26,-41.41,-13.52,-121.1,94.41,8.946,62.05,-63.41};
357 t=std::sqrt(2.0*eki);
358 G4double tSq = t*t;
359 f0=(b[0]+b[1]*a+b[2]*a*a)*t+(b[3]+b[4]*a+b[5]*a*a)*(tSq)+(b[6]+b[7]*a+b[8]*a*a)*(tSq*t)+
360 (b[9]+b[10]*a+b[11]*a*a)*(tSq*tSq);
361 return f0;
362}
363
364G4double G4PenelopeGammaConversion::GetScreeningRadius(G4double Z)
365{
366 char* path = getenv("G4LEDATA");
367 if (!path)
368 {
369 G4String excep = "G4PenelopeGammaConversion - G4LEDATA environment variable not set!";
370 G4Exception(excep);
371 }
372 G4String pathString(path);
373 G4String pathFile = pathString + "/penelope/pp-pen.dat";
374 std::ifstream file(pathFile);
375 std::filebuf* lsdp = file.rdbuf();
376
377 if (!(lsdp->is_open()))
378 {
379 G4String excep = "G4PenelopeGammaConversion - data file " + pathFile + "not found!";
380 G4Exception(excep);
381 }
382 G4int k;
383 G4double a1,a2;
384 while(!file.eof()) {
385 file >> k >> a1 >> a2;
386 if ((G4double) k == Z)
387 {
388 return a1;
389 }
390 }
391 G4String excep = "G4PenelopeGammaConversion - Screening Radius for not found in the data file";
392 G4Exception(excep);
393 return 0;
394}
Note: See TracBrowser for help on using the repository browser.