| 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: G4VMscModel.hh 71484 2013-06-17 08:12:51Z gcosmo $
>> 27 //
26 // ------------------------------------------- 28 // -------------------------------------------------------------------
27 // 29 //
28 // GEANT4 Class header file 30 // GEANT4 Class header file
29 // 31 //
30 // 32 //
31 // File name: G4VMscModel 33 // File name: G4VMscModel
32 // 34 //
33 // Author: Vladimir Ivanchenko 35 // Author: Vladimir Ivanchenko
34 // 36 //
35 // Creation date: 07.03.2008 37 // Creation date: 07.03.2008
36 // 38 //
37 // Modifications: 39 // Modifications:
38 // 07.04.2009 V.Ivanchenko moved msc methods f 40 // 07.04.2009 V.Ivanchenko moved msc methods from G4VEmModel to G4VMscModel
39 // 26.03.2012 V.Ivanchenko added transport x-s 41 // 26.03.2012 V.Ivanchenko added transport x-section pointer and Get?Set methods
40 // 42 //
41 // Class Description: 43 // Class Description:
42 // 44 //
43 // General interface to msc models 45 // General interface to msc models
44 46
45 // ------------------------------------------- 47 // -------------------------------------------------------------------
46 // 48 //
47 #ifndef G4VMscModel_h 49 #ifndef G4VMscModel_h
48 #define G4VMscModel_h 1 50 #define G4VMscModel_h 1
49 51
50 #include <CLHEP/Units/SystemOfUnits.h> 52 #include <CLHEP/Units/SystemOfUnits.h>
51 53
52 #include "G4VEmModel.hh" 54 #include "G4VEmModel.hh"
53 #include "G4MscStepLimitType.hh" 55 #include "G4MscStepLimitType.hh"
54 #include "globals.hh" 56 #include "globals.hh"
55 #include "G4ThreeVector.hh" 57 #include "G4ThreeVector.hh"
56 #include "G4Track.hh" 58 #include "G4Track.hh"
57 #include "G4SafetyHelper.hh" 59 #include "G4SafetyHelper.hh"
>> 60 #include "G4VEnergyLossProcess.hh"
>> 61 #include "G4LossTableManager.hh"
58 #include "G4PhysicsTable.hh" 62 #include "G4PhysicsTable.hh"
59 #include "G4ThreeVector.hh" 63 #include "G4ThreeVector.hh"
60 #include <vector> 64 #include <vector>
61 65
62 class G4ParticleChangeForMSC; 66 class G4ParticleChangeForMSC;
63 class G4ParticleDefinition; <<
64 class G4VEnergyLossProcess; <<
65 67
66 class G4VMscModel : public G4VEmModel 68 class G4VMscModel : public G4VEmModel
67 { 69 {
68 70
69 public: 71 public:
70 72
71 explicit G4VMscModel(const G4String& nam); << 73 G4VMscModel(const G4String& nam);
72 74
73 ~G4VMscModel() override; << 75 virtual ~G4VMscModel();
74 76
75 virtual G4double ComputeTruePathLengthLimit( 77 virtual G4double ComputeTruePathLengthLimit(const G4Track& track,
76 G4double& stepLimit) = 0; << 78 G4double& stepLimit);
77 79
78 virtual G4double ComputeGeomPathLength(G4dou << 80 virtual G4double ComputeGeomPathLength(G4double truePathLength);
79 81
80 virtual G4double ComputeTrueStepLength(G4dou << 82 virtual G4double ComputeTrueStepLength(G4double geomPathLength);
81 83
82 virtual G4ThreeVector& SampleScattering(cons 84 virtual G4ThreeVector& SampleScattering(const G4ThreeVector&,
83 G4double safety) = 0; << 85 G4double safety);
84 <<
85 void InitialiseParameters(const G4ParticleDe <<
86 <<
87 void DumpParameters(std::ostream& out) const <<
88 86
89 // empty method 87 // empty method
90 void SampleSecondaries(std::vector<G4Dynamic << 88 virtual void SampleSecondaries(std::vector<G4DynamicParticle*>*,
91 const G4MaterialCutsCouple*, << 89 const G4MaterialCutsCouple*,
92 const G4DynamicParticle*, << 90 const G4DynamicParticle*,
93 G4double tmin, G4double tmax) override; << 91 G4double tmin,
>> 92 G4double tmax);
94 93
95 //========================================== 94 //================================================================
96 // Set parameters of multiple scattering mo 95 // Set parameters of multiple scattering models
97 //========================================== 96 //================================================================
98 97
99 inline void SetStepLimitType(G4MscStepLimitT 98 inline void SetStepLimitType(G4MscStepLimitType);
100 99
101 inline void SetLateralDisplasmentFlag(G4bool 100 inline void SetLateralDisplasmentFlag(G4bool val);
102 101
103 inline void SetRangeFactor(G4double); 102 inline void SetRangeFactor(G4double);
104 103
105 inline void SetGeomFactor(G4double); 104 inline void SetGeomFactor(G4double);
106 105
107 inline void SetSkin(G4double); 106 inline void SetSkin(G4double);
108 107
109 inline void SetLambdaLimit(G4double); <<
110 <<
111 inline void SetSafetyFactor(G4double); <<
112 <<
113 inline void SetSampleZ(G4bool); 108 inline void SetSampleZ(G4bool);
114 109
115 //========================================== 110 //================================================================
116 // Get/Set access to Physics Tables 111 // Get/Set access to Physics Tables
117 //========================================== 112 //================================================================
118 113
119 inline G4VEnergyLossProcess* GetIonisation() 114 inline G4VEnergyLossProcess* GetIonisation() const;
120 115
121 inline void SetIonisation(G4VEnergyLossProce 116 inline void SetIonisation(G4VEnergyLossProcess*,
122 const G4ParticleDefinition* part); 117 const G4ParticleDefinition* part);
123 118
124 //========================================== 119 //================================================================
125 // Run time methods 120 // Run time methods
126 //========================================== 121 //================================================================
127 122
128 protected: 123 protected:
129 124
130 // initialisation of the ParticleChange for 125 // initialisation of the ParticleChange for the model
131 // initialisation of interface with geometry 126 // initialisation of interface with geometry and ionisation
132 G4ParticleChangeForMSC* 127 G4ParticleChangeForMSC*
133 GetParticleChangeForMSC(const G4ParticleDefi << 128 GetParticleChangeForMSC(const G4ParticleDefinition* p = 0);
134 129
135 // convert true length to geometry length 130 // convert true length to geometry length
136 inline G4double ConvertTrueToGeom(G4double& 131 inline G4double ConvertTrueToGeom(G4double& tLength, G4double& gLength);
137 132
138 // should be set before initialisation <<
139 inline void SetUseSplineForMSC(G4bool val); <<
140 <<
141 public: 133 public:
142 134
143 // compute safety 135 // compute safety
144 inline G4double ComputeSafety(const G4ThreeV << 136 inline G4double ComputeSafety(const G4ThreeVector& position, G4double limit);
145 G4double limit= DBL_MAX); <<
146 137
147 // compute linear distance to a geometry bou 138 // compute linear distance to a geometry boundary
148 inline G4double ComputeGeomLimit(const G4Tra 139 inline G4double ComputeGeomLimit(const G4Track&, G4double& presafety,
149 G4double limit); 140 G4double limit);
150 141
151 G4double GetDEDX(const G4ParticleDefinition* << 142 inline G4double GetDEDX(const G4ParticleDefinition* part,
152 G4double kineticEner << 143 G4double kineticEnergy,
153 const G4MaterialCuts << 144 const G4MaterialCutsCouple* couple);
154 <<
155 G4double GetDEDX(const G4ParticleDefinition* <<
156 G4double kineticEner <<
157 const G4MaterialCuts <<
158 G4double logKineticE <<
159 <<
160 G4double GetRange(const G4ParticleDefinition <<
161 G4double kineticEne <<
162 const G4MaterialCut <<
163 145
164 G4double GetRange(const G4ParticleDefinition << 146 inline G4double GetRange(const G4ParticleDefinition* part,
165 G4double kineticEne 147 G4double kineticEnergy,
166 const G4MaterialCut << 148 const G4MaterialCutsCouple* couple);
167 G4double logKinetic <<
168 149
169 G4double GetEnergy(const G4ParticleDefinitio << 150 inline G4double GetEnergy(const G4ParticleDefinition* part,
170 G4double range, 151 G4double range,
171 const G4MaterialCutsCouple* couple); 152 const G4MaterialCutsCouple* couple);
172 153
173 // G4MaterialCutsCouple should be defined be 154 // G4MaterialCutsCouple should be defined before call to this method
174 inline 155 inline
175 G4double GetTransportMeanFreePath(const G4Pa 156 G4double GetTransportMeanFreePath(const G4ParticleDefinition* part,
176 G4double k << 157 G4double kinEnergy);
177 158
178 inline << 159 private:
179 G4double GetTransportMeanFreePath(const G4Pa <<
180 G4double k <<
181 G4double l <<
182 160
183 // hide assignment operator 161 // hide assignment operator
184 G4VMscModel & operator=(const G4VMscModel & << 162 G4VMscModel & operator=(const G4VMscModel &right);
185 G4VMscModel(const G4VMscModel&) = delete; << 163 G4VMscModel(const G4VMscModel&);
186 <<
187 private: <<
188 164
189 G4SafetyHelper* safetyHelper = nullptr; << 165 G4SafetyHelper* safetyHelper;
190 G4VEnergyLossProcess* ionisation = nullptr; << 166 G4VEnergyLossProcess* ionisation;
191 const G4ParticleDefinition* currentPart = nu << 167 const G4ParticleDefinition* currentPart;
192 << 168 G4LossTableManager* man;
193 G4double dedx = 0.0; << 169
194 G4double localtkin = 0.0; << 170 G4double dedx;
195 G4double localrange = DBL_MAX; << 171 G4double localtkin;
>> 172 G4double localrange;
196 173
197 protected: 174 protected:
198 175
199 G4double facrange = 0.04; << 176 G4double facrange;
200 G4double facgeom = 2.5; << 177 G4double facgeom;
201 G4double facsafety = 0.6; << 178 G4double facsafety;
202 G4double skin = 1.0; << 179 G4double skin;
203 G4double dtrl = 0.05; << 180 G4double dtrl;
204 G4double lambdalimit; 181 G4double lambdalimit;
205 G4double geomMin; 182 G4double geomMin;
206 G4double geomMax; 183 G4double geomMax;
207 184
208 G4ThreeVector fDisplacement; << 185 G4ThreeVector fDisplacement;
209 G4MscStepLimitType steppingAlgorithm; 186 G4MscStepLimitType steppingAlgorithm;
210 187
211 G4bool samplez = false; << 188 G4bool samplez;
212 G4bool latDisplasment = true; << 189 G4bool latDisplasment;
213 190
214 private: <<
215 <<
216 G4bool useSpline = true; <<
217 }; 191 };
218 192
219 //....oooOO0OOooo........oooOO0OOooo........oo 193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
220 //....oooOO0OOooo........oooOO0OOooo........oo 194 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
221 195
222 inline void G4VMscModel::SetLateralDisplasment 196 inline void G4VMscModel::SetLateralDisplasmentFlag(G4bool val)
223 { 197 {
224 if(!IsLocked()) { latDisplasment = val; } << 198 latDisplasment = val;
225 } 199 }
226 200
227 //....oooOO0OOooo........oooOO0OOooo........oo 201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
228 202
229 inline void G4VMscModel::SetSkin(G4double val) 203 inline void G4VMscModel::SetSkin(G4double val)
230 { 204 {
231 if(!IsLocked()) { skin = val; } << 205 skin = val;
232 } 206 }
233 207
234 //....oooOO0OOooo........oooOO0OOooo........oo 208 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
235 209
236 inline void G4VMscModel::SetRangeFactor(G4doub 210 inline void G4VMscModel::SetRangeFactor(G4double val)
237 { 211 {
238 if(!IsLocked()) { facrange = val; } << 212 facrange = val;
239 } 213 }
240 214
241 //....oooOO0OOooo........oooOO0OOooo........oo 215 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
242 216
243 inline void G4VMscModel::SetGeomFactor(G4doubl 217 inline void G4VMscModel::SetGeomFactor(G4double val)
244 { 218 {
245 if(!IsLocked()) { facgeom = val; } << 219 facgeom = val;
246 } <<
247 <<
248 //....oooOO0OOooo........oooOO0OOooo........oo <<
249 <<
250 inline void G4VMscModel::SetLambdaLimit(G4doub <<
251 { <<
252 if(!IsLocked()) { lambdalimit = val; } <<
253 } <<
254 <<
255 //....oooOO0OOooo........oooOO0OOooo........oo <<
256 <<
257 inline void G4VMscModel::SetSafetyFactor(G4dou <<
258 { <<
259 if(!IsLocked()) { facsafety = val; } <<
260 } 220 }
261 221
262 //....oooOO0OOooo........oooOO0OOooo........oo 222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
263 223
264 inline void G4VMscModel::SetStepLimitType(G4Ms 224 inline void G4VMscModel::SetStepLimitType(G4MscStepLimitType val)
265 { 225 {
266 if(!IsLocked()) { steppingAlgorithm = val; } << 226 steppingAlgorithm = val;
267 } 227 }
268 228
269 //....oooOO0OOooo........oooOO0OOooo........oo 229 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
270 230
271 inline void G4VMscModel::SetSampleZ(G4bool val 231 inline void G4VMscModel::SetSampleZ(G4bool val)
272 { 232 {
273 if(!IsLocked()) { samplez = val; } << 233 samplez = val;
274 } 234 }
275 235
276 //....oooOO0OOooo........oooOO0OOooo........oo 236 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
277 237
278 inline G4double G4VMscModel::ComputeSafety(con 238 inline G4double G4VMscModel::ComputeSafety(const G4ThreeVector& position,
279 G4double limit) << 239 G4double)
280 { 240 {
281 return safetyHelper->ComputeSafety(position << 241 return safetyHelper->ComputeSafety(position);
282 } 242 }
283 243
284 //....oooOO0OOooo........oooOO0OOooo........oo 244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
285 245
286 inline G4double G4VMscModel::ConvertTrueToGeom 246 inline G4double G4VMscModel::ConvertTrueToGeom(G4double& tlength,
287 G4double& glength) 247 G4double& glength)
288 { 248 {
289 glength = ComputeGeomPathLength(tlength); 249 glength = ComputeGeomPathLength(tlength);
290 // should return true length 250 // should return true length
291 return tlength; 251 return tlength;
292 } 252 }
293 253
294 //....oooOO0OOooo........oooOO0OOooo........oo 254 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
295 255
296 inline G4double G4VMscModel::ComputeGeomLimit( 256 inline G4double G4VMscModel::ComputeGeomLimit(const G4Track& track,
297 G4double& presafety, 257 G4double& presafety,
298 G4double limit) 258 G4double limit)
299 { 259 {
300 return safetyHelper->CheckNextStep( << 260 G4double res = geomMax;
>> 261 if(track.GetVolume() != safetyHelper->GetWorldVolume()) {
>> 262 res = safetyHelper->CheckNextStep(
301 track.GetStep()->GetPreStepPoint()-> 263 track.GetStep()->GetPreStepPoint()->GetPosition(),
302 track.GetMomentumDirection(), 264 track.GetMomentumDirection(),
303 limit, presafety); 265 limit, presafety);
>> 266 }
>> 267 return res;
304 } 268 }
305 269
306 //....oooOO0OOooo........oooOO0OOooo........oo 270 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
307 271
308 inline G4VEnergyLossProcess* G4VMscModel::GetI << 272 inline G4double
>> 273 G4VMscModel::GetDEDX(const G4ParticleDefinition* part,
>> 274 G4double kinEnergy, const G4MaterialCutsCouple* couple)
309 { 275 {
310 return ionisation; << 276 G4double x;
>> 277 if(ionisation) { x = ionisation->GetDEDX(kinEnergy, couple); }
>> 278 else {
>> 279 G4double q = part->GetPDGCharge()/CLHEP::eplus;
>> 280 x = dedx*q*q;
>> 281 }
>> 282 return x;
311 } 283 }
312 284
313 //....oooOO0OOooo........oooOO0OOooo........oo 285 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
314 286
315 inline void G4VMscModel::SetIonisation(G4VEner << 287 inline G4double
316 const G4ParticleDefinition* par << 288 G4VMscModel::GetRange(const G4ParticleDefinition* part,
>> 289 G4double kinEnergy, const G4MaterialCutsCouple* couple)
317 { 290 {
318 ionisation = p; << 291 //G4cout << "G4VMscModel::GetRange E(MeV)= " << kinEnergy << " " << ionisation
319 currentPart = part; << 292 // << " " << part->GetParticleName()
>> 293 // << G4endl;
>> 294 localtkin = kinEnergy;
>> 295 if(ionisation) {
>> 296 localrange = ionisation->GetRangeForLoss(kinEnergy, couple);
>> 297 } else {
>> 298 G4double q = part->GetPDGCharge()/CLHEP::eplus;
>> 299 localrange = kinEnergy/(dedx*q*q*couple->GetMaterial()->GetDensity());
>> 300 }
>> 301 //G4cout << "R(mm)= " << localrange << " " << ionisation << G4endl;
>> 302 return localrange;
320 } 303 }
321 304
322 //....oooOO0OOooo........oooOO0OOooo........oo 305 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
323 306
324 inline G4double 307 inline G4double
325 G4VMscModel::GetTransportMeanFreePath(const G4 << 308 G4VMscModel::GetEnergy(const G4ParticleDefinition* part,
326 G4double << 309 G4double range, const G4MaterialCutsCouple* couple)
327 { 310 {
328 G4double x; << 311 G4double e;
329 if (nullptr != xSectionTable) { << 312 //G4cout << "G4VMscModel::GetEnergy R(mm)= " << range << " " << ionisation
330 x = pFactor*(*xSectionTable)[basedCoupleIn << 313 // << " Rlocal(mm)= " << localrange << " Elocal(MeV)= " << localtkin
331 } else { << 314 // << G4endl;
332 x = pFactor*CrossSectionPerVolume(pBaseMat << 315 if(ionisation) { e = ionisation->GetKineticEnergy(range, couple); }
>> 316 else {
>> 317 e = localtkin;
>> 318 if(localrange > range) {
>> 319 G4double q = part->GetPDGCharge()/CLHEP::eplus;
>> 320 e -= (localrange - range)*dedx*q*q*couple->GetMaterial()->GetDensity();
>> 321 }
333 } 322 }
334 return (x > 0.0) ? 1.0/x : DBL_MAX; << 323 return e;
335 } 324 }
336 325
337 //....oooOO0OOooo........oooOO0OOooo........oo 326 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
338 327
339 inline G4double << 328 inline G4VEnergyLossProcess* G4VMscModel::GetIonisation() const
340 G4VMscModel::GetTransportMeanFreePath(const G4 <<
341 G4double <<
342 { 329 {
343 G4double x; << 330 return ionisation;
344 if (nullptr != xSectionTable) { <<
345 x = pFactor*(*xSectionTable)[basedCoupleIn <<
346 } else { <<
347 x = pFactor*CrossSectionPerVolume(pBaseMat <<
348 } <<
349 return (x > 0.0) ? 1.0/x : DBL_MAX; <<
350 } 331 }
351 332
352 //....oooOO0OOooo........oooOO0OOooo........oo 333 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
353 334
354 inline void G4VMscModel::SetUseSplineForMSC(G4 << 335 inline void G4VMscModel::SetIonisation(G4VEnergyLossProcess* p,
>> 336 const G4ParticleDefinition* part)
355 { 337 {
356 useSpline = val; << 338 ionisation = p;
>> 339 currentPart = part;
>> 340 }
>> 341
>> 342 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 343
>> 344 inline G4double
>> 345 G4VMscModel::GetTransportMeanFreePath(const G4ParticleDefinition* part,
>> 346 G4double ekin)
>> 347 {
>> 348 G4double x;
>> 349 if(xSectionTable) {
>> 350 G4int idx = CurrentCouple()->GetIndex();
>> 351 x = (*xSectionTable)[(*theDensityIdx)[idx]]->Value(ekin, idxTable)
>> 352 *(*theDensityFactor)[idx]/(ekin*ekin);
>> 353 } else {
>> 354 x = CrossSectionPerVolume(CurrentCouple()->GetMaterial(), part, ekin,
>> 355 0.0, DBL_MAX);
>> 356 }
>> 357 if(0.0 >= x) { x = DBL_MAX; }
>> 358 else { x = 1.0/x; }
>> 359 return x;
357 } 360 }
358 361
359 //....oooOO0OOooo........oooOO0OOooo........oo 362 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
360 363
361 #endif 364 #endif
>> 365
362 366