source: trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeGammaConversionModel.cc@ 1005

Last change on this file since 1005 was 991, checked in by garnier, 17 years ago

update
File size: 17.9 KB
RevLine 
[968]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// $Id: G4PenelopeGammaConversionModel.cc,v 1.2 2008/12/04 14:09:36 pandola Exp $
[991]27// GEANT4 tag $Name: geant4-09-02 $
[968]28//
29// Author: Luciano Pandola
30//
31// History:
32// --------
33// 06 Oct 2008 L Pandola Migration from process to model
34//
35
36#include "G4PenelopeGammaConversionModel.hh"
37#include "G4ParticleDefinition.hh"
38#include "G4MaterialCutsCouple.hh"
39#include "G4ProductionCutsTable.hh"
40#include "G4DynamicParticle.hh"
41#include "G4Element.hh"
42#include "G4Gamma.hh"
43#include "G4Electron.hh"
44#include "G4Positron.hh"
45#include "G4CrossSectionHandler.hh"
46
47//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48
49
50G4PenelopeGammaConversionModel::G4PenelopeGammaConversionModel(const G4ParticleDefinition*,
51 const G4String& nam)
52 :G4VEmModel(nam),fTheScreeningRadii(0),crossSectionHandler(0),isInitialised(false)
53{
54 fIntrinsicLowEnergyLimit = 2.0*electron_mass_c2;
55 fIntrinsicHighEnergyLimit = 100.0*GeV;
56 fSmallEnergy = 1.1*MeV;
57
58 SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
59 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
60 //
61 verboseLevel= 0;
62 // Verbosity scale:
63 // 0 = nothing
64 // 1 = warning for energy non-conservation
65 // 2 = details of energy budget
66 // 3 = calculation of cross sections, file openings, sampling of atoms
67 // 4 = entering in methods
68}
69
70//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71
72G4PenelopeGammaConversionModel::~G4PenelopeGammaConversionModel()
73{
74 if (crossSectionHandler) delete crossSectionHandler;
75 if (fTheScreeningRadii) delete fTheScreeningRadii;
76}
77
78//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
79
80void G4PenelopeGammaConversionModel::Initialise(const G4ParticleDefinition*,
81 const G4DataVector& )
82{
83 if (verboseLevel > 3)
84 G4cout << "Calling G4PenelopeGammaConversionModel::Initialise()" << G4endl;
85
86 //Delete the old cross section handler, if necessary
87 if (crossSectionHandler)
88 {
89 crossSectionHandler->Clear();
90 delete crossSectionHandler;
91 }
92
93 //Check energy limits
94 if (LowEnergyLimit() < fIntrinsicLowEnergyLimit)
95 {
96 G4cout << "G4PenelopeGammaConversionModel: low energy limit increased from " <<
97 LowEnergyLimit()/eV << " eV to " << fIntrinsicLowEnergyLimit/eV << " eV" << G4endl;
98 SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
99 }
100 if (HighEnergyLimit() > fIntrinsicHighEnergyLimit)
101 {
102 G4cout << "G4PenelopeGammaConversionModel: high energy limit decreased from " <<
103 HighEnergyLimit()/GeV << " GeV to " << fIntrinsicHighEnergyLimit/GeV << " GeV" << G4endl;
104 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
105 }
106
107 //Re-initialize cross section handler
108 crossSectionHandler = new G4CrossSectionHandler();
109 crossSectionHandler->Initialise(0,LowEnergyLimit(),HighEnergyLimit(),400);
110 crossSectionHandler->Clear();
111 G4String crossSectionFile = "penelope/pp-cs-pen-";
112 crossSectionHandler->LoadData(crossSectionFile);
113 //This is used to retrieve cross section values later on
114 crossSectionHandler->BuildMeanFreePathForMaterials();
115
116 if (verboseLevel > 2)
117 G4cout << "Loaded cross section files for PenelopeGammaConversion" << G4endl;
118
119 G4cout << "Penelope Gamma Conversion model is initialized " << G4endl
120 << "Energy range: "
121 << LowEnergyLimit() / MeV << " MeV - "
122 << HighEnergyLimit() / GeV << " GeV"
123 << G4endl;
124
125 if(isInitialised) return;
126
127 if(pParticleChange)
128 fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
129 else
130 fParticleChange = new G4ParticleChangeForGamma();
131 isInitialised = true;
132}
133
134//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
135
136G4double G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom(
137 const G4ParticleDefinition*,
138 G4double energy,
139 G4double Z, G4double,
140 G4double, G4double)
141{
142 //
143 // Penelope model.
144 // Cross section (including triplet production) read from database and managed
145 // through the G4CrossSectionHandler utility. Cross section data are from
146 // M.J. Berger and J.H. Hubbel (XCOM), Report NBSIR 887-3598
147 //
148
149 if (verboseLevel > 3)
150 G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopePhotoElectricModel" << G4endl;
151
152 G4int iZ = (G4int) Z;
153 if (!crossSectionHandler)
154 {
155 G4cout << "G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom" << G4endl;
156 G4cout << "The cross section handler is not correctly initialized" << G4endl;
157 G4Exception();
158 }
159 G4double cs = crossSectionHandler->FindValue(iZ,energy);
160
161 if (verboseLevel > 2)
162 G4cout << "Gamma conversion cross section at " << energy/MeV << " MeV for Z=" << Z <<
163 " = " << cs/barn << " barn" << G4endl;
164 return cs;
165}
166
167//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
168
169void G4PenelopeGammaConversionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
170 const G4MaterialCutsCouple* couple,
171 const G4DynamicParticle* aDynamicGamma,
172 G4double,
173 G4double)
174{
175 //
176 // Penelope model.
177 // Final state is sampled according to the Bethe-Heitler model with Coulomb
178 // corrections, according to the semi-empirical model of
179 // J. Baro' et al., Radiat. Phys. Chem. 44 (1994) 531.
180 //
181 // The model uses the high energy Coulomb correction from
182 // H. Davies et al., Phys. Rev. 93 (1954) 788
183 // and atomic screening radii tabulated from
184 // J.H. Hubbel et al., J. Phys. Chem. Ref. Data 9 (1980) 1023
185 // for Z= 1 to 92. This managed in this model by the method
186 // GetScreeningRadius().
187 //
188 if (verboseLevel > 3)
189 G4cout << "Calling SamplingSecondaries() of G4PenelopeGammaConversionModel" << G4endl;
190
191 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
192
193 if (photonEnergy <= LowEnergyLimit())
194 {
195 fParticleChange->ProposeTrackStatus(fStopAndKill);
196 fParticleChange->SetProposedKineticEnergy(0.);
197 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
198 return ;
199 }
200
201 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
202
203 G4double eps ;
204 G4double eki = electron_mass_c2 / photonEnergy ;
205
206 // Do it fast if photon energy < 1.1 MeV
207 if (photonEnergy < fSmallEnergy )
208 {
209 eps = eki + (1-2*eki) * G4UniformRand();
210 }
211 else
212 {
213 // Select randomly one element in the current material
214 if (verboseLevel > 2)
215 G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl;
216 //use crossSectionHandler instead of G4EmElementSelector because in this case
217 //the dimension of the table is equal to the dimension of the database
218 //(less interpolation errors)
219 G4int Z_int = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
220 if (verboseLevel > 2)
221 G4cout << "Selected Z = " << Z_int << G4endl;
222
223 //Low energy and Coulomb corrections
224 G4double Z=(G4double) Z_int;
225 G4double ZAlpha = Z*fine_structure_const;
226 G4double ScreenRadius = GetScreeningRadius(Z);
227 G4double funct1=0,g0=0;
228 G4double g1min=0,g2min=0;
229 funct1 = 4.0*std::log(ScreenRadius);
230 g0 = funct1-4*CoulombCorrection(ZAlpha)+LowEnergyCorrection(ZAlpha,eki);
231 G4double bmin = 2*eki*ScreenRadius;
232 std::vector<G4double> ScreenFunctionValues = ScreenFunction(bmin);
233 if (ScreenFunctionValues.size() != 2)
234 {
235 G4cout << "G4PenelopeGammaConversionModel::SampleSecondaries" << G4endl;
236 G4cout << "ScreenFunction did not return 2 values! Something wrong! " << G4endl;
237 G4Exception();
238 }
239 g1min=g0+ScreenFunctionValues[0];
240 g2min=g0+ScreenFunctionValues[1];
241 G4double xr,a1,p1;
242 xr=0.5-eki;
243 a1=(2.0/3.0)*g1min*xr*xr;
244 p1=a1/(a1+g2min);
245
246 //Random sampling of eps
247 G4double rand1,rand2,rand3,b;
248 G4double g1;
249
250 do{
251 rand1 = G4UniformRand();
252 if (rand1 < p1) {
253 rand2 = 2.0*G4UniformRand()-1.0;
254 if (rand2 < 0) {
255 eps = 0.5 - xr*std::pow(std::abs(rand2),(1./3.));
256 }
257 else
258 {
259 eps = 0.5 + xr*std::pow(rand2,(1./3.));
260 }
261 b = (eki*ScreenRadius)/(2*eps*(1.0-eps));
262 std::vector<G4double> ScreenFunctionSampling = ScreenFunction(b);
263 g1 = g0+ScreenFunctionSampling[0];
264 if (g1 < 0) g1=0;
265 rand3 = G4UniformRand()*g1min;
266 }
267 else
268 {
269 eps = eki+2.0*xr*G4UniformRand();
270 b = (eki*ScreenRadius)/(2*eps*(1.0-eps));
271 std::vector<G4double> ScreenFunctionSampling = ScreenFunction(b);
272 g1 = g0+ScreenFunctionSampling[1];
273 if (g1 < 0) g1=0;
274 rand3 = G4UniformRand()*g2min;
275 }
276 } while (rand3>g1);
277 } //End of eps sampling
278
279 G4double electronTotEnergy = eps*photonEnergy;
280 G4double positronTotEnergy = (1.0-eps)*photonEnergy;
281
282 // Scattered electron (positron) angles. ( Z - axis along the parent photon)
283
284 //electron kinematics
285 G4double costheta_el,costheta_po;
286 G4double phi_el,phi_po;
287 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
288 costheta_el = G4UniformRand()*2.0-1.0;
289 G4double kk = std::sqrt(electronKineEnergy*(electronKineEnergy+2.*electron_mass_c2));
290 costheta_el = (costheta_el*electronTotEnergy+kk)/(electronTotEnergy+costheta_el*kk);
291 phi_el = twopi * G4UniformRand() ;
292 G4double dirX_el = std::sqrt(1.-costheta_el*costheta_el) * std::cos(phi_el);
293 G4double dirY_el = std::sqrt(1.-costheta_el*costheta_el) * std::sin(phi_el);
294 G4double dirZ_el = costheta_el;
295
296 //positron kinematics
297 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
298 costheta_po = G4UniformRand()*2.0-1.0;
299 kk = std::sqrt(positronKineEnergy*(positronKineEnergy+2.*electron_mass_c2));
300 costheta_po = (costheta_po*positronTotEnergy+kk)/(positronTotEnergy+costheta_po*kk);
301 phi_po = twopi * G4UniformRand() ;
302 G4double dirX_po = std::sqrt(1.-costheta_po*costheta_po) * std::cos(phi_po);
303 G4double dirY_po = std::sqrt(1.-costheta_po*costheta_po) * std::sin(phi_po);
304 G4double dirZ_po = costheta_po;
305
306 // Kinematics of the created pair:
307 // the electron and positron are assumed to have a symetric angular
308 // distribution with respect to the Z axis along the parent photon
309 G4double localEnergyDeposit = 0. ;
310
311 //Generate explicitely the electron in the pair, only if it is > threshold
312 const G4ProductionCutsTable* theCoupleTable=
313 G4ProductionCutsTable::GetProductionCutsTable();
314 size_t indx = couple->GetIndex();
315 G4double cutE = (*(theCoupleTable->GetEnergyCutsVector(1)))[indx];
316 //G4double cutP = (*(theCoupleTable->GetEnergyCutsVector(2)))[indx];
317
318 if (electronKineEnergy > cutE)
319 {
320 G4ThreeVector electronDirection ( dirX_el, dirY_el, dirZ_el);
321 electronDirection.rotateUz(photonDirection);
322 G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
323 electronDirection,
324 electronKineEnergy);
325 fvect->push_back(electron);
326 }
327 else
328 {
329 localEnergyDeposit += electronKineEnergy;
330 electronKineEnergy = 0;
331 }
332
333 //Generate the positron. Real particle in any case, because it will annihilate. If below
334 //threshold, produce it at rest
335 if (positronKineEnergy < cutE)
336 {
337 localEnergyDeposit += positronKineEnergy;
338 positronKineEnergy = 0; //produce it at rest
339 }
340 G4ThreeVector positronDirection(dirX_po,dirY_po,dirZ_po);
341 positronDirection.rotateUz(photonDirection);
342 G4DynamicParticle* positron = new G4DynamicParticle(G4Positron::Positron(),
343 positronDirection, positronKineEnergy);
344 fvect->push_back(positron);
345
346 //Update the status of the primary gamma (kill it)
347 fParticleChange->SetProposedKineticEnergy(0.);
348 fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
349 fParticleChange->ProposeTrackStatus(fStopAndKill);
350
351
352 if (verboseLevel > 1)
353 {
354 G4cout << "-----------------------------------------------------------" << G4endl;
355 G4cout << "Energy balance from G4PenelopeGammaConversion" << G4endl;
356 G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl;
357 G4cout << "-----------------------------------------------------------" << G4endl;
358 if (electronKineEnergy)
359 G4cout << "Electron (explicitely produced) " << electronKineEnergy/keV << " keV" << G4endl;
360 if (positronKineEnergy)
361 G4cout << "Positron (not at rest) " << positronKineEnergy/keV << " keV" << G4endl;
362 G4cout << "Rest masses of e+/- " << 2.0*electron_mass_c2/keV << " keV" << G4endl;
363 if (localEnergyDeposit)
364 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
365 G4cout << "Total final state: " << (electronKineEnergy+positronKineEnergy+
366 localEnergyDeposit+2.0*electron_mass_c2)/keV <<
367 " keV" << G4endl;
368 G4cout << "-----------------------------------------------------------" << G4endl;
369 }
370 if (verboseLevel > 0)
371 {
372 G4double energyDiff = std::fabs(electronKineEnergy+positronKineEnergy+
373 localEnergyDeposit+2.0*electron_mass_c2-photonEnergy);
374 if (energyDiff > 0.05*keV)
375 G4cout << "Warning from G4PenelopeGammaConversion: problem with energy conservation: " <<
376 (electronKineEnergy+positronKineEnergy+
377 localEnergyDeposit+2.0*electron_mass_c2)/keV << " keV (final) vs. " <<
378 photonEnergy/keV << " keV (initial)" << G4endl;
379 }
380
381}
382
383//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
384
385std::vector<G4double> G4PenelopeGammaConversionModel::ScreenFunction(G4double b)
386{
387 std::vector<G4double> result;
388 result.clear();
389 G4double bsquare=b*b;
390 G4double a0,f1,f2;
391 f1=2.0-2*std::log(1+bsquare);
392 f2=f1-(2.0/3.0);
393 if (b < 1.0e-10)
394 {
395 f1=f1-twopi*b;
396 }
397 else
398 {
399 a0 = 4*b*std::atan(1.0/b);
400 f1 = f1 - a0;
401 f2 = f2+2*bsquare*(4.0-a0-3*std::log((1+bsquare)/bsquare));
402 }
403 result.push_back(0.5*(3*f1-f2));
404 result.push_back(0.25*(3*f1+f2));
405 return result;
406}
407
408//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
409
410G4double G4PenelopeGammaConversionModel::GetScreeningRadius(G4double Z)
411{
412 G4double result = 0;
413 G4bool foundElement = false;
414 G4int iZ = (G4int) Z;
415 if (!fTheScreeningRadii)
416 fTheScreeningRadii = new std::map<G4int,G4double>;
417
418 if (fTheScreeningRadii->count(iZ))
419 {
420 //The element is already loaded: just return it
421 result = fTheScreeningRadii->find(iZ)->second;
422 return result;
423 }
424 else //retrieve all from file
425 {
426 char* path = getenv("G4LEDATA");
427 if (!path)
428 {
429 G4String excep = "G4PenelopeGammaConversionModel - G4LEDATA environment variable not set!";
430 G4Exception(excep);
431 }
432 G4String pathString(path);
433 G4String pathFile = pathString + "/penelope/pp-pen.dat";
434 std::ifstream file(pathFile);
435
436 if (!(file.is_open()))
437 {
438 G4String excep = "G4PenelopeGammaConversionModel - data file " + pathFile + "not found!";
439 G4Exception(excep);
440 }
441 G4int k;
442 G4double a1,a2;
443 while(!file.eof()) {
444 file >> k >> a1 >> a2;
445 fTheScreeningRadii->insert(std::make_pair(k,a1));
446 if ((G4double) k == Z)
447 {
448 result = a1;
449 foundElement = true;
450 }
451 }
452 file.close();
453 if (verboseLevel > 2)
454 G4cout << "Read file pp-pen.dat" << G4endl;
455 if (foundElement)
456 return result;
457 else
458 {
459 G4String excep = "G4PenelopeGammaConversionModel - Screening Radius for not found in the data file";
460 G4Exception(excep);
461 return 0;
462 }
463 }
464}
465
466//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
467
468G4double G4PenelopeGammaConversionModel::CoulombCorrection(G4double a)
469{
470 G4double fc=0;
471 G4double b[7] = {0.202059,-0.03693,0.00835,-0.00201,0.00049,-0.00012,0.00003};
472 G4double aSquared = a*a;
473 G4double aFourth = aSquared*aSquared;
474 G4double aEighth = aFourth*aFourth;
475
476 fc = ((1.0/(1.0+a*a))+b[0]+b[1]*aSquared+b[2]*aFourth+b[3]*(aSquared*aFourth)+
477 b[4]*aEighth+b[5]*(aEighth*aSquared)+b[6]*(aEighth*aFourth));
478 fc=aSquared*fc;
479 return fc;
480}
481
482//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
483
484G4double G4PenelopeGammaConversionModel::LowEnergyCorrection(G4double a,G4double eki)
485{
486 G4double f0=0,t=0;
487 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};
488 t=std::sqrt(2.0*eki);
489 G4double tSq = t*t;
490 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)+
491 (b[9]+b[10]*a+b[11]*a*a)*(tSq*tSq);
492 return f0;
493
494}
Note: See TracBrowser for help on using the repository browser.