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

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

geant4 tag 9.4
File size: 13.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: G4LivermoreBremsstrahlungModel.cc,v 1.8 2010/12/02 16:07:05 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-04-ref-00 $
28//
29// Author: Luciano Pandola
30//
31// History:
32// --------
33// 03 Mar 2009   L Pandola    Migration from process to model
34// 12 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//                  - added MinEnergyCut method
38//                  - do not change track status
39//                  - do not initialize element selectors
40//                  - use cut value from the interface
41//                  - fixed bug in sampling of angles between keV and MeV
42// 19 May 2009   L Pandola    Explicitely set to zero pointers deleted in
43//                            Initialise(), since they might be checked later on
44//
45
46#include "G4LivermoreBremsstrahlungModel.hh"
47#include "G4ParticleDefinition.hh"
48#include "G4MaterialCutsCouple.hh"
49
50#include "G4DynamicParticle.hh"
51#include "G4Element.hh"
52#include "G4Gamma.hh"
53#include "G4Electron.hh"
54#include "G4SemiLogInterpolation.hh"
55//
56#include "G4VBremAngularDistribution.hh"
57#include "G4ModifiedTsai.hh"
58#include "G4Generator2BS.hh"
59#include "G4Generator2BN.hh"
60//
61#include "G4BremsstrahlungCrossSectionHandler.hh"
62//
63#include "G4VEnergySpectrum.hh"
64#include "G4eBremsstrahlungSpectrum.hh"
65
66//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
67
68
69G4LivermoreBremsstrahlungModel::G4LivermoreBremsstrahlungModel(const G4ParticleDefinition*,
70                                                               const G4String& nam)
71  :G4VEmModel(nam),isInitialised(false),crossSectionHandler(0),
72   energySpectrum(0)
73{
74  fIntrinsicLowEnergyLimit = 10.0*eV;
75  fIntrinsicHighEnergyLimit = 100.0*GeV;
76  fNBinEnergyLoss = 360;
77  //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
78  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
79  //
80  verboseLevel = 0;
81  //
82  generatorName = "tsai";
83  angularDistribution = new G4ModifiedTsai("TsaiGenerator"); //default generator
84  //
85  TsaiAngularDistribution = new G4ModifiedTsai("TsaiGenerator");
86  //
87}
88
89//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90
91G4LivermoreBremsstrahlungModel::~G4LivermoreBremsstrahlungModel()
92{
93  if (crossSectionHandler) delete crossSectionHandler;
94  if (energySpectrum) delete energySpectrum;
95  energyBins.clear();
96  delete angularDistribution;
97  delete TsaiAngularDistribution;
98}
99
100//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101
102void G4LivermoreBremsstrahlungModel::Initialise(const G4ParticleDefinition* particle,
103                                                const G4DataVector& cuts)
104{
105  //Check that the Livermore Bremsstrahlung is NOT attached to e+
106  if (particle != G4Electron::Electron())
107    {
108      G4cout << "ERROR: Livermore Bremsstrahlung Model is applicable only to electrons" 
109             << G4endl;
110      G4cout << "It cannot be registered to " << particle->GetParticleName() << G4endl;
111      G4Exception();
112    }
113  //Prepare energy spectrum
114  if (energySpectrum) 
115    {
116      delete energySpectrum;
117      energySpectrum = 0;
118    }
119
120  energyBins.clear();
121  for(size_t i=0; i<15; i++) 
122    {
123      G4double x = 0.1*((G4double)i);
124      if(i == 0)  x = 0.01;
125      if(i == 10) x = 0.95;
126      if(i == 11) x = 0.97;
127      if(i == 12) x = 0.99;
128      if(i == 13) x = 0.995;
129      if(i == 14) x = 1.0;
130      energyBins.push_back(x);
131    }
132  const G4String dataName("/brem/br-sp.dat");
133  energySpectrum = new G4eBremsstrahlungSpectrum(energyBins,dataName);
134 
135  if (verboseLevel > 0)
136    G4cout << "G4eBremsstrahlungSpectrum is initialized" << G4endl;
137
138  //Initialize cross section handler
139  if (crossSectionHandler) 
140    {
141      delete crossSectionHandler;
142      crossSectionHandler = 0;
143    }
144  G4VDataSetAlgorithm* interpolation = 0;//new G4SemiLogInterpolation();
145  crossSectionHandler = new G4BremsstrahlungCrossSectionHandler(energySpectrum,interpolation);
146  crossSectionHandler->Initialise(0,LowEnergyLimit(),HighEnergyLimit(),
147                                  fNBinEnergyLoss);
148  crossSectionHandler->Clear();
149  crossSectionHandler->LoadShellData("brem/br-cs-");
150  //This is used to retrieve cross section values later on
151  crossSectionHandler->BuildMeanFreePathForMaterials(&cuts);
152 
153  if (verboseLevel > 0)
154    {
155      G4cout << "Livermore Bremsstrahlung model is initialized " << G4endl
156             << "Energy range: "
157             << LowEnergyLimit() / keV << " keV - "
158             << HighEnergyLimit() / GeV << " GeV"
159             << G4endl;
160    }
161
162  if (verboseLevel > 1)
163    {
164      G4cout << "Cross section data: " << G4endl; 
165      crossSectionHandler->PrintData();
166      G4cout << "Parameters: " << G4endl;
167      energySpectrum->PrintData();
168    }
169
170  if(isInitialised) return;
171  fParticleChange = GetParticleChangeForLoss();
172  isInitialised = true; 
173}
174
175//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
176
177G4double G4LivermoreBremsstrahlungModel::MinEnergyCut(const G4ParticleDefinition*,
178                                                      const G4MaterialCutsCouple*)
179{
180  return 250.*eV;
181}
182
183//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
184
185G4double
186G4LivermoreBremsstrahlungModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
187                                                           G4double energy,
188                                                           G4double Z, G4double,
189                                                           G4double cutEnergy, 
190                                                           G4double)
191{
192  G4int iZ = (G4int) Z;
193  if (!crossSectionHandler)
194    {
195      G4cout << "G4LivermoreBremsstrahlungModel::ComputeCrossSectionPerAtom" << G4endl;
196      G4cout << "The cross section handler is not correctly initialized" << G4endl;
197      G4Exception();
198      return 0;
199    }
200 
201  //The cut is already included in the crossSectionHandler
202  G4double cs = 
203    crossSectionHandler->GetCrossSectionAboveThresholdForElement(energy,cutEnergy,iZ);
204
205  if (verboseLevel > 1)
206    {
207      G4cout << "G4LivermoreBremsstrahlungModel " << G4endl;
208      G4cout << "Cross section for gamma emission > " << cutEnergy/keV << " keV at " <<
209        energy/keV << " keV and Z = " << iZ << " --> " << cs/barn << " barn" << G4endl;
210    }
211  return cs;
212}
213
214
215//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
216
217G4double G4LivermoreBremsstrahlungModel::ComputeDEDXPerVolume(const G4Material* material,
218                                           const G4ParticleDefinition* ,
219                                           G4double kineticEnergy,
220                                           G4double cutEnergy)
221{
222  G4double sPower = 0.0;
223
224  const G4ElementVector* theElementVector = material->GetElementVector();
225  size_t NumberOfElements = material->GetNumberOfElements() ;
226  const G4double* theAtomicNumDensityVector =
227                    material->GetAtomicNumDensityVector();
228
229  // loop for elements in the material
230  for (size_t iel=0; iel<NumberOfElements; iel++ ) 
231    {
232      G4int iZ = (G4int)((*theElementVector)[iel]->GetZ());
233      G4double e = energySpectrum->AverageEnergy(iZ, 0.0,cutEnergy,
234                                                 kineticEnergy);
235      G4double cs= crossSectionHandler->FindValue(iZ,kineticEnergy);
236      sPower   += e * cs * theAtomicNumDensityVector[iel];
237    }
238
239  if (verboseLevel > 2)
240    {
241      G4cout << "G4LivermoreBremsstrahlungModel " << G4endl;
242      G4cout << "Stopping power < " << cutEnergy/keV << " keV at " << 
243        kineticEnergy/keV << " keV = " << sPower/(keV/mm) << " keV/mm" << G4endl;
244    }
245   
246  return sPower;
247}
248
249//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
250
251void G4LivermoreBremsstrahlungModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
252                                              const G4MaterialCutsCouple* couple,
253                                              const G4DynamicParticle* aDynamicParticle,
254                                              G4double energyCut,
255                                              G4double)
256{
257 
258  G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
259
260  // this is neede for pathalogical cases of no ionisation
261  if (kineticEnergy <= fIntrinsicLowEnergyLimit)
262    {
263      fParticleChange->SetProposedKineticEnergy(0.);
264      fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy);
265      return;
266    }
267
268  //Sample material
269  G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy);
270
271  //Sample gamma energy
272  G4double tGamma = energySpectrum->SampleEnergy(Z, energyCut, kineticEnergy, kineticEnergy);
273  G4double totalEnergy = kineticEnergy + electron_mass_c2;
274  G4double finalEnergy = kineticEnergy - tGamma; // electron final energy 
275  G4double theta = 0;
276
277  if (tGamma == 0.) //nothing happens
278    return;
279
280  //Sample gamma direction
281  //Use alternative algorithms, if it is the case.
282  if((kineticEnergy < MeV && kineticEnergy > keV))
283    { 
284      theta = angularDistribution->PolarAngle(kineticEnergy,finalEnergy,Z);
285    }
286  else
287    //Otherwise, use tsai
288    {
289      theta = TsaiAngularDistribution->PolarAngle(kineticEnergy,finalEnergy,Z);
290    } 
291
292  G4double phi   = twopi * G4UniformRand();
293  G4double dirZ  = std::cos(theta);
294  G4double sinTheta  = std::sqrt(1. - dirZ*dirZ);
295  G4double dirX  = sinTheta*std::cos(phi);
296  G4double dirY  = sinTheta*std::sin(phi);
297
298  G4ThreeVector gammaDirection (dirX, dirY, dirZ);
299  G4ThreeVector electronDirection = aDynamicParticle->GetMomentumDirection();
300
301  //Update the incident particle
302  gammaDirection.rotateUz(electronDirection);   
303   
304  if (finalEnergy < 0.) 
305    {
306      // Kinematic problem
307      tGamma = kineticEnergy;
308      fParticleChange->SetProposedKineticEnergy(0.);
309    }
310  else
311    {
312      G4double momentum = std::sqrt((totalEnergy + electron_mass_c2)*kineticEnergy);
313      G4double finalX = momentum*electronDirection.x() - tGamma*gammaDirection.x();
314      G4double finalY = momentum*electronDirection.y() - tGamma*gammaDirection.y();
315      G4double finalZ = momentum*electronDirection.z() - tGamma*gammaDirection.z();
316      G4double norm = 1./std::sqrt(finalX*finalX + finalY*finalY + finalZ*finalZ);
317     
318      fParticleChange->ProposeMomentumDirection(finalX*norm, finalY*norm, finalZ*norm);
319      fParticleChange->SetProposedKineticEnergy(finalEnergy);
320    }
321
322  //Generate the bremsstrahlung gamma
323  G4DynamicParticle* aGamma= new G4DynamicParticle (G4Gamma::Gamma(),
324                                                    gammaDirection, tGamma);
325  fvect->push_back(aGamma);
326
327  if (verboseLevel > 1)
328    {
329      G4cout << "-----------------------------------------------------------" << G4endl;
330      G4cout << "Energy balance from G4LivermoreBremsstrahlung" << G4endl;
331      G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
332      G4cout << "-----------------------------------------------------------" << G4endl;
333      G4cout << "Outgoing primary energy: " << finalEnergy/keV << " keV" << G4endl;
334      G4cout << "Gamma ray " << tGamma/keV << " keV" << G4endl;
335      G4cout << "Total final state: " << (finalEnergy+tGamma)/keV << " keV" << G4endl;
336      G4cout << "-----------------------------------------------------------" << G4endl;
337    }
338  if (verboseLevel > 0)
339    {
340      G4double energyDiff = std::fabs(finalEnergy+tGamma-kineticEnergy);
341      if (energyDiff > 0.05*keV)
342        G4cout << "G4LivermoreBremsstrahlung WARNING: problem with energy conservation: " 
343               << (finalEnergy+tGamma)/keV << " keV (final) vs. " 
344               << kineticEnergy/keV << " keV (initial)" << G4endl;
345    }
346  return;
347}
348
349//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
350
351void 
352G4LivermoreBremsstrahlungModel::SetAngularGenerator(G4VBremAngularDistribution* distribution)
353{
354  if(angularDistribution == distribution) return;
355  if(angularDistribution) delete angularDistribution;
356  angularDistribution = distribution;
357  angularDistribution->PrintGeneratorInformation();
358}
359
360//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
361
362void G4LivermoreBremsstrahlungModel::SetAngularGenerator(const G4String& name)
363{
364  if(name == generatorName) return;
365  if (name == "tsai") 
366    {
367      delete angularDistribution;
368      angularDistribution = new G4ModifiedTsai("TsaiGenerator");
369      generatorName = name;
370    }
371  else if (name == "2bn")
372    {
373      delete angularDistribution;
374      angularDistribution = new G4Generator2BN("2BNGenerator");
375      generatorName = name;
376    }
377  else if (name == "2bs")
378    {
379      delete angularDistribution;
380      angularDistribution = new G4Generator2BS("2BSGenerator");
381      generatorName = name;
382    }
383  else
384    {
385      G4cout << "### G4LivermoreBremsstrahlungModel::SetAngularGenerator WARNING:"
386             << " generator <" << name << "> is not known" << G4endl;
387      return; 
388
389    }
390
391  angularDistribution->PrintGeneratorInformation();
392}
Note: See TracBrowser for help on using the repository browser.