| Geant4 Cross Reference |
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 /// \file electromagnetic/TestEm7/include/G4ScreenedNuclearRecoil.hh
27 /// \brief Definition of the G4ScreenedNuclearRecoil class
28 //
29 //
30 //
31 // G4ScreenedNuclearRecoil.hh,v 1.24 2008/05/01 19:58:59 marcus Exp
32 // GEANT4 tag
33 //
34 //
35 //
36 // Class Description
37 // Process for screened electromagnetic nuclear elastic scattering;
38 // Physics comes from:
39 // Marcus H. Mendenhall and Robert A. Weller,
40 // "Algorithms for the rapid computation of classical cross sections
41 // for screened Coulomb collisions "
42 // Nuclear Instruments and Methods in Physics Research B58 (1991) 11-17
43 // The only input required is a screening function phi(r/a) which is the ratio
44 // of the actual interatomic potential for two atoms with atomic numbers Z1 and
45 // Z2,
46 // to the unscreened potential Z1*Z2*e^2/r where e^2 is elm_coupling in Geant4
47 // units
48 // the actual screening tables are computed externally in a python module
49 // "screened_scattering.py"
50 // to allow very specific screening functions to be added if desired, without
51 // messing with the insides of this code.
52 //
53 // First version, April 2004, Marcus H. Mendenhall, Vanderbilt University
54 // May 1, 2008 -- Added code to allow process to have zero cross section above
55 // max energy, to coordinate with G4MSC. -- mhm
56 //
57 // Class Description - End
58
59 #ifndef G4ScreenedNuclearRecoil_h
60 #define G4ScreenedNuclearRecoil_h 1
61
62 #include "c2_function.hh"
63
64 #include "G4ParticleChange.hh"
65 #include "G4PhysicalConstants.hh"
66 #include "G4SystemOfUnits.hh"
67 #include "G4VDiscreteProcess.hh"
68 #include "globals.hh"
69
70 #include <map>
71 #include <vector>
72
73 class G4VNIELPartition;
74
75 typedef c2_const_ptr<G4double> G4_c2_const_ptr;
76 typedef c2_ptr<G4double> G4_c2_ptr;
77 typedef c2_function<G4double> G4_c2_function;
78
79 typedef struct G4ScreeningTables
80 {
81 G4double z1, z2, m1, m2, au, emin;
82 G4_c2_const_ptr EMphiData;
83 } G4ScreeningTables;
84
85 // A class for loading ScreenedCoulombCrossSections
86 class G4ScreenedCoulombCrossSectionInfo
87 {
88 public:
89 G4ScreenedCoulombCrossSectionInfo() {}
90 ~G4ScreenedCoulombCrossSectionInfo() {}
91
92 static const char* CVSHeaderVers()
93 {
94 return "G4ScreenedNuclearRecoil.hh,v 1.24 2008/05/01 19:58:59 marcus Exp GEANT4 tag ";
95 }
96 static const char* CVSFileVers();
97 };
98
99 // A class for loading ScreenedCoulombCrossSections
100 class G4ScreenedCoulombCrossSection : public G4ScreenedCoulombCrossSectionInfo
101 {
102 public:
103 G4ScreenedCoulombCrossSection() : verbosity(1) {}
104 G4ScreenedCoulombCrossSection(const G4ScreenedCoulombCrossSection& src)
105 : G4ScreenedCoulombCrossSectionInfo(), verbosity(src.verbosity)
106 {}
107 virtual ~G4ScreenedCoulombCrossSection();
108
109 typedef std::map<G4int, G4ScreeningTables> ScreeningMap;
110
111 // a local, fast-access mapping of a particle's Z to its full definition
112 typedef std::map<G4int, class G4ParticleDefinition*> ParticleCache;
113
114 // LoadData is called by G4ScreenedNuclearRecoil::GetMeanFreePath
115 // It loads the data tables, builds the elemental cross-section tables.
116 virtual void LoadData(G4String screeningKey, G4int z1, G4double m1, G4double recoilCutoff) = 0;
117
118 // BuildMFPTables is called by G4ScreenedNuclearRecoil::GetMeanFreePath
119 // to build the MFP tables for each material
120 void BuildMFPTables(void); // scan the MaterialsTable and construct MFP
121 // tables
122
123 virtual G4ScreenedCoulombCrossSection* create() = 0;
124 // a 'virtual constructor' which clones the class
125 const G4ScreeningTables* GetScreening(G4int Z) { return &(screeningData[Z]); }
126 void SetVerbosity(G4int v) { verbosity = v; }
127
128 // this process needs element selection weighted only by number density
129 G4ParticleDefinition* SelectRandomUnweightedTarget(const G4MaterialCutsCouple* couple);
130
131 enum
132 {
133 nMassMapElements = 116
134 };
135
136 G4double standardmass(G4int z1) { return z1 <= nMassMapElements ? massmap[z1] : 2.5 * z1; }
137
138 // get the mean-free-path table for the indexed material
139 const G4_c2_function* operator[](G4int materialIndex)
140 {
141 return MFPTables.find(materialIndex) != MFPTables.end() ? &(MFPTables[materialIndex].get())
142 : (G4_c2_function*)0;
143 }
144
145 protected:
146 ScreeningMap screeningData; // screening tables for each element
147 ParticleCache targetMap;
148 G4int verbosity;
149 std::map<G4int, G4_c2_const_ptr> sigmaMap;
150 // total cross section for each element
151 std::map<G4int, G4_c2_const_ptr> MFPTables; // MFP for each material
152
153 private:
154 static const G4double massmap[nMassMapElements + 1];
155 };
156
157 typedef struct G4CoulombKinematicsInfo
158 {
159 G4double impactParameter;
160 G4ScreenedCoulombCrossSection* crossSection;
161 G4double a1, a2, sinTheta, cosTheta, sinZeta, cosZeta, eRecoil;
162 G4ParticleDefinition* recoilIon;
163 const G4Material* targetMaterial;
164 } G4CoulombKinematicsInfo;
165
166 class G4ScreenedCollisionStage
167 {
168 public:
169 virtual void DoCollisionStep(class G4ScreenedNuclearRecoil* master, const class G4Track& aTrack,
170 const class G4Step& aStep) = 0;
171 virtual ~G4ScreenedCollisionStage() {}
172 };
173
174 class G4ScreenedCoulombClassicalKinematics : public G4ScreenedCoulombCrossSectionInfo,
175 public G4ScreenedCollisionStage
176 {
177 public:
178 G4ScreenedCoulombClassicalKinematics();
179 virtual void DoCollisionStep(class G4ScreenedNuclearRecoil* master, const class G4Track& aTrack,
180 const class G4Step& aStep);
181
182 G4bool DoScreeningComputation(class G4ScreenedNuclearRecoil* master,
183 const G4ScreeningTables* screen, G4double eps, G4double beta);
184
185 virtual ~G4ScreenedCoulombClassicalKinematics() {}
186
187 protected:
188 // the c2_functions we need to do the work.
189 c2_const_plugin_function_p<G4double>& phifunc;
190 c2_linear_p<G4double>& xovereps;
191 G4_c2_ptr diff;
192 };
193
194 class G4SingleScatter : public G4ScreenedCoulombCrossSectionInfo, public G4ScreenedCollisionStage
195 {
196 public:
197 G4SingleScatter() {}
198 virtual void DoCollisionStep(class G4ScreenedNuclearRecoil* master, const class G4Track& aTrack,
199 const class G4Step& aStep);
200 virtual ~G4SingleScatter() {}
201 };
202
203 /**
204 \brief A process which handles screened Coulomb collisions between nuclei
205
206 */
207
208 class G4ScreenedNuclearRecoil : public G4ScreenedCoulombCrossSectionInfo, public G4VDiscreteProcess
209 {
210 public:
211 friend class G4ScreenedCollisionStage;
212
213 /// \brief Construct the process and set some physics parameters for it.
214 /// \param processName the name to assign the process
215 /// \param ScreeningKey the name of a screening function to use.
216 /// The default functions are "zbl" (recommended for soft scattering),
217 /// "lj" (recommended for backscattering) and "mol" (Moliere potential)
218 /// \param GenerateRecoils if frue, ions struck by primary are converted
219 /// into new moving particles.
220 /// If false, energy is deposited, but no new moving ions are created.
221 /// \param RecoilCutoff energy below which no new moving particles will be
222 /// created, even if a GenerateRecoils is true.
223 /// Also, a moving primary particle will be stopped if its energy falls
224 /// below this limit.
225 /// \param PhysicsCutoff the energy transfer to which screening tables are
226 /// calucalted.
227 /// There is no really
228 /// compelling reason to change it from the 10.0 eV default.
229 /// However, see the paper on running this
230 /// in thin targets for further discussion, and its interaction
231 /// with SetMFPScaling()
232
233 G4ScreenedNuclearRecoil(const G4String& processName = "ScreenedElastic",
234 const G4String& ScreeningKey = "zbl", G4bool GenerateRecoils = 1,
235 G4double RecoilCutoff = 100.0 * eV, G4double PhysicsCutoff = 10.0 * eV);
236
237 /// \brief destructor
238 virtual ~G4ScreenedNuclearRecoil();
239
240 /// \brief used internally by Geant4 machinery
241 virtual G4double GetMeanFreePath(const G4Track&, G4double, G4ForceCondition*);
242
243 /// \brief used internally by Geant4 machinery
244 virtual G4VParticleChange* PostStepDoIt(const G4Track& aTrack, const G4Step& aStep);
245 /// \brief test if a prticle of type \a aParticleType can use this process
246 /// \param aParticleType the particle to test
247 virtual G4bool IsApplicable(const G4ParticleDefinition& aParticleType);
248 /// \brief Build physics tables in advance. Not Implemented.
249 /// \param aParticleType the type of particle to build tables for
250 virtual void BuildPhysicsTable(const G4ParticleDefinition& aParticleType);
251 /// \brief Export physics tables for persistency. Not Implemented.
252 /// \param aParticleType the type of particle to build tables for
253 virtual void DumpPhysicsTable(const G4ParticleDefinition& aParticleType);
254 /// \brief deterine if the moving particle is within the strong force
255 /// range of the selected nucleus
256 /// \param A the nucleon number of the beam
257 /// \param A1 the nucleon number of the target
258 /// \param apsis the distance of closest approach
259
260 virtual G4bool CheckNuclearCollision(G4double A, G4double A1,
261 G4double apsis); // return true if hard collision
262
263 virtual G4ScreenedCoulombCrossSection* GetNewCrossSectionHandler(void);
264
265 /// \brief Get non-ionizing energy loss for last step
266 G4double GetNIEL() const { return NIEL; }
267
268 /// \brief clear precomputed screening tables
269 void ResetTables();
270 // clear all data tables to allow changing energy cutoff, materials, etc.
271
272 /// \brief set the upper energy beyond which this process has no
273 /// cross section
274 ///
275 /// This funciton is used to coordinate this process with G4MSC.
276 /// Typically, G4MSC should
277 /// not be allowed to operate in a range which overlaps that of this
278 /// process. The criterion which is most reasonable
279 /// is that the transition should be somewhere in the modestly
280 /// relativistic regime (500 MeV/u for example).
281 /// param energy energy per nucleon for the cutoff
282
283 void SetMaxEnergyForScattering(G4double energy) { processMaxEnergy = energy; }
284
285 /// \brief find out what screening function we are using
286 std::string GetScreeningKey() const { return screeningKey; }
287
288 /// \brief enable or disable all energy deposition by this process
289 /// \param flag if true, enable deposition of energy (the default).
290 /// If false, disable deposition.
291
292 void AllowEnergyDeposition(G4bool flag) { registerDepositedEnergy = flag; }
293
294 /// \brief get flag indicating whether deposition is enabled
295 G4bool GetAllowEnergyDeposition() const { return registerDepositedEnergy; }
296
297 /// \brief enable or disable the generation of recoils.
298 /// If recoils are disabled, the energy they would have received is just
299 /// deposited.
300 /// param flag if true, create recoil ions in cases in which the energy
301 /// is above the recoilCutoff.
302 /// If false, just deposit the energy.
303
304 void EnableRecoils(G4bool flag) { generateRecoils = flag; }
305
306 /// \brief find out if generation of recoils is enabled.
307 G4bool GetEnableRecoils() const { return generateRecoils; }
308
309 /// \brief set the mean free path scaling as specified
310 /// \param scale the factor by which the default MFP will be scaled.
311 /// Set to less than 1 for very thin films, typically,
312 /// to sample multiple scattering,
313 /// or to greater than 1 for quick simulations with a very long flight path
314
315 void SetMFPScaling(G4double scale) { MFPScale = scale; }
316
317 /// \brief get the MFPScaling parameter
318 G4double GetMFPScaling() const { return MFPScale; }
319
320 /// \brief enable or disable whether this process will skip collisions
321 /// which are close enough they need hadronic phsyics.
322 /// Default is true (skip close collisions).
323 /// Disabling this results in excess nuclear stopping power.
324 /// \param flag true results in hard collisions being skipped.
325 /// false allows hard collisions.
326
327 void AvoidNuclearReactions(G4bool flag) { avoidReactions = flag; }
328
329 /// \brief get the flag indicating whether hadronic collisions are ignored.
330 G4bool GetAvoidNuclearReactions() const { return avoidReactions; }
331
332 /// \brief set the minimum energy (per nucleon) at which recoils can
333 /// be generated,
334 /// and the energy (per nucleon) below which all ions are stopped.
335 /// \param energy energy per nucleon
336
337 void SetRecoilCutoff(G4double energy) { recoilCutoff = energy; }
338
339 /// \brief get the recoil cutoff
340 G4double GetRecoilCutoff() const { return recoilCutoff; }
341
342 /// \brief set the energy to which screening tables are computed.
343 /// Typically, this is 10 eV or so, and not often changed.
344 /// \param energy the cutoff energy
345
346 void SetPhysicsCutoff(G4double energy)
347 {
348 physicsCutoff = energy;
349 ResetTables();
350 }
351
352 /// \brief get the physics cutoff energy.
353 G4double GetPhysicsCutoff() const { return physicsCutoff; }
354
355 /// \brief set the pointer to a class for paritioning energy into NIEL
356 /// \brief part the pointer to the class.
357
358 void SetNIELPartitionFunction(const G4VNIELPartition* part);
359
360 /// \brief set the cross section boost to provide faster computation of
361 /// backscattering
362 /// \param fraction the fraction of particles to have their cross section
363 /// boosted.
364 /// \param HardeningFactor the factor by which to boost the scattering
365 /// cross section.
366
367 void SetCrossSectionHardening(G4double fraction, G4double HardeningFactor)
368 {
369 hardeningFraction = fraction;
370 hardeningFactor = HardeningFactor;
371 }
372
373 /// \brief get the fraction of particles which will have boosted scattering
374 G4double GetHardeningFraction() const { return hardeningFraction; }
375
376 /// \brief get the boost factor in use.
377 G4double GetHardeningFactor() const { return hardeningFactor; }
378
379 /// \brief the the interaciton length used in the last scattering.
380 G4double GetCurrentInteractionLength() const { return currentInteractionLength; }
381
382 /// \brief set a function to compute screening tables,
383 /// if the user needs non-standard behavior.
384 /// \param cs a class which constructs the screening tables.
385
386 void SetExternalCrossSectionHandler(G4ScreenedCoulombCrossSection* cs)
387 {
388 externalCrossSectionConstructor = cs;
389 }
390
391 /// \brief get the verbosity.
392 G4int GetVerboseLevel() const { return verboseLevel; }
393
394 std::map<G4int, G4ScreenedCoulombCrossSection*>& GetCrossSectionHandlers()
395 {
396 return crossSectionHandlers;
397 }
398
399 void ClearStages(void);
400
401 void AddStage(G4ScreenedCollisionStage* stage) { collisionStages.push_back(stage); }
402
403 G4CoulombKinematicsInfo& GetKinematics() { return kinematics; }
404
405 void SetValidCollision(G4bool flag) { validCollision = flag; }
406
407 G4bool GetValidCollision() const { return validCollision; }
408
409 /// \brief get the pointer to our ParticleChange object.
410 /// for internal use, primarily.
411 class G4ParticleChange& GetParticleChange()
412 {
413 return static_cast<G4ParticleChange&>(*pParticleChange);
414 }
415
416 /// \brief take the given energy, and use the material information
417 /// to partition it into NIEL and ionizing energy.
418
419 void DepositEnergy(G4int z1, G4double a1, const G4Material* material, G4double energy);
420
421 protected:
422 /// \brief the energy per nucleon above which the MFP is constant
423 G4double highEnergyLimit;
424
425 /// \brief the energy per nucleon below which the MFP is zero
426 G4double lowEnergyLimit;
427
428 /// \brief the energy per nucleon beyond which the cross section is zero,
429 /// to cross over to G4MSC
430 G4double processMaxEnergy;
431 G4String screeningKey;
432 G4bool generateRecoils, avoidReactions;
433 G4double recoilCutoff, physicsCutoff;
434 G4bool registerDepositedEnergy;
435 G4double IonizingLoss, NIEL;
436 G4double MFPScale;
437 G4double hardeningFraction, hardeningFactor;
438
439 G4ScreenedCoulombCrossSection* externalCrossSectionConstructor;
440 std::vector<G4ScreenedCollisionStage*> collisionStages;
441
442 std::map<G4int, G4ScreenedCoulombCrossSection*> crossSectionHandlers;
443
444 G4bool validCollision;
445 G4CoulombKinematicsInfo kinematics;
446 const G4VNIELPartition* NIELPartitionFunction;
447 };
448
449 // A customized G4CrossSectionHandler which gets its data from
450 // an external program
451
452 class G4NativeScreenedCoulombCrossSection : public G4ScreenedCoulombCrossSection
453 {
454 public:
455 G4NativeScreenedCoulombCrossSection();
456
457 G4NativeScreenedCoulombCrossSection(const G4NativeScreenedCoulombCrossSection& src)
458 : G4ScreenedCoulombCrossSection(src), phiMap(src.phiMap)
459 {}
460
461 G4NativeScreenedCoulombCrossSection(const G4ScreenedCoulombCrossSection& src)
462 : G4ScreenedCoulombCrossSection(src)
463 {}
464
465 virtual ~G4NativeScreenedCoulombCrossSection();
466
467 virtual void LoadData(G4String screeningKey, G4int z1, G4double m1, G4double recoilCutoff);
468
469 virtual G4ScreenedCoulombCrossSection* create()
470 {
471 return new G4NativeScreenedCoulombCrossSection(*this);
472 }
473
474 // get a list of available keys
475 std::vector<G4String> GetScreeningKeys() const;
476
477 typedef G4_c2_function& (*ScreeningFunc)(G4int z1, G4int z2, size_t nPoints, G4double rMax,
478 G4double* au);
479
480 void AddScreeningFunction(G4String name, ScreeningFunc fn) { phiMap[name] = fn; }
481
482 private:
483 // this is a map used to look up screening function generators
484 std::map<std::string, ScreeningFunc> phiMap;
485 };
486
487 G4_c2_function& ZBLScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double* auval);
488
489 G4_c2_function& MoliereScreening(G4int z1, G4int z2, size_t npoints, G4double rMax,
490 G4double* auval);
491
492 G4_c2_function& LJScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double* auval);
493
494 G4_c2_function& LJZBLScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double* auval);
495
496 #endif
497