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

Last change on this file since 1199 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
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.6 2009/06/11 15:47:08 mantero Exp $
27// GEANT4 tag $Name: geant4-09-03-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.