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

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

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