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