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

Last change on this file since 1337 was 1337, checked in by garnier, 14 years ago

tag geant4.9.4 beta 1 + modifs locales
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.7 2009/06/11 15:47:08 mantero Exp $
27// GEANT4 tag $Name: geant4-09-04-beta-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.