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

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

geant4 tag 9.4
File size: 20.0 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.8 2010/11/25 09:44:05 pandola Exp $
27// GEANT4 tag $Name: geant4-09-04-ref-00 $
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  G4VEMDataSet* emdata = 
186    crossSectionHandler->BuildMeanFreePathForMaterials();
187  //The method BuildMeanFreePathForMaterials() is required here only to force
188  //the building of an internal table: the output pointer can be deleted
189  delete emdata;   
190
191  if (verboseLevel > 2)
192    G4cout << "Loaded cross section files for PenelopeBremsstrahlungModel" << G4endl;
193 
194  if (verboseLevel > 0) {
195    G4cout << "Penelope Bremsstrahlung model is initialized " << G4endl
196           << "Energy range: "
197           << LowEnergyLimit() / keV << " keV - "
198           << HighEnergyLimit() / GeV << " GeV"
199           << G4endl;
200  } 
201
202  //This has to be invoked AFTER the crossSectionHandler has been created,
203  //because it makes use of ComputeCrossSectionPerAtom()
204  InitialiseElementSelectors(particle,cuts);
205
206  if(isInitialised) return;
207  fParticleChange = GetParticleChangeForLoss();
208  isInitialised = true;
209}
210
211//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
212
213G4double G4PenelopeBremsstrahlungModel::MinEnergyCut(const G4ParticleDefinition*,
214                                                     const G4MaterialCutsCouple*)
215{
216  return 250.*eV;
217}
218
219//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
220
221G4double G4PenelopeBremsstrahlungModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
222                                                                   G4double kinEnergy,
223                                                                   G4double Z,
224                                                                   G4double,
225                                                                   G4double cutEnergy,
226                                                                   G4double)
227{
228  // Penelope model to calculate cross section for hard bremsstrahlung emission
229  // (gamma energy > cutEnergy).
230  //
231  // The total bremsstrahlung cross section is read from database, following data
232  // reported in the EEDL library
233  //  D.E.Cullen et al., Report UCRL-50400 (Lawrence Livermore National Laboratory) (1989)
234  // The probability to have photon emission above a given threshold is calculated
235  // analytically using the differential cross section model dSigma/dW = F(x)/x, where
236  // W is the outgoing photon energy and x = W/E is the ratio of the photon energy to the
237  // incident energy. The function F(x) is tabulated (for all elements) using 32 points in x
238  // ranging from 1e-12 to 1. Data are derived from
239  //  S.M.Seltzer and M.J.Berger, At.Data Nucl.Data Tables 35,345 (1986)
240  // Differential cross sections for electrons and positrons dSigma/dW are assumed to scale
241  // with a function S(Z,E) which does not depend on W; therefore, only overall cross section
242  // changes but not the shape of the photon energy spectrum.
243  //
244
245  if (verboseLevel > 3)
246    G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopeBremsstrahlungModel" << G4endl;
247 
248  G4int iZ = (G4int) Z;
249
250  // VI - not needed in run time
251  // if (!crossSectionHandler)
252  //  {
253  //    G4cout << "G4PenelopeBremsstrahlungModel::ComputeCrossSectionPerAtom" << G4endl;
254  //    G4cout << "The cross section handler is not correctly initialized" << G4endl;
255  //    G4Exception();
256  //  }
257  G4double totalCs = crossSectionHandler->FindValue(iZ,kinEnergy);
258  G4double cs = totalCs * energySpectrum->Probability(iZ,cutEnergy,kinEnergy,kinEnergy);
259
260  if (verboseLevel > 2)
261    {
262      G4cout << "Bremsstrahlung cross section at " << kinEnergy/MeV << " MeV for Z=" << Z <<
263        " and energy > " << cutEnergy/keV << " keV --> " << cs/barn << " barn" << G4endl;
264      G4cout << "Total bremsstrahlung cross section at " << kinEnergy/MeV << " MeV for Z=" << 
265      Z << " --> " << totalCs/barn << " barn" << G4endl;
266    }
267  return cs;
268
269}
270
271//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
272G4double
273G4PenelopeBremsstrahlungModel::ComputeDEDXPerVolume(const G4Material* theMaterial,
274                                                    const G4ParticleDefinition* theParticle,
275                                                    G4double kineticEnergy,
276                                                    G4double cutEnergy)
277{
278  // Penelope model to calculate the stopping power (in [Energy]/[Length]) for soft
279  // bremsstrahlung emission (gamma energy < cutEnergy).
280  //
281  // The actual calculation is performed by the helper class
282  // G4PenelopeBremsstrahlungContinuous and its method CalculateStopping(). Notice:
283  // CalculateStopping() gives the stopping cross section, namely the first momentum of
284  // dSigma/dW, restricted to W < cut (W = gamma energy) This is dimensionally:
285  //  [Energy]*[Surface]
286  // The calculation is performed by interpolation (in E = incident energy and
287  // x=W/E) from the tabulated data derived from
288  //  M.J.Berger and S.M.Seltzer, Report NBSIR 82-2550 (National Bureau of Standards) (1982);
289  // for electrons.
290  // For positrons, dSigma/dW are assumed to scale with a function S(Z,E) with respect to electrons.
291  // An analytical approximation for the scaling function S(Z,E) is given in
292  //  L.Kim et al., Phys.Rev.A 33,3002 (1986)
293  //
294  if (!stoppingPowerData)
295    stoppingPowerData = new std::map<std::pair<G4int,G4double>,
296      G4PenelopeBremsstrahlungContinuous*>;
297 
298  const G4ElementVector* theElementVector = theMaterial->GetElementVector();
299  const G4double* theAtomicNumDensityVector = theMaterial->GetAtomicNumDensityVector();
300
301  G4double sPower = 0.0;
302
303  //Loop on the elements of the material
304  for (size_t iel=0;iel<theMaterial->GetNumberOfElements();iel++)
305    {
306      G4int iZ = (G4int) ((*theElementVector)[iel]->GetZ());
307      G4PenelopeBremsstrahlungContinuous* theContinuousCalculator = 
308        GetStoppingPowerData(iZ,cutEnergy,theParticle);
309      sPower += theContinuousCalculator->CalculateStopping(kineticEnergy)*
310        theAtomicNumDensityVector[iel];
311    }
312 
313   if (verboseLevel > 2)
314    {
315      G4cout << "Bremsstrahlung stopping power at " << kineticEnergy/MeV
316             << " MeV for material " << theMaterial->GetName() 
317             << " and energy < " << cutEnergy/keV << " keV --> " 
318             << sPower/(keV/mm) << " keV/mm" << G4endl;
319    }
320
321  return sPower;
322}
323
324//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
325
326void G4PenelopeBremsstrahlungModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
327                                                      const G4MaterialCutsCouple* couple,
328                                                      const G4DynamicParticle* aDynamicParticle,
329                                                      G4double cutG,G4double)
330{
331  // Penelope model to sample the final state for hard bremsstrahlung emission
332  // (gamma energy < cutEnergy).
333  // The energy distributionof the emitted photons is sampled according to the F(x)
334  // function tabulated in the database from
335  //  S.M.Seltzer and M.J.Berger, At.Data Nucl.Data Tables 35,345 (1986)
336  // The database contains the function F(x) (32 points) for 57 energies of the
337  // incident electron between 1 keV and 100 GeV. For other primary energies,
338  // logarithmic interpolation is used to obtain the values of the function F(x).
339  // The double differential cross section dSigma/(dW dOmega), with W=photon energy,
340  // is described as
341  //  dSigma/(dW dOmega) = dSigma/dW * p(Z,E,x,cosTheta)
342  // where the shape function p depends on atomic number, incident energy and x=W/E.
343  // Numerical values of the shape function, calculated by partial-waves methods, have been
344  // reported in
345  //  L.Kissel et al., At.Data Nucl.Data.Tab. 28,381 (1983);
346  // 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
347  // function p(Z,E,x,cosTheta) is approximated by a Lorentz-dipole function as reported in
348  //  Penelope - A Code System for Monte Carlo Simulation of Electron and Photon Transport,
349  //  Workshop Proceedings Issy-les-Moulineaux, France, 5-7 November 2001, AEN-NEA;
350  // The analytical function contains two adjustable parameters that are obtained by fitting
351  // the complete solution from
352  //  L.Kissel et al., At.Data Nucl.Data.Tab. 28,381 (1983);
353  // This allows the evaluation of p(Z,E,x,cosTheta) for any choice of Z, E and x. The actual
354  // sampling of cos(theta) is performed in the helper class
355  //  G4PenelopeBremsstrahlungAngular, method ExtractCosTheta()
356  // Energy and direction of the primary particle are updated according to energy-momentum
357  // conservation. For positrons, it is sampled the same final state as for electrons.
358  //
359  if (verboseLevel > 3)
360    G4cout << "Calling SampleSecondaries() of G4PenelopeBremsstrahlungModel" << G4endl;
361
362  G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
363  const G4ParticleDefinition* theParticle = aDynamicParticle->GetDefinition();
364
365  if (kineticEnergy <= fIntrinsicLowEnergyLimit)
366   {
367     fParticleChange->SetProposedKineticEnergy(0.);
368     fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy);
369     return ;
370   }
371 
372  G4ParticleMomentum particleDirection0 = aDynamicParticle->GetMomentumDirection();
373  //This is the momentum
374  G4ThreeVector initialMomentum =  aDynamicParticle->GetMomentum();
375 
376  //One can use Vladimir's selector!
377  if (verboseLevel > 2)
378    G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl;
379  // atom can be selected effitiantly if element selectors are initialised
380  const G4Element* anElement = SelectRandomAtom(couple,theParticle,kineticEnergy);
381  G4int iZ = (G4int) anElement->GetZ();
382  if (verboseLevel > 2)
383    G4cout << "Selected " << anElement->GetName() << G4endl;
384  //
385
386  //Sample gamma's energy according to the spectrum
387  G4double gammaEnergy = energySpectrum->SampleEnergy(iZ,cutG,kineticEnergy,kineticEnergy);
388 
389  //Now sample cosTheta for the Gamma
390  G4double cosThetaPrimary = GetAngularDataForZ(iZ)->ExtractCosTheta(kineticEnergy,gammaEnergy);
391 
392  G4double residualPrimaryEnergy = kineticEnergy-gammaEnergy;
393  if (residualPrimaryEnergy < 0)
394    {
395      //Ok we have a problem, all energy goes with the gamma
396      gammaEnergy += residualPrimaryEnergy;
397      residualPrimaryEnergy = 0.0;
398    }
399
400  //Get primary kinematics
401  G4double sinTheta = std::sqrt(1. - cosThetaPrimary*cosThetaPrimary);
402  G4double phi  = twopi * G4UniformRand(); 
403  G4ThreeVector gammaDirection1(sinTheta* std::cos(phi),
404                                sinTheta* std::sin(phi),
405                                cosThetaPrimary);
406 
407  gammaDirection1.rotateUz(particleDirection0);
408 
409  //Produce final state according to momentum conservation
410  G4ThreeVector particleDirection1 = initialMomentum - gammaEnergy*gammaDirection1;
411  particleDirection1 = particleDirection1.unit(); //normalize   
412
413  //Update the primary particle
414  if (residualPrimaryEnergy > 0.)
415    {
416      fParticleChange->ProposeMomentumDirection(particleDirection1);
417      fParticleChange->SetProposedKineticEnergy(residualPrimaryEnergy);
418    }
419  else
420    {
421      fParticleChange->SetProposedKineticEnergy(0.);
422    }
423
424  //Now produce the photon
425  G4DynamicParticle* theGamma = new G4DynamicParticle(G4Gamma::Gamma(),
426                                                      gammaDirection1,
427                                                      gammaEnergy);
428  fvect->push_back(theGamma);
429
430  if (verboseLevel > 1)
431    {
432      G4cout << "-----------------------------------------------------------" << G4endl;
433      G4cout << "Energy balance from G4PenelopeBremsstrahlung" << G4endl;
434      G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
435      G4cout << "-----------------------------------------------------------" << G4endl;
436      G4cout << "Outgoing primary energy: " << residualPrimaryEnergy/keV << " keV" << G4endl;
437      G4cout << "Bremsstrahlung photon " << gammaEnergy/keV << " keV" << G4endl;
438      G4cout << "Total final state: " << (residualPrimaryEnergy+gammaEnergy)/keV
439             << " keV" << G4endl;
440      G4cout << "-----------------------------------------------------------" << G4endl;
441    }
442  if (verboseLevel > 0)
443    {
444      G4double energyDiff = std::fabs(residualPrimaryEnergy+gammaEnergy-kineticEnergy);
445      if (energyDiff > 0.05*keV)
446        G4cout << "Warning from G4PenelopeBremsstrahlung: problem with energy conservation: " <<
447          (residualPrimaryEnergy+gammaEnergy)/keV <<
448          " keV (final) vs. " <<
449          kineticEnergy/keV << " keV (initial)" << G4endl;
450    }
451}
452
453//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
454
455G4PenelopeBremsstrahlungAngular* G4PenelopeBremsstrahlungModel::GetAngularDataForZ(G4int iZ)
456{ 
457  if (!angularData)
458    angularData = new std::map<G4int,G4PenelopeBremsstrahlungAngular*>;
459
460  if (angularData->count(iZ)) //the material already exists
461    return angularData->find(iZ)->second;
462
463  //Otherwise create a new object, store it and return it
464  G4PenelopeBremsstrahlungAngular* theAngular = new G4PenelopeBremsstrahlungAngular(iZ);
465  angularData->insert(std::make_pair(iZ,theAngular));
466
467  if (angularData->count(iZ)) //the material should exist now
468    return angularData->find(iZ)->second;
469  else
470    {
471      G4Exception("Problem in G4PenelopeBremsstrahlungModel::GetAngularDataForZ()");
472      return 0;
473    }
474}
475
476//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
477
478G4PenelopeBremsstrahlungContinuous* 
479G4PenelopeBremsstrahlungModel::GetStoppingPowerData(G4int iZ,G4double energyCut,
480                                                    const G4ParticleDefinition* 
481                                                    theParticle)
482{ 
483  if (!stoppingPowerData)
484    stoppingPowerData = new std::map<std::pair<G4int,G4double>,G4PenelopeBremsstrahlungContinuous*>;
485
486  std::pair<G4int,G4double> theKey = std::make_pair(iZ,energyCut);
487
488  if (stoppingPowerData->count(theKey)) //the material already exists
489    return stoppingPowerData->find(theKey)->second;
490
491  //Otherwise create a new object, store it and return it
492  G4String theParticleName = theParticle->GetParticleName();
493  G4PenelopeBremsstrahlungContinuous* theContinuous = new 
494    G4PenelopeBremsstrahlungContinuous(iZ,energyCut,LowEnergyLimit(),HighEnergyLimit(),theParticleName);
495  stoppingPowerData->insert(std::make_pair(theKey,theContinuous));
496
497  if (stoppingPowerData->count(theKey)) //the material should exist now
498    return stoppingPowerData->find(theKey)->second;
499  else
500    {
501      G4Exception("Problem in G4PenelopeBremsstrahlungModel::GetStoppingPowerData()");
502      return 0;
503    }
504}
Note: See TracBrowser for help on using the repository browser.