| 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: G4UrbanMscModel.hh,v 1.33 2008/03/10 10:39:21 vnivanch Exp $
>> 27 // GEANT4 tag $Name: geant4-09-02-patch-03 $
26 // 28 //
27 // ------------------------------------------- 29 // -------------------------------------------------------------------
28 // 30 //
29 // 31 //
30 // GEANT4 Class header file 32 // GEANT4 Class header file
31 // 33 //
32 // 34 //
33 // File name: G4UrbanMscModel 35 // File name: G4UrbanMscModel
34 // 36 //
35 // Author: Laszlo Urban 37 // Author: Laszlo Urban
36 // 38 //
37 // Creation date: 19.02.2013 << 39 // Creation date: 01.03.2001
38 // 40 //
39 // Created from G4UrbanMscModel96 << 41 // Modifications:
>> 42 //
>> 43 // 27-03-03 Move model part from G4MultipleScattering (V.Ivanchenko)
>> 44 // 27-03-03 Rename (V.Ivanchenko)
>> 45 //
>> 46 // 05-08-03 angle distribution has been modified (L.Urban)
>> 47 // 26-11-03 new data member currentRange (L.Urban)
>> 48 // 01-03-04 changes in data members + signature changed in SampleCosineTheta
>> 49 // 11-03-04 changes in data members (L.Urban)
>> 50 // 23-04-04 changes in data members and in signature of SampleCosineTheta
>> 51 // (L.Urban)
>> 52 // 17-08-04 name of data member facxsi changed to factail (L.Urban)
>> 53 // 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
>> 54 // 15-04-05 optimize internal interface - add SampleSecondaries method (V.Ivanchenko)
>> 55 // 11-08-05 computation of lateral correlation added (L.Urban)
>> 56 // 02-10-05 nuclear size correction computation removed, the correction
>> 57 // included in the (theoretical) tabulated values (L.Urban)
>> 58 // 16-02-06 data members b and xsi have been removed (L.Urban)
>> 59 // 17-02-06 Save table of transport cross sections not mfp (V.Ivanchenko)
>> 60 // 07-03-06 Create G4UrbanMscModel and move there step limit calculation (V.Ivanchenko)
>> 61 // 10-05-06 SetMscStepLimitation at initialisation (V.Ivantchenko)
>> 62 // 11-05-06 name of data member safety changed to presafety, some new data
>> 63 // members added (frscaling1,frscaling2,tlimitminfix,nstepmax)
>> 64 // (L.Urban)
>> 65 // 13-10-06 data member factail removed, data member tkinlimit changed
>> 66 // to lambdalimit,
>> 67 // new data members tgeom,tnow,skin,skindepth,Zeff,geomlimit
>> 68 // G4double GeomLimit(const G4Track& track) changed to
>> 69 // void GeomLimit(const G4Track& track) (L.Urban)
>> 70 // 20-10-06 parameter theta0 now computed in the (public)
>> 71 // function ComputeTheta0,
>> 72 // single scattering modified allowing not small
>> 73 // angles as well (L.Urban)
>> 74 // 31-01-07 code cleaning (L.Urban)
>> 75 // 06-02-07 Move SetMscStepLimitation method into the source (VI)
>> 76 // 15-02-07 new data member : smallstep (L.Urban)
>> 77 // 10-04-07 remove navigator, smallstep, tnow (V.Ivanchenko)
40 // 78 //
41 // New parametrization for theta0 <<
42 // Correction for very small step length <<
43 // 79 //
44 // Class Description: 80 // Class Description:
45 // 81 //
46 // Implementation of the model of multiple sca 82 // Implementation of the model of multiple scattering based on
47 // H.W.Lewis Phys Rev 78 (1950) 526 and L.Urba 83 // H.W.Lewis Phys Rev 78 (1950) 526 and L.Urban model
48 84
49 // ------------------------------------------- 85 // -------------------------------------------------------------------
50 // 86 //
51 87
52 #ifndef G4UrbanMscModel_h 88 #ifndef G4UrbanMscModel_h
53 #define G4UrbanMscModel_h 1 89 #define G4UrbanMscModel_h 1
54 90
55 //....oooOO0OOooo........oooOO0OOooo........oo 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 92
57 #include <CLHEP/Units/SystemOfUnits.h> <<
58 <<
59 #include "G4VMscModel.hh" 93 #include "G4VMscModel.hh"
>> 94 #include "G4PhysicsTable.hh"
60 #include "G4MscStepLimitType.hh" 95 #include "G4MscStepLimitType.hh"
61 #include "G4Log.hh" <<
62 #include "G4Exp.hh" <<
63 96
64 class G4ParticleChangeForMSC; 97 class G4ParticleChangeForMSC;
65 class G4SafetyHelper; 98 class G4SafetyHelper;
>> 99 class G4LossTableManager;
66 100
67 //....oooOO0OOooo........oooOO0OOooo........oo 101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
68 102
69 class G4UrbanMscModel : public G4VMscModel 103 class G4UrbanMscModel : public G4VMscModel
70 { 104 {
71 105
72 public: 106 public:
73 107
74 explicit G4UrbanMscModel(const G4String& nam << 108 G4UrbanMscModel(const G4String& nam = "UrbanMscUni");
75 <<
76 ~G4UrbanMscModel() override; <<
77 109
78 void Initialise(const G4ParticleDefinition*, << 110 virtual ~G4UrbanMscModel();
79 const G4DataVector&) override; <<
80 111
81 void StartTracking(G4Track*) override; << 112 void Initialise(const G4ParticleDefinition*, const G4DataVector&);
>> 113
>> 114 G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition* particle,
>> 115 G4double KineticEnergy,
>> 116 G4double AtomicNumber,
>> 117 G4double AtomicWeight=0.,
>> 118 G4double cut =0.,
>> 119 G4double emax=DBL_MAX);
>> 120
>> 121 void SampleSecondaries(std::vector<G4DynamicParticle*>*,
>> 122 const G4MaterialCutsCouple*,
>> 123 const G4DynamicParticle*,
>> 124 G4double,
>> 125 G4double);
82 126
83 G4double << 127 void SampleScattering(const G4DynamicParticle*,
84 ComputeCrossSectionPerAtom(const G4ParticleD << 128 G4double safety);
85 G4double KineticEnergy, <<
86 G4double AtomicNumber, <<
87 G4double AtomicWeight=0., <<
88 G4double cut =0., <<
89 G4double emax=DBL_MAX) override; <<
90 <<
91 G4ThreeVector& SampleScattering(const G4Thre <<
92 G4double safety) override; <<
93 129
94 G4double ComputeTruePathLengthLimit(const G4 130 G4double ComputeTruePathLengthLimit(const G4Track& track,
95 G4double& currentMinimalSt << 131 G4PhysicsTable* theLambdaTable,
96 << 132 G4double currentMinimalStep);
97 G4double ComputeGeomPathLength(G4double true <<
98 <<
99 G4double ComputeTrueStepLength(G4double geom <<
100 133
101 G4double ComputeTheta0(G4double truePathLeng << 134 G4double ComputeGeomPathLength(G4double truePathLength);
102 135
103 inline void SetDisplacementAlgorithm96(const << 136 G4double ComputeTrueStepLength(G4double geomStepLength);
104 137
105 inline void SetPositronCorrection(const G4bo << 138 G4double ComputeTheta0(G4double truePathLength,
106 << 139 G4double KineticEnergy);
107 // hide assignment operator <<
108 G4UrbanMscModel & operator=(const G4UrbanMs <<
109 G4UrbanMscModel(const G4UrbanMscModel&) = d <<
110 140
111 private: 141 private:
112 142
113 G4double SampleCosineTheta(G4double trueStep 143 G4double SampleCosineTheta(G4double trueStepLength, G4double KineticEnergy);
114 144
115 void SampleDisplacement(G4double sinTheta, G << 145 G4double SampleDisplacement();
116 <<
117 void SampleDisplacementNew(G4double sinTheta <<
118 146
119 void InitialiseModelCache(); << 147 G4double LatCorrelation();
120 148
121 inline void SetParticle(const G4ParticleDefi << 149 void GeomLimit(const G4Track& track);
122 150
123 inline G4double Randomizetlimit(); << 151 inline G4double GetLambda(G4double kinEnergy);
124 <<
125 inline G4double SimpleScattering(); <<
126 152
127 inline G4double ComputeStepmin(); << 153 inline void SetParticle(const G4ParticleDefinition*);
128 154
129 inline G4double ComputeTlimitmin(); << 155 // hide assignment operator
>> 156 G4UrbanMscModel & operator=(const G4UrbanMscModel &right);
>> 157 G4UrbanMscModel(const G4UrbanMscModel&);
130 158
131 CLHEP::HepRandomEngine* rndmEngineMod; << 159 const G4ParticleDefinition* particle;
>> 160 G4ParticleChangeForMSC* fParticleChange;
132 161
133 const G4ParticleDefinition* particle = nullp << 162 G4SafetyHelper* safetyHelper;
134 const G4ParticleDefinition* positron; << 163 G4PhysicsTable* theLambdaTable;
135 G4ParticleChangeForMSC* fParticleChange = nu << 164 const G4MaterialCutsCouple* couple;
136 const G4MaterialCutsCouple* couple = nullptr << 165 G4LossTableManager* theManager;
137 166
138 G4double mass; 167 G4double mass;
139 G4double charge,chargeSquare; << 168 G4double charge;
140 G4double masslimite,fr; << 169
>> 170 G4double masslimite,masslimitmu;
141 171
142 G4double taubig; 172 G4double taubig;
143 G4double tausmall; 173 G4double tausmall;
144 G4double taulim; 174 G4double taulim;
145 G4double currentTau; 175 G4double currentTau;
>> 176 G4double frscaling1,frscaling2;
146 G4double tlimit; 177 G4double tlimit;
147 G4double tlimitmin; 178 G4double tlimitmin;
148 G4double tlimitminfix,tlimitminfix2; << 179 G4double tlimitminfix;
149 G4double tgeom; <<
150 180
>> 181 G4double nstepmax;
151 G4double geombig; 182 G4double geombig;
152 G4double geommin; 183 G4double geommin;
153 G4double geomlimit; 184 G4double geomlimit;
154 G4double skindepth; 185 G4double skindepth;
155 G4double smallstep; 186 G4double smallstep;
156 187
157 G4double presafety; 188 G4double presafety;
158 189
159 G4double lambda0; 190 G4double lambda0;
160 G4double lambdaeff; 191 G4double lambdaeff;
161 G4double tPathLength; 192 G4double tPathLength;
162 G4double zPathLength; 193 G4double zPathLength;
163 G4double par1,par2,par3; 194 G4double par1,par2,par3;
164 195
165 G4double stepmin; 196 G4double stepmin;
166 197
167 G4double currentKinEnergy; 198 G4double currentKinEnergy;
168 G4double currentLogKinEnergy; <<
169 G4double currentRange; 199 G4double currentRange;
170 G4double rangeinit; <<
171 G4double currentRadLength; 200 G4double currentRadLength;
172 201
173 G4double drr,finalr; << 202 G4double Zeff;
>> 203
>> 204 G4int currentMaterialIndex;
>> 205
>> 206 G4bool isInitialized;
>> 207 G4bool inside;
>> 208 G4bool insideskin;
174 209
175 G4double tlow; <<
176 G4double invmev; <<
177 G4double xmeanth = 0.0; <<
178 G4double x2meanth = 1./3.; <<
179 G4double rndmarray[2]; <<
180 <<
181 struct mscData { <<
182 G4double Z23, sqrtZ, factmin; <<
183 G4double coeffth1, coeffth2; <<
184 G4double coeffc1, coeffc2, coeffc3, coeffc <<
185 G4double stepmina, stepminb; <<
186 G4double doverra, doverrb; <<
187 G4double posa, posb, posc, posd, pose; <<
188 }; <<
189 static std::vector<mscData*> msc; <<
190 <<
191 // index of G4MaterialCutsCouple <<
192 G4int idx = 0; <<
193 <<
194 G4bool firstStep = true; <<
195 G4bool insideskin = false; <<
196 <<
197 G4bool latDisplasmentbackup = false; <<
198 G4bool dispAlg96 = true; <<
199 G4bool fPosiCorrection = true; <<
200 G4bool isFirstInstance = false; <<
201 }; 210 };
202 211
203 //....oooOO0OOooo........oooOO0OOooo........oo 212 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
204 //....oooOO0OOooo........oooOO0OOooo........oo <<
205 213
206 inline 214 inline
207 void G4UrbanMscModel::SetParticle(const G4Part << 215 G4double G4UrbanMscModel::GetLambda(G4double e)
208 { 216 {
209 if (p != particle) { << 217 G4double x;
210 particle = p; << 218 if(theLambdaTable) {
211 mass = p->GetPDGMass(); << 219 G4bool b;
212 charge = p->GetPDGCharge()/CLHEP::eplus; << 220 x = ((*theLambdaTable)[currentMaterialIndex])->GetValue(e, b);
213 chargeSquare = charge*charge; << 221 } else {
>> 222 x = CrossSection(couple,particle,e);
214 } 223 }
>> 224 if(x > DBL_MIN) x = 1./x;
>> 225 else x = DBL_MAX;
>> 226 return x;
215 } 227 }
216 228
217 //....oooOO0OOooo........oooOO0OOooo........oo 229 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
218 230
219 inline G4double G4UrbanMscModel::Randomizetlim << 231 inline
>> 232 void G4UrbanMscModel::SetParticle(const G4ParticleDefinition* p)
220 { 233 {
221 G4double res = tlimitmin; << 234 if (p != particle) {
222 if(tlimit > tlimitmin) << 235 particle = p;
223 { << 236 mass = p->GetPDGMass();
224 res = G4RandGauss::shoot(rndmEngineMod,tli << 237 charge = p->GetPDGCharge()/eplus;
225 res = std::max(res, tlimitmin); <<
226 } 238 }
227 return res; <<
228 } <<
229 <<
230 //....oooOO0OOooo........oooOO0OOooo........oo <<
231 <<
232 inline G4double G4UrbanMscModel::SimpleScatter <<
233 { <<
234 // 'large angle scattering' <<
235 // 2 model functions with correct xmean and <<
236 const G4double a = (2.*xmeanth+9.*x2meanth-3 <<
237 const G4double prob = (a+2.)*xmeanth/a; <<
238 <<
239 // sampling <<
240 rndmEngineMod->flatArray(2, rndmarray); <<
241 return (rndmarray[1] < prob) ? <<
242 -1.+2.*G4Exp(G4Log(rndmarray[0])/(a+1.)) : <<
243 } <<
244 <<
245 //....oooOO0OOooo........oooOO0OOooo........oo <<
246 <<
247 inline G4double G4UrbanMscModel::ComputeStepmi <<
248 { <<
249 // define stepmin using estimation of the ra <<
250 // of lambda_elastic/lambda_transport <<
251 const G4double rat = currentKinEnergy*invmev <<
252 return lambda0*msc[idx]->factmin/ <<
253 (0.002 + rat*(msc[idx]->stepmina + msc[idx <<
254 } <<
255 <<
256 //....oooOO0OOooo........oooOO0OOooo........oo <<
257 <<
258 inline G4double G4UrbanMscModel::ComputeTlimit <<
259 { <<
260 G4double x = (particle == positron) ? <<
261 0.7*msc[idx]->sqrtZ*stepmin : 0.87*msc[idx <<
262 if(currentKinEnergy < tlow) { x *= 0.5*(1.+c <<
263 return std::max(x, tlimitminfix); <<
264 } <<
265 <<
266 //....oooOO0OOooo........oooOO0OOooo........oo <<
267 <<
268 inline void G4UrbanMscModel::SetDisplacementAl <<
269 { <<
270 dispAlg96 = val; <<
271 } <<
272 <<
273 //....oooOO0OOooo........oooOO0OOooo........oo <<
274 <<
275 inline void G4UrbanMscModel::SetPositronCorrec <<
276 { <<
277 fPosiCorrection = val; <<
278 } 239 }
279 240
280 //....oooOO0OOooo........oooOO0OOooo........oo 241 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
281 242
282 #endif 243 #endif
283 244
284 245