Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/include/G4eDPWAElasticDCS.hh

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 // GEANT4 Class header file
 30 //
 31 //
 32 // File name:     G4eDPWAElasticDCS
 33 //
 34 // Author:        Mihaly Novak
 35 //
 36 // Creation date: 02.07.2020
 37 //
 38 // Modifications:
 39 //
 40 // Class Description:
 41 //
 42 // Contains numerical Differential Cross Sections (DCS) for e-/e+ Coulomb
 43 // scattering computed by Dirac Partial Wave Analysis (DPWA) [1]:
 44 // - electrostatic interaction, with a local exchange correction in the case of
 45 //   electrons (using Dirac-Fock e- densities; finite nuclear size with Fermi
 46 //   charge distribution; exchange potential with Furness and McCarthy for e-)[2]
 47 // - correlation-polarization (projectiles cause the polarization of the charge
 48 //   cloud of the target atom and the induced dipole moment acts back on the
 49 //   projectile) was accounted by using Local-Density Approximation (LDA) [2]
 50 // - absorption: not included since it's an inelastic channel [2] (the cor-
 51 //   responding excitations needs to be modelled by a separate, independent,
 52 //   inelastic model).
 53 // Using the above mentioned DPWA computation with a free atom approximation
 54 // might lead to questionable results below few hundred [eV] where possible
 55 // solid state or bounding effects might start to affect the potential.
 56 // Nevertheless, the lower energy was set to 10 eV in order to provide(at least)
 57 // some model even at low energies (with this caution). The highest projectile
 58 // kinetic energy is 100 [MeV].
 59 //
 60 // The class provides interface methods for elastic, first-, second-transport
 61 // cross section computations as well as for sampling cosine of polar angular
 62 // deflections. These interface methods are also available for resricted cross
 63 // section computations and angular deflection sampling.
 64 //
 65 // References:
 66 //
 67 // [1] Salvat, F., Jablonski, A. and Powell, C.J., 2005. ELSEPA—Dirac partial-
 68 //     wave calculation of elastic scattering of electrons and positrons by
 69 //     atoms, positive ions and molecules. Computer physics communications,
 70 //     165(2), pp.157-190.
 71 // [2] Salvat, F., 2003. Optical-model potential for electron and positron
 72 //     elastic scattering by atoms. Physical Review A, 68(1), p.012708.
 73 // [3] Benedito, E., Fernández-Varea, J.M. and Salvat, F.,2001. Mixed simulation
 74 //     of the multiple elastic scattering of electrons and positrons using
 75 //     partial-wave differential cross-sections. Nuclear Instruments and Methods
 76 //     in Physics Research Section B: Beam Interactions with Materials and Atoms,
 77 //     174(1-2), pp.91-110.
 78 //
 79 // -------------------------------------------------------------------
 80 
 81 #ifndef G4eDPWAElasticDCS_h
 82 #define G4eDPWAElasticDCS_h 1
 83 
 84 
 85 #include <vector>
 86 #include <fstream>
 87 #include <iomanip>
 88 #include <sstream>
 89 
 90 #include "globals.hh"
 91 #include "G4String.hh"
 92 #include "G4Physics2DVector.hh"
 93 
 94 
 95 #include "G4MaterialTable.hh"
 96 #include "G4Material.hh"
 97 #include "G4Element.hh"
 98 #include "G4MaterialCutsCouple.hh"
 99 #include "G4ProductionCutsTable.hh"
100 
101 
102 #include "G4Log.hh"
103 #include "G4Exp.hh"
104 
105 
106 class G4eDPWAElasticDCS {
107 
108 public:
109 
110   // CTR:
111   // - iselectron   : data for e- (for e+ otherwise)
112   // - isrestricted : sampling of angular deflection on restricted interavl is
113   //                  required (i.e. in case of mixed-simulation models)
114   G4eDPWAElasticDCS(G4bool iselectron=true, G4bool isrestricted=false);
115 
116    // DTR
117  ~G4eDPWAElasticDCS();
118 
119   // initialise for a given 'iz' atomic number:
120   //  - nothing happens if it has already been initialised for that Z.
121   void InitialiseForZ(std::size_t iz);
122 
123   // Computes the elastic, first and second cross sections for the given kinetic
124   // energy and target atom.
125   // Cross sections are zero ff ekin is below/above the kinetic energy grid
126   void  ComputeCSPerAtom(G4int iz, G4double ekin, G4double& elcs, G4double& tr1cs,
127                         G4double& tr2cs, G4double mumin=0.0, G4double mumax=1.0);
128 
129 
130   // samples cos(theta) i.e. cosine of the polar angle of scattering in elastic
131   // interaction (Coulomb scattering) of the projectile (e- or e+ depending on
132   // fIsElectron) with kinetic energy of exp('lekin'), target atom with atomic
133   // muber of 'iz'. See the 'SampleCosineThetaRestricted' for obtain samples on
134   // a restricted inteval.
135   G4double SampleCosineTheta(std::size_t iz, G4double lekin, G4double r1,
136                              G4double r2, G4double r3);
137 
138   // samples cos(theta) i.e. cosine of the polar angle of scattering in elastic
139   // interaction (Coulomb scattering) of the projectile (e- or e+ depending on
140   // fIsElectron) with kinetic energy of exp('lekin'), target atom with atomic
141   // muber of 'iz'.
142   // The cosine theta will be in the [costMin, costMax] interval where costMin
143   // corresponds to a maximum allowed polar scattering angle thetaMax while
144   // costMin corresponds to minimum allowed polar scatterin angle thetaMin.
145   // See the 'SampleCosineTheta' for obtain samples on the entire [-1,1] range.
146   G4double SampleCosineThetaRestricted(std::size_t iz, G4double lekin,
147                                        G4double r1, G4double r2,
148                                        G4double costMax, G4double costMin);
149 
150   // interpolate scattering power correction form table buit at init.
151   G4double ComputeScatteringPowerCorrection(const G4MaterialCutsCouple *matcut,
152                                             G4double ekin);
153 
154   // build scattering power correction table at init.
155   void     InitSCPCorrection(G4double lowEnergyLimit, G4double highEnergyLimit);
156 
157 
158 private:
159 
160   // data structure to store one sampling table: combined Alias + RatIn
161   // NOTE: when Alias is used, sampling on a resctricted interval is not possible
162   //       However, Alias makes possible faster sampling. Alias is used in case
163   //       of single scattering model while it's not used in case of mixed-model
164   //       when restricted interval sampling is needed. This is controlled by
165   //       the fIsRestrictedSamplingRequired flag (false by default).
166   struct OneSamplingTable {
167     OneSamplingTable () = default;
168     void SetSize(std::size_t nx, G4bool useAlias)  {
169       fN = nx;
170       // Alias
171       if (useAlias) {
172         fW.resize(nx);
173         fI.resize(nx);
174       }
175       // Ratin
176       fCum.resize(nx);
177       fA.resize(nx);
178       fB.resize(nx);
179     }
180 
181     // members
182     std::size_t           fN;            // # data points
183     G4double              fScreenParA;   // the screening parameter
184     std::vector<G4double> fW;
185     std::vector<G4double> fCum;
186     std::vector<G4double> fA;
187     std::vector<G4double> fB;
188     std::vector<G4int>    fI;
189   };
190 
191 
192   // loads the kinetic energy and theta grids for the DCS data (first init step)
193   // should be called only by the master
194   void LoadGrid();
195 
196   // load DCS data for a given Z
197   void LoadDCSForZ(G4int iz);
198 
199   // loads sampling table for the given Z over the enrgy grid
200   void BuildSmplingTableForZ(G4int iz);
201 
202   G4double SampleMu(std::size_t izet, std::size_t ie, G4double r1, G4double r2);
203 
204   G4double FindCumValue(G4double u, const OneSamplingTable& stable,
205                         const std::vector<G4double>& uvect);
206 
207   // muMin and muMax : no checks on these
208   G4double SampleMu(std::size_t izet, std::size_t ie, G4double r1, G4double muMin,
209                     G4double muMax);
210 
211   // set the DCS data directory path
212   const G4String& FindDirectoryPath();
213 
214   // uncompress one data file into the input string stream
215   void ReadCompressedFile(G4String fname, std::istringstream &iss);
216 
217   // compute Molier material dependent parameters
218   void ComputeMParams(const G4Material* mat, G4double& theBc, G4double& theXc2);
219 
220 
221 // members
222 private:
223 
224   // indicates if the object is for mixed-simulation (single scatterin otherwise)
225   G4bool                           fIsRestrictedSamplingRequired;
226   // indicates if the object is for e- (for e+ otherwise)
227   G4bool                           fIsElectron;
228   // indicates if the ekin, mu grids has already been loaded (only once)
229   static G4bool                    gIsGridLoaded;
230   // data directory
231   static G4String                  gDataDirectory;
232   // max atomic number (Z) for which DCS has been computed (103)
233   static constexpr std::size_t     gMaxZ = 103;
234   // energy and theta grid(s) relaed variables: loaded from gridinfo by LoadGrid
235   static std::size_t               gNumEnergies;
236   static std::size_t               gIndxEnergyLim;// the energy index just above 2 [keV]
237   static std::size_t               gNumThetas1;   // used for e- below 2 [keV]
238   static std::size_t               gNumThetas2;   // used for e+ and for e- bove 2 [keV]
239   static std::vector<G4double>     gTheEnergies;  // log-kinetic energy grid
240   static std::vector<G4double>     gTheMus1;      // mu(theta) = 0.5[1-cos(theta)]
241   static std::vector<G4double>     gTheMus2;
242   static std::vector<G4double>     gTheU1;        // u(mu; A'=0.01) = (A'+1)mu/(mu+A')
243   static std::vector<G4double>     gTheU2;
244   static G4double                  gLogMinEkin;   // log(gTheEnergies[0])
245   static G4double                  gInvDelLogEkin;// 1./log(gTheEnergies[i+1]/gTheEnergies[i])
246   // abscissas and weights of an 8 point Gauss-Legendre quadrature
247   // for numerical integration on [0,1]
248   static const G4double            gXGL[8];
249   static const G4double            gWGL[8];
250   //
251   std::vector<G4Physics2DVector*>  fDCS;    // log(DCS) data per Z
252   std::vector<G4Physics2DVector*>  fDCSLow; // only for e- E < 2keV
253   // sampling tables: only one of the followings will be utilized
254   std::vector< std::vector<OneSamplingTable>* > fSamplingTables;
255   //
256   // scattering power correction: to account sub-threshold inelastic deflections
257   const G4int                      fNumSPCEbinPerDec = 3;
258   struct SCPCorrection {
259     G4bool   fIsUse;               //
260     G4double fPrCut;               // sec. e- production cut energy
261     G4double fLEmin;               // log min energy
262     G4double fILDel;               // inverse log delta kinetic energy
263     std::vector<G4double> fVSCPC;  // scattering power correction vector
264   };
265   std::vector<SCPCorrection*>      fSCPCPerMatCuts;
266 
267 
268 };
269 
270 #endif
271