source: trunk/source/processes/electromagnetic/lowenergy/src/G4Penelope08GammaConversionModel.cc@ 1316

Last change on this file since 1316 was 1316, checked in by garnier, 15 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 21.4 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// $Id: G4Penelope08GammaConversionModel.cc,v 1.2 2010/04/23 14:49:46 pandola Exp $
27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
28//
29// Author: Luciano Pandola
30//
31// History:
32// --------
33// 13 Jan 2010 L Pandola First implementation (updated to Penelope08)
34//
35
36#include "G4Penelope08GammaConversionModel.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 "G4PhysicsFreeVector.hh"
46
47//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48
49
50G4Penelope08GammaConversionModel::G4Penelope08GammaConversionModel(const G4ParticleDefinition*,
51 const G4String& nam)
52 :G4VEmModel(nam),logAtomicCrossSection(0),fEffectiveCharge(0),fMaterialInvScreeningRadius(0),
53 fScreeningFunction(0),isInitialised(false)
54{
55 fIntrinsicLowEnergyLimit = 2.0*electron_mass_c2;
56 fIntrinsicHighEnergyLimit = 100.0*GeV;
57 fSmallEnergy = 1.1*MeV;
58 InitializeScreeningRadii();
59
60 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
61 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
62 //
63 verboseLevel= 0;
64 // Verbosity scale:
65 // 0 = nothing
66 // 1 = warning for energy non-conservation
67 // 2 = details of energy budget
68 // 3 = calculation of cross sections, file openings, sampling of atoms
69 // 4 = entering in methods
70}
71
72//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73
74G4Penelope08GammaConversionModel::~G4Penelope08GammaConversionModel()
75{
76 std::map <const G4int,G4PhysicsFreeVector*>::iterator i;
77 if (logAtomicCrossSection)
78 {
79 for (i=logAtomicCrossSection->begin();i != logAtomicCrossSection->end();i++)
80 if (i->second) delete i->second;
81 delete logAtomicCrossSection;
82 }
83 if (fEffectiveCharge)
84 delete fEffectiveCharge;
85 if (fMaterialInvScreeningRadius)
86 delete fMaterialInvScreeningRadius;
87 if (fScreeningFunction)
88 delete fScreeningFunction;
89
90}
91
92
93//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
94
95void G4Penelope08GammaConversionModel::Initialise(const G4ParticleDefinition*,
96 const G4DataVector&)
97{
98 if (verboseLevel > 3)
99 G4cout << "Calling G4Penelope08GammaConversionModel::Initialise()" << G4endl;
100
101 // logAtomicCrossSection is created only once, since it is never cleared
102 if (!logAtomicCrossSection)
103 logAtomicCrossSection = new std::map<const G4int,G4PhysicsFreeVector*>;
104
105 //delete old material data...
106 if (fEffectiveCharge)
107 {
108 delete fEffectiveCharge;
109 fEffectiveCharge = 0;
110 }
111 if (fMaterialInvScreeningRadius)
112 {
113 delete fMaterialInvScreeningRadius;
114 fMaterialInvScreeningRadius = 0;
115 }
116 if (fScreeningFunction)
117 {
118 delete fScreeningFunction;
119 fScreeningFunction = 0;
120 }
121 //and create new ones
122 fEffectiveCharge = new std::map<const G4Material*,G4double>;
123 fMaterialInvScreeningRadius = new std::map<const G4Material*,G4double>;
124 fScreeningFunction = new std::map<const G4Material*,std::pair<G4double,G4double> >;
125
126 if (verboseLevel > 0) {
127 G4cout << "Penelope Gamma Conversion model is initialized " << G4endl
128 << "Energy range: "
129 << LowEnergyLimit() / MeV << " MeV - "
130 << HighEnergyLimit() / GeV << " GeV"
131 << G4endl;
132 }
133
134 if(isInitialised) return;
135 fParticleChange = GetParticleChangeForGamma();
136 isInitialised = true;
137}
138
139//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
140
141G4double G4Penelope08GammaConversionModel::ComputeCrossSectionPerAtom(
142 const G4ParticleDefinition*,
143 G4double energy,
144 G4double Z, G4double,
145 G4double, G4double)
146{
147 //
148 // Penelope model.
149 // Cross section (including triplet production) read from database and managed
150 // through the G4CrossSectionHandler utility. Cross section data are from
151 // M.J. Berger and J.H. Hubbel (XCOM), Report NBSIR 887-3598
152 //
153
154 if (energy < fIntrinsicLowEnergyLimit)
155 return 0;
156
157 G4int iZ = (G4int) Z;
158
159 //read data files
160 if (!logAtomicCrossSection->count(iZ))
161 ReadDataFile(iZ);
162 //now it should be ok
163 if (!logAtomicCrossSection->count(iZ))
164 {
165 G4cout << "Problem in G4Penelope08GammaConversion::ComputeCrossSectionPerAtom"
166 << G4endl;
167 G4Exception();
168 }
169
170 G4double cs = 0;
171 G4double logene = log(energy);
172 G4PhysicsFreeVector* theVec = logAtomicCrossSection->find(iZ)->second;
173
174 G4double logXS = theVec->Value(logene);
175 cs = exp(logXS);
176
177 if (verboseLevel > 2)
178 G4cout << "Gamma conversion cross section at " << energy/MeV << " MeV for Z=" << Z <<
179 " = " << cs/barn << " barn" << G4endl;
180 return cs;
181}
182
183//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
184
185void
186G4Penelope08GammaConversionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
187 const G4MaterialCutsCouple* couple,
188 const G4DynamicParticle* aDynamicGamma,
189 G4double,
190 G4double)
191{
192 //
193 // Penelope model.
194 // Final state is sampled according to the Bethe-Heitler model with Coulomb
195 // corrections, according to the semi-empirical model of
196 // J. Baro' et al., Radiat. Phys. Chem. 44 (1994) 531.
197 //
198 // The model uses the high energy Coulomb correction from
199 // H. Davies et al., Phys. Rev. 93 (1954) 788
200 // and atomic screening radii tabulated from
201 // J.H. Hubbel et al., J. Phys. Chem. Ref. Data 9 (1980) 1023
202 // for Z= 1 to 92.
203 //
204 if (verboseLevel > 3)
205 G4cout << "Calling SamplingSecondaries() of G4Penelope08GammaConversionModel" << G4endl;
206
207 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
208
209 // Always kill primary
210 fParticleChange->ProposeTrackStatus(fStopAndKill);
211 fParticleChange->SetProposedKineticEnergy(0.);
212
213 if (photonEnergy <= fIntrinsicLowEnergyLimit)
214 {
215 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
216 return ;
217 }
218
219 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
220 const G4Material* mat = couple->GetMaterial();
221
222 //check if material data are available
223 if (!fEffectiveCharge->count(mat))
224 InitializeScreeningFunctions(mat);
225 if (!fEffectiveCharge->count(mat))
226 {
227 G4cout << "Problem in G4Penelope08GammaConversion::SampleSecondaries()" << G4endl;
228 G4cout << "Unable to allocate the EffectiveCharge data" << G4endl;
229 G4Exception();
230 }
231
232 // eps is the fraction of the photon energy assigned to e- (including rest mass)
233 G4double eps = 0;
234 G4double eki = electron_mass_c2/photonEnergy;
235
236 //Do it fast for photon energy < 1.1 MeV (close to threshold)
237 if (photonEnergy < fSmallEnergy)
238 eps = eki + (1.0-2.0*eki)*G4UniformRand();
239 else
240 {
241 //Complete calculation
242 G4double effC = fEffectiveCharge->find(mat)->second;
243 G4double alz = effC*fine_structure_const;
244 G4double T = sqrt(2.0*eki);
245 G4double F00=(-1.774-1.210e1*alz+1.118e1*alz*alz)*T
246 +(8.523+7.326e1*alz-4.441e1*alz*alz)*T*T
247 -(1.352e1+1.211e2*alz-9.641e1*alz*alz)*T*T*T
248 +(8.946+6.205e1*alz-6.341e1*alz*alz)*T*T*T*T;
249
250 G4double F0b = fScreeningFunction->find(mat)->second.second;
251 G4double g0 = F0b + F00;
252 G4double invRad = fMaterialInvScreeningRadius->find(mat)->second;
253 G4double bmin = 4.0*eki/invRad;
254 std::pair<G4double,G4double> scree = GetScreeningFunctions(bmin);
255 G4double g1 = scree.first;
256 G4double g2 = scree.second;
257 G4double g1min = g1+g0;
258 G4double g2min = g2+g0;
259 G4double xr = 0.5-eki;
260 G4double a1 = 2.*g1min*xr*xr/3.;
261 G4double p1 = a1/(a1+g2min);
262
263 G4bool loopAgain = false;
264 //Random sampling of eps
265 do{
266 loopAgain = false;
267 if (G4UniformRand() <= p1)
268 {
269 G4double ru2m1 = 2.0*G4UniformRand()-1.0;
270 if (ru2m1 < 0)
271 eps = 0.5-xr*pow(std::abs(ru2m1),1./3.);
272 else
273 eps = 0.5+xr*pow(ru2m1,1./3.);
274 G4double B = eki/(invRad*eps*(1.0-eps));
275 scree = GetScreeningFunctions(B);
276 g1 = scree.first;
277 g1 = std::max(g1+g0,0.);
278 if (G4UniformRand()*g1min > g1)
279 loopAgain = true;
280 }
281 else
282 {
283 eps = eki+2.0*xr*G4UniformRand();
284 G4double B = eki/(invRad*eps*(1.0-eps));
285 scree = GetScreeningFunctions(B);
286 g2 = scree.second;
287 g2 = std::max(g2+g0,0.);
288 if (G4UniformRand()*g2min > g2)
289 loopAgain = true;
290 }
291 }while(loopAgain);
292
293 }
294 if (verboseLevel > 4)
295 G4cout << "Sampled eps = " << eps << G4endl;
296
297 G4double electronTotEnergy = eps*photonEnergy;
298 G4double positronTotEnergy = (1.0-eps)*photonEnergy;
299
300 // Scattered electron (positron) angles. ( Z - axis along the parent photon)
301
302 //electron kinematics
303 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
304 G4double costheta_el = G4UniformRand()*2.0-1.0;
305 G4double kk = std::sqrt(electronKineEnergy*(electronKineEnergy+2.*electron_mass_c2));
306 costheta_el = (costheta_el*electronTotEnergy+kk)/(electronTotEnergy+costheta_el*kk);
307 G4double phi_el = twopi * G4UniformRand() ;
308 G4double dirX_el = std::sqrt(1.-costheta_el*costheta_el) * std::cos(phi_el);
309 G4double dirY_el = std::sqrt(1.-costheta_el*costheta_el) * std::sin(phi_el);
310 G4double dirZ_el = costheta_el;
311
312 //positron kinematics
313 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
314 G4double costheta_po = G4UniformRand()*2.0-1.0;
315 kk = std::sqrt(positronKineEnergy*(positronKineEnergy+2.*electron_mass_c2));
316 costheta_po = (costheta_po*positronTotEnergy+kk)/(positronTotEnergy+costheta_po*kk);
317 G4double phi_po = twopi * G4UniformRand() ;
318 G4double dirX_po = std::sqrt(1.-costheta_po*costheta_po) * std::cos(phi_po);
319 G4double dirY_po = std::sqrt(1.-costheta_po*costheta_po) * std::sin(phi_po);
320 G4double dirZ_po = costheta_po;
321
322 // Kinematics of the created pair:
323 // the electron and positron are assumed to have a symetric angular
324 // distribution with respect to the Z axis along the parent photon
325 G4double localEnergyDeposit = 0. ;
326
327 if (electronKineEnergy > 0.0)
328 {
329 G4ThreeVector electronDirection ( dirX_el, dirY_el, dirZ_el);
330 electronDirection.rotateUz(photonDirection);
331 G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
332 electronDirection,
333 electronKineEnergy);
334 fvect->push_back(electron);
335 }
336 else
337 {
338 localEnergyDeposit += electronKineEnergy;
339 electronKineEnergy = 0;
340 }
341
342 //Generate the positron. Real particle in any case, because it will annihilate. If below
343 //threshold, produce it at rest
344 if (positronKineEnergy < 0.0)
345 {
346 localEnergyDeposit += positronKineEnergy;
347 positronKineEnergy = 0; //produce it at rest
348 }
349 G4ThreeVector positronDirection(dirX_po,dirY_po,dirZ_po);
350 positronDirection.rotateUz(photonDirection);
351 G4DynamicParticle* positron = new G4DynamicParticle(G4Positron::Positron(),
352 positronDirection, positronKineEnergy);
353 fvect->push_back(positron);
354
355 //Add rest of energy to the local energy deposit
356 fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
357
358 if (verboseLevel > 1)
359 {
360 G4cout << "-----------------------------------------------------------" << G4endl;
361 G4cout << "Energy balance from G4Penelope08GammaConversion" << G4endl;
362 G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl;
363 G4cout << "-----------------------------------------------------------" << G4endl;
364 if (electronKineEnergy)
365 G4cout << "Electron (explicitely produced) " << electronKineEnergy/keV << " keV"
366 << G4endl;
367 if (positronKineEnergy)
368 G4cout << "Positron (not at rest) " << positronKineEnergy/keV << " keV" << G4endl;
369 G4cout << "Rest masses of e+/- " << 2.0*electron_mass_c2/keV << " keV" << G4endl;
370 if (localEnergyDeposit)
371 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
372 G4cout << "Total final state: " << (electronKineEnergy+positronKineEnergy+
373 localEnergyDeposit+2.0*electron_mass_c2)/keV <<
374 " keV" << G4endl;
375 G4cout << "-----------------------------------------------------------" << G4endl;
376 }
377 if (verboseLevel > 0)
378 {
379 G4double energyDiff = std::fabs(electronKineEnergy+positronKineEnergy+
380 localEnergyDeposit+2.0*electron_mass_c2-photonEnergy);
381 if (energyDiff > 0.05*keV)
382 G4cout << "Warning from G4Penelope08GammaConversion: problem with energy conservation: "
383 << (electronKineEnergy+positronKineEnergy+
384 localEnergyDeposit+2.0*electron_mass_c2)/keV
385 << " keV (final) vs. " << photonEnergy/keV << " keV (initial)" << G4endl;
386 }
387}
388
389//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
390
391void G4Penelope08GammaConversionModel::ReadDataFile(const G4int Z)
392{
393 if (verboseLevel > 2)
394 {
395 G4cout << "G4Penelope08GammaConversionModel::ReadDataFile()" << G4endl;
396 G4cout << "Going to read Gamma Conversion data files for Z=" << Z << G4endl;
397 }
398
399 char* path = getenv("G4LEDATA");
400 if (!path)
401 {
402 G4String excep =
403 "G4Penelope08GammaConversionModel - G4LEDATA environment variable not set!";
404 G4Exception(excep);
405 }
406
407 /*
408 Read the cross section file
409 */
410 std::ostringstream ost;
411 if (Z>9)
412 ost << path << "/penelope/pairproduction/pdgpp" << Z << ".p08";
413 else
414 ost << path << "/penelope/pairproduction/pdgpp0" << Z << ".p08";
415 std::ifstream file(ost.str().c_str());
416 if (!file.is_open())
417 {
418 G4String excep = "G4Penelope08GammaConversionModel - data file " +
419 G4String(ost.str()) + " not found!";
420 G4Exception(excep);
421 }
422
423 //I have to know in advance how many points are in the data list
424 //to initialize the G4PhysicsFreeVector()
425 size_t ndata=0;
426 G4String line;
427 while( getline(file, line) )
428 ndata++;
429 ndata -= 1; //remove one header line
430 //G4cout << "Found: " << ndata << " lines" << G4endl;
431
432 file.clear();
433 file.close();
434 file.open(ost.str().c_str());
435 G4int readZ =0;
436 file >> readZ;
437
438 if (verboseLevel > 3)
439 G4cout << "Element Z=" << Z << G4endl;
440
441 //check the right file is opened.
442 if (readZ != Z)
443 {
444 G4cout << "G4Penelope08GammaConversionModel::ReadDataFile()" << G4endl;
445 G4cout << "Corrupted data file for Z=" << Z << G4endl;
446 G4Exception();
447 }
448
449 G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector(ndata);
450 G4double ene=0,xs=0;
451 for (size_t i=0;i<ndata;i++)
452 {
453 file >> ene >> xs;
454 //dimensional quantities
455 ene *= eV;
456 xs *= barn;
457 if (xs < 1e-40*cm2) //protection against log(0)
458 xs = 1e-40*cm2;
459 theVec->PutValue(i,log(ene),log(xs));
460 }
461 file.close();
462
463 if (!logAtomicCrossSection)
464 {
465 G4cout << "G4Penelope08RayleighModel::ReadDataFile()" << G4endl;
466 G4cout << "Problem with allocation of logAtomicCrossSection data table " << G4endl;
467 G4Exception();
468 }
469 logAtomicCrossSection->insert(std::make_pair(Z,theVec));
470
471 return;
472
473}
474
475//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
476
477void G4Penelope08GammaConversionModel::InitializeScreeningRadii()
478{
479 G4double temp[99] = {1.2281e+02,7.3167e+01,6.9228e+01,6.7301e+01,6.4696e+01,
480 6.1228e+01,5.7524e+01,5.4033e+01,5.0787e+01,4.7851e+01,4.6373e+01,
481 4.5401e+01,4.4503e+01,4.3815e+01,4.3074e+01,4.2321e+01,4.1586e+01,
482 4.0953e+01,4.0524e+01,4.0256e+01,3.9756e+01,3.9144e+01,3.8462e+01,
483 3.7778e+01,3.7174e+01,3.6663e+01,3.5986e+01,3.5317e+01,3.4688e+01,
484 3.4197e+01,3.3786e+01,3.3422e+01,3.3068e+01,3.2740e+01,3.2438e+01,
485 3.2143e+01,3.1884e+01,3.1622e+01,3.1438e+01,3.1142e+01,3.0950e+01,
486 3.0758e+01,3.0561e+01,3.0285e+01,3.0097e+01,2.9832e+01,2.9581e+01,
487 2.9411e+01,2.9247e+01,2.9085e+01,2.8930e+01,2.8721e+01,2.8580e+01,
488 2.8442e+01,2.8312e+01,2.8139e+01,2.7973e+01,2.7819e+01,2.7675e+01,
489 2.7496e+01,2.7285e+01,2.7093e+01,2.6911e+01,2.6705e+01,2.6516e+01,
490 2.6304e+01,2.6108e+01,2.5929e+01,2.5730e+01,2.5577e+01,2.5403e+01,
491 2.5245e+01,2.5100e+01,2.4941e+01,2.4790e+01,2.4655e+01,2.4506e+01,
492 2.4391e+01,2.4262e+01,2.4145e+01,2.4039e+01,2.3922e+01,2.3813e+01,
493 2.3712e+01,2.3621e+01,2.3523e+01,2.3430e+01,2.3331e+01,2.3238e+01,
494 2.3139e+01,2.3048e+01,2.2967e+01,2.2833e+01,2.2694e+01,2.2624e+01,
495 2.2545e+01,2.2446e+01,2.2358e+01,2.2264e+01};
496
497 //copy temporary vector in class data member
498 for (G4int i=0;i<99;i++)
499 fAtomicScreeningRadius[i] = temp[i];
500}
501
502//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
503
504void G4Penelope08GammaConversionModel::InitializeScreeningFunctions(const G4Material* material)
505{
506 // This is subroutine GPPa0 of Penelope
507 //
508 // 1) calculate the effective Z for the purpose
509 //
510 G4double zeff = 0;
511 G4int intZ = 0;
512 G4int nElements = material->GetNumberOfElements();
513 const G4ElementVector* elementVector = material->GetElementVector();
514
515 //avoid calculations if only one building element!
516 if (nElements == 1)
517 {
518 zeff = (*elementVector)[0]->GetZ();
519 intZ = (G4int) zeff;
520 }
521 else // many elements...let's do the calculation
522 {
523 const G4double* fractionVector = material->GetVecNbOfAtomsPerVolume();
524
525 G4double atot = 0;
526 for (G4int i=0;i<nElements;i++)
527 {
528 G4double Zelement = (*elementVector)[i]->GetZ();
529 G4double Aelement = (*elementVector)[i]->GetA();
530 atot += Aelement*fractionVector[i];
531 zeff += Zelement*Aelement*fractionVector[i]; //average with the number of nuclei
532 }
533 atot /= material->GetTotNbOfAtomsPerVolume();
534 zeff /= (material->GetTotNbOfAtomsPerVolume()*atot);
535
536 intZ = (G4int) (zeff+0.25);
537 if (intZ <= 0)
538 intZ = 1;
539 if (intZ > 99)
540 intZ = 99;
541 }
542
543 if (fEffectiveCharge)
544 fEffectiveCharge->insert(std::make_pair(material,zeff));
545
546 //
547 // 2) Calculate Coulomb Correction
548 //
549 G4double alz = fine_structure_const*zeff;
550 G4double alzSquared = alz*alz;
551 G4double fc = alzSquared*(0.202059-alzSquared*
552 (0.03693-alzSquared*
553 (0.00835-alzSquared*(0.00201-alzSquared*
554 (0.00049-alzSquared*
555 (0.00012-alzSquared*0.00003)))))
556 +1.0/(alzSquared+1.0));
557 //
558 // 3) Screening functions and low-energy corrections
559 //
560 G4double matRadius = 2.0/ fAtomicScreeningRadius[intZ-1];
561 if (fMaterialInvScreeningRadius)
562 fMaterialInvScreeningRadius->insert(std::make_pair(material,matRadius));
563
564 std::pair<G4double,G4double> myPair(0,0);
565 G4double f0a = 4.0*log(fAtomicScreeningRadius[intZ-1]);
566 G4double f0b = f0a - 4.0*fc;
567 myPair.first = f0a;
568 myPair.second = f0b;
569
570 if (fScreeningFunction)
571 fScreeningFunction->insert(std::make_pair(material,myPair));
572
573 if (verboseLevel > 2)
574 {
575 G4cout << "Average Z for material " << material->GetName() << " = " <<
576 zeff << G4endl;
577 G4cout << "Effective radius for material " << material->GetName() << " = " <<
578 fAtomicScreeningRadius[intZ-1] << " m_e*c/hbar --> BCB = " <<
579 matRadius << G4endl;
580 G4cout << "Screening parameters F0 for material " << material->GetName() << " = " <<
581 f0a << "," << f0b << G4endl;
582 }
583 return;
584}
585
586//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
587
588std::pair<G4double,G4double>
589G4Penelope08GammaConversionModel::GetScreeningFunctions(G4double B)
590{
591 // This is subroutine SCHIFF of Penelope
592 //
593 // Screening Functions F1(B) and F2(B) in the Bethe-Heitler differential cross
594 // section for pair production
595 //
596 std::pair<G4double,G4double> result(0.,0.);
597 G4double BSquared = B*B;
598 G4double f1 = 2.0-2.0*log(1.0+BSquared);
599 G4double f2 = f1 - 6.66666666e-1; // (-2/3)
600 if (B < 1.0e-10)
601 f1 = f1-twopi*B;
602 else
603 {
604 G4double a0 = 4.0*B*std::atan(1./B);
605 f1 = f1 - a0;
606 f2 += 2.0*BSquared*(4.0-a0-3.0*log((1.0+BSquared)/BSquared));
607 }
608 G4double g1 = 0.5*(3.0*f1-f2);
609 G4double g2 = 0.25*(3.0*f1+f2);
610
611 result.first = g1;
612 result.second = g2;
613
614 return result;
615}
Note: See TracBrowser for help on using the repository browser.