Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LowEnergyBremsstrahlung.cc
Go to the documentation of this file.
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$
27 // GEANT4 tag $Name: geant4-09-01-ref-00 $
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 
70 #include "G4RDModifiedTsai.hh"
71 #include "G4RDGenerator2BS.hh"
72 #include "G4RDGenerator2BN.hh"
73 #include "G4RDVDataSetAlgorithm.hh"
75 #include "G4RDVEMDataSet.hh"
76 #include "G4EnergyLossTables.hh"
77 #include "G4PhysicalConstants.hh"
78 #include "G4SystemOfUnits.hh"
79 #include "G4UnitsTable.hh"
80 #include "G4Electron.hh"
81 #include "G4Gamma.hh"
82 #include "G4ProductionCutsTable.hh"
83 
84 
86  : G4eLowEnergyLoss(nam),
87  crossSectionHandler(0),
88  theMeanFreePath(0),
89  energySpectrum(0)
90 {
91  cutForPhotons = 0.;
92  verboseLevel = 0;
93  generatorName = "tsai";
94  angularDistribution = new G4RDModifiedTsai("TsaiGenerator"); // default generator
95 // angularDistribution->PrintGeneratorInformation();
96  TsaiAngularDistribution = new G4RDModifiedTsai("TsaiGenerator");
97 }
98 
99 /*
100 G4LowEnergyBremsstrahlung::G4LowEnergyBremsstrahlung(const G4String& nam, G4RDVBremAngularDistribution* distribution)
101  : G4eLowEnergyLoss(nam),
102  crossSectionHandler(0),
103  theMeanFreePath(0),
104  energySpectrum(0),
105  angularDistribution(distribution)
106 {
107  cutForPhotons = 0.;
108  verboseLevel = 0;
109 
110  angularDistribution->PrintGeneratorInformation();
111 
112  TsaiAngularDistribution = new G4RDModifiedTsai("TsaiGenerator");
113 }
114 */
115 
117 {
118  if(crossSectionHandler) delete crossSectionHandler;
119  if(energySpectrum) delete energySpectrum;
120  if(theMeanFreePath) delete theMeanFreePath;
121  delete angularDistribution;
122  delete TsaiAngularDistribution;
123  energyBins.clear();
124 }
125 
126 
128 {
129  if(verboseLevel > 0) {
130  G4cout << "G4LowEnergyBremsstrahlung::BuildPhysicsTable start"
131  << G4endl;
132  }
133 
134  cutForSecondaryPhotons.clear();
135 
136  // Create and fill BremsstrahlungParameters once
137  if( energySpectrum != 0 ) delete energySpectrum;
138  energyBins.clear();
139  for(size_t i=0; i<15; i++) {
140  G4double x = 0.1*((G4double)i);
141  if(i == 0) x = 0.01;
142  if(i == 10) x = 0.95;
143  if(i == 11) x = 0.97;
144  if(i == 12) x = 0.99;
145  if(i == 13) x = 0.995;
146  if(i == 14) x = 1.0;
147  energyBins.push_back(x);
148  }
149  const G4String dataName("/brem/br-sp.dat");
150  energySpectrum = new G4RDeBremsstrahlungSpectrum(energyBins,dataName);
151 
152  if(verboseLevel > 0) {
153  G4cout << "G4LowEnergyBremsstrahlungSpectrum is initialized"
154  << G4endl;
155  }
156 
157  // Create and fill G4RDCrossSectionHandler once
158 
159  if( crossSectionHandler != 0 ) delete crossSectionHandler;
160  G4RDVDataSetAlgorithm* interpolation = new G4RDLogLogInterpolation();
161  G4double lowKineticEnergy = GetLowerBoundEloss();
162  G4double highKineticEnergy = GetUpperBoundEloss();
163  G4int totBin = GetNbinEloss();
164  crossSectionHandler = new G4RDBremsstrahlungCrossSectionHandler(energySpectrum, interpolation);
165  crossSectionHandler->Initialise(0,lowKineticEnergy, highKineticEnergy, totBin);
166  crossSectionHandler->LoadShellData("brem/br-cs-");
167 
168  if (verboseLevel > 0) {
170  << " is created; Cross section data: "
171  << G4endl;
172  crossSectionHandler->PrintData();
173  G4cout << "Parameters: "
174  << G4endl;
175  energySpectrum->PrintData();
176  }
177 
178  // Build loss table for Bremsstrahlung
179 
180  BuildLossTable(aParticleType);
181 
182  if(verboseLevel > 0) {
183  G4cout << "The loss table is built"
184  << G4endl;
185  }
186 
187  if (&aParticleType==G4Electron::Electron()) {
188 
189  RecorderOfElectronProcess[CounterOfElectronProcess] = (*this).theLossTable;
192 
193  } else {
194 
195  RecorderOfPositronProcess[CounterOfPositronProcess] = (*this).theLossTable;
197  }
198 
199  // Build mean free path data using cut values
200 
201  if( theMeanFreePath != 0 ) delete theMeanFreePath;
202  theMeanFreePath = crossSectionHandler->
203  BuildMeanFreePathForMaterials(&cutForSecondaryPhotons);
204 
205  if(verboseLevel > 0) {
206  G4cout << "The MeanFreePath table is built"
207  << G4endl;
208  }
209 
210  // Build common DEDX table for all ionisation processes
211 
212  BuildDEDXTable(aParticleType);
213 
214  if(verboseLevel > 0) {
215  G4cout << "G4LowEnergyBremsstrahlung::BuildPhysicsTable end"
216  << G4endl;
217  }
218 
219 }
220 
221 
222 void G4LowEnergyBremsstrahlung::BuildLossTable(const G4ParticleDefinition& )
223 {
224  // Build table for energy loss due to soft brems
225  // the tables are built for *MATERIALS* binning is taken from LowEnergyLoss
226 
227  G4double lowKineticEnergy = GetLowerBoundEloss();
228  G4double highKineticEnergy = GetUpperBoundEloss();
229  size_t totBin = GetNbinEloss();
230 
231  // create table
232 
233  if (theLossTable) {
235  delete theLossTable;
236  }
237  const G4ProductionCutsTable* theCoupleTable=
239  size_t numOfCouples = theCoupleTable->GetTableSize();
240  theLossTable = new G4PhysicsTable(numOfCouples);
241 
242  // Clean up the vector of cuts
243  cutForSecondaryPhotons.clear();
244 
245  // Loop for materials
246 
247  for (size_t j=0; j<numOfCouples; j++) {
248 
249  // create physics vector and fill it
250  G4PhysicsLogVector* aVector = new G4PhysicsLogVector(lowKineticEnergy,
251  highKineticEnergy,
252  totBin);
253 
254  // get material parameters needed for the energy loss calculation
255  const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
256  const G4Material* material= couple->GetMaterial();
257 
258  // the cut cannot be below lowest limit
259  G4double tCut = (*(theCoupleTable->GetEnergyCutsVector(0)))[j];
260  tCut = std::min(highKineticEnergy, tCut);
261  cutForSecondaryPhotons.push_back(tCut);
262 
263  const G4ElementVector* theElementVector = material->GetElementVector();
264  size_t NumberOfElements = material->GetNumberOfElements() ;
265  const G4double* theAtomicNumDensityVector =
266  material->GetAtomicNumDensityVector();
267  if(verboseLevel > 1) {
268  G4cout << "Energy loss for material # " << j
269  << " tCut(keV)= " << tCut/keV
270  << G4endl;
271  }
272 
273  // now comes the loop for the kinetic energy values
274  for (size_t i = 0; i<totBin; i++) {
275 
276  G4double lowEdgeEnergy = aVector->GetLowEdgeEnergy(i);
277  G4double ionloss = 0.;
278 
279  // loop for elements in the material
280  for (size_t iel=0; iel<NumberOfElements; iel++ ) {
281  G4int Z = (G4int)((*theElementVector)[iel]->GetZ());
282  G4double e = energySpectrum->AverageEnergy(Z, 0.0, tCut, lowEdgeEnergy);
283  G4double cs= crossSectionHandler->FindValue(Z, lowEdgeEnergy);
284  ionloss += e * cs * theAtomicNumDensityVector[iel];
285  if(verboseLevel > 1) {
286  G4cout << "Z= " << Z
287  << "; tCut(keV)= " << tCut/keV
288  << "; E(keV)= " << lowEdgeEnergy/keV
289  << "; Eav(keV)= " << e/keV
290  << "; cs= " << cs
291  << "; loss= " << ionloss
292  << G4endl;
293  }
294  }
295  aVector->PutValue(i,ionloss);
296  }
297  theLossTable->insert(aVector);
298  }
299 }
300 
301 
303  const G4Step& step)
304 {
306 
307  const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
308  G4double kineticEnergy = track.GetKineticEnergy();
309  G4int index = couple->GetIndex();
310  G4double tCut = cutForSecondaryPhotons[index];
311 
312  // Control limits
313  if(tCut >= kineticEnergy)
314  return G4VContinuousDiscreteProcess::PostStepDoIt(track, step);
315 
316  G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy);
317 
318  G4double tGamma = energySpectrum->SampleEnergy(Z, tCut, kineticEnergy, kineticEnergy);
319  G4double totalEnergy = kineticEnergy + electron_mass_c2;
320  G4double finalEnergy = kineticEnergy - tGamma; // electron/positron final energy
321  G4double theta = 0;
322 
323  if((kineticEnergy < 1*MeV && kineticEnergy > 1*keV && generatorName == "2bn")){
324  theta = angularDistribution->PolarAngle(kineticEnergy,finalEnergy,Z);
325  }else{
326  theta = TsaiAngularDistribution->PolarAngle(kineticEnergy,finalEnergy,Z);
327  }
328 
329  G4double phi = twopi * G4UniformRand();
330  G4double dirZ = std::cos(theta);
331  G4double sinTheta = std::sqrt(1. - dirZ*dirZ);
332  G4double dirX = sinTheta*std::cos(phi);
333  G4double dirY = sinTheta*std::sin(phi);
334 
335  G4ThreeVector gammaDirection (dirX, dirY, dirZ);
336  G4ThreeVector electronDirection = track.GetMomentumDirection();
337 
338  //
339  // Update the incident particle
340  //
341  gammaDirection.rotateUz(electronDirection);
342 
343  // Kinematic problem
344  if (finalEnergy < 0.) {
345  tGamma += finalEnergy;
346  finalEnergy = 0.0;
347  }
348 
349  G4double momentum = std::sqrt((totalEnergy + electron_mass_c2)*kineticEnergy);
350 
351  G4double finalX = momentum*electronDirection.x() - tGamma*gammaDirection.x();
352  G4double finalY = momentum*electronDirection.y() - tGamma*gammaDirection.y();
353  G4double finalZ = momentum*electronDirection.z() - tGamma*gammaDirection.z();
354 
356  G4double norm = 1./std::sqrt(finalX*finalX + finalY*finalY + finalZ*finalZ);
357  aParticleChange.ProposeMomentumDirection(finalX*norm, finalY*norm, finalZ*norm);
358  aParticleChange.ProposeEnergy( finalEnergy );
359 
360  // create G4DynamicParticle object for the gamma
362  gammaDirection, tGamma);
364 
365  return G4VContinuousDiscreteProcess::PostStepDoIt(track, step);
366 }
367 
368 
370 {
371  G4String comments = "Total cross sections from EEDL database.";
372  comments += "\n Gamma energy sampled from a parameterised formula.";
373  comments += "\n Implementation of the continuous dE/dx part.";
374  comments += "\n At present it can be used for electrons ";
375  comments += "in the energy range [250eV,100GeV].";
376  comments += "\n The process must work with G4LowEnergyIonisation.";
377 
378  G4cout << G4endl << GetProcessName() << ": " << comments << G4endl;
379 }
380 
382 {
383  return ( (&particle == G4Electron::Electron()) );
384 }
385 
386 
388  G4double,
389  G4ForceCondition* cond)
390 {
391  *cond = NotForced;
392  G4int index = (track.GetMaterialCutsCouple())->GetIndex();
393  const G4RDVEMDataSet* data = theMeanFreePath->GetComponent(index);
394  G4double meanFreePath = data->FindValue(track.GetKineticEnergy());
395  return meanFreePath;
396 }
397 
399 {
400  cutForPhotons = cut;
401 }
402 
404 {
405  angularDistribution = distribution;
406  angularDistribution->PrintGeneratorInformation();
407 }
408 
410 {
411  if (name == "tsai")
412  {
413  delete angularDistribution;
414  angularDistribution = new G4RDModifiedTsai("TsaiGenerator");
415  generatorName = name;
416  }
417  else if (name == "2bn")
418  {
419  delete angularDistribution;
420  angularDistribution = new G4RDGenerator2BN("2BNGenerator");
421  generatorName = name;
422  }
423  else if (name == "2bs")
424  {
425  delete angularDistribution;
426  angularDistribution = new G4RDGenerator2BS("2BSGenerator");
427  generatorName = name;
428  }
429  else
430  {
431  G4Exception("G4LowEnergyBremsstrahlung::SetAngularGenerator()",
432  "InvalidSetup", FatalException, "Generator does not exist!");
433  }
434 
435  angularDistribution->PrintGeneratorInformation();
436 }
437 
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
G4int verboseLevel
Definition: G4VProcess.hh:368
std::vector< G4Element * > G4ElementVector
virtual G4double PolarAngle(const G4double initial_energy, const G4double final_energy, const G4int Z)=0
void insert(G4PhysicsVector *)
double x() const
static G4int CounterOfElectronProcess
void Initialise(G4RDVDataSetAlgorithm *interpolation=0, G4double minE=250 *CLHEP::eV, G4double maxE=100 *CLHEP::GeV, G4int numberOfBins=200, G4double unitE=CLHEP::MeV, G4double unitData=CLHEP::barn, G4int minZ=1, G4int maxZ=99)
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
void LoadShellData(const G4String &dataFile)
const XML_Char * name
virtual const G4RDVEMDataSet * GetComponent(G4int componentId) const =0
static G4double GetLowerBoundEloss()
G4double FindValue(G4int Z, G4double e) const
G4double GetLowEdgeEnergy(size_t binNumber) const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:188
int G4int
Definition: G4Types.hh:78
double z() const
string material
Definition: eplot.py:19
G4double GetKineticEnergy() const
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
void PutValue(size_t index, G4double theValue)
void BuildDEDXTable(const G4ParticleDefinition &aParticleType)
float electron_mass_c2
Definition: hepunit.py:274
Definition: G4Step.hh:76
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
virtual void PrintData() const =0
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
const G4double * GetAtomicNumDensityVector() const
Definition: G4Material.hh:214
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4PhysicsTable ** RecorderOfPositronProcess
virtual void Initialize(const G4Track &)
static G4ProductionCutsTable * GetProductionCutsTable()
static G4double GetUpperBoundEloss()
virtual G4double AverageEnergy(G4int Z, G4double minKineticEnergy, G4double maxKineticEnergy, G4double kineticEnergy, G4int shell=0, const G4ParticleDefinition *pd=0) const =0
static G4PhysicsTable ** RecorderOfElectronProcess
const G4ThreeVector & GetMomentumDirection() const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
void SetNumberOfSecondaries(G4int totSecondaries)
static G4int GetNbinEloss()
void SetAngularGenerator(G4RDVBremAngularDistribution *distribution)
G4LowEnergyBremsstrahlung(const G4String &processName="LowEnBrem")
double y() const
G4int SelectRandomAtom(const G4MaterialCutsCouple *couple, G4double e) const
void ProposeEnergy(G4double finalEnergy)
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void AddSecondary(G4Track *aSecondary)
void BuildPhysicsTable(const G4ParticleDefinition &particleType)
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
size_t GetNumberOfElements() const
Definition: G4Material.hh:184
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
G4PhysicsTable * theLossTable
static G4int CounterOfPositronProcess
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
double G4double
Definition: G4Types.hh:76
G4double GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
G4ForceCondition
G4VParticleChange * PostStepDoIt(const G4Track &track, const G4Step &step)
virtual G4double SampleEnergy(G4int Z, G4double minKineticEnergy, G4double maxKineticEnergy, G4double kineticEnergy, G4int shell=0, const G4ParticleDefinition *pd=0) const =0
const XML_Char const XML_Char * data
void clearAndDestroy()
const G4Material * GetMaterial() const
G4bool IsApplicable(const G4ParticleDefinition &)
virtual void PrintGeneratorInformation() const =0