| 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: G4WentzelVIModel.hh,v 1.7 2008/08/04 08:49:09 vnivanch Exp $
>> 27 // GEANT4 tag $Name: geant4-09-02 $
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: G4WentzelVIModel 35 // File name: G4WentzelVIModel
34 // 36 //
35 // Author: V.Ivanchenko 37 // Author: V.Ivanchenko
36 // 38 //
37 // Creation date: 09.04.2008 from G4MuMscModel 39 // Creation date: 09.04.2008 from G4MuMscModel
38 // 40 //
39 // Modifications: 41 // Modifications:
40 // 27-05-2010 V.Ivanchenko added G4WentzelOKan << 42 //
41 // compute cross sections and sam <<
42 // 43 //
43 // Class Description: 44 // Class Description:
44 // 45 //
45 // Implementation of the model of multiple sca 46 // Implementation of the model of multiple scattering based on
46 // G.Wentzel, Z. Phys. 40 (1927) 590. 47 // G.Wentzel, Z. Phys. 40 (1927) 590.
47 // H.W.Lewis, Phys Rev 78 (1950) 526. 48 // H.W.Lewis, Phys Rev 78 (1950) 526.
48 // J.M. Fernandez-Varea et al., NIM B73 (1993) 49 // J.M. Fernandez-Varea et al., NIM B73 (1993) 447.
49 // L.Urban, CERN-OPEN-2006-077. 50 // L.Urban, CERN-OPEN-2006-077.
50 51
51 // ------------------------------------------- 52 // -------------------------------------------------------------------
52 // 53 //
53 54
54 #ifndef G4WentzelVIModel_h 55 #ifndef G4WentzelVIModel_h
55 #define G4WentzelVIModel_h 1 56 #define G4WentzelVIModel_h 1
56 57
57 //....oooOO0OOooo........oooOO0OOooo........oo 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58 59
59 #include "G4VMscModel.hh" 60 #include "G4VMscModel.hh"
>> 61 #include "G4PhysicsTable.hh"
>> 62 #include "G4MscStepLimitType.hh"
60 #include "G4MaterialCutsCouple.hh" 63 #include "G4MaterialCutsCouple.hh"
61 #include "G4WentzelOKandVIxSection.hh" << 64 #include "G4NistManager.hh"
>> 65
>> 66 class G4LossTableManager;
>> 67 class G4ParticleChangeForMSC;
>> 68 class G4SafetyHelper;
>> 69 class G4ParticleDefinition;
62 70
63 //....oooOO0OOooo........oooOO0OOooo........oo 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
64 72
65 class G4WentzelVIModel : public G4VMscModel 73 class G4WentzelVIModel : public G4VMscModel
66 { 74 {
67 75
68 public: 76 public:
69 77
70 explicit G4WentzelVIModel(G4bool comb=true, << 78 G4WentzelVIModel(const G4String& nam = "WentzelVIUni");
71 <<
72 ~G4WentzelVIModel() override; <<
73 <<
74 void Initialise(const G4ParticleDefinition*, <<
75 79
76 void InitialiseLocal(const G4ParticleDefinit << 80 virtual ~G4WentzelVIModel();
77 G4VEmModel* masterModel) override; <<
78 81
79 void StartTracking(G4Track*) override; << 82 void Initialise(const G4ParticleDefinition*, const G4DataVector&);
80 83
81 G4double ComputeCrossSectionPerAtom(const G4 84 G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
82 G4double KineticEnergy, 85 G4double KineticEnergy,
83 G4double AtomicNumber, 86 G4double AtomicNumber,
84 G4double AtomicWeight=0., 87 G4double AtomicWeight=0.,
85 G4double cut = DBL_MAX, 88 G4double cut = DBL_MAX,
86 G4double emax= DBL_MAX) override << 89 G4double emax= DBL_MAX);
87 90
88 G4ThreeVector& SampleScattering(const G4Thre << 91 void SampleScattering(const G4DynamicParticle*, G4double safety);
89 G4double safety) override; <<
90 92
91 G4double << 93 void SampleSecondaries(std::vector<G4DynamicParticle*>*,
92 ComputeTruePathLengthLimit(const G4Track& tr << 94 const G4MaterialCutsCouple*,
93 G4double& currentMinimalStep) overr << 95 const G4DynamicParticle*,
>> 96 G4double,
>> 97 G4double);
94 98
95 G4double ComputeGeomPathLength(G4double true << 99 G4double ComputeTruePathLengthLimit(const G4Track& track,
>> 100 G4PhysicsTable* theLambdaTable,
>> 101 G4double currentMinimalStep);
96 102
97 G4double ComputeTrueStepLength(G4double geom << 103 G4double ComputeGeomPathLength(G4double truePathLength);
98 104
99 // defines low energy limit on energy transf << 105 G4double ComputeTrueStepLength(G4double geomStepLength);
100 void SetFixedCut(G4double); <<
101 106
102 // low energy limit on energy transfer to at << 107 private:
103 G4double GetFixedCut() const; <<
104 <<
105 // access to cross section class <<
106 void SetWVICrossSection(G4WentzelOKandVIxSec <<
107 <<
108 G4WentzelOKandVIxSection* GetWVICrossSection <<
109 108
110 void SetUseSecondMoment(G4bool); << 109 G4double ComputeTransportXSectionPerVolume();
111 110
112 G4bool UseSecondMoment() const; << 111 G4double ComputeXSectionPerVolume();
113 112
114 G4PhysicsTable* GetSecondMomentTable(); << 113 void ComputeMaxElectronScattering(G4double cut);
115 114
116 G4double SecondMoment(const G4ParticleDefini << 115 inline G4double GetLambda(G4double kinEnergy);
117 const G4MaterialCutsCouple*, <<
118 G4double kineticEnergy); <<
119 116
120 void SetSingleScatteringFactor(G4double); << 117 inline void SetupParticle(const G4ParticleDefinition*);
121 118
122 void DefineMaterial(const G4MaterialCutsCoup << 119 inline void SetupKinematic(G4double kinEnergy, G4double cut);
>> 120
>> 121 inline void SetupTarget(G4double Z, G4double kinEnergy);
123 122
124 G4WentzelVIModel & operator=(const G4Wentzel << 123 inline void DefineMaterial(const G4MaterialCutsCouple*);
125 G4WentzelVIModel(const G4WentzelVIModel&) = <<
126 124
127 protected: << 125 // hide assignment operator
>> 126 G4WentzelVIModel & operator=(const G4WentzelVIModel &right);
>> 127 G4WentzelVIModel(const G4WentzelVIModel&);
128 128
129 G4double ComputeTransportXSectionPerVolume(G << 129 const G4ParticleDefinition* theProton;
>> 130 const G4ParticleDefinition* theElectron;
>> 131 const G4ParticleDefinition* thePositron;
130 132
131 inline void SetupParticle(const G4ParticleDe << 133 G4ParticleChangeForMSC* fParticleChange;
132 134
133 private: << 135 G4SafetyHelper* safetyHelper;
>> 136 G4PhysicsTable* theLambdaTable;
>> 137 G4PhysicsTable* theLambda2Table;
>> 138 G4LossTableManager* theManager;
>> 139 const G4DataVector* currentCuts;
134 140
135 G4double ComputeSecondMoment(const G4Particl << 141 G4NistManager* fNistManager;
136 G4double kineticEnergy); <<
137 142
138 protected: << 143 G4double numlimit;
>> 144 G4double tlimitminfix;
>> 145 G4double invsqrt12;
139 146
140 G4WentzelOKandVIxSection* wokvi; << 147 // cash
141 const G4MaterialCutsCouple* currentCouple = << 148 G4double preKinEnergy;
142 const G4Material* currentMaterial = nullptr; << 149 G4double ecut;
143 << 150 G4double lambda0;
144 const G4ParticleDefinition* particle = nullp << 151 G4double tPathLength;
145 G4ParticleChangeForMSC* fParticleChange = nu << 152 G4double zPathLength;
146 const G4DataVector* currentCuts = nullptr; << 153 G4double lambdaeff;
147 G4PhysicsTable* fSecondMoments = nullptr; << 154 G4double currentRange;
>> 155 G4double par1;
>> 156 G4double par2;
>> 157 G4double par3;
148 158
149 G4double lowEnergyLimit; << 159 G4double xtsec;
150 G4double tlimitminfix; << 160 std::vector<G4double> xsecn;
151 G4double ssFactor = 1.05; << 161 std::vector<G4double> prob;
152 G4double invssFactor = 1.0; << 162 G4int nelments;
153 163
154 // cache kinematics << 164 G4int nbins;
155 G4double preKinEnergy = 0.0; << 165 G4int nwarnings;
156 G4double tPathLength = 0.0; << 166 G4int nwarnlimit;
157 G4double zPathLength = 0.0; <<
158 G4double lambdaeff = 0.0; <<
159 G4double currentRange = 0.0; <<
160 G4double cosTetMaxNuc = 0.0; <<
161 167
162 G4double fixedCut = -1.0; << 168 G4int currentMaterialIndex;
163 169
164 // cache kinematics << 170 const G4MaterialCutsCouple* currentCouple;
165 G4double effKinEnergy = 0.0; << 171 const G4Material* currentMaterial;
166 172
167 // single scattering parameters 173 // single scattering parameters
168 G4double cosThetaMin = 1.0; << 174 G4double coeff;
169 G4double cosThetaMax = -1.0; << 175 G4double constn;
170 G4double xtsec = 0.0; << 176 G4double cosThetaMin;
171 << 177 G4double cosThetaMax;
172 G4int currentMaterialIndex = 0; << 178 G4double cosTetMaxNuc;
173 size_t idx2 = 0; << 179 G4double cosTetMaxNuc2;
>> 180 G4double cosTetMaxElec;
>> 181 G4double cosTetMaxElec2;
>> 182 G4double q2Limit;
>> 183 G4double alpha2;
>> 184 G4double a0;
>> 185
>> 186 // projectile
>> 187 const G4ParticleDefinition* particle;
>> 188
>> 189 G4double chargeSquare;
>> 190 G4double spin;
>> 191 G4double mass;
>> 192 G4double tkin;
>> 193 G4double mom2;
>> 194 G4double invbeta2;
>> 195 G4double etag;
>> 196 G4double lowEnergyLimit;
174 197
175 // data for single scattering mode << 198 // target
176 G4int nelments = 0; << 199 G4double targetZ;
>> 200 G4double screenZ;
>> 201 G4double formfactA;
>> 202 G4double FF[100];
177 203
178 // flags 204 // flags
179 G4bool singleScatteringMode; << 205 G4bool isInitialized;
180 G4bool isCombined; << 206 G4bool inside;
181 G4bool useSecondMoment; <<
182 <<
183 std::vector<G4double> xsecn; <<
184 std::vector<G4double> prob; <<
185 }; 207 };
186 208
187 //....oooOO0OOooo........oooOO0OOooo........oo 209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
188 //....oooOO0OOooo........oooOO0OOooo........oo 210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
189 211
190 inline void G4WentzelVIModel::SetupParticle(co << 212 inline
191 { << 213 void G4WentzelVIModel::DefineMaterial(const G4MaterialCutsCouple* cup)
192 // Initialise mass and charge << 214 {
193 if(p != particle) { << 215 if(cup != currentCouple) {
194 particle = p; << 216 currentCouple = cup;
195 wokvi->SetupParticle(p); << 217 currentMaterial = cup->GetMaterial();
>> 218 currentMaterialIndex = currentCouple->GetIndex();
196 } 219 }
197 } 220 }
198 221
199 //....oooOO0OOooo........oooOO0OOooo........oo << 222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
200 <<
201 inline void G4WentzelVIModel::SetFixedCut(G4do <<
202 { <<
203 fixedCut = val; <<
204 } <<
205 <<
206 //....oooOO0OOooo........oooOO0OOooo........oo <<
207 <<
208 inline G4double G4WentzelVIModel::GetFixedCut( <<
209 { <<
210 return fixedCut; <<
211 } <<
212 <<
213 //....oooOO0OOooo........oooOO0OOooo........oo <<
214 223
215 inline void G4WentzelVIModel::SetWVICrossSecti << 224 inline
>> 225 G4double G4WentzelVIModel::GetLambda(G4double e)
216 { 226 {
217 if(ptr != wokvi) { << 227 G4double x;
218 delete wokvi; << 228 if(theLambdaTable) {
219 wokvi = ptr; << 229 G4bool b;
>> 230 x = ((*theLambdaTable)[currentMaterialIndex])->GetValue(e, b);
>> 231 } else {
>> 232 x = CrossSection(currentCouple,particle,e,
>> 233 (*currentCuts)[currentMaterialIndex]);
220 } 234 }
>> 235 if(x > DBL_MIN) x = 1./x;
>> 236 else x = DBL_MAX;
>> 237 return x;
221 } 238 }
222 239
223 //....oooOO0OOooo........oooOO0OOooo........oo 240 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
224 241
225 inline G4WentzelOKandVIxSection* G4WentzelVIMo << 242 inline
226 { << 243 void G4WentzelVIModel::SetupParticle(const G4ParticleDefinition* p)
227 return wokvi; <<
228 } <<
229 <<
230 //....oooOO0OOooo........oooOO0OOooo........oo <<
231 <<
232 inline void G4WentzelVIModel::SetUseSecondMome <<
233 { <<
234 useSecondMoment = val; <<
235 } <<
236 <<
237 //....oooOO0OOooo........oooOO0OOooo........oo <<
238 <<
239 inline G4bool G4WentzelVIModel::UseSecondMomen <<
240 { 244 {
241 return useSecondMoment; << 245 // Initialise mass and charge
>> 246 if(p != particle) {
>> 247 particle = p;
>> 248 mass = particle->GetPDGMass();
>> 249 spin = particle->GetPDGSpin();
>> 250 G4double q = particle->GetPDGCharge()/eplus;
>> 251 chargeSquare = q*q;
>> 252 tkin = 0.0;
>> 253 lowEnergyLimit = keV*mass/electron_mass_c2;
>> 254 }
242 } 255 }
243 256
244 //....oooOO0OOooo........oooOO0OOooo........oo 257 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
245 258
246 inline G4PhysicsTable* G4WentzelVIModel::GetSe << 259 inline void G4WentzelVIModel::SetupKinematic(G4double ekin, G4double cut)
247 { 260 {
248 return fSecondMoments; << 261 if(ekin != tkin || ecut != cut) {
>> 262 tkin = ekin;
>> 263 mom2 = tkin*(tkin + 2.0*mass);
>> 264 invbeta2 = 1.0 + mass*mass/mom2;
>> 265 cosTetMaxNuc = cosThetaMax;
>> 266 if(ekin <= 10.*cut && mass < MeV) {
>> 267 cosTetMaxNuc = ekin*(cosThetaMax + 1.0)/(10.*cut) - 1.0;
>> 268 }
>> 269 ComputeMaxElectronScattering(cut);
>> 270 }
249 } 271 }
250 272
251 //....oooOO0OOooo........oooOO0OOooo........oo 273 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
252 << 274
253 inline G4double << 275 inline void G4WentzelVIModel::SetupTarget(G4double Z, G4double e)
254 G4WentzelVIModel::SecondMoment(const G4Particl <<
255 const G4MaterialCutsCouple* coupl <<
256 G4double ekin) <<
257 { 276 {
258 G4double x = 0.0; << 277 if(Z != targetZ || e != etag) {
259 if(useSecondMoment) { << 278 etag = e;
260 DefineMaterial(couple); << 279 targetZ = Z;
261 x = (fSecondMoments) ? << 280 G4int iz= G4int(Z);
262 (*fSecondMoments)[(*theDensityIdx)[curre << 281 if(iz > 99) iz = 99;
263 *(*theDensityFactor)[currentMaterialInde << 282 G4double x = fNistManager->GetZ13(iz);
264 : ComputeSecondMoment(part, ekin); << 283 screenZ = a0*x*x/mom2;
265 } << 284 if(iz > 1) screenZ *=(1.13 + 3.76*invbeta2*Z*Z*chargeSquare*alpha2);
266 return x; << 285 // screenZ = a0*x*x*(1.13 + 3.76*Z*Z*chargeSquare*alpha2)/mom2;
267 } << 286 // A.V. Butkevich et al., NIM A 488 (2002) 282
>> 287 formfactA = FF[iz];
>> 288 if(formfactA == 0.0) {
>> 289 x = fNistManager->GetA27(iz);
>> 290 formfactA = constn*x*x;
>> 291 FF[iz] = formfactA;
>> 292 }
>> 293 formfactA *= mom2;
>> 294 cosTetMaxNuc2 = cosTetMaxNuc;
>> 295 /*
>> 296 G4double ee = 10.*eV*Z;
>> 297 if(1 == iz) ee *= 2.0;
>> 298 G4double z = std::min(cosTetMaxElec, 1.0 - std::max(ecut,ee)*amu_c2
>> 299 *fNistManager->GetAtomicMassAmu(iz)/mom2);
>> 300 cosTetMaxElec2 = std::max(cosTetMaxNuc2, z);
>> 301 */
>> 302 cosTetMaxElec2 = cosTetMaxElec;
>> 303 }
>> 304 }
268 305
269 //....oooOO0OOooo........oooOO0OOooo........oo 306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
270 307
271 #endif 308 #endif
272 309
273 310