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