| 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 // $Id: G4GoudsmitSaundersonMscModel.hh 106953 2017-10-31 08:30:13Z gcosmo $
26 // 27 //
27 // ------------------------------------------- 28 // ----------------------------------------------------------------------------
28 // 29 //
29 // GEANT4 Class header file 30 // GEANT4 Class header file
30 // 31 //
31 // File name: G4GoudsmitSaundersonMscModel 32 // File name: G4GoudsmitSaundersonMscModel
32 // 33 //
33 // Author: Mihaly Novak / (Omrane Kadri 34 // Author: Mihaly Novak / (Omrane Kadri)
34 // 35 //
35 // Creation date: 20.02.2009 36 // Creation date: 20.02.2009
36 // 37 //
37 // Modifications: 38 // Modifications:
38 // 04.03.2009 V.Ivanchenko cleanup and format 39 // 04.03.2009 V.Ivanchenko cleanup and format according to Geant4 EM style
39 // 12.05.2010 O.Kadri: adding Qn1 and Qn12 as 40 // 12.05.2010 O.Kadri: adding Qn1 and Qn12 as private doubles
40 // 18.05.2015 M. Novak provide PLERIMINARYY ve 41 // 18.05.2015 M. Novak provide PLERIMINARYY version of updated class.
41 // All algorithms of the class were 42 // All algorithms of the class were revised and updated, new methods added.
42 // A new version of Kawrakow-Bielaj 43 // A new version of Kawrakow-Bielajew Goudsmit-Saunderson MSC model
43 // based on the screened Rutherford 44 // based on the screened Rutherford DCS for elastic scattering of
44 // electrons/positrons has been int 45 // electrons/positrons has been introduced[1,2]. The corresponding MSC
45 // angular distributions over a 2D 46 // angular distributions over a 2D parameter grid have been recomputed
46 // and the CDFs are now stored in a 47 // and the CDFs are now stored in a variable transformed (smooth) form[2,3]
47 // together with the corresponding 48 // together with the corresponding rational interpolation parameters.
48 // These angular distributions are 49 // These angular distributions are handled by the new
49 // G4GoudsmitSaundersonTable class 50 // G4GoudsmitSaundersonTable class that is responsible to sample if
50 // it was no, single, few or multip 51 // it was no, single, few or multiple scattering case and delivers the
51 // angular deflection (i.e. cos(the 52 // angular deflection (i.e. cos(theta) and sin(theta)).
52 // Two screening options are provid 53 // Two screening options are provided:
53 // - if fgIsUsePWATotalXsecData=TR 54 // - if fgIsUsePWATotalXsecData=TRUE i.e. SetOptionPWAScreening(TRUE)
54 // was called before initialisat 55 // was called before initialisation: screening parameter value A is
55 // determined such that the firs 56 // determined such that the first transport coefficient G1(A)
56 // computed according to the scr 57 // computed according to the screened Rutherford DCS for elastic
57 // scattering will reproduce the 58 // scattering will reproduce the one computed from the PWA elastic
58 // and first transport mean free 59 // and first transport mean free paths[4].
59 // - if fgIsUsePWATotalXsecData=FA 60 // - if fgIsUsePWATotalXsecData=FALSE i.e. default value or
60 // SetOptionPWAScreening(FALSE) 61 // SetOptionPWAScreening(FALSE) was called before initialisation:
61 // screening parameter value A i 62 // screening parameter value A is computed according to Moliere's
62 // formula (by using material de 63 // formula (by using material dependent parameters \chi_cc2 and b_c
63 // precomputed for each material 64 // precomputed for each material used at initialization in
64 // G4GoudsmitSaundersonTable) [3 65 // G4GoudsmitSaundersonTable) [3]
65 // Elastic and first trasport mean 66 // Elastic and first trasport mean free paths are used consistently.
66 // The new version is self-consiste 67 // The new version is self-consistent, several times faster, more
67 // robust and accurate compared to 68 // robust and accurate compared to the earlier version.
68 // Spin effects as well as a more a 69 // Spin effects as well as a more accurate energy loss correction and
69 // computations of Lewis moments wi 70 // computations of Lewis moments will be implemented later on.
70 // 02.09.2015 M. Novak: first version of new s 71 // 02.09.2015 M. Novak: first version of new step limit is provided.
71 // fUseSafetyPlus corresponds to Ur 72 // fUseSafetyPlus corresponds to Urban fUseSafety (default)
72 // fUseDistanceToBoundary correspon 73 // fUseDistanceToBoundary corresponds to Urban fUseDistanceToBoundary
73 // fUseSafety corresponds to EGSnr 74 // fUseSafety corresponds to EGSnrc error-free stepping algorithm
74 // Range factor can be significantl 75 // Range factor can be significantly higher at each case than in Urban.
75 // 23.08.2017 M. Novak: added corrections to a 76 // 23.08.2017 M. Novak: added corrections to account spin effects (Mott-correction).
76 // It can be activated by setting t 77 // It can be activated by setting the fIsMottCorrection flag to be true
77 // before initialization using the 78 // before initialization using the SetOptionMottCorrection() public method.
78 // The fMottCorrection member is re 79 // The fMottCorrection member is responsible to handle pre-computed Mott
79 // correction (rejection) functions 80 // correction (rejection) functions obtained by numerically computing
80 // Goudsmit-Saunderson agnular dist 81 // Goudsmit-Saunderson agnular distributions based on a DCS accounting spin
81 // effects and screening correction 82 // effects and screening corrections. The DCS used to compute the accurate
82 // GS angular distributions is: DCS 83 // GS angular distributions is: DCS_{cor} = DCS_{SR}x[ DCS_{R}/DCS_{Mott}] where :
83 // # DCS_{SR} is the relativisti 84 // # DCS_{SR} is the relativistic Screened-Rutherford DCS (first Born approximate
84 // solution of the Klein-Gordo 85 // solution of the Klein-Gordon i.e. relativistic Schrodinger equation =>
85 // scattering of spinless e- o 86 // scattering of spinless e- on exponentially screened Coulomb potential)
86 // note: the default (without 87 // note: the default (without using Mott-correction) GS angular distributions
87 // are based on this DCS_{SR} 88 // are based on this DCS_{SR} with Moliere's screening parameter!
88 // # DCS_{R} is the Rutherford D 89 // # DCS_{R} is the Rutherford DCS which is the same as above but without
89 // screening 90 // screening
90 // # DCS_{Mott} is the Mott DCS 91 // # DCS_{Mott} is the Mott DCS i.e. solution of the Dirac equation with a bare
91 // Coulomb potential i.e. scat 92 // Coulomb potential i.e. scattering of particles with spin (e- or e+) on a
92 // point-like unscreened Coulo 93 // point-like unscreened Coulomb potential
93 // # moreover, the screening par 94 // # moreover, the screening parameter of the DCS_{cor} was determined such that
94 // the DCS_{cor} with this cor 95 // the DCS_{cor} with this corrected screening parameter reproduce the first
95 // transport cross sections ob 96 // transport cross sections obtained from the corresponding most accurate DCS
96 // (i.e. from elsepa [4]) 97 // (i.e. from elsepa [4])
97 // Unlike the default GS, the Mott- 98 // Unlike the default GS, the Mott-corrected angular distributions are particle type
98 // (different for e- and e+ <= the 99 // (different for e- and e+ <= the DCS_{Mott} and the screening correction) and target
99 // (Z and material) dependent. 100 // (Z and material) dependent.
100 // 02.02.2018 M. Novak: implemented CrossSecti <<
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 virtual ~G4GoudsmitSaundersonMscModel();
145 145
146 void Initialise(const G4ParticleDefinition*, << 146 virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&);
147 147
148 void InitialiseLocal(const G4ParticleDefinit << 148 virtual void InitialiseLocal(const G4ParticleDefinition* p, G4VEmModel* masterModel);
149 149
150 G4ThreeVector& SampleScattering(const G4Thre <<
151 150
152 G4double ComputeTruePathLengthLimit(const G4 << 151 virtual G4ThreeVector& SampleScattering(const G4ThreeVector&, G4double safety);
153 152
154 G4double ComputeGeomPathLength(G4double true << 153 virtual G4double ComputeTruePathLengthLimit(const G4Track& track, G4double& currentMinimalStep);
155 154
156 G4double ComputeTrueStepLength(G4double geom << 155 virtual G4double ComputeGeomPathLength(G4double truePathLength);
157 156
158 // method to compute first transport cross s << 157 virtual G4double ComputeTrueStepLength(G4double geomStepLength);
159 // method is used only for testing and not d <<
160 G4double CrossSectionPerVolume(const G4Mater <<
161 158
162 void StartTracking(G4Track*) override; << 159 void StartTracking(G4Track*);
163 160
164 void SampleMSC(); 161 void SampleMSC();
165 162
166 G4double GetTransportMeanFreePath(const G4Pa 163 G4double GetTransportMeanFreePath(const G4ParticleDefinition*, G4double);
167 164
168 void SetOptionPWACorrection(G4bool opt) { 165 void SetOptionPWACorrection(G4bool opt) { fIsUsePWACorrection = opt; }
169 166
170 G4bool GetOptionPWACorrection() const { 167 G4bool GetOptionPWACorrection() const { return fIsUsePWACorrection; }
171 168
172 void SetOptionMottCorrection(G4bool opt) { 169 void SetOptionMottCorrection(G4bool opt) { fIsUseMottCorrection = opt; }
173 170
174 G4bool GetOptionMottCorrection() const { 171 G4bool GetOptionMottCorrection() const { return fIsUseMottCorrection; }
175 172
176 G4GoudsmitSaundersonTable* GetGSTable() 173 G4GoudsmitSaundersonTable* GetGSTable() { return fGSTable; }
177 174
178 G4GSPWACorrections* GetPWACorrection( 175 G4GSPWACorrections* GetPWACorrection() { return fPWACorrection; }
179 176
180 // hide assignment operator <<
181 G4GoudsmitSaundersonMscModel & operator=(con <<
182 G4GoudsmitSaundersonMscModel(const G4Goudsm <<
183 <<
184 private: 177 private:
185 inline void SetParticle(const G4Particle 178 inline void SetParticle(const G4ParticleDefinition* p);
186 179
187 inline G4double GetLambda(G4double); 180 inline G4double GetLambda(G4double);
188 181
>> 182 // hide assignment operator
>> 183 G4GoudsmitSaundersonMscModel & operator=(const G4GoudsmitSaundersonMscModel &right);
>> 184 G4GoudsmitSaundersonMscModel(const G4GoudsmitSaundersonMscModel&);
189 G4double GetTransportMeanFreePathOnly(const 185 G4double GetTransportMeanFreePathOnly(const G4ParticleDefinition*,G4double);
190 186
191 inline G4double Randomizetlimit(); 187 inline G4double Randomizetlimit();
192 188
193 private: 189 private:
194 CLHEP::HepRandomEngine* rndmEngineMod; 190 CLHEP::HepRandomEngine* rndmEngineMod;
195 // 191 //
196 G4double currentKinEnergy; 192 G4double currentKinEnergy;
197 G4double currentRange; 193 G4double currentRange;
198 // 194 //
199 G4double fr; 195 G4double fr;
200 G4double rangeinit; 196 G4double rangeinit;
201 G4double geombig; 197 G4double geombig;
202 G4double geomlimit; 198 G4double geomlimit;
>> 199 G4double lambdalimit;
203 G4double tlimit; 200 G4double tlimit;
204 G4double tgeom; 201 G4double tgeom;
205 // 202 //
206 G4double par1; 203 G4double par1;
207 G4double par2; 204 G4double par2;
208 G4double par3; 205 G4double par3;
209 G4double tlimitminfix2; 206 G4double tlimitminfix2;
210 G4double tausmall; 207 G4double tausmall;
211 G4double mass; 208 G4double mass;
212 G4double taulim; 209 G4double taulim;
213 // 210 //
214 // 211 //
215 G4double presafety; 212 G4double presafety;
216 G4double fZeff; 213 G4double fZeff;
217 // 214 //
218 G4int charge; 215 G4int charge;
219 G4int currentMaterialIndex; 216 G4int currentMaterialIndex;
220 // 217 //
221 G4bool firstStep; 218 G4bool firstStep;
222 // 219 //
223 G4LossTableManager* theManager; 220 G4LossTableManager* theManager;
224 const G4ParticleDefinition* particle; 221 const G4ParticleDefinition* particle;
225 G4ParticleChangeForMSC* fParticleChange; 222 G4ParticleChangeForMSC* fParticleChange;
226 const G4MaterialCutsCouple* currentCouple; 223 const G4MaterialCutsCouple* currentCouple;
227 224
228 G4GoudsmitSaundersonTable* fGSTable; 225 G4GoudsmitSaundersonTable* fGSTable;
229 G4GSPWACorrections* fPWACorrection; 226 G4GSPWACorrections* fPWACorrection;
230 227
231 G4bool fIsUsePWACorrection; 228 G4bool fIsUsePWACorrection;
232 G4bool fIsUseMottCorrection; 229 G4bool fIsUseMottCorrection;
233 // 230 //
234 G4double fLambda0; // elastic mean free path 231 G4double fLambda0; // elastic mean free path
235 G4double fLambda1; // first transport mean f 232 G4double fLambda1; // first transport mean free path
236 G4double fScrA; // screening parameter 233 G4double fScrA; // screening parameter
237 G4double fG1; // first transport coef. 234 G4double fG1; // first transport coef.
238 // in case of Mott-correction 235 // in case of Mott-correction
239 G4double fMCtoScrA; 236 G4double fMCtoScrA;
240 G4double fMCtoQ1; 237 G4double fMCtoQ1;
241 G4double fMCtoG2PerG1; 238 G4double fMCtoG2PerG1;
242 // 239 //
243 G4double fTheTrueStepLenght; 240 G4double fTheTrueStepLenght;
244 G4double fTheTransportDistance; 241 G4double fTheTransportDistance;
245 G4double fTheZPathLenght; 242 G4double fTheZPathLenght;
246 // 243 //
247 G4ThreeVector fTheDisplacementVector; 244 G4ThreeVector fTheDisplacementVector;
248 G4ThreeVector fTheNewDirection; 245 G4ThreeVector fTheNewDirection;
249 // 246 //
250 G4bool fIsEndedUpOnBoundary; // step ended 247 G4bool fIsEndedUpOnBoundary; // step ended up on boundary i.e. transportation is the winer
251 G4bool fIsMultipleSacettring; 248 G4bool fIsMultipleSacettring;
252 G4bool fIsSingleScattering; 249 G4bool fIsSingleScattering;
253 G4bool fIsEverythingWasDone; 250 G4bool fIsEverythingWasDone;
254 G4bool fIsNoScatteringInMSC; 251 G4bool fIsNoScatteringInMSC;
255 G4bool fIsNoDisplace; 252 G4bool fIsNoDisplace;
256 G4bool fIsInsideSkin; 253 G4bool fIsInsideSkin;
257 G4bool fIsWasOnBoundary; 254 G4bool fIsWasOnBoundary;
258 G4bool fIsFirstRealStep; 255 G4bool fIsFirstRealStep;
259 // 256 //
260 static G4bool gIsUseAccurate; 257 static G4bool gIsUseAccurate;
261 static G4bool gIsOptimizationOn; 258 static G4bool gIsOptimizationOn;
262 }; 259 };
263 260
264 ////////////////////////////////////////////// 261 ////////////////////////////////////////////////////////////////////////////////
265 inline 262 inline
266 void G4GoudsmitSaundersonMscModel::SetParticle 263 void G4GoudsmitSaundersonMscModel::SetParticle(const G4ParticleDefinition* p)
267 { 264 {
268 if (p != particle) { 265 if (p != particle) {
269 particle = p; 266 particle = p;
270 charge = (G4int)(p->GetPDGCharge()/CLHEP:: 267 charge = (G4int)(p->GetPDGCharge()/CLHEP::eplus);
271 mass = p->GetPDGMass(); 268 mass = p->GetPDGMass();
272 } 269 }
273 } 270 }
274 271
275 272
276 //....oooOO0OOooo........oooOO0OOooo........oo 273 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
277 inline 274 inline
278 G4double G4GoudsmitSaundersonMscModel::Randomi 275 G4double G4GoudsmitSaundersonMscModel::Randomizetlimit()
279 { 276 {
280 G4double temptlimit; << 277 G4double temptlimit = tlimit;
281 do { 278 do {
282 temptlimit = G4RandGauss::shoot(rndmE 279 temptlimit = G4RandGauss::shoot(rndmEngineMod,tlimit,0.1*tlimit);
283 } while ( (temptlimit<0.) || (temptlimi 280 } while ( (temptlimit<0.) || (temptlimit>2.*tlimit));
284 281
285 return temptlimit; 282 return temptlimit;
286 } 283 }
287 284
288 285
289 286
290 #endif 287 #endif
291 288