source: trunk/source/processes/electromagnetic/lowenergy/src/G4LowEnergyBremsstrahlung.cc@ 1201

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

update CVS release candidate geant4.9.3.01
File size: 15.5 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: G4LowEnergyBremsstrahlung.cc,v 1.74 2009/06/11 15:47:08 mantero Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
28//
29// --------------------------------------------------------------
30//
31// File name: G4LowEnergyBremsstrahlung
32//
33// Author: Alessandra Forti, Vladimir Ivanchenko
34//
35// Creation date: March 1999
36//
37// Modifications:
38// 18.04.2000 V.L.
39// - First implementation of continuous energy loss.
40// 17.02.2000 Veronique Lefebure
41// - correct bug : the gamma energy was not deposited when the gamma was
42// not produced when its energy was < cutForLowEnergySecondaryPhotons
43//
44// Added Livermore data table construction methods A. Forti
45// Modified BuildMeanFreePath to read new data tables A. Forti
46// Modified PostStepDoIt to insert sampling with with EEDL data A. Forti
47// Added SelectRandomAtom A. Forti
48// Added map of the elements A. Forti
49// 20.09.00 update printout V.Ivanchenko
50// 24.04.01 V.Ivanchenko remove RogueWave
51// 29.09.2001 V.Ivanchenko: major revision based on design iteration
52// 10.10.2001 MGP Revision to improve code quality and consistency with design
53// 18.10.2001 MGP Revision to improve code quality
54// 28.10.2001 VI Update printout
55// 29.11.2001 VI New parametrisation
56// 30.07.2002 VI Fix in restricted energy loss
57// 21.01.2003 VI Cut per region
58// 21.02.2003 V.Ivanchenko Energy bins for spectrum are defined here
59// 28.02.03 V.Ivanchenko Filename is defined in the constructor
60// 24.03.2003 P.Rodrigues Changes to accommodate new angular generators
61// 20.05.2003 MGP Removed memory leak related to angularDistribution
62// 06.11.2003 MGP Improved user interface to select angular distribution model
63//
64// --------------------------------------------------------------
65
66#include "G4LowEnergyBremsstrahlung.hh"
67#include "G4eBremsstrahlungSpectrum.hh"
68#include "G4BremsstrahlungCrossSectionHandler.hh"
69#include "G4VBremAngularDistribution.hh"
70#include "G4ModifiedTsai.hh"
71#include "G4Generator2BS.hh"
72#include "G4Generator2BN.hh"
73#include "G4VDataSetAlgorithm.hh"
74#include "G4LogLogInterpolation.hh"
75#include "G4VEMDataSet.hh"
76#include "G4EnergyLossTables.hh"
77#include "G4UnitsTable.hh"
78#include "G4Electron.hh"
79#include "G4Gamma.hh"
80#include "G4ProductionCutsTable.hh"
81
82
83G4LowEnergyBremsstrahlung::G4LowEnergyBremsstrahlung(const G4String& nam)
84 : G4eLowEnergyLoss(nam),
85 crossSectionHandler(0),
86 theMeanFreePath(0),
87 energySpectrum(0)
88{
89 cutForPhotons = 0.;
90 verboseLevel = 0;
91 generatorName = "tsai";
92 angularDistribution = new G4ModifiedTsai("TsaiGenerator"); // default generator
93// angularDistribution->PrintGeneratorInformation();
94 TsaiAngularDistribution = new G4ModifiedTsai("TsaiGenerator");
95
96
97 G4cout << G4endl;
98 G4cout << "*******************************************************************************" << G4endl;
99 G4cout << "*******************************************************************************" << G4endl;
100 G4cout << " The class G4LowEnergyBremsstrahlung is NOT SUPPORTED ANYMORE. " << G4endl;
101 G4cout << " It will be REMOVED with the next major release of Geant4. " << G4endl;
102 G4cout << " Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
103 G4cout << "*******************************************************************************" << G4endl;
104 G4cout << "*******************************************************************************" << G4endl;
105 G4cout << G4endl;
106}
107
108/*
109G4LowEnergyBremsstrahlung::G4LowEnergyBremsstrahlung(const G4String& nam, G4VBremAngularDistribution* distribution)
110 : G4eLowEnergyLoss(nam),
111 crossSectionHandler(0),
112 theMeanFreePath(0),
113 energySpectrum(0),
114 angularDistribution(distribution)
115{
116 cutForPhotons = 0.;
117 verboseLevel = 0;
118
119 angularDistribution->PrintGeneratorInformation();
120
121 TsaiAngularDistribution = new G4ModifiedTsai("TsaiGenerator");
122}
123*/
124
125G4LowEnergyBremsstrahlung::~G4LowEnergyBremsstrahlung()
126{
127 if(crossSectionHandler) delete crossSectionHandler;
128 if(energySpectrum) delete energySpectrum;
129 if(theMeanFreePath) delete theMeanFreePath;
130 delete angularDistribution;
131 delete TsaiAngularDistribution;
132 energyBins.clear();
133}
134
135
136void G4LowEnergyBremsstrahlung::BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
137{
138 if(verboseLevel > 0) {
139 G4cout << "G4LowEnergyBremsstrahlung::BuildPhysicsTable start"
140 << G4endl;
141 }
142
143 cutForSecondaryPhotons.clear();
144
145 // Create and fill BremsstrahlungParameters once
146 if( energySpectrum != 0 ) delete energySpectrum;
147 energyBins.clear();
148 for(size_t i=0; i<15; i++) {
149 G4double x = 0.1*((G4double)i);
150 if(i == 0) x = 0.01;
151 if(i == 10) x = 0.95;
152 if(i == 11) x = 0.97;
153 if(i == 12) x = 0.99;
154 if(i == 13) x = 0.995;
155 if(i == 14) x = 1.0;
156 energyBins.push_back(x);
157 }
158 const G4String dataName("/brem/br-sp.dat");
159 energySpectrum = new G4eBremsstrahlungSpectrum(energyBins,dataName);
160
161 if(verboseLevel > 0) {
162 G4cout << "G4LowEnergyBremsstrahlungSpectrum is initialized"
163 << G4endl;
164 }
165
166 // Create and fill G4CrossSectionHandler once
167
168 if( crossSectionHandler != 0 ) delete crossSectionHandler;
169 G4VDataSetAlgorithm* interpolation = new G4LogLogInterpolation();
170 G4double lowKineticEnergy = GetLowerBoundEloss();
171 G4double highKineticEnergy = GetUpperBoundEloss();
172 G4int totBin = GetNbinEloss();
173 crossSectionHandler = new G4BremsstrahlungCrossSectionHandler(energySpectrum, interpolation);
174 crossSectionHandler->Initialise(0,lowKineticEnergy, highKineticEnergy, totBin);
175 crossSectionHandler->LoadShellData("brem/br-cs-");
176
177 if (verboseLevel > 0) {
178 G4cout << GetProcessName()
179 << " is created; Cross section data: "
180 << G4endl;
181 crossSectionHandler->PrintData();
182 G4cout << "Parameters: "
183 << G4endl;
184 energySpectrum->PrintData();
185 }
186
187 // Build loss table for Bremsstrahlung
188
189 BuildLossTable(aParticleType);
190
191 if(verboseLevel > 0) {
192 G4cout << "The loss table is built"
193 << G4endl;
194 }
195
196 if (&aParticleType==G4Electron::Electron()) {
197
198 RecorderOfElectronProcess[CounterOfElectronProcess] = (*this).theLossTable;
199 CounterOfElectronProcess++;
200 PrintInfoDefinition();
201
202 } else {
203
204 RecorderOfPositronProcess[CounterOfPositronProcess] = (*this).theLossTable;
205 CounterOfPositronProcess++;
206 }
207
208 // Build mean free path data using cut values
209
210 if( theMeanFreePath != 0 ) delete theMeanFreePath;
211 theMeanFreePath = crossSectionHandler->
212 BuildMeanFreePathForMaterials(&cutForSecondaryPhotons);
213
214 if(verboseLevel > 0) {
215 G4cout << "The MeanFreePath table is built"
216 << G4endl;
217 }
218
219 // Build common DEDX table for all ionisation processes
220
221 BuildDEDXTable(aParticleType);
222
223 if(verboseLevel > 0) {
224 G4cout << "G4LowEnergyBremsstrahlung::BuildPhysicsTable end"
225 << G4endl;
226 }
227
228}
229
230
231void G4LowEnergyBremsstrahlung::BuildLossTable(const G4ParticleDefinition& )
232{
233 // Build table for energy loss due to soft brems
234 // the tables are built for *MATERIALS* binning is taken from LowEnergyLoss
235
236 G4double lowKineticEnergy = GetLowerBoundEloss();
237 G4double highKineticEnergy = GetUpperBoundEloss();
238 size_t totBin = GetNbinEloss();
239
240 // create table
241
242 if (theLossTable) {
243 theLossTable->clearAndDestroy();
244 delete theLossTable;
245 }
246 const G4ProductionCutsTable* theCoupleTable=
247 G4ProductionCutsTable::GetProductionCutsTable();
248 size_t numOfCouples = theCoupleTable->GetTableSize();
249 theLossTable = new G4PhysicsTable(numOfCouples);
250
251 // Clean up the vector of cuts
252 cutForSecondaryPhotons.clear();
253
254 // Loop for materials
255
256 for (size_t j=0; j<numOfCouples; j++) {
257
258 // create physics vector and fill it
259 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(lowKineticEnergy,
260 highKineticEnergy,
261 totBin);
262
263 // get material parameters needed for the energy loss calculation
264 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
265 const G4Material* material= couple->GetMaterial();
266
267 // the cut cannot be below lowest limit
268 G4double tCut = (*(theCoupleTable->GetEnergyCutsVector(0)))[j];
269 tCut = std::min(highKineticEnergy, tCut);
270 cutForSecondaryPhotons.push_back(tCut);
271
272 const G4ElementVector* theElementVector = material->GetElementVector();
273 size_t NumberOfElements = material->GetNumberOfElements() ;
274 const G4double* theAtomicNumDensityVector =
275 material->GetAtomicNumDensityVector();
276 if(verboseLevel > 1) {
277 G4cout << "Energy loss for material # " << j
278 << " tCut(keV)= " << tCut/keV
279 << G4endl;
280 }
281
282 // now comes the loop for the kinetic energy values
283 for (size_t i = 0; i<totBin; i++) {
284
285 G4double lowEdgeEnergy = aVector->GetLowEdgeEnergy(i);
286 G4double ionloss = 0.;
287
288 // loop for elements in the material
289 for (size_t iel=0; iel<NumberOfElements; iel++ ) {
290 G4int Z = (G4int)((*theElementVector)[iel]->GetZ());
291 G4double e = energySpectrum->AverageEnergy(Z, 0.0, tCut, lowEdgeEnergy);
292 G4double cs= crossSectionHandler->FindValue(Z, lowEdgeEnergy);
293 ionloss += e * cs * theAtomicNumDensityVector[iel];
294 if(verboseLevel > 1) {
295 G4cout << "Z= " << Z
296 << "; tCut(keV)= " << tCut/keV
297 << "; E(keV)= " << lowEdgeEnergy/keV
298 << "; Eav(keV)= " << e/keV
299 << "; cs= " << cs
300 << "; loss= " << ionloss
301 << G4endl;
302 }
303 }
304 aVector->PutValue(i,ionloss);
305 }
306 theLossTable->insert(aVector);
307 }
308}
309
310
311G4VParticleChange* G4LowEnergyBremsstrahlung::PostStepDoIt(const G4Track& track,
312 const G4Step& step)
313{
314 aParticleChange.Initialize(track);
315
316 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
317 G4double kineticEnergy = track.GetKineticEnergy();
318 G4int index = couple->GetIndex();
319 G4double tCut = cutForSecondaryPhotons[index];
320
321 // Control limits
322 if(tCut >= kineticEnergy)
323 return G4VContinuousDiscreteProcess::PostStepDoIt(track, step);
324
325 G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy);
326
327 G4double tGamma = energySpectrum->SampleEnergy(Z, tCut, kineticEnergy, kineticEnergy);
328 G4double totalEnergy = kineticEnergy + electron_mass_c2;
329 G4double finalEnergy = kineticEnergy - tGamma; // electron/positron final energy
330 G4double theta = 0;
331
332 if((kineticEnergy < 1*MeV && kineticEnergy > 1*keV && generatorName == "2bn")){
333 theta = angularDistribution->PolarAngle(kineticEnergy,finalEnergy,Z);
334 }else{
335 theta = TsaiAngularDistribution->PolarAngle(kineticEnergy,finalEnergy,Z);
336 }
337
338 G4double phi = twopi * G4UniformRand();
339 G4double dirZ = std::cos(theta);
340 G4double sinTheta = std::sqrt(1. - dirZ*dirZ);
341 G4double dirX = sinTheta*std::cos(phi);
342 G4double dirY = sinTheta*std::sin(phi);
343
344 G4ThreeVector gammaDirection (dirX, dirY, dirZ);
345 G4ThreeVector electronDirection = track.GetMomentumDirection();
346
347 //
348 // Update the incident particle
349 //
350 gammaDirection.rotateUz(electronDirection);
351
352 // Kinematic problem
353 if (finalEnergy < 0.) {
354 tGamma += finalEnergy;
355 finalEnergy = 0.0;
356 }
357
358 G4double momentum = std::sqrt((totalEnergy + electron_mass_c2)*kineticEnergy);
359
360 G4double finalX = momentum*electronDirection.x() - tGamma*gammaDirection.x();
361 G4double finalY = momentum*electronDirection.y() - tGamma*gammaDirection.y();
362 G4double finalZ = momentum*electronDirection.z() - tGamma*gammaDirection.z();
363
364 aParticleChange.SetNumberOfSecondaries(1);
365 G4double norm = 1./std::sqrt(finalX*finalX + finalY*finalY + finalZ*finalZ);
366 aParticleChange.ProposeMomentumDirection(finalX*norm, finalY*norm, finalZ*norm);
367 aParticleChange.ProposeEnergy( finalEnergy );
368
369 // create G4DynamicParticle object for the gamma
370 G4DynamicParticle* aGamma= new G4DynamicParticle (G4Gamma::Gamma(),
371 gammaDirection, tGamma);
372 aParticleChange.AddSecondary(aGamma);
373
374 return G4VContinuousDiscreteProcess::PostStepDoIt(track, step);
375}
376
377
378void G4LowEnergyBremsstrahlung::PrintInfoDefinition()
379{
380 G4String comments = "Total cross sections from EEDL database.";
381 comments += "\n Gamma energy sampled from a parameterised formula.";
382 comments += "\n Implementation of the continuous dE/dx part.";
383 comments += "\n At present it can be used for electrons ";
384 comments += "in the energy range [250eV,100GeV].";
385 comments += "\n The process must work with G4LowEnergyIonisation.";
386
387 G4cout << G4endl << GetProcessName() << ": " << comments << G4endl;
388}
389
390G4bool G4LowEnergyBremsstrahlung::IsApplicable(const G4ParticleDefinition& particle)
391{
392 return ( (&particle == G4Electron::Electron()) );
393}
394
395
396G4double G4LowEnergyBremsstrahlung::GetMeanFreePath(const G4Track& track,
397 G4double,
398 G4ForceCondition* cond)
399{
400 *cond = NotForced;
401 G4int index = (track.GetMaterialCutsCouple())->GetIndex();
402 const G4VEMDataSet* data = theMeanFreePath->GetComponent(index);
403 G4double meanFreePath = data->FindValue(track.GetKineticEnergy());
404 return meanFreePath;
405}
406
407void G4LowEnergyBremsstrahlung::SetCutForLowEnSecPhotons(G4double cut)
408{
409 cutForPhotons = cut;
410}
411
412void G4LowEnergyBremsstrahlung::SetAngularGenerator(G4VBremAngularDistribution* distribution)
413{
414 angularDistribution = distribution;
415 angularDistribution->PrintGeneratorInformation();
416}
417
418void G4LowEnergyBremsstrahlung::SetAngularGenerator(const G4String& name)
419{
420 if (name == "tsai")
421 {
422 delete angularDistribution;
423 angularDistribution = new G4ModifiedTsai("TsaiGenerator");
424 generatorName = name;
425 }
426 else if (name == "2bn")
427 {
428 delete angularDistribution;
429 angularDistribution = new G4Generator2BN("2BNGenerator");
430 generatorName = name;
431 }
432 else if (name == "2bs")
433 {
434 delete angularDistribution;
435 angularDistribution = new G4Generator2BS("2BSGenerator");
436 generatorName = name;
437 }
438 else
439 {
440 G4Exception("G4LowEnergyBremsstrahlung::SetAngularGenerator - generator does not exist");
441 }
442
443 angularDistribution->PrintGeneratorInformation();
444}
445
Note: See TracBrowser for help on using the repository browser.