| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // ------------------------------------------- 27 // ----------------------------------------------------------------------------
28 // 28 //
29 // GEANT4 Class header file 29 // GEANT4 Class header file
30 // 30 //
31 // File name: G4GoudsmitSaundersonMscModel 31 // File name: G4GoudsmitSaundersonMscModel
32 // 32 //
33 // Author: Mihaly Novak / (Omrane Kadri 33 // Author: Mihaly Novak / (Omrane Kadri)
34 // 34 //
35 // Creation date: 20.02.2009 35 // Creation date: 20.02.2009
36 // 36 //
37 // Modifications: 37 // Modifications:
38 // 04.03.2009 V.Ivanchenko cleanup and format 38 // 04.03.2009 V.Ivanchenko cleanup and format according to Geant4 EM style
39 // 12.05.2010 O.Kadri: adding Qn1 and Qn12 as 39 // 12.05.2010 O.Kadri: adding Qn1 and Qn12 as private doubles
40 // 18.05.2015 M. Novak provide PLERIMINARYY ve 40 // 18.05.2015 M. Novak provide PLERIMINARYY version of updated class.
41 // All algorithms of the class were 41 // All algorithms of the class were revised and updated, new methods added.
42 // A new version of Kawrakow-Bielaj 42 // A new version of Kawrakow-Bielajew Goudsmit-Saunderson MSC model
43 // based on the screened Rutherford 43 // based on the screened Rutherford DCS for elastic scattering of
44 // electrons/positrons has been int 44 // electrons/positrons has been introduced[1,2]. The corresponding MSC
45 // angular distributions over a 2D 45 // angular distributions over a 2D parameter grid have been recomputed
46 // and the CDFs are now stored in a 46 // and the CDFs are now stored in a variable transformed (smooth) form[2,3]
47 // together with the corresponding 47 // together with the corresponding rational interpolation parameters.
48 // These angular distributions are 48 // These angular distributions are handled by the new
49 // G4GoudsmitSaundersonTable class 49 // G4GoudsmitSaundersonTable class that is responsible to sample if
50 // it was no, single, few or multip 50 // it was no, single, few or multiple scattering case and delivers the
51 // angular deflection (i.e. cos(the 51 // angular deflection (i.e. cos(theta) and sin(theta)).
52 // Two screening options are provid 52 // Two screening options are provided:
53 // - if fgIsUsePWATotalXsecData=TR 53 // - if fgIsUsePWATotalXsecData=TRUE i.e. SetOptionPWAScreening(TRUE)
54 // was called before initialisat 54 // was called before initialisation: screening parameter value A is
55 // determined such that the firs 55 // determined such that the first transport coefficient G1(A)
56 // computed according to the scr 56 // computed according to the screened Rutherford DCS for elastic
57 // scattering will reproduce the 57 // scattering will reproduce the one computed from the PWA elastic
58 // and first transport mean free 58 // and first transport mean free paths[4].
59 // - if fgIsUsePWATotalXsecData=FA 59 // - if fgIsUsePWATotalXsecData=FALSE i.e. default value or
60 // SetOptionPWAScreening(FALSE) 60 // SetOptionPWAScreening(FALSE) was called before initialisation:
61 // screening parameter value A i 61 // screening parameter value A is computed according to Moliere's
62 // formula (by using material de 62 // formula (by using material dependent parameters \chi_cc2 and b_c
63 // precomputed for each material 63 // precomputed for each material used at initialization in
64 // G4GoudsmitSaundersonTable) [3 64 // G4GoudsmitSaundersonTable) [3]
65 // Elastic and first trasport mean 65 // Elastic and first trasport mean free paths are used consistently.
66 // The new version is self-consiste 66 // The new version is self-consistent, several times faster, more
67 // robust and accurate compared to 67 // robust and accurate compared to the earlier version.
68 // Spin effects as well as a more a 68 // Spin effects as well as a more accurate energy loss correction and
69 // computations of Lewis moments wi 69 // computations of Lewis moments will be implemented later on.
70 // 02.09.2015 M. Novak: first version of new s 70 // 02.09.2015 M. Novak: first version of new step limit is provided.
71 // fUseSafetyPlus corresponds to Ur 71 // fUseSafetyPlus corresponds to Urban fUseSafety (default)
72 // fUseDistanceToBoundary correspon 72 // fUseDistanceToBoundary corresponds to Urban fUseDistanceToBoundary
73 // fUseSafety corresponds to EGSnr 73 // fUseSafety corresponds to EGSnrc error-free stepping algorithm
74 // Range factor can be significantl 74 // Range factor can be significantly higher at each case than in Urban.
75 // 23.08.2017 M. Novak: added corrections to a 75 // 23.08.2017 M. Novak: added corrections to account spin effects (Mott-correction).
76 // It can be activated by setting t 76 // It can be activated by setting the fIsMottCorrection flag to be true
77 // before initialization using the 77 // before initialization using the SetOptionMottCorrection() public method.
78 // The fMottCorrection member is re 78 // The fMottCorrection member is responsible to handle pre-computed Mott
79 // correction (rejection) functions 79 // correction (rejection) functions obtained by numerically computing
80 // Goudsmit-Saunderson agnular dist 80 // Goudsmit-Saunderson agnular distributions based on a DCS accounting spin
81 // effects and screening correction 81 // effects and screening corrections. The DCS used to compute the accurate
82 // GS angular distributions is: DCS 82 // GS angular distributions is: DCS_{cor} = DCS_{SR}x[ DCS_{R}/DCS_{Mott}] where :
83 // # DCS_{SR} is the relativisti 83 // # DCS_{SR} is the relativistic Screened-Rutherford DCS (first Born approximate
84 // solution of the Klein-Gordo 84 // solution of the Klein-Gordon i.e. relativistic Schrodinger equation =>
85 // scattering of spinless e- o 85 // scattering of spinless e- on exponentially screened Coulomb potential)
86 // note: the default (without 86 // note: the default (without using Mott-correction) GS angular distributions
87 // are based on this DCS_{SR} 87 // are based on this DCS_{SR} with Moliere's screening parameter!
88 // # DCS_{R} is the Rutherford D 88 // # DCS_{R} is the Rutherford DCS which is the same as above but without
89 // screening 89 // screening
90 // # DCS_{Mott} is the Mott DCS 90 // # DCS_{Mott} is the Mott DCS i.e. solution of the Dirac equation with a bare
91 // Coulomb potential i.e. scat 91 // Coulomb potential i.e. scattering of particles with spin (e- or e+) on a
92 // point-like unscreened Coulo 92 // point-like unscreened Coulomb potential
93 // # moreover, the screening par 93 // # moreover, the screening parameter of the DCS_{cor} was determined such that
94 // the DCS_{cor} with this cor 94 // the DCS_{cor} with this corrected screening parameter reproduce the first
95 // transport cross sections ob 95 // transport cross sections obtained from the corresponding most accurate DCS
96 // (i.e. from elsepa [4]) 96 // (i.e. from elsepa [4])
97 // Unlike the default GS, the Mott- 97 // Unlike the default GS, the Mott-corrected angular distributions are particle type
98 // (different for e- and e+ <= the 98 // (different for e- and e+ <= the DCS_{Mott} and the screening correction) and target
99 // (Z and material) dependent. 99 // (Z and material) dependent.
100 // 02.02.2018 M. Novak: implemented CrossSecti 100 // 02.02.2018 M. Novak: implemented CrossSectionPerVolume interface method (used only for testing)
101 // 101 //
102 // Class description: 102 // Class description:
103 // Kawrakow-Bielajew Goudsmit-Saunderson MSC 103 // Kawrakow-Bielajew Goudsmit-Saunderson MSC model based on the screened Rutherford DCS
104 // for elastic scattering of e-/e+. Option, 104 // for elastic scattering of e-/e+. Option, to include (Mott) correction (see above), is
105 // also available now (SetOptionMottCorrecti 105 // also available now (SetOptionMottCorrection(true)). An EGSnrc like error-free stepping
106 // algorithm (UseSafety) is available beyond 106 // algorithm (UseSafety) is available beyond the usual Geant4 step limitation algorithms
107 // and true to geomerty and geometry to true 107 // and true to geomerty and geometry to true step length computations that were adopted
108 // from the Urban model[5]. The most accurat 108 // from the Urban model[5]. The most accurate setting: error-free stepping (UseSafety)
109 // with Mott-correction (SetOptionMottCorrec 109 // with Mott-correction (SetOptionMottCorrection(true)).
110 // 110 //
111 // References: 111 // References:
112 // [1] A.F.Bielajew, NIMB 111 (1996) 195-208 112 // [1] A.F.Bielajew, NIMB 111 (1996) 195-208
113 // [2] I.Kawrakow, A.F.Bielajew, NIMB 134(19 113 // [2] I.Kawrakow, A.F.Bielajew, NIMB 134(1998) 325-336
114 // [3] I.Kawrakow, E.Mainegra-Hing, D.W.O.Ro 114 // [3] I.Kawrakow, E.Mainegra-Hing, D.W.O.Rogers, F.Tessier,B.R.B.Walters, NRCC
115 // Report PIRS-701 (2013) 115 // Report PIRS-701 (2013)
116 // [4] F.Salvat, A.Jablonski, C.J. Powell, C 116 // [4] F.Salvat, A.Jablonski, C.J. Powell, CPC 165(2005) 157-190
117 // [5] L.Urban, Preprint CERN-OPEN-2006-077 117 // [5] L.Urban, Preprint CERN-OPEN-2006-077 (2006)
118 // 118 //
119 // ------------------------------------------- 119 // -----------------------------------------------------------------------------
120 120
121 #ifndef G4GoudsmitSaundersonMscModel_h 121 #ifndef G4GoudsmitSaundersonMscModel_h
122 #define G4GoudsmitSaundersonMscModel_h 1 122 #define G4GoudsmitSaundersonMscModel_h 1
123 123
124 #include <CLHEP/Units/SystemOfUnits.h> 124 #include <CLHEP/Units/SystemOfUnits.h>
125 125
126 #include "G4VMscModel.hh" 126 #include "G4VMscModel.hh"
127 #include "G4PhysicsTable.hh" 127 #include "G4PhysicsTable.hh"
128 #include "G4MaterialCutsCouple.hh" 128 #include "G4MaterialCutsCouple.hh"
129 #include "globals.hh" 129 #include "globals.hh"
130 130
131 131
132 class G4DataVector; 132 class G4DataVector;
133 class G4ParticleChangeForMSC; 133 class G4ParticleChangeForMSC;
134 class G4LossTableManager; 134 class G4LossTableManager;
135 class G4GoudsmitSaundersonTable; 135 class G4GoudsmitSaundersonTable;
136 class G4GSPWACorrections; 136 class G4GSPWACorrections;
137 137
138 class G4GoudsmitSaundersonMscModel : public G4 138 class G4GoudsmitSaundersonMscModel : public G4VMscModel
139 { 139 {
140 public: 140 public:
141 141
142 G4GoudsmitSaundersonMscModel(const G4String& 142 G4GoudsmitSaundersonMscModel(const G4String& nam = "GoudsmitSaunderson");
143 143
144 ~G4GoudsmitSaundersonMscModel() override; 144 ~G4GoudsmitSaundersonMscModel() override;
145 145
146 void Initialise(const G4ParticleDefinition*, 146 void Initialise(const G4ParticleDefinition*, const G4DataVector&) override;
147 147
148 void InitialiseLocal(const G4ParticleDefinit 148 void InitialiseLocal(const G4ParticleDefinition* p, G4VEmModel* masterModel) override;
149 149
150 G4ThreeVector& SampleScattering(const G4Thre 150 G4ThreeVector& SampleScattering(const G4ThreeVector&, G4double safety) override;
151 151
152 G4double ComputeTruePathLengthLimit(const G4 152 G4double ComputeTruePathLengthLimit(const G4Track& track, G4double& currentMinimalStep) override;
153 153
154 G4double ComputeGeomPathLength(G4double true 154 G4double ComputeGeomPathLength(G4double truePathLength) override;
155 155
156 G4double ComputeTrueStepLength(G4double geom 156 G4double ComputeTrueStepLength(G4double geomStepLength) override;
157 157
158 // method to compute first transport cross s 158 // method to compute first transport cross section per Volume (i.e. macroscropic first transport cross section; this
159 // method is used only for testing and not d 159 // method is used only for testing and not during a normal simulation)
160 G4double CrossSectionPerVolume(const G4Mater 160 G4double CrossSectionPerVolume(const G4Material*, const G4ParticleDefinition*, G4double kineticEnergy, G4double cutEnergy = 0.0, G4double maxEnergy = DBL_MAX) override;
161 161
162 void StartTracking(G4Track*) override; 162 void StartTracking(G4Track*) override;
163 163
164 void SampleMSC(); 164 void SampleMSC();
165 165
166 G4double GetTransportMeanFreePath(const G4Pa 166 G4double GetTransportMeanFreePath(const G4ParticleDefinition*, G4double);
167 167
168 void SetOptionPWACorrection(G4bool opt) { 168 void SetOptionPWACorrection(G4bool opt) { fIsUsePWACorrection = opt; }
169 169
170 G4bool GetOptionPWACorrection() const { 170 G4bool GetOptionPWACorrection() const { return fIsUsePWACorrection; }
171 171
172 void SetOptionMottCorrection(G4bool opt) { 172 void SetOptionMottCorrection(G4bool opt) { fIsUseMottCorrection = opt; }
173 173
174 G4bool GetOptionMottCorrection() const { 174 G4bool GetOptionMottCorrection() const { return fIsUseMottCorrection; }
175 175
176 G4GoudsmitSaundersonTable* GetGSTable() 176 G4GoudsmitSaundersonTable* GetGSTable() { return fGSTable; }
177 177
178 G4GSPWACorrections* GetPWACorrection( 178 G4GSPWACorrections* GetPWACorrection() { return fPWACorrection; }
179 179
180 // hide assignment operator 180 // hide assignment operator
181 G4GoudsmitSaundersonMscModel & operator=(con 181 G4GoudsmitSaundersonMscModel & operator=(const G4GoudsmitSaundersonMscModel &right) = delete;
182 G4GoudsmitSaundersonMscModel(const G4Goudsm 182 G4GoudsmitSaundersonMscModel(const G4GoudsmitSaundersonMscModel&) = delete;
183 183
184 private: 184 private:
185 inline void SetParticle(const G4Particle 185 inline void SetParticle(const G4ParticleDefinition* p);
186 186
187 inline G4double GetLambda(G4double); 187 inline G4double GetLambda(G4double);
188 188
189 G4double GetTransportMeanFreePathOnly(const 189 G4double GetTransportMeanFreePathOnly(const G4ParticleDefinition*,G4double);
190 190
191 inline G4double Randomizetlimit(); 191 inline G4double Randomizetlimit();
192 192
193 private: 193 private:
194 CLHEP::HepRandomEngine* rndmEngineMod; 194 CLHEP::HepRandomEngine* rndmEngineMod;
195 // 195 //
196 G4double currentKinEnergy; 196 G4double currentKinEnergy;
197 G4double currentRange; 197 G4double currentRange;
198 // 198 //
199 G4double fr; 199 G4double fr;
200 G4double rangeinit; 200 G4double rangeinit;
201 G4double geombig; 201 G4double geombig;
202 G4double geomlimit; 202 G4double geomlimit;
203 G4double tlimit; 203 G4double tlimit;
204 G4double tgeom; 204 G4double tgeom;
205 // 205 //
206 G4double par1; 206 G4double par1;
207 G4double par2; 207 G4double par2;
208 G4double par3; 208 G4double par3;
209 G4double tlimitminfix2; 209 G4double tlimitminfix2;
210 G4double tausmall; 210 G4double tausmall;
211 G4double mass; 211 G4double mass;
212 G4double taulim; 212 G4double taulim;
213 // 213 //
214 // 214 //
215 G4double presafety; 215 G4double presafety;
216 G4double fZeff; 216 G4double fZeff;
217 // 217 //
218 G4int charge; 218 G4int charge;
219 G4int currentMaterialIndex; 219 G4int currentMaterialIndex;
220 // 220 //
221 G4bool firstStep; 221 G4bool firstStep;
222 // 222 //
223 G4LossTableManager* theManager; 223 G4LossTableManager* theManager;
224 const G4ParticleDefinition* particle; 224 const G4ParticleDefinition* particle;
225 G4ParticleChangeForMSC* fParticleChange; 225 G4ParticleChangeForMSC* fParticleChange;
226 const G4MaterialCutsCouple* currentCouple; 226 const G4MaterialCutsCouple* currentCouple;
227 227
228 G4GoudsmitSaundersonTable* fGSTable; 228 G4GoudsmitSaundersonTable* fGSTable;
229 G4GSPWACorrections* fPWACorrection; 229 G4GSPWACorrections* fPWACorrection;
230 230
231 G4bool fIsUsePWACorrection; 231 G4bool fIsUsePWACorrection;
232 G4bool fIsUseMottCorrection; 232 G4bool fIsUseMottCorrection;
233 // 233 //
234 G4double fLambda0; // elastic mean free path 234 G4double fLambda0; // elastic mean free path
235 G4double fLambda1; // first transport mean f 235 G4double fLambda1; // first transport mean free path
236 G4double fScrA; // screening parameter 236 G4double fScrA; // screening parameter
237 G4double fG1; // first transport coef. 237 G4double fG1; // first transport coef.
238 // in case of Mott-correction 238 // in case of Mott-correction
239 G4double fMCtoScrA; 239 G4double fMCtoScrA;
240 G4double fMCtoQ1; 240 G4double fMCtoQ1;
241 G4double fMCtoG2PerG1; 241 G4double fMCtoG2PerG1;
242 // 242 //
243 G4double fTheTrueStepLenght; 243 G4double fTheTrueStepLenght;
244 G4double fTheTransportDistance; 244 G4double fTheTransportDistance;
245 G4double fTheZPathLenght; 245 G4double fTheZPathLenght;
246 // 246 //
247 G4ThreeVector fTheDisplacementVector; 247 G4ThreeVector fTheDisplacementVector;
248 G4ThreeVector fTheNewDirection; 248 G4ThreeVector fTheNewDirection;
249 // 249 //
250 G4bool fIsEndedUpOnBoundary; // step ended 250 G4bool fIsEndedUpOnBoundary; // step ended up on boundary i.e. transportation is the winer
251 G4bool fIsMultipleSacettring; 251 G4bool fIsMultipleSacettring;
252 G4bool fIsSingleScattering; 252 G4bool fIsSingleScattering;
253 G4bool fIsEverythingWasDone; 253 G4bool fIsEverythingWasDone;
254 G4bool fIsNoScatteringInMSC; 254 G4bool fIsNoScatteringInMSC;
255 G4bool fIsNoDisplace; 255 G4bool fIsNoDisplace;
256 G4bool fIsInsideSkin; 256 G4bool fIsInsideSkin;
257 G4bool fIsWasOnBoundary; 257 G4bool fIsWasOnBoundary;
258 G4bool fIsFirstRealStep; 258 G4bool fIsFirstRealStep;
259 // 259 //
260 static G4bool gIsUseAccurate; 260 static G4bool gIsUseAccurate;
261 static G4bool gIsOptimizationOn; 261 static G4bool gIsOptimizationOn;
262 }; 262 };
263 263
264 ////////////////////////////////////////////// 264 ////////////////////////////////////////////////////////////////////////////////
265 inline 265 inline
266 void G4GoudsmitSaundersonMscModel::SetParticle 266 void G4GoudsmitSaundersonMscModel::SetParticle(const G4ParticleDefinition* p)
267 { 267 {
268 if (p != particle) { 268 if (p != particle) {
269 particle = p; 269 particle = p;
270 charge = (G4int)(p->GetPDGCharge()/CLHEP:: 270 charge = (G4int)(p->GetPDGCharge()/CLHEP::eplus);
271 mass = p->GetPDGMass(); 271 mass = p->GetPDGMass();
272 } 272 }
273 } 273 }
274 274
275 275
276 //....oooOO0OOooo........oooOO0OOooo........oo 276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
277 inline 277 inline
278 G4double G4GoudsmitSaundersonMscModel::Randomi 278 G4double G4GoudsmitSaundersonMscModel::Randomizetlimit()
279 { 279 {
280 G4double temptlimit; 280 G4double temptlimit;
281 do { 281 do {
282 temptlimit = G4RandGauss::shoot(rndmE 282 temptlimit = G4RandGauss::shoot(rndmEngineMod,tlimit,0.1*tlimit);
283 } while ( (temptlimit<0.) || (temptlimi 283 } while ( (temptlimit<0.) || (temptlimit>2.*tlimit));
284 284
285 return temptlimit; 285 return temptlimit;
286 } 286 }
287 287
288 288
289 289
290 #endif 290 #endif
291 291