source: trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeBremsstrahlungModel.cc@ 1142

Last change on this file since 1142 was 1058, checked in by garnier, 17 years ago

file release beta
File size: 19.8 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: G4PenelopeBremsstrahlungModel.cc,v 1.5 2009/05/14 10:56:09 pandola Exp $
27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
28//
29// Author: Luciano Pandola
30// --------
31// 05 Dec 2008 L Pandola Migration from process to model
32// 25 Mar 2008 L Pandola Fixed .unit() call
33// 16 Apr 2009 V Ivanchenko Cleanup initialisation and generation of secondaries:
34// - apply internal high-energy limit only in constructor
35// - do not apply low-energy limit (default is 0)
36// - added MinEnergyCut method
37// - do not change track status
38// 14 May 2009 L Pandola Explicitely set to zero pointers deleted in
39// Initialise(), since they are checked later on
40//
41#include "G4PenelopeBremsstrahlungModel.hh"
42#include "G4PenelopeBremsstrahlungContinuous.hh"
43#include "G4PenelopeBremsstrahlungAngular.hh"
44#include "G4eBremsstrahlungSpectrum.hh"
45#include "G4CrossSectionHandler.hh"
46#include "G4VEMDataSet.hh"
47#include "G4DataVector.hh"
48#include "G4Positron.hh"
49#include "G4Electron.hh"
50#include "G4Gamma.hh"
51#include "G4MaterialCutsCouple.hh"
52#include "G4LogLogInterpolation.hh"
53
54//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
55
56G4PenelopeBremsstrahlungModel::G4PenelopeBremsstrahlungModel(const G4ParticleDefinition*,
57 const G4String& nam)
58 :G4VEmModel(nam),isInitialised(false),energySpectrum(0),
59 angularData(0),stoppingPowerData(0),crossSectionHandler(0)
60{
61 fIntrinsicLowEnergyLimit = 100.0*eV;
62 fIntrinsicHighEnergyLimit = 100.0*GeV;
63 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
64 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
65 //
66
67 verboseLevel= 0;
68
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 //These vectors do not change when materials or cut change.
77 //Therefore I can read it at the constructor
78 angularData = new std::map<G4int,G4PenelopeBremsstrahlungAngular*>;
79
80 //These data do not depend on materials and cuts.
81 G4DataVector eBins;
82
83 eBins.push_back(1.0e-12);
84 eBins.push_back(0.05);
85 eBins.push_back(0.075);
86 eBins.push_back(0.1);
87 eBins.push_back(0.125);
88 eBins.push_back(0.15);
89 eBins.push_back(0.2);
90 eBins.push_back(0.25);
91 eBins.push_back(0.3);
92 eBins.push_back(0.35);
93 eBins.push_back(0.40);
94 eBins.push_back(0.45);
95 eBins.push_back(0.50);
96 eBins.push_back(0.55);
97 eBins.push_back(0.60);
98 eBins.push_back(0.65);
99 eBins.push_back(0.70);
100 eBins.push_back(0.75);
101 eBins.push_back(0.80);
102 eBins.push_back(0.85);
103 eBins.push_back(0.90);
104 eBins.push_back(0.925);
105 eBins.push_back(0.95);
106 eBins.push_back(0.97);
107 eBins.push_back(0.99);
108 eBins.push_back(0.995);
109 eBins.push_back(0.999);
110 eBins.push_back(0.9995);
111 eBins.push_back(0.9999);
112 eBins.push_back(0.99995);
113 eBins.push_back(0.99999);
114 eBins.push_back(1.0);
115
116 const G4String dataName("/penelope/br-sp-pen.dat");
117 energySpectrum = new G4eBremsstrahlungSpectrum(eBins,dataName);
118}
119
120//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
121
122G4PenelopeBremsstrahlungModel::~G4PenelopeBremsstrahlungModel()
123{
124 if (crossSectionHandler)
125 delete crossSectionHandler;
126
127 if (energySpectrum)
128 delete energySpectrum;
129
130 if (angularData)
131 {
132 std::map <G4int,G4PenelopeBremsstrahlungAngular*>::iterator i;
133 for (i=angularData->begin();i != angularData->end();i++)
134 if (i->second) delete i->second;
135 delete angularData;
136 }
137
138 if (stoppingPowerData)
139 {
140 std::map <std::pair<G4int,G4double>,G4PenelopeBremsstrahlungContinuous*>::iterator j;
141 for (j=stoppingPowerData->begin();j != stoppingPowerData->end();j++)
142 if (j->second) delete j->second;
143 delete stoppingPowerData;
144 }
145}
146
147//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
148
149void G4PenelopeBremsstrahlungModel::Initialise(const G4ParticleDefinition* particle,
150 const G4DataVector& cuts)
151{
152 if (verboseLevel > 3)
153 G4cout << "Calling G4PenelopeBremsstrahlungModel::Initialise()" << G4endl;
154
155 // Delete everything, but angular data (do not depend on cuts)
156 if (crossSectionHandler)
157 {
158 crossSectionHandler->Clear();
159 delete crossSectionHandler;
160 crossSectionHandler = 0;
161 }
162
163 if (stoppingPowerData)
164 {
165 std::map <std::pair<G4int,G4double>,G4PenelopeBremsstrahlungContinuous*>::iterator j;
166 for (j=stoppingPowerData->begin();j != stoppingPowerData->end();j++)
167 if (j->second)
168 {
169 delete j->second;
170 j->second = 0;
171 }
172 delete stoppingPowerData;
173 stoppingPowerData = 0;
174 }
175
176 crossSectionHandler = new G4CrossSectionHandler();
177 crossSectionHandler->Clear();
178 //
179 if (particle==G4Electron::Electron())
180 crossSectionHandler->LoadData("brem/br-cs-");
181 else
182 crossSectionHandler->LoadData("penelope/br-cs-pos-"); //cross section for positrons
183
184 //This is used to retrieve cross section values later on
185 crossSectionHandler->BuildMeanFreePathForMaterials();
186
187 if (verboseLevel > 2)
188 G4cout << "Loaded cross section files for PenelopeBremsstrahlungModel" << G4endl;
189
190 if (verboseLevel > 0) {
191 G4cout << "Penelope Bremsstrahlung model is initialized " << G4endl
192 << "Energy range: "
193 << LowEnergyLimit() / keV << " keV - "
194 << HighEnergyLimit() / GeV << " GeV"
195 << G4endl;
196 }
197
198 //This has to be invoked AFTER the crossSectionHandler has been created,
199 //because it makes use of ComputeCrossSectionPerAtom()
200 InitialiseElementSelectors(particle,cuts);
201
202 if(isInitialised) return;
203 fParticleChange = GetParticleChangeForLoss();
204 isInitialised = true;
205}
206
207//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
208
209G4double G4PenelopeBremsstrahlungModel::MinEnergyCut(const G4ParticleDefinition*,
210 const G4MaterialCutsCouple*)
211{
212 return 250.*eV;
213}
214
215//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
216
217G4double G4PenelopeBremsstrahlungModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
218 G4double kinEnergy,
219 G4double Z,
220 G4double,
221 G4double cutEnergy,
222 G4double)
223{
224 // Penelope model to calculate cross section for hard bremsstrahlung emission
225 // (gamma energy > cutEnergy).
226 //
227 // The total bremsstrahlung cross section is read from database, following data
228 // reported in the EEDL library
229 // D.E.Cullen et al., Report UCRL-50400 (Lawrence Livermore National Laboratory) (1989)
230 // The probability to have photon emission above a given threshold is calculated
231 // analytically using the differential cross section model dSigma/dW = F(x)/x, where
232 // W is the outgoing photon energy and x = W/E is the ratio of the photon energy to the
233 // incident energy. The function F(x) is tabulated (for all elements) using 32 points in x
234 // ranging from 1e-12 to 1. Data are derived from
235 // S.M.Seltzer and M.J.Berger, At.Data Nucl.Data Tables 35,345 (1986)
236 // Differential cross sections for electrons and positrons dSigma/dW are assumed to scale
237 // with a function S(Z,E) which does not depend on W; therefore, only overall cross section
238 // changes but not the shape of the photon energy spectrum.
239 //
240
241 if (verboseLevel > 3)
242 G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopeBremsstrahlungModel" << G4endl;
243
244 G4int iZ = (G4int) Z;
245
246 // VI - not needed in run time
247 // if (!crossSectionHandler)
248 // {
249 // G4cout << "G4PenelopeBremsstrahlungModel::ComputeCrossSectionPerAtom" << G4endl;
250 // G4cout << "The cross section handler is not correctly initialized" << G4endl;
251 // G4Exception();
252 // }
253 G4double totalCs = crossSectionHandler->FindValue(iZ,kinEnergy);
254 G4double cs = totalCs * energySpectrum->Probability(iZ,cutEnergy,kinEnergy,kinEnergy);
255
256 if (verboseLevel > 2)
257 {
258 G4cout << "Bremsstrahlung cross section at " << kinEnergy/MeV << " MeV for Z=" << Z <<
259 " and energy > " << cutEnergy/keV << " keV --> " << cs/barn << " barn" << G4endl;
260 G4cout << "Total bremsstrahlung cross section at " << kinEnergy/MeV << " MeV for Z=" <<
261 Z << " --> " << totalCs/barn << " barn" << G4endl;
262 }
263 return cs;
264
265}
266
267//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
268G4double
269G4PenelopeBremsstrahlungModel::ComputeDEDXPerVolume(const G4Material* theMaterial,
270 const G4ParticleDefinition* theParticle,
271 G4double kineticEnergy,
272 G4double cutEnergy)
273{
274 // Penelope model to calculate the stopping power (in [Energy]/[Length]) for soft
275 // bremsstrahlung emission (gamma energy < cutEnergy).
276 //
277 // The actual calculation is performed by the helper class
278 // G4PenelopeBremsstrahlungContinuous and its method CalculateStopping(). Notice:
279 // CalculateStopping() gives the stopping cross section, namely the first momentum of
280 // dSigma/dW, restricted to W < cut (W = gamma energy) This is dimensionally:
281 // [Energy]*[Surface]
282 // The calculation is performed by interpolation (in E = incident energy and
283 // x=W/E) from the tabulated data derived from
284 // M.J.Berger and S.M.Seltzer, Report NBSIR 82-2550 (National Bureau of Standards) (1982);
285 // for electrons.
286 // For positrons, dSigma/dW are assumed to scale with a function S(Z,E) with respect to electrons.
287 // An analytical approximation for the scaling function S(Z,E) is given in
288 // L.Kim et al., Phys.Rev.A 33,3002 (1986)
289 //
290 if (!stoppingPowerData)
291 stoppingPowerData = new std::map<std::pair<G4int,G4double>,
292 G4PenelopeBremsstrahlungContinuous*>;
293
294 const G4ElementVector* theElementVector = theMaterial->GetElementVector();
295 const G4double* theAtomicNumDensityVector = theMaterial->GetAtomicNumDensityVector();
296
297 G4double sPower = 0.0;
298
299 //Loop on the elements of the material
300 for (size_t iel=0;iel<theMaterial->GetNumberOfElements();iel++)
301 {
302 G4int iZ = (G4int) ((*theElementVector)[iel]->GetZ());
303 G4PenelopeBremsstrahlungContinuous* theContinuousCalculator =
304 GetStoppingPowerData(iZ,cutEnergy,theParticle);
305 sPower += theContinuousCalculator->CalculateStopping(kineticEnergy)*
306 theAtomicNumDensityVector[iel];
307 }
308
309 if (verboseLevel > 2)
310 {
311 G4cout << "Bremsstrahlung stopping power at " << kineticEnergy/MeV
312 << " MeV for material " << theMaterial->GetName()
313 << " and energy < " << cutEnergy/keV << " keV --> "
314 << sPower/(keV/mm) << " keV/mm" << G4endl;
315 }
316
317 return sPower;
318}
319
320//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
321
322void G4PenelopeBremsstrahlungModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
323 const G4MaterialCutsCouple* couple,
324 const G4DynamicParticle* aDynamicParticle,
325 G4double cutG,G4double)
326{
327 // Penelope model to sample the final state for hard bremsstrahlung emission
328 // (gamma energy < cutEnergy).
329 // The energy distributionof the emitted photons is sampled according to the F(x)
330 // function tabulated in the database from
331 // S.M.Seltzer and M.J.Berger, At.Data Nucl.Data Tables 35,345 (1986)
332 // The database contains the function F(x) (32 points) for 57 energies of the
333 // incident electron between 1 keV and 100 GeV. For other primary energies,
334 // logarithmic interpolation is used to obtain the values of the function F(x).
335 // The double differential cross section dSigma/(dW dOmega), with W=photon energy,
336 // is described as
337 // dSigma/(dW dOmega) = dSigma/dW * p(Z,E,x,cosTheta)
338 // where the shape function p depends on atomic number, incident energy and x=W/E.
339 // Numerical values of the shape function, calculated by partial-waves methods, have been
340 // reported in
341 // L.Kissel et al., At.Data Nucl.Data.Tab. 28,381 (1983);
342 // for Z=2,8,13,47,79 and 92; E=1,5,10,50,100 and 500 keV; x=0,0.6,0.8 and 0.95. The
343 // function p(Z,E,x,cosTheta) is approximated by a Lorentz-dipole function as reported in
344 // Penelope - A Code System for Monte Carlo Simulation of Electron and Photon Transport,
345 // Workshop Proceedings Issy-les-Moulineaux, France, 5-7 November 2001, AEN-NEA;
346 // The analytical function contains two adjustable parameters that are obtained by fitting
347 // the complete solution from
348 // L.Kissel et al., At.Data Nucl.Data.Tab. 28,381 (1983);
349 // This allows the evaluation of p(Z,E,x,cosTheta) for any choice of Z, E and x. The actual
350 // sampling of cos(theta) is performed in the helper class
351 // G4PenelopeBremsstrahlungAngular, method ExtractCosTheta()
352 // Energy and direction of the primary particle are updated according to energy-momentum
353 // conservation. For positrons, it is sampled the same final state as for electrons.
354 //
355 if (verboseLevel > 3)
356 G4cout << "Calling SampleSecondaries() of G4PenelopeBremsstrahlungModel" << G4endl;
357
358 G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
359 const G4ParticleDefinition* theParticle = aDynamicParticle->GetDefinition();
360
361 if (kineticEnergy <= fIntrinsicLowEnergyLimit)
362 {
363 fParticleChange->SetProposedKineticEnergy(0.);
364 fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy);
365 return ;
366 }
367
368 G4ParticleMomentum particleDirection0 = aDynamicParticle->GetMomentumDirection();
369 //This is the momentum
370 G4ThreeVector initialMomentum = aDynamicParticle->GetMomentum();
371
372 //One can use Vladimir's selector!
373 if (verboseLevel > 2)
374 G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl;
375 // atom can be selected effitiantly if element selectors are initialised
376 const G4Element* anElement = SelectRandomAtom(couple,theParticle,kineticEnergy);
377 G4int iZ = (G4int) anElement->GetZ();
378 if (verboseLevel > 2)
379 G4cout << "Selected " << anElement->GetName() << G4endl;
380 //
381
382 //Sample gamma's energy according to the spectrum
383 G4double gammaEnergy = energySpectrum->SampleEnergy(iZ,cutG,kineticEnergy,kineticEnergy);
384
385 //Now sample cosTheta for the Gamma
386 G4double cosThetaPrimary = GetAngularDataForZ(iZ)->ExtractCosTheta(kineticEnergy,gammaEnergy);
387
388 G4double residualPrimaryEnergy = kineticEnergy-gammaEnergy;
389 if (residualPrimaryEnergy < 0)
390 {
391 //Ok we have a problem, all energy goes with the gamma
392 gammaEnergy += residualPrimaryEnergy;
393 residualPrimaryEnergy = 0.0;
394 }
395
396 //Get primary kinematics
397 G4double sinTheta = std::sqrt(1. - cosThetaPrimary*cosThetaPrimary);
398 G4double phi = twopi * G4UniformRand();
399 G4ThreeVector gammaDirection1(sinTheta* std::cos(phi),
400 sinTheta* std::sin(phi),
401 cosThetaPrimary);
402
403 gammaDirection1.rotateUz(particleDirection0);
404
405 //Produce final state according to momentum conservation
406 G4ThreeVector particleDirection1 = initialMomentum - gammaEnergy*gammaDirection1;
407 particleDirection1 = particleDirection1.unit(); //normalize
408
409 //Update the primary particle
410 if (residualPrimaryEnergy > 0.)
411 {
412 fParticleChange->ProposeMomentumDirection(particleDirection1);
413 fParticleChange->SetProposedKineticEnergy(residualPrimaryEnergy);
414 }
415 else
416 {
417 fParticleChange->SetProposedKineticEnergy(0.);
418 }
419
420 //Now produce the photon
421 G4DynamicParticle* theGamma = new G4DynamicParticle(G4Gamma::Gamma(),
422 gammaDirection1,
423 gammaEnergy);
424 fvect->push_back(theGamma);
425
426 if (verboseLevel > 1)
427 {
428 G4cout << "-----------------------------------------------------------" << G4endl;
429 G4cout << "Energy balance from G4PenelopeBremsstrahlung" << G4endl;
430 G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
431 G4cout << "-----------------------------------------------------------" << G4endl;
432 G4cout << "Outgoing primary energy: " << residualPrimaryEnergy/keV << " keV" << G4endl;
433 G4cout << "Bremsstrahlung photon " << gammaEnergy/keV << " keV" << G4endl;
434 G4cout << "Total final state: " << (residualPrimaryEnergy+gammaEnergy)/keV
435 << " keV" << G4endl;
436 G4cout << "-----------------------------------------------------------" << G4endl;
437 }
438 if (verboseLevel > 0)
439 {
440 G4double energyDiff = std::fabs(residualPrimaryEnergy+gammaEnergy-kineticEnergy);
441 if (energyDiff > 0.05*keV)
442 G4cout << "Warning from G4PenelopeBremsstrahlung: problem with energy conservation: " <<
443 (residualPrimaryEnergy+gammaEnergy)/keV <<
444 " keV (final) vs. " <<
445 kineticEnergy/keV << " keV (initial)" << G4endl;
446 }
447}
448
449//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
450
451G4PenelopeBremsstrahlungAngular* G4PenelopeBremsstrahlungModel::GetAngularDataForZ(G4int iZ)
452{
453 if (!angularData)
454 angularData = new std::map<G4int,G4PenelopeBremsstrahlungAngular*>;
455
456 if (angularData->count(iZ)) //the material already exists
457 return angularData->find(iZ)->second;
458
459 //Otherwise create a new object, store it and return it
460 G4PenelopeBremsstrahlungAngular* theAngular = new G4PenelopeBremsstrahlungAngular(iZ);
461 angularData->insert(std::make_pair(iZ,theAngular));
462
463 if (angularData->count(iZ)) //the material should exist now
464 return angularData->find(iZ)->second;
465 else
466 {
467 G4Exception("Problem in G4PenelopeBremsstrahlungModel::GetAngularDataForZ()");
468 return 0;
469 }
470}
471
472//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
473
474G4PenelopeBremsstrahlungContinuous*
475G4PenelopeBremsstrahlungModel::GetStoppingPowerData(G4int iZ,G4double energyCut,
476 const G4ParticleDefinition*
477 theParticle)
478{
479 if (!stoppingPowerData)
480 stoppingPowerData = new std::map<std::pair<G4int,G4double>,G4PenelopeBremsstrahlungContinuous*>;
481
482 std::pair<G4int,G4double> theKey = std::make_pair(iZ,energyCut);
483
484 if (stoppingPowerData->count(theKey)) //the material already exists
485 return stoppingPowerData->find(theKey)->second;
486
487 //Otherwise create a new object, store it and return it
488 G4String theParticleName = theParticle->GetParticleName();
489 G4PenelopeBremsstrahlungContinuous* theContinuous = new
490 G4PenelopeBremsstrahlungContinuous(iZ,energyCut,LowEnergyLimit(),HighEnergyLimit(),theParticleName);
491 stoppingPowerData->insert(std::make_pair(theKey,theContinuous));
492
493 if (stoppingPowerData->count(theKey)) //the material should exist now
494 return stoppingPowerData->find(theKey)->second;
495 else
496 {
497 G4Exception("Problem in G4PenelopeBremsstrahlungModel::GetStoppingPowerData()");
498 return 0;
499 }
500}
Note: See TracBrowser for help on using the repository browser.