| 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 name: G4UrbanAdjointMscModel
27 //
28 // Author: Laszlo Urban
29 //
30 // Class Description:
31 // Implementation of the model of multiple scattering based on
32 // H.W.Lewis Phys Rev 78 (1950) 526 and L.Urban model
33 // -------------------------------------------------------------------
34
35 #ifndef G4UrbanAdjointMscModel_h
36 #define G4UrbanAdjointMscModel_h 1
37
38 #include "G4Electron.hh"
39 #include "G4Exp.hh"
40 #include "G4Log.hh"
41 #include "G4MscStepLimitType.hh"
42 #include "G4VMscModel.hh"
43
44 class G4LossTableManager;
45 class G4MaterialCutsCouple;
46 class G4ParticleChangeForMSC;
47 class G4ParticleDefinition;
48 class G4SafetyHelper;
49
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
51
52 class G4UrbanAdjointMscModel : public G4VMscModel
53 {
54 public:
55 explicit G4UrbanAdjointMscModel(const G4String& nam = "UrbanMsc");
56
57 ~G4UrbanAdjointMscModel() override;
58
59 void Initialise(const G4ParticleDefinition*, const G4DataVector&) override;
60
61 void StartTracking(G4Track*) override;
62
63 G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition* particle,
64 G4double KineticEnergy,
65 G4double AtomicNumber,
66 G4double AtomicWeight = 0.,
67 G4double cut = 0.,
68 G4double emax = DBL_MAX) override;
69
70 G4ThreeVector& SampleScattering(const G4ThreeVector&,
71 G4double safety) override;
72
73 G4double ComputeTruePathLengthLimit(const G4Track& track,
74 G4double& currentMinimalStep) override;
75
76 G4double ComputeGeomPathLength(G4double truePathLength) override;
77
78 G4double ComputeTrueStepLength(G4double geomStepLength) override;
79
80 G4double ComputeTheta0(G4double truePathLength, G4double KineticEnergy);
81
82 inline void SetNewDisplacementFlag(G4bool);
83
84 G4UrbanAdjointMscModel& operator=(const G4UrbanAdjointMscModel& right) =
85 delete;
86 G4UrbanAdjointMscModel(const G4UrbanAdjointMscModel&) = delete;
87
88 private:
89 G4double SampleCosineTheta(G4double trueStepLength, G4double KineticEnergy);
90
91 void SampleDisplacement(G4double sinTheta, G4double phi);
92
93 void SampleDisplacementNew(G4double sinTheta, G4double phi);
94
95 inline void SetParticle(const G4ParticleDefinition*);
96
97 inline void UpdateCache();
98
99 inline G4double Randomizetlimit();
100
101 inline G4double SimpleScattering(G4double xmeanth, G4double x2meanth);
102
103 CLHEP::HepRandomEngine* rndmEngineMod;
104
105 const G4ParticleDefinition* particle;
106 const G4ParticleDefinition* positron;
107 G4ParticleChangeForMSC* fParticleChange;
108
109 const G4MaterialCutsCouple* couple;
110 G4LossTableManager* theManager;
111
112 G4double mass;
113 G4double charge, ChargeSquare;
114 G4double masslimite, lambdalimit, fr;
115
116 G4double taubig;
117 G4double tausmall;
118 G4double taulim;
119 G4double currentTau;
120 G4double tlimit;
121 G4double tlimitmin;
122 G4double tlimitminfix, tlimitminfix2;
123 G4double tgeom;
124
125 G4double geombig;
126 G4double geommin;
127 G4double geomlimit;
128 G4double skindepth;
129 G4double smallstep;
130
131 G4double presafety;
132
133 G4double lambda0;
134 G4double lambdaeff;
135 G4double tPathLength;
136 G4double zPathLength;
137 G4double par1, par2, par3;
138
139 G4double stepmin;
140
141 G4double currentKinEnergy;
142 G4double currentRange;
143 G4double rangeinit;
144 G4double currentRadLength;
145
146 G4double Zold;
147 G4double Zeff, Z2, Z23, lnZ;
148 G4double coeffth1, coeffth2;
149 G4double coeffc1, coeffc2, coeffc3, coeffc4;
150
151 G4double rangecut;
152 G4double drr, finalr;
153
154 G4int currentMaterialIndex;
155
156 G4bool firstStep;
157 G4bool insideskin;
158
159 G4bool latDisplasmentbackup;
160 G4bool displacementFlag;
161 };
162
163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
164 inline void G4UrbanAdjointMscModel::SetNewDisplacementFlag(G4bool val)
165 {
166 displacementFlag = val;
167 }
168
169 inline void G4UrbanAdjointMscModel::SetParticle(const G4ParticleDefinition* p)
170 {
171 const G4ParticleDefinition* p1 = p;
172
173 if(p->GetParticleName() == "adj_e-")
174 p1 = G4Electron::Electron();
175
176 if(p1 != particle)
177 {
178 particle = p1;
179 mass = p1->GetPDGMass();
180 charge = p1->GetPDGCharge() / CLHEP::eplus;
181 ChargeSquare = charge * charge;
182 }
183 }
184
185 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
186 inline G4double G4UrbanAdjointMscModel::Randomizetlimit()
187 {
188 G4double temptlimit = tlimit;
189 if(tlimit > tlimitmin)
190 {
191 G4double delta = tlimit - tlimitmin;
192 do
193 {
194 temptlimit = G4RandGauss::shoot(rndmEngineMod, tlimit, 0.1 * delta);
195 // Loop checking, 10-Apr-2016, Laszlo Urban
196 } while((temptlimit < tlimit - delta) || (temptlimit > tlimit + delta));
197 }
198 else
199 {
200 temptlimit = tlimitmin;
201 }
202
203 return temptlimit;
204 }
205
206 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
207 inline void G4UrbanAdjointMscModel::UpdateCache()
208 {
209 lnZ = G4Log(Zeff);
210 // correction in theta0 formula
211 G4double w = G4Exp(lnZ / 6.);
212 G4double facz = 0.990395 + w * (-0.168386 + w * 0.093286);
213 coeffth1 = facz * (1. - 8.7780e-2 / Zeff);
214 coeffth2 = facz * (4.0780e-2 + 1.7315e-4 * Zeff);
215
216 // tail parameters
217 G4double Z13 = w * w;
218 coeffc1 = 2.3785 - Z13 * (4.1981e-1 - Z13 * 6.3100e-2);
219 coeffc2 = 4.7526e-1 + Z13 * (1.7694 - Z13 * 3.3885e-1);
220 coeffc3 = 2.3683e-1 - Z13 * (1.8111 - Z13 * 3.2774e-1);
221 coeffc4 = 1.7888e-2 + Z13 * (1.9659e-2 - Z13 * 2.6664e-3);
222
223 Z2 = Zeff * Zeff;
224 Z23 = Z13 * Z13;
225
226 Zold = Zeff;
227 }
228
229 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
230 inline G4double G4UrbanAdjointMscModel::SimpleScattering(G4double xmeanth,
231 G4double x2meanth)
232 {
233 // 'large angle scattering'
234 // 2 model functions with correct xmean and x2mean
235 G4double a =
236 (2. * xmeanth + 9. * x2meanth - 3.) / (2. * xmeanth - 3. * x2meanth + 1.);
237 G4double prob = (a + 2.) * xmeanth / a;
238
239 // sampling
240 G4double cth = 1.;
241 if(rndmEngineMod->flat() < prob)
242 {
243 cth = -1. + 2. * G4Exp(G4Log(rndmEngineMod->flat()) / (a + 1.));
244 }
245 else
246 {
247 cth = -1. + 2. * rndmEngineMod->flat();
248 }
249 return cth;
250 }
251
252 #endif
253