| 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::extrapolatorMutex = 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 currentParticle = nullptr;
74 emax = 100.*CLHEP::TeV; << 74 currentMaterial = nullptr;
>> 75
>> 76 linLossLimit = 0.01;
>> 77 emin = 1.*MeV;
>> 78 emax = 10.*TeV;
>> 79 nbins = 70;
>> 80
>> 81 nmat = index = 0;
>> 82
>> 83 mass = charge2 = electronDensity = radLength = bg2 = beta2
>> 84 = kineticEnergy = tmax = 0.0;
>> 85 gam = 1.0;
>> 86
>> 87 idxDedxElectron = idxDedxPositron = idxDedxMuon = idxDedxProton
>> 88 = idxRangeElectron = idxRangePositron = idxRangeMuon = idxRangeProton
>> 89 = idxInvRangeElectron = idxInvRangePositron = idxInvRangeMuon
>> 90 = idxInvRangeProton = idxMscElectron = 0;
>> 91
>> 92 electron = positron = proton = muonPlus = muonMinus = nullptr;
75 } 93 }
76 94
77 //....oooOO0OOooo........oooOO0OOooo........oo 95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78 96
79 G4EnergyLossForExtrapolator::~G4EnergyLossForE 97 G4EnergyLossForExtrapolator::~G4EnergyLossForExtrapolator()
80 { 98 {
81 if(isMaster) { << 99 if(tables) {
82 delete tables; << 100 #ifdef G4MULTITHREADED
83 tables = nullptr; << 101 G4MUTEXLOCK(&extrapolatorMutex);
>> 102 if (tables) {
>> 103 #endif
>> 104 delete tables;
>> 105 tables = nullptr;
>> 106 #ifdef G4MULTITHREADED
>> 107 }
>> 108 G4MUTEXUNLOCK(&extrapolatorMutex);
>> 109 #endif
84 } 110 }
85 } 111 }
86 112
87 //....oooOO0OOooo........oooOO0OOooo........oo 113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88 114
89 G4double 115 G4double
90 G4EnergyLossForExtrapolator::EnergyAfterStep(G 116 G4EnergyLossForExtrapolator::EnergyAfterStep(G4double kinEnergy,
91 G4double stepLength, 117 G4double stepLength,
92 const G4Material* mat, 118 const G4Material* mat,
93 const G4ParticleDefinition* par 119 const G4ParticleDefinition* part)
94 { 120 {
>> 121 if(0 == nmat) { Initialisation(); }
95 G4double kinEnergyFinal = kinEnergy; 122 G4double kinEnergyFinal = kinEnergy;
96 if(SetupKinematics(part, mat, kinEnergy)) { 123 if(SetupKinematics(part, mat, kinEnergy)) {
97 G4double step = TrueStepLength(kinEnergy,s 124 G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
98 G4double r = ComputeRange(kinEnergy,part, << 125 G4double r = ComputeRange(kinEnergy,part);
99 if(r <= step) { 126 if(r <= step) {
100 kinEnergyFinal = 0.0; 127 kinEnergyFinal = 0.0;
101 } else if(step < linLossLimit*r) { 128 } else if(step < linLossLimit*r) {
102 kinEnergyFinal -= step*ComputeDEDX(kinEn << 129 kinEnergyFinal -= step*ComputeDEDX(kinEnergy,part);
103 } else { 130 } else {
104 G4double r1 = r - step; 131 G4double r1 = r - step;
105 kinEnergyFinal = ComputeEnergy(r1,part,m << 132 kinEnergyFinal = ComputeEnergy(r1,part);
106 } 133 }
107 } 134 }
108 return kinEnergyFinal; 135 return kinEnergyFinal;
109 } 136 }
110 137
111 //....oooOO0OOooo........oooOO0OOooo........oo 138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
112 139
113 G4double 140 G4double
114 G4EnergyLossForExtrapolator::EnergyBeforeStep( 141 G4EnergyLossForExtrapolator::EnergyBeforeStep(G4double kinEnergy,
115 G4double stepLength, 142 G4double stepLength,
116 const G4Material* mat, 143 const G4Material* mat,
117 const G4ParticleDefinition* pa 144 const G4ParticleDefinition* part)
118 { 145 {
119 //G4cout << "G4EnergyLossForExtrapolator::En 146 //G4cout << "G4EnergyLossForExtrapolator::EnergyBeforeStep" << G4endl;
>> 147 if(0 == nmat) { Initialisation(); }
120 G4double kinEnergyFinal = kinEnergy; 148 G4double kinEnergyFinal = kinEnergy;
121 149
122 if(SetupKinematics(part, mat, kinEnergy)) { 150 if(SetupKinematics(part, mat, kinEnergy)) {
123 G4double step = TrueStepLength(kinEnergy,s 151 G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
124 G4double r = ComputeRange(kinEnergy,part,m << 152 G4double r = ComputeRange(kinEnergy,part);
125 153
126 if(step < linLossLimit*r) { 154 if(step < linLossLimit*r) {
127 kinEnergyFinal += step*ComputeDEDX(kinEn << 155 kinEnergyFinal += step*ComputeDEDX(kinEnergy,part);
128 } else { 156 } else {
129 G4double r1 = r + step; 157 G4double r1 = r + step;
130 kinEnergyFinal = ComputeEnergy(r1,part,m << 158 kinEnergyFinal = ComputeEnergy(r1,part);
131 } 159 }
132 } 160 }
133 return kinEnergyFinal; 161 return kinEnergyFinal;
134 } 162 }
135 163
136 //....oooOO0OOooo........oooOO0OOooo........oo 164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
137 165
138 G4double 166 G4double
139 G4EnergyLossForExtrapolator::TrueStepLength(G4 167 G4EnergyLossForExtrapolator::TrueStepLength(G4double kinEnergy,
140 G4double stepLength, 168 G4double stepLength,
141 const G4Material* mat, 169 const G4Material* mat,
142 const G4ParticleDefinition* part 170 const G4ParticleDefinition* part)
143 { 171 {
144 G4double res = stepLength; 172 G4double res = stepLength;
>> 173 if(0 == nmat) { Initialisation(); }
145 //G4cout << "## G4EnergyLossForExtrapolator: 174 //G4cout << "## G4EnergyLossForExtrapolator::TrueStepLength L= " << res
146 // << " " << part->GetParticleName() << 175 // << " " << part->GetParticleName() << G4endl;
147 if(SetupKinematics(part, mat, kinEnergy)) { 176 if(SetupKinematics(part, mat, kinEnergy)) {
148 if(part == electron || part == positron) { 177 if(part == electron || part == positron) {
149 const G4double x = stepLength* << 178 G4double x = stepLength*
150 ComputeValue(kinEnergy, GetPhysicsTable(fMsc << 179 ComputeValue(kinEnergy, GetPhysicsTable(fMscElectron), idxMscElectron);
151 //G4cout << " x= " << x << G4endl; 180 //G4cout << " x= " << x << G4endl;
152 if(x < 0.2) { res *= (1.0 + 0.5* 181 if(x < 0.2) { res *= (1.0 + 0.5*x + x*x/3.0); }
153 else if(x < 0.9999) { res = -G4Log(1.0 - 182 else if(x < 0.9999) { res = -G4Log(1.0 - x)*stepLength/x; }
154 else { res = ComputeRange(kinEnergy, par << 183 else { res = ComputeRange(kinEnergy,part); }
155 } else { 184 } else {
156 res = ComputeTrueStep(mat,part,kinEnergy 185 res = ComputeTrueStep(mat,part,kinEnergy,stepLength);
157 } 186 }
158 } 187 }
159 return res; 188 return res;
160 } 189 }
161 190
162 //....oooOO0OOooo........oooOO0OOooo........oo 191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
163 192
164 G4bool 193 G4bool
165 G4EnergyLossForExtrapolator::SetupKinematics(c 194 G4EnergyLossForExtrapolator::SetupKinematics(const G4ParticleDefinition* part,
166 const G4Material* mat, 195 const G4Material* mat,
167 G4double kinEnergy) 196 G4double kinEnergy)
168 { 197 {
169 if(mat->GetNumberOfMaterials() != nmat) { In << 198 if(0 == nmat) { Initialisation(); }
170 if(nullptr == part || nullptr == mat || kinE << 199 if(!part || !mat || kinEnergy < keV) { return false; }
171 { return false; } << 200 G4bool flag = false;
172 if(part != currentParticle) { 201 if(part != currentParticle) {
>> 202 flag = true;
173 currentParticle = part; 203 currentParticle = part;
>> 204 mass = part->GetPDGMass();
174 G4double q = part->GetPDGCharge()/eplus; 205 G4double q = part->GetPDGCharge()/eplus;
175 charge2 = q*q; 206 charge2 = q*q;
176 } 207 }
177 if(mat != currentMaterial) { 208 if(mat != currentMaterial) {
178 size_t i = mat->GetIndex(); << 209 G4int i = mat->GetIndex();
179 if(i >= nmat) { 210 if(i >= nmat) {
180 G4cout << "### G4EnergyLossForExtrapolat << 211 G4cout << "### G4EnergyLossForExtrapolator WARNING:index i= "
181 << i << " above number of materials " < << 212 << i << " is out of table - NO extrapolation" << G4endl;
182 return false; <<
183 } else { 213 } else {
>> 214 flag = true;
184 currentMaterial = mat; 215 currentMaterial = mat;
185 electronDensity = mat->GetElectronDensit 216 electronDensity = mat->GetElectronDensity();
186 radLength = mat->GetRadlen(); 217 radLength = mat->GetRadlen();
>> 218 index = i;
187 } 219 }
188 } 220 }
189 if(kinEnergy != kineticEnergy) { << 221 if(flag || kinEnergy != kineticEnergy) {
190 kineticEnergy = kinEnergy; 222 kineticEnergy = kinEnergy;
191 G4double mass = part->GetPDGMass(); <<
192 G4double tau = kinEnergy/mass; 223 G4double tau = kinEnergy/mass;
193 224
194 gam = tau + 1.0; 225 gam = tau + 1.0;
195 bg2 = tau * (tau + 2.0); 226 bg2 = tau * (tau + 2.0);
196 beta2 = bg2/(gam*gam); 227 beta2 = bg2/(gam*gam);
197 tmax = kinEnergy; 228 tmax = kinEnergy;
198 if(part == electron) tmax *= 0.5; 229 if(part == electron) tmax *= 0.5;
199 else if(part != positron) { 230 else if(part != positron) {
200 G4double r = CLHEP::electron_mass_c2/mas << 231 G4double r = electron_mass_c2/mass;
201 tmax = 2.0*bg2*CLHEP::electron_mass_c2/( << 232 tmax = 2.0*bg2*electron_mass_c2/(1.0 + 2.0*gam*r + r*r);
202 } 233 }
203 tmax = std::min(tmax, maxEnergyTransfer); << 234 if(tmax > maxEnergyTransfer) { tmax = maxEnergyTransfer; }
204 } 235 }
205 return true; 236 return true;
206 } 237 }
207 238
208 //....oooOO0OOooo........oooOO0OOooo........oo 239 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
209 240
210 const G4ParticleDefinition* 241 const G4ParticleDefinition*
211 G4EnergyLossForExtrapolator::FindParticle(cons 242 G4EnergyLossForExtrapolator::FindParticle(const G4String& name)
212 { 243 {
213 currentParticle = G4ParticleTable::GetPartic << 244 const G4ParticleDefinition* p = nullptr;
214 if(nullptr == currentParticle) { << 245 if(name != currentParticleName) {
215 G4cout << "### G4EnergyLossForExtrapolator << 246 p = G4ParticleTable::GetParticleTable()->FindParticle(name);
216 << "FindParticle() fails to find " << nam << 247 if(!p) {
>> 248 G4cout << "### G4EnergyLossForExtrapolator WARNING: "
>> 249 << "FindParticle() fails to find "
>> 250 << name << G4endl;
>> 251 }
>> 252 } else {
>> 253 p = currentParticle;
217 } 254 }
218 return currentParticle; << 255 return p;
219 } 256 }
220 257
221 //....oooOO0OOooo........oooOO0OOooo........oo 258 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
222 259
223 G4double 260 G4double
224 G4EnergyLossForExtrapolator::ComputeDEDX(G4dou 261 G4EnergyLossForExtrapolator::ComputeDEDX(G4double ekin,
225 const G4ParticleDefinition* part, << 262 const G4ParticleDefinition* part)
226 const <<
227 { 263 {
228 if(mat->GetNumberOfMaterials() != nmat) { In <<
229 G4double x = 0.0; 264 G4double x = 0.0;
230 if(part == electron) { 265 if(part == electron) {
231 x = ComputeValue(ekin, GetPhysicsTable(fDe << 266 x = ComputeValue(ekin, GetPhysicsTable(fDedxElectron), idxDedxElectron);
232 } else if(part == positron) { 267 } else if(part == positron) {
233 x = ComputeValue(ekin, GetPhysicsTable(fDe << 268 x = ComputeValue(ekin, GetPhysicsTable(fDedxPositron), idxDedxPositron);
234 } else if(part == muonPlus || part == muonMi 269 } else if(part == muonPlus || part == muonMinus) {
235 x = ComputeValue(ekin, GetPhysicsTable(fDe << 270 x = ComputeValue(ekin, GetPhysicsTable(fDedxMuon), idxDedxMuon);
236 } else { 271 } else {
237 G4double e = ekin*CLHEP::proton_mass_c2/pa << 272 G4double e = ekin*proton_mass_c2/mass;
238 G4double q = part->GetPDGCharge()/CLHEP::e << 273 x = ComputeValue(e, GetPhysicsTable(fDedxProton), idxDedxProton)*charge2;
239 x = ComputeValue(e, GetPhysicsTable(fDedxP <<
240 } 274 }
241 return x; 275 return x;
242 } 276 }
243 277
244 //....oooOO0OOooo........oooOO0OOooo........oo 278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
245 279
246 G4double 280 G4double
247 G4EnergyLossForExtrapolator::ComputeRange(G4do 281 G4EnergyLossForExtrapolator::ComputeRange(G4double ekin,
248 const G4ParticleDefinition* part, << 282 const G4ParticleDefinition* part)
249 const G4Material* mat) <<
250 { 283 {
251 if(mat->GetNumberOfMaterials() != nmat) { In <<
252 G4double x = 0.0; 284 G4double x = 0.0;
253 if(part == electron) { 285 if(part == electron) {
254 x = ComputeValue(ekin, GetPhysicsTable(fRa << 286 x = ComputeValue(ekin, GetPhysicsTable(fRangeElectron), idxRangeElectron);
255 } else if(part == positron) { 287 } else if(part == positron) {
256 x = ComputeValue(ekin, GetPhysicsTable(fRa << 288 x = ComputeValue(ekin, GetPhysicsTable(fRangePositron), idxRangePositron);
257 } else if(part == muonPlus || part == muonMi 289 } else if(part == muonPlus || part == muonMinus) {
258 x = ComputeValue(ekin, GetPhysicsTable(fRa << 290 x = ComputeValue(ekin, GetPhysicsTable(fRangeMuon), idxRangeMuon);
259 } else { 291 } else {
260 G4double massratio = CLHEP::proton_mass_c2 << 292 G4double massratio = proton_mass_c2/mass;
261 G4double e = ekin*massratio; 293 G4double e = ekin*massratio;
262 G4double q = part->GetPDGCharge()/CLHEP::e << 294 x = ComputeValue(e, GetPhysicsTable(fRangeProton), idxRangeProton)
263 x = ComputeValue(e, GetPhysicsTable(fRange << 295 /(charge2*massratio);
264 /(q*q*massratio); <<
265 } 296 }
266 return x; 297 return x;
267 } 298 }
268 299
269 //....oooOO0OOooo........oooOO0OOooo........oo 300 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
270 301
271 G4double 302 G4double
272 G4EnergyLossForExtrapolator::ComputeEnergy(G4d 303 G4EnergyLossForExtrapolator::ComputeEnergy(G4double range,
273 const G4ParticleDefinition* part, << 304 const G4ParticleDefinition* part)
274 const G4Material* mat) <<
275 { 305 {
276 if(mat->GetNumberOfMaterials() != nmat) { In <<
277 G4double x = 0.0; 306 G4double x = 0.0;
278 if(part == electron) { 307 if(part == electron) {
279 x = ComputeValue(range,GetPhysicsTable(fIn << 308 x = ComputeValue(range, GetPhysicsTable(fInvRangeElectron),
>> 309 idxInvRangeElectron);
280 } else if(part == positron) { 310 } else if(part == positron) {
281 x = ComputeValue(range,GetPhysicsTable(fIn << 311 x = ComputeValue(range, GetPhysicsTable(fInvRangePositron),
>> 312 idxInvRangePositron);
282 } else if(part == muonPlus || part == muonMi 313 } else if(part == muonPlus || part == muonMinus) {
283 x = ComputeValue(range, GetPhysicsTable(fI << 314 x = ComputeValue(range, GetPhysicsTable(fInvRangeMuon), idxInvRangeMuon);
284 } else { 315 } else {
285 G4double massratio = CLHEP::proton_mass_c2 << 316 G4double massratio = proton_mass_c2/mass;
286 G4double q = part->GetPDGCharge()/CLHEP::e << 317 G4double r = range*massratio*charge2;
287 G4double r = range*massratio*q*q; << 318 x = ComputeValue(r, GetPhysicsTable(fInvRangeProton),
288 x = ComputeValue(r, GetPhysicsTable(fInvRa << 319 idxInvRangeProton)/massratio;
289 } 320 }
290 return x; 321 return x;
291 } 322 }
292 323
293 //....oooOO0OOooo........oooOO0OOooo........oo 324 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
294 <<
295 G4double <<
296 G4EnergyLossForExtrapolator::EnergyDispersion( <<
297 G4double stepLength, <<
298 const G4Material* mat, <<
299 const G4ParticleDefinition* pa <<
300 { <<
301 G4double sig2 = 0.0; <<
302 if(SetupKinematics(part, mat, kinEnergy)) { <<
303 G4double step = ComputeTrueStep(mat,part,k <<
304 sig2 = (1.0/beta2 - 0.5) <<
305 *CLHEP::twopi_mc2_rcl2*tmax*step*electro <<
306 } <<
307 return sig2; <<
308 } <<
309 <<
310 //....oooOO0OOooo........oooOO0OOooo........oo <<
311 <<
312 G4double G4EnergyLossForExtrapolator::AverageS <<
313 G4double kinEnergy, <<
314 G4double stepLength, <<
315 const G4Material* mat, <<
316 const G4ParticleDefinition* part) <<
317 { <<
318 G4double theta = 0.0; <<
319 if(SetupKinematics(part, mat, kinEnergy)) { <<
320 G4double t = stepLength/radLength; <<
321 G4double y = std::max(0.001, t); <<
322 theta = 19.23*CLHEP::MeV*std::sqrt(charge2 <<
323 /(beta2*gam*part->GetPDGMass()); <<
324 } <<
325 return theta; <<
326 } <<
327 <<
328 //....oooOO0OOooo........oooOO0OOooo........oo <<
329 325
330 void G4EnergyLossForExtrapolator::Initialisati 326 void G4EnergyLossForExtrapolator::Initialisation()
331 { 327 {
332 if(verbose>0) { << 328 if(verbose>1) {
333 G4cout << "### G4EnergyLossForExtrapolator << 329 G4cout << "### G4EnergyLossForExtrapolator::Initialisation" << G4endl;
334 << tables << G4endl; <<
335 } 330 }
>> 331 currentParticle = nullptr;
>> 332 currentMaterial = nullptr;
>> 333 kineticEnergy = 0.0;
>> 334
336 electron = G4Electron::Electron(); 335 electron = G4Electron::Electron();
337 positron = G4Positron::Positron(); 336 positron = G4Positron::Positron();
338 proton = G4Proton::Proton(); 337 proton = G4Proton::Proton();
339 muonPlus = G4MuonPlus::MuonPlus(); 338 muonPlus = G4MuonPlus::MuonPlus();
340 muonMinus= G4MuonMinus::MuonMinus(); 339 muonMinus= G4MuonMinus::MuonMinus();
341 340
342 // initialisation for the 1st run << 341 nmat = G4Material::GetNumberOfMaterials();
343 if(nullptr == tables) { << 342 currentParticleName = "";
>> 343 BuildTables();
>> 344 }
>> 345
>> 346 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 347
>> 348 void G4EnergyLossForExtrapolator::BuildTables()
>> 349 {
>> 350 if(!tables) {
344 #ifdef G4MULTITHREADED 351 #ifdef G4MULTITHREADED
345 G4MUTEXLOCK(&extrMutex); << 352 G4MUTEXLOCK(&extrapolatorMutex);
346 if(nullptr == tables) { << 353 if (!tables) {
347 #endif 354 #endif
348 isMaster = true; <<
349 tables = new G4TablesForExtrapolator(ver <<
350 tables->Initialisation(); <<
351 nmat = G4Material::GetNumberOfMaterials( <<
352 if(verbose > 0) { 355 if(verbose > 0) {
353 G4cout << "### G4EnergyLossForExtrapol << 356 G4cout << "### G4EnergyLossForExtrapolator::BuildTables for "
354 << nmat << " materials Nbins= " << 357 << nmat << " materials Nbins= " << nbins
355 << nbins << " Emin(MeV)= " << e << 358 << " Emin(MeV)= " << emin << " Emax(MeV)= " << emax << G4endl;
356 << G4endl; <<
357 } 359 }
>> 360 tables = new G4TablesForExtrapolator(verbose, nbins, emin, emax);
358 #ifdef G4MULTITHREADED 361 #ifdef G4MULTITHREADED
359 } 362 }
360 G4MUTEXUNLOCK(&extrMutex); << 363 G4MUTEXUNLOCK(&extrapolatorMutex);
361 #endif <<
362 } <<
363 <<
364 // initialisation for the next run <<
365 if(isMaster && G4Material::GetNumberOfMateri <<
366 #ifdef G4MULTITHREADED <<
367 G4MUTEXLOCK(&extrMutex); <<
368 #endif <<
369 tables->Initialisation(); <<
370 #ifdef G4MULTITHREADED <<
371 G4MUTEXUNLOCK(&extrMutex); <<
372 #endif 364 #endif
373 } 365 }
374 nmat = G4Material::GetNumberOfMaterials(); <<
375 } 366 }
376 367
377 //....oooOO0OOooo........oooOO0OOooo........oo 368 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 369
378 370