| 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 // 26 //
27 //-------------------------------------------- 27 //---------------------------------------------------------------------------
28 // 28 //
29 // ClassName: G4EnergyLossForExtrapolator 29 // ClassName: G4EnergyLossForExtrapolator
30 // 30 //
31 // Description: This class provide calculatio 31 // Description: This class provide calculation of energy loss, fluctuation,
32 // and msc angle 32 // and msc angle
33 // 33 //
34 // Author: 09.12.04 V.Ivanchenko 34 // Author: 09.12.04 V.Ivanchenko
35 // 35 //
36 // Modification: 36 // Modification:
37 // 08-04-05 Rename Propogator -> Extrapolator 37 // 08-04-05 Rename Propogator -> Extrapolator (V.Ivanchenko)
38 // 16-03-06 Add muon tables and fix bug in uni 38 // 16-03-06 Add muon tables and fix bug in units (V.Ivanchenko)
39 // 21-03-06 Add verbosity defined in the const 39 // 21-03-06 Add verbosity defined in the constructor and Initialisation
40 // start only when first public metho 40 // start only when first public method is called (V.Ivanchenko)
41 // 03-05-06 Remove unused pointer G4Material* 41 // 03-05-06 Remove unused pointer G4Material* from number of methods (VI)
42 // 12-05-06 SEt linLossLimit=0.001 (VI) 42 // 12-05-06 SEt linLossLimit=0.001 (VI)
43 // 43 //
44 //-------------------------------------------- 44 //----------------------------------------------------------------------------
45 // 45 //
46 46
47 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48 48
49 #include "G4EnergyLossForExtrapolator.hh" 49 #include "G4EnergyLossForExtrapolator.hh"
50 #include "G4PhysicalConstants.hh" 50 #include "G4PhysicalConstants.hh"
51 #include "G4SystemOfUnits.hh" 51 #include "G4SystemOfUnits.hh"
52 #include "G4ParticleDefinition.hh" 52 #include "G4ParticleDefinition.hh"
53 #include "G4Material.hh" 53 #include "G4Material.hh"
54 #include "G4MaterialCutsCouple.hh" 54 #include "G4MaterialCutsCouple.hh"
55 #include "G4Electron.hh" 55 #include "G4Electron.hh"
56 #include "G4Positron.hh" 56 #include "G4Positron.hh"
57 #include "G4Proton.hh" 57 #include "G4Proton.hh"
58 #include "G4MuonPlus.hh" 58 #include "G4MuonPlus.hh"
59 #include "G4MuonMinus.hh" 59 #include "G4MuonMinus.hh"
60 #include "G4ParticleTable.hh" 60 #include "G4ParticleTable.hh"
61 61
62 //....oooOO0OOooo........oooOO0OOooo........oo 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63 63
64 #ifdef G4MULTITHREADED 64 #ifdef G4MULTITHREADED
65 G4Mutex G4EnergyLossForExtrapolator::extrMutex 65 G4Mutex G4EnergyLossForExtrapolator::extrMutex = G4MUTEX_INITIALIZER;
66 #endif 66 #endif
67 67
68 G4TablesForExtrapolator* G4EnergyLossForExtrap 68 G4TablesForExtrapolator* G4EnergyLossForExtrapolator::tables = nullptr;
69 69
70 G4EnergyLossForExtrapolator::G4EnergyLossForEx 70 G4EnergyLossForExtrapolator::G4EnergyLossForExtrapolator(G4int verb)
71 : maxEnergyTransfer(DBL_MAX), verbose(verb) 71 : maxEnergyTransfer(DBL_MAX), verbose(verb)
72 { 72 {
73 emin = 1.*CLHEP::MeV; 73 emin = 1.*CLHEP::MeV;
74 emax = 100.*CLHEP::TeV; 74 emax = 100.*CLHEP::TeV;
75 } 75 }
76 76
77 //....oooOO0OOooo........oooOO0OOooo........oo 77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78 78
79 G4EnergyLossForExtrapolator::~G4EnergyLossForE 79 G4EnergyLossForExtrapolator::~G4EnergyLossForExtrapolator()
80 { 80 {
81 if(isMaster) { 81 if(isMaster) {
82 delete tables; 82 delete tables;
83 tables = nullptr; 83 tables = nullptr;
84 } 84 }
85 } 85 }
86 86
87 //....oooOO0OOooo........oooOO0OOooo........oo 87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88 88
89 G4double 89 G4double
90 G4EnergyLossForExtrapolator::EnergyAfterStep(G 90 G4EnergyLossForExtrapolator::EnergyAfterStep(G4double kinEnergy,
91 G4double stepLength, 91 G4double stepLength,
92 const G4Material* mat, 92 const G4Material* mat,
93 const G4ParticleDefinition* par 93 const G4ParticleDefinition* part)
94 { 94 {
95 G4double kinEnergyFinal = kinEnergy; 95 G4double kinEnergyFinal = kinEnergy;
96 if(SetupKinematics(part, mat, kinEnergy)) { 96 if(SetupKinematics(part, mat, kinEnergy)) {
97 G4double step = TrueStepLength(kinEnergy,s 97 G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
98 G4double r = ComputeRange(kinEnergy,part, 98 G4double r = ComputeRange(kinEnergy,part,mat);
99 if(r <= step) { 99 if(r <= step) {
100 kinEnergyFinal = 0.0; 100 kinEnergyFinal = 0.0;
101 } else if(step < linLossLimit*r) { 101 } else if(step < linLossLimit*r) {
102 kinEnergyFinal -= step*ComputeDEDX(kinEn 102 kinEnergyFinal -= step*ComputeDEDX(kinEnergy,part,mat);
103 } else { 103 } else {
104 G4double r1 = r - step; 104 G4double r1 = r - step;
105 kinEnergyFinal = ComputeEnergy(r1,part,m 105 kinEnergyFinal = ComputeEnergy(r1,part,mat);
106 } 106 }
107 } 107 }
108 return kinEnergyFinal; 108 return kinEnergyFinal;
109 } 109 }
110 110
111 //....oooOO0OOooo........oooOO0OOooo........oo 111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
112 112
113 G4double 113 G4double
114 G4EnergyLossForExtrapolator::EnergyBeforeStep( 114 G4EnergyLossForExtrapolator::EnergyBeforeStep(G4double kinEnergy,
115 G4double stepLength, 115 G4double stepLength,
116 const G4Material* mat, 116 const G4Material* mat,
117 const G4ParticleDefinition* pa 117 const G4ParticleDefinition* part)
118 { 118 {
119 //G4cout << "G4EnergyLossForExtrapolator::En 119 //G4cout << "G4EnergyLossForExtrapolator::EnergyBeforeStep" << G4endl;
120 G4double kinEnergyFinal = kinEnergy; 120 G4double kinEnergyFinal = kinEnergy;
121 121
122 if(SetupKinematics(part, mat, kinEnergy)) { 122 if(SetupKinematics(part, mat, kinEnergy)) {
123 G4double step = TrueStepLength(kinEnergy,s 123 G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
124 G4double r = ComputeRange(kinEnergy,part,m 124 G4double r = ComputeRange(kinEnergy,part,mat);
125 125
126 if(step < linLossLimit*r) { 126 if(step < linLossLimit*r) {
127 kinEnergyFinal += step*ComputeDEDX(kinEn 127 kinEnergyFinal += step*ComputeDEDX(kinEnergy,part,mat);
128 } else { 128 } else {
129 G4double r1 = r + step; 129 G4double r1 = r + step;
130 kinEnergyFinal = ComputeEnergy(r1,part,m 130 kinEnergyFinal = ComputeEnergy(r1,part,mat);
131 } 131 }
132 } 132 }
133 return kinEnergyFinal; 133 return kinEnergyFinal;
134 } 134 }
135 135
136 //....oooOO0OOooo........oooOO0OOooo........oo 136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
137 137
138 G4double 138 G4double
139 G4EnergyLossForExtrapolator::TrueStepLength(G4 139 G4EnergyLossForExtrapolator::TrueStepLength(G4double kinEnergy,
140 G4double stepLength, 140 G4double stepLength,
141 const G4Material* mat, 141 const G4Material* mat,
142 const G4ParticleDefinition* part 142 const G4ParticleDefinition* part)
143 { 143 {
144 G4double res = stepLength; 144 G4double res = stepLength;
145 //G4cout << "## G4EnergyLossForExtrapolator: 145 //G4cout << "## G4EnergyLossForExtrapolator::TrueStepLength L= " << res
146 // << " " << part->GetParticleName() << 146 // << " " << part->GetParticleName() << G4endl;
147 if(SetupKinematics(part, mat, kinEnergy)) { 147 if(SetupKinematics(part, mat, kinEnergy)) {
148 if(part == electron || part == positron) { 148 if(part == electron || part == positron) {
149 const G4double x = stepLength* 149 const G4double x = stepLength*
150 ComputeValue(kinEnergy, GetPhysicsTable(fMsc 150 ComputeValue(kinEnergy, GetPhysicsTable(fMscElectron), mat->GetIndex());
151 //G4cout << " x= " << x << G4endl; 151 //G4cout << " x= " << x << G4endl;
152 if(x < 0.2) { res *= (1.0 + 0.5* 152 if(x < 0.2) { res *= (1.0 + 0.5*x + x*x/3.0); }
153 else if(x < 0.9999) { res = -G4Log(1.0 - 153 else if(x < 0.9999) { res = -G4Log(1.0 - x)*stepLength/x; }
154 else { res = ComputeRange(kinEnergy, par 154 else { res = ComputeRange(kinEnergy, part, mat); }
155 } else { 155 } else {
156 res = ComputeTrueStep(mat,part,kinEnergy 156 res = ComputeTrueStep(mat,part,kinEnergy,stepLength);
157 } 157 }
158 } 158 }
159 return res; 159 return res;
160 } 160 }
161 161
162 //....oooOO0OOooo........oooOO0OOooo........oo 162 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
163 163
164 G4bool 164 G4bool
165 G4EnergyLossForExtrapolator::SetupKinematics(c 165 G4EnergyLossForExtrapolator::SetupKinematics(const G4ParticleDefinition* part,
166 const G4Material* mat, 166 const G4Material* mat,
167 G4double kinEnergy) 167 G4double kinEnergy)
168 { 168 {
169 if(mat->GetNumberOfMaterials() != nmat) { In 169 if(mat->GetNumberOfMaterials() != nmat) { Initialisation(); }
170 if(nullptr == part || nullptr == mat || kinE 170 if(nullptr == part || nullptr == mat || kinEnergy < CLHEP::keV)
171 { return false; } 171 { return false; }
172 if(part != currentParticle) { 172 if(part != currentParticle) {
173 currentParticle = part; 173 currentParticle = part;
174 G4double q = part->GetPDGCharge()/eplus; 174 G4double q = part->GetPDGCharge()/eplus;
175 charge2 = q*q; 175 charge2 = q*q;
176 } 176 }
177 if(mat != currentMaterial) { 177 if(mat != currentMaterial) {
178 size_t i = mat->GetIndex(); 178 size_t i = mat->GetIndex();
179 if(i >= nmat) { 179 if(i >= nmat) {
180 G4cout << "### G4EnergyLossForExtrapolat 180 G4cout << "### G4EnergyLossForExtrapolator WARNING: material index i= "
181 << i << " above number of materials " < 181 << i << " above number of materials " << nmat << G4endl;
182 return false; 182 return false;
183 } else { 183 } else {
184 currentMaterial = mat; 184 currentMaterial = mat;
185 electronDensity = mat->GetElectronDensit 185 electronDensity = mat->GetElectronDensity();
186 radLength = mat->GetRadlen(); 186 radLength = mat->GetRadlen();
187 } 187 }
188 } 188 }
189 if(kinEnergy != kineticEnergy) { 189 if(kinEnergy != kineticEnergy) {
190 kineticEnergy = kinEnergy; 190 kineticEnergy = kinEnergy;
191 G4double mass = part->GetPDGMass(); 191 G4double mass = part->GetPDGMass();
192 G4double tau = kinEnergy/mass; 192 G4double tau = kinEnergy/mass;
193 193
194 gam = tau + 1.0; 194 gam = tau + 1.0;
195 bg2 = tau * (tau + 2.0); 195 bg2 = tau * (tau + 2.0);
196 beta2 = bg2/(gam*gam); 196 beta2 = bg2/(gam*gam);
197 tmax = kinEnergy; 197 tmax = kinEnergy;
198 if(part == electron) tmax *= 0.5; 198 if(part == electron) tmax *= 0.5;
199 else if(part != positron) { 199 else if(part != positron) {
200 G4double r = CLHEP::electron_mass_c2/mas 200 G4double r = CLHEP::electron_mass_c2/mass;
201 tmax = 2.0*bg2*CLHEP::electron_mass_c2/( 201 tmax = 2.0*bg2*CLHEP::electron_mass_c2/(1.0 + 2.0*gam*r + r*r);
202 } 202 }
203 tmax = std::min(tmax, maxEnergyTransfer); 203 tmax = std::min(tmax, maxEnergyTransfer);
204 } 204 }
205 return true; 205 return true;
206 } 206 }
207 207
208 //....oooOO0OOooo........oooOO0OOooo........oo 208 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
209 209
210 const G4ParticleDefinition* 210 const G4ParticleDefinition*
211 G4EnergyLossForExtrapolator::FindParticle(cons 211 G4EnergyLossForExtrapolator::FindParticle(const G4String& name)
212 { 212 {
213 currentParticle = G4ParticleTable::GetPartic 213 currentParticle = G4ParticleTable::GetParticleTable()->FindParticle(name);
214 if(nullptr == currentParticle) { 214 if(nullptr == currentParticle) {
215 G4cout << "### G4EnergyLossForExtrapolator 215 G4cout << "### G4EnergyLossForExtrapolator WARNING: "
216 << "FindParticle() fails to find " << nam 216 << "FindParticle() fails to find " << name << G4endl;
217 } 217 }
218 return currentParticle; 218 return currentParticle;
219 } 219 }
220 220
221 //....oooOO0OOooo........oooOO0OOooo........oo 221 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
222 222
223 G4double 223 G4double
224 G4EnergyLossForExtrapolator::ComputeDEDX(G4dou 224 G4EnergyLossForExtrapolator::ComputeDEDX(G4double ekin,
225 const G4ParticleDefinition* part, 225 const G4ParticleDefinition* part,
226 const 226 const G4Material* mat)
227 { 227 {
228 if(mat->GetNumberOfMaterials() != nmat) { In 228 if(mat->GetNumberOfMaterials() != nmat) { Initialisation(); }
229 G4double x = 0.0; 229 G4double x = 0.0;
230 if(part == electron) { 230 if(part == electron) {
231 x = ComputeValue(ekin, GetPhysicsTable(fDe 231 x = ComputeValue(ekin, GetPhysicsTable(fDedxElectron), mat->GetIndex());
232 } else if(part == positron) { 232 } else if(part == positron) {
233 x = ComputeValue(ekin, GetPhysicsTable(fDe 233 x = ComputeValue(ekin, GetPhysicsTable(fDedxPositron), mat->GetIndex());
234 } else if(part == muonPlus || part == muonMi 234 } else if(part == muonPlus || part == muonMinus) {
235 x = ComputeValue(ekin, GetPhysicsTable(fDe 235 x = ComputeValue(ekin, GetPhysicsTable(fDedxMuon), mat->GetIndex());
236 } else { 236 } else {
237 G4double e = ekin*CLHEP::proton_mass_c2/pa 237 G4double e = ekin*CLHEP::proton_mass_c2/part->GetPDGMass();
238 G4double q = part->GetPDGCharge()/CLHEP::e 238 G4double q = part->GetPDGCharge()/CLHEP::eplus;
239 x = ComputeValue(e, GetPhysicsTable(fDedxP 239 x = ComputeValue(e, GetPhysicsTable(fDedxProton), mat->GetIndex())*q*q;
240 } 240 }
241 return x; 241 return x;
242 } 242 }
243 243
244 //....oooOO0OOooo........oooOO0OOooo........oo 244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
245 245
246 G4double 246 G4double
247 G4EnergyLossForExtrapolator::ComputeRange(G4do 247 G4EnergyLossForExtrapolator::ComputeRange(G4double ekin,
248 const G4ParticleDefinition* part, 248 const G4ParticleDefinition* part,
249 const G4Material* mat) 249 const G4Material* mat)
250 { 250 {
251 if(mat->GetNumberOfMaterials() != nmat) { In 251 if(mat->GetNumberOfMaterials() != nmat) { Initialisation(); }
252 G4double x = 0.0; 252 G4double x = 0.0;
253 if(part == electron) { 253 if(part == electron) {
254 x = ComputeValue(ekin, GetPhysicsTable(fRa 254 x = ComputeValue(ekin, GetPhysicsTable(fRangeElectron), mat->GetIndex());
255 } else if(part == positron) { 255 } else if(part == positron) {
256 x = ComputeValue(ekin, GetPhysicsTable(fRa 256 x = ComputeValue(ekin, GetPhysicsTable(fRangePositron), mat->GetIndex());
257 } else if(part == muonPlus || part == muonMi 257 } else if(part == muonPlus || part == muonMinus) {
258 x = ComputeValue(ekin, GetPhysicsTable(fRa 258 x = ComputeValue(ekin, GetPhysicsTable(fRangeMuon), mat->GetIndex());
259 } else { 259 } else {
260 G4double massratio = CLHEP::proton_mass_c2 260 G4double massratio = CLHEP::proton_mass_c2/part->GetPDGMass();
261 G4double e = ekin*massratio; 261 G4double e = ekin*massratio;
262 G4double q = part->GetPDGCharge()/CLHEP::e 262 G4double q = part->GetPDGCharge()/CLHEP::eplus;
263 x = ComputeValue(e, GetPhysicsTable(fRange 263 x = ComputeValue(e, GetPhysicsTable(fRangeProton), mat->GetIndex())
264 /(q*q*massratio); 264 /(q*q*massratio);
265 } 265 }
266 return x; 266 return x;
267 } 267 }
268 268
269 //....oooOO0OOooo........oooOO0OOooo........oo 269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
270 270
271 G4double 271 G4double
272 G4EnergyLossForExtrapolator::ComputeEnergy(G4d 272 G4EnergyLossForExtrapolator::ComputeEnergy(G4double range,
273 const G4ParticleDefinition* part, 273 const G4ParticleDefinition* part,
274 const G4Material* mat) 274 const G4Material* mat)
275 { 275 {
276 if(mat->GetNumberOfMaterials() != nmat) { In 276 if(mat->GetNumberOfMaterials() != nmat) { Initialisation(); }
277 G4double x = 0.0; 277 G4double x = 0.0;
278 if(part == electron) { 278 if(part == electron) {
279 x = ComputeValue(range,GetPhysicsTable(fIn 279 x = ComputeValue(range,GetPhysicsTable(fInvRangeElectron),mat->GetIndex());
280 } else if(part == positron) { 280 } else if(part == positron) {
281 x = ComputeValue(range,GetPhysicsTable(fIn 281 x = ComputeValue(range,GetPhysicsTable(fInvRangePositron),mat->GetIndex());
282 } else if(part == muonPlus || part == muonMi 282 } else if(part == muonPlus || part == muonMinus) {
283 x = ComputeValue(range, GetPhysicsTable(fI 283 x = ComputeValue(range, GetPhysicsTable(fInvRangeMuon), mat->GetIndex());
284 } else { 284 } else {
285 G4double massratio = CLHEP::proton_mass_c2 285 G4double massratio = CLHEP::proton_mass_c2/part->GetPDGMass();
286 G4double q = part->GetPDGCharge()/CLHEP::e 286 G4double q = part->GetPDGCharge()/CLHEP::eplus;
287 G4double r = range*massratio*q*q; 287 G4double r = range*massratio*q*q;
288 x = ComputeValue(r, GetPhysicsTable(fInvRa 288 x = ComputeValue(r, GetPhysicsTable(fInvRangeProton), mat->GetIndex())/massratio;
289 } 289 }
290 return x; 290 return x;
291 } 291 }
292 292
293 //....oooOO0OOooo........oooOO0OOooo........oo 293 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
294 294
295 G4double 295 G4double
296 G4EnergyLossForExtrapolator::EnergyDispersion( 296 G4EnergyLossForExtrapolator::EnergyDispersion(G4double kinEnergy,
297 G4double stepLength, 297 G4double stepLength,
298 const G4Material* mat, 298 const G4Material* mat,
299 const G4ParticleDefinition* pa 299 const G4ParticleDefinition* part)
300 { 300 {
301 G4double sig2 = 0.0; 301 G4double sig2 = 0.0;
302 if(SetupKinematics(part, mat, kinEnergy)) { 302 if(SetupKinematics(part, mat, kinEnergy)) {
303 G4double step = ComputeTrueStep(mat,part,k 303 G4double step = ComputeTrueStep(mat,part,kinEnergy,stepLength);
304 sig2 = (1.0/beta2 - 0.5) 304 sig2 = (1.0/beta2 - 0.5)
305 *CLHEP::twopi_mc2_rcl2*tmax*step*electro 305 *CLHEP::twopi_mc2_rcl2*tmax*step*electronDensity*charge2;
306 } 306 }
307 return sig2; 307 return sig2;
308 } 308 }
309 309
310 //....oooOO0OOooo........oooOO0OOooo........oo 310 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
311 311
312 G4double G4EnergyLossForExtrapolator::AverageS 312 G4double G4EnergyLossForExtrapolator::AverageScatteringAngle(
313 G4double kinEnergy, 313 G4double kinEnergy,
314 G4double stepLength, 314 G4double stepLength,
315 const G4Material* mat, 315 const G4Material* mat,
316 const G4ParticleDefinition* part) 316 const G4ParticleDefinition* part)
317 { 317 {
318 G4double theta = 0.0; 318 G4double theta = 0.0;
319 if(SetupKinematics(part, mat, kinEnergy)) { 319 if(SetupKinematics(part, mat, kinEnergy)) {
320 G4double t = stepLength/radLength; 320 G4double t = stepLength/radLength;
321 G4double y = std::max(0.001, t); 321 G4double y = std::max(0.001, t);
322 theta = 19.23*CLHEP::MeV*std::sqrt(charge2 322 theta = 19.23*CLHEP::MeV*std::sqrt(charge2*t)*(1.0 + 0.038*G4Log(y))
323 /(beta2*gam*part->GetPDGMass()); 323 /(beta2*gam*part->GetPDGMass());
324 } 324 }
325 return theta; 325 return theta;
326 } 326 }
327 327
328 //....oooOO0OOooo........oooOO0OOooo........oo 328 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
329 329
330 void G4EnergyLossForExtrapolator::Initialisati 330 void G4EnergyLossForExtrapolator::Initialisation()
331 { 331 {
332 if(verbose>0) { 332 if(verbose>0) {
333 G4cout << "### G4EnergyLossForExtrapolator 333 G4cout << "### G4EnergyLossForExtrapolator::Initialisation tables= "
334 << tables << G4endl; 334 << tables << G4endl;
335 } 335 }
336 electron = G4Electron::Electron(); 336 electron = G4Electron::Electron();
337 positron = G4Positron::Positron(); 337 positron = G4Positron::Positron();
338 proton = G4Proton::Proton(); 338 proton = G4Proton::Proton();
339 muonPlus = G4MuonPlus::MuonPlus(); 339 muonPlus = G4MuonPlus::MuonPlus();
340 muonMinus= G4MuonMinus::MuonMinus(); 340 muonMinus= G4MuonMinus::MuonMinus();
341 341
342 // initialisation for the 1st run 342 // initialisation for the 1st run
343 if(nullptr == tables) { 343 if(nullptr == tables) {
344 #ifdef G4MULTITHREADED 344 #ifdef G4MULTITHREADED
345 G4MUTEXLOCK(&extrMutex); 345 G4MUTEXLOCK(&extrMutex);
346 if(nullptr == tables) { 346 if(nullptr == tables) {
347 #endif 347 #endif
348 isMaster = true; 348 isMaster = true;
349 tables = new G4TablesForExtrapolator(ver 349 tables = new G4TablesForExtrapolator(verbose, nbins, emin, emax);
350 tables->Initialisation(); 350 tables->Initialisation();
351 nmat = G4Material::GetNumberOfMaterials( 351 nmat = G4Material::GetNumberOfMaterials();
352 if(verbose > 0) { 352 if(verbose > 0) {
353 G4cout << "### G4EnergyLossForExtrapol 353 G4cout << "### G4EnergyLossForExtrapolator::BuildTables for "
354 << nmat << " materials Nbins= " 354 << nmat << " materials Nbins= "
355 << nbins << " Emin(MeV)= " << e 355 << nbins << " Emin(MeV)= " << emin << " Emax(MeV)= " << emax
356 << G4endl; 356 << G4endl;
357 } 357 }
358 #ifdef G4MULTITHREADED 358 #ifdef G4MULTITHREADED
359 } 359 }
360 G4MUTEXUNLOCK(&extrMutex); 360 G4MUTEXUNLOCK(&extrMutex);
361 #endif 361 #endif
362 } 362 }
363 363
364 // initialisation for the next run 364 // initialisation for the next run
365 if(isMaster && G4Material::GetNumberOfMateri 365 if(isMaster && G4Material::GetNumberOfMaterials() != nmat) {
366 #ifdef G4MULTITHREADED 366 #ifdef G4MULTITHREADED
367 G4MUTEXLOCK(&extrMutex); 367 G4MUTEXLOCK(&extrMutex);
368 #endif 368 #endif
369 tables->Initialisation(); 369 tables->Initialisation();
370 #ifdef G4MULTITHREADED 370 #ifdef G4MULTITHREADED
371 G4MUTEXUNLOCK(&extrMutex); 371 G4MUTEXUNLOCK(&extrMutex);
372 #endif 372 #endif
373 } 373 }
374 nmat = G4Material::GetNumberOfMaterials(); 374 nmat = G4Material::GetNumberOfMaterials();
375 } 375 }
376 376
377 //....oooOO0OOooo........oooOO0OOooo........oo 377 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
378 378