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