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

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

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