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

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

file release beta
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.4 2009/05/19 14:57:01 pandola Exp $
27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
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 = 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    }
199 
200  //The cut is already included in the crossSectionHandler
201  G4double cs = 
202    crossSectionHandler->GetCrossSectionAboveThresholdForElement(energy,cutEnergy,iZ);
203
204  if (verboseLevel > 1)
205    {
206      G4cout << "G4LivermoreBremsstrahlungModel " << G4endl;
207      G4cout << "Cross section for gamma emission > " << cutEnergy/keV << " keV at " <<
208        energy/keV << " keV and Z = " << iZ << " --> " << cs/barn << " barn" << G4endl;
209    }
210  return cs;
211}
212
213
214//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
215
216G4double G4LivermoreBremsstrahlungModel::ComputeDEDXPerVolume(const G4Material* material,
217                                           const G4ParticleDefinition* ,
218                                           G4double kineticEnergy,
219                                           G4double cutEnergy)
220{
221  G4double sPower = 0.0;
222
223  const G4ElementVector* theElementVector = material->GetElementVector();
224  size_t NumberOfElements = material->GetNumberOfElements() ;
225  const G4double* theAtomicNumDensityVector =
226                    material->GetAtomicNumDensityVector();
227
228  // loop for elements in the material
229  for (size_t iel=0; iel<NumberOfElements; iel++ ) 
230    {
231      G4int iZ = (G4int)((*theElementVector)[iel]->GetZ());
232      G4double e = energySpectrum->AverageEnergy(iZ, 0.0,cutEnergy,
233                                                 kineticEnergy);
234      G4double cs= crossSectionHandler->FindValue(iZ,kineticEnergy);
235      sPower   += e * cs * theAtomicNumDensityVector[iel];
236    }
237
238  if (verboseLevel > 2)
239    {
240      G4cout << "G4LivermoreBremsstrahlungModel " << G4endl;
241      G4cout << "Stopping power < " << cutEnergy/keV << " keV at " << 
242        kineticEnergy/keV << " keV = " << sPower/(keV/mm) << " keV/mm" << G4endl;
243    }
244   
245  return sPower;
246}
247
248//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
249
250void G4LivermoreBremsstrahlungModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
251                                              const G4MaterialCutsCouple* couple,
252                                              const G4DynamicParticle* aDynamicParticle,
253                                              G4double energyCut,
254                                              G4double)
255{
256 
257  G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
258
259  // this is neede for pathalogical cases of no ionisation
260  if (kineticEnergy <= fIntrinsicLowEnergyLimit)
261    {
262      fParticleChange->SetProposedKineticEnergy(0.);
263      fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy);
264      return;
265    }
266
267  //Sample material
268  G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy);
269
270  //Sample gamma energy
271  G4double tGamma = energySpectrum->SampleEnergy(Z, energyCut, kineticEnergy, kineticEnergy);
272  G4double totalEnergy = kineticEnergy + electron_mass_c2;
273  G4double finalEnergy = kineticEnergy - tGamma; // electron final energy 
274  G4double theta = 0;
275
276  if (tGamma == 0.) //nothing happens
277    return;
278
279  //Sample gamma direction
280  //Use alternative algorithms, if it is the case.
281  if((kineticEnergy < MeV && kineticEnergy > keV))
282    { 
283      theta = angularDistribution->PolarAngle(kineticEnergy,finalEnergy,Z);
284    }
285  else
286    //Otherwise, use tsai
287    {
288      theta = TsaiAngularDistribution->PolarAngle(kineticEnergy,finalEnergy,Z);
289    } 
290
291  G4double phi   = twopi * G4UniformRand();
292  G4double dirZ  = std::cos(theta);
293  G4double sinTheta  = std::sqrt(1. - dirZ*dirZ);
294  G4double dirX  = sinTheta*std::cos(phi);
295  G4double dirY  = sinTheta*std::sin(phi);
296
297  G4ThreeVector gammaDirection (dirX, dirY, dirZ);
298  G4ThreeVector electronDirection = aDynamicParticle->GetMomentumDirection();
299
300  //Update the incident particle
301  gammaDirection.rotateUz(electronDirection);   
302   
303  if (finalEnergy < 0.) 
304    {
305      // Kinematic problem
306      tGamma = kineticEnergy;
307      fParticleChange->SetProposedKineticEnergy(0.);
308    }
309  else
310    {
311      G4double momentum = std::sqrt((totalEnergy + electron_mass_c2)*kineticEnergy);
312      G4double finalX = momentum*electronDirection.x() - tGamma*gammaDirection.x();
313      G4double finalY = momentum*electronDirection.y() - tGamma*gammaDirection.y();
314      G4double finalZ = momentum*electronDirection.z() - tGamma*gammaDirection.z();
315      G4double norm = 1./std::sqrt(finalX*finalX + finalY*finalY + finalZ*finalZ);
316     
317      fParticleChange->ProposeMomentumDirection(finalX*norm, finalY*norm, finalZ*norm);
318      fParticleChange->SetProposedKineticEnergy(finalEnergy);
319    }
320
321  //Generate the bremsstrahlung gamma
322  G4DynamicParticle* aGamma= new G4DynamicParticle (G4Gamma::Gamma(),
323                                                    gammaDirection, tGamma);
324  fvect->push_back(aGamma);
325
326  if (verboseLevel > 1)
327    {
328      G4cout << "-----------------------------------------------------------" << G4endl;
329      G4cout << "Energy balance from G4LivermoreBremsstrahlung" << G4endl;
330      G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
331      G4cout << "-----------------------------------------------------------" << G4endl;
332      G4cout << "Outgoing primary energy: " << finalEnergy/keV << " keV" << G4endl;
333      G4cout << "Gamma ray " << tGamma/keV << " keV" << G4endl;
334      G4cout << "Total final state: " << (finalEnergy+tGamma)/keV << " keV" << G4endl;
335      G4cout << "-----------------------------------------------------------" << G4endl;
336    }
337  if (verboseLevel > 0)
338    {
339      G4double energyDiff = std::fabs(finalEnergy+tGamma-kineticEnergy);
340      if (energyDiff > 0.05*keV)
341        G4cout << "G4LivermoreBremsstrahlung WARNING: problem with energy conservation: " 
342               << (finalEnergy+tGamma)/keV << " keV (final) vs. " 
343               << kineticEnergy/keV << " keV (initial)" << G4endl;
344    }
345  return;
346}
347
348//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
349
350void 
351G4LivermoreBremsstrahlungModel::SetAngularGenerator(G4VBremAngularDistribution* distribution)
352{
353  if(angularDistribution == distribution) return;
354  if(angularDistribution) delete angularDistribution;
355  angularDistribution = distribution;
356  angularDistribution->PrintGeneratorInformation();
357}
358
359//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
360
361void G4LivermoreBremsstrahlungModel::SetAngularGenerator(const G4String& name)
362{
363  if(name == generatorName) return;
364  if (name == "tsai") 
365    {
366      delete angularDistribution;
367      angularDistribution = new G4ModifiedTsai("TsaiGenerator");
368      generatorName = name;
369    }
370  else if (name == "2bn")
371    {
372      delete angularDistribution;
373      angularDistribution = new G4Generator2BN("2BNGenerator");
374      generatorName = name;
375    }
376  else if (name == "2bs")
377    {
378      delete angularDistribution;
379      angularDistribution = new G4Generator2BS("2BSGenerator");
380      generatorName = name;
381    }
382  else
383    {
384      G4cout << "### G4LivermoreBremsstrahlungModel::SetAngularGenerator WARNING:"
385             << " generator <" << name << "> is not known" << G4endl;
386      return; 
387
388    }
389
390  angularDistribution->PrintGeneratorInformation();
391}
Note: See TracBrowser for help on using the repository browser.