| 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: G4EnergyLossForExtrapolator.cc,v 1.19 2009/07/09 17:04:55 vnivanch Exp $
>> 27 // GEANT4 tag $Name: geant4-09-03-patch-01 $
26 // 28 //
27 //-------------------------------------------- 29 //---------------------------------------------------------------------------
28 // 30 //
29 // ClassName: G4EnergyLossForExtrapolator 31 // ClassName: G4EnergyLossForExtrapolator
30 // 32 //
31 // Description: This class provide calculatio 33 // Description: This class provide calculation of energy loss, fluctuation,
32 // and msc angle 34 // and msc angle
33 // 35 //
34 // Author: 09.12.04 V.Ivanchenko 36 // Author: 09.12.04 V.Ivanchenko
35 // 37 //
36 // Modification: 38 // Modification:
37 // 08-04-05 Rename Propogator -> Extrapolator 39 // 08-04-05 Rename Propogator -> Extrapolator (V.Ivanchenko)
38 // 16-03-06 Add muon tables and fix bug in uni 40 // 16-03-06 Add muon tables and fix bug in units (V.Ivanchenko)
39 // 21-03-06 Add verbosity defined in the const 41 // 21-03-06 Add verbosity defined in the constructor and Initialisation
40 // start only when first public metho 42 // start only when first public method is called (V.Ivanchenko)
41 // 03-05-06 Remove unused pointer G4Material* 43 // 03-05-06 Remove unused pointer G4Material* from number of methods (VI)
42 // 12-05-06 SEt linLossLimit=0.001 (VI) 44 // 12-05-06 SEt linLossLimit=0.001 (VI)
43 // 45 //
44 //-------------------------------------------- 46 //----------------------------------------------------------------------------
45 // 47 //
46 48
47 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48 50
49 #include "G4EnergyLossForExtrapolator.hh" 51 #include "G4EnergyLossForExtrapolator.hh"
50 #include "G4PhysicalConstants.hh" << 52 #include "G4PhysicsLogVector.hh"
51 #include "G4SystemOfUnits.hh" <<
52 #include "G4ParticleDefinition.hh" 53 #include "G4ParticleDefinition.hh"
53 #include "G4Material.hh" 54 #include "G4Material.hh"
54 #include "G4MaterialCutsCouple.hh" 55 #include "G4MaterialCutsCouple.hh"
55 #include "G4Electron.hh" 56 #include "G4Electron.hh"
56 #include "G4Positron.hh" 57 #include "G4Positron.hh"
57 #include "G4Proton.hh" 58 #include "G4Proton.hh"
58 #include "G4MuonPlus.hh" 59 #include "G4MuonPlus.hh"
59 #include "G4MuonMinus.hh" 60 #include "G4MuonMinus.hh"
60 #include "G4ParticleTable.hh" 61 #include "G4ParticleTable.hh"
>> 62 #include "G4LossTableBuilder.hh"
>> 63 #include "G4MollerBhabhaModel.hh"
>> 64 #include "G4BetheBlochModel.hh"
>> 65 #include "G4eBremsstrahlungModel.hh"
>> 66 #include "G4MuPairProductionModel.hh"
>> 67 #include "G4MuBremsstrahlungModel.hh"
>> 68 #include "G4ProductionCuts.hh"
>> 69 #include "G4LossTableManager.hh"
>> 70 #include "G4WentzelVIModel.hh"
61 71
62 //....oooOO0OOooo........oooOO0OOooo........oo 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63 73
64 #ifdef G4MULTITHREADED <<
65 G4Mutex G4EnergyLossForExtrapolator::extrMutex <<
66 #endif <<
67 <<
68 G4TablesForExtrapolator* G4EnergyLossForExtrap <<
69 <<
70 G4EnergyLossForExtrapolator::G4EnergyLossForEx 74 G4EnergyLossForExtrapolator::G4EnergyLossForExtrapolator(G4int verb)
71 : maxEnergyTransfer(DBL_MAX), verbose(verb) << 75 :maxEnergyTransfer(DBL_MAX),verbose(verb),isInitialised(false)
72 { << 76 {}
73 emin = 1.*CLHEP::MeV; <<
74 emax = 100.*CLHEP::TeV; <<
75 } <<
76 77
77 //....oooOO0OOooo........oooOO0OOooo........oo 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78 79
79 G4EnergyLossForExtrapolator::~G4EnergyLossForE << 80 G4EnergyLossForExtrapolator:: ~G4EnergyLossForExtrapolator()
80 { 81 {
81 if(isMaster) { << 82 for(G4int i=0; i<nmat; i++) {delete couples[i];}
82 delete tables; << 83 delete dedxElectron;
83 tables = nullptr; << 84 delete dedxPositron;
84 } << 85 delete dedxProton;
>> 86 delete rangeElectron;
>> 87 delete rangePositron;
>> 88 delete rangeProton;
>> 89 delete invRangeElectron;
>> 90 delete invRangePositron;
>> 91 delete invRangeProton;
>> 92 delete mscElectron;
>> 93 delete cuts;
85 } 94 }
86 95
87 //....oooOO0OOooo........oooOO0OOooo........oo 96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88 97
89 G4double << 98 G4double G4EnergyLossForExtrapolator::EnergyAfterStep(G4double kinEnergy,
90 G4EnergyLossForExtrapolator::EnergyAfterStep(G << 99 G4double stepLength,
91 G4double stepLength, << 100 const G4Material* mat,
92 const G4Material* mat, << 101 const G4ParticleDefinition* part)
93 const G4ParticleDefinition* par <<
94 { 102 {
>> 103 if(!isInitialised) Initialisation();
95 G4double kinEnergyFinal = kinEnergy; 104 G4double kinEnergyFinal = kinEnergy;
96 if(SetupKinematics(part, mat, kinEnergy)) { 105 if(SetupKinematics(part, mat, kinEnergy)) {
97 G4double step = TrueStepLength(kinEnergy,s 106 G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
98 G4double r = ComputeRange(kinEnergy,part, << 107 G4double r = ComputeRange(kinEnergy,part);
99 if(r <= step) { 108 if(r <= step) {
100 kinEnergyFinal = 0.0; 109 kinEnergyFinal = 0.0;
101 } else if(step < linLossLimit*r) { 110 } else if(step < linLossLimit*r) {
102 kinEnergyFinal -= step*ComputeDEDX(kinEn << 111 kinEnergyFinal -= step*ComputeDEDX(kinEnergy,part);
103 } else { 112 } else {
104 G4double r1 = r - step; 113 G4double r1 = r - step;
105 kinEnergyFinal = ComputeEnergy(r1,part,m << 114 kinEnergyFinal = ComputeEnergy(r1,part);
106 } 115 }
107 } 116 }
108 return kinEnergyFinal; 117 return kinEnergyFinal;
109 } 118 }
110 119
111 //....oooOO0OOooo........oooOO0OOooo........oo 120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
112 121
113 G4double << 122 G4double G4EnergyLossForExtrapolator::EnergyBeforeStep(G4double kinEnergy,
114 G4EnergyLossForExtrapolator::EnergyBeforeStep( << 123 G4double stepLength,
115 G4double stepLength, << 124 const G4Material* mat,
116 const G4Material* mat, << 125 const G4ParticleDefinition* part)
117 const G4ParticleDefinition* pa <<
118 { 126 {
119 //G4cout << "G4EnergyLossForExtrapolator::En << 127 if(!isInitialised) Initialisation();
120 G4double kinEnergyFinal = kinEnergy; 128 G4double kinEnergyFinal = kinEnergy;
121 129
122 if(SetupKinematics(part, mat, kinEnergy)) { 130 if(SetupKinematics(part, mat, kinEnergy)) {
123 G4double step = TrueStepLength(kinEnergy,s 131 G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
124 G4double r = ComputeRange(kinEnergy,part,m << 132 G4double r = ComputeRange(kinEnergy,part);
125 133
126 if(step < linLossLimit*r) { 134 if(step < linLossLimit*r) {
127 kinEnergyFinal += step*ComputeDEDX(kinEn << 135 kinEnergyFinal += step*ComputeDEDX(kinEnergy,part);
128 } else { 136 } else {
129 G4double r1 = r + step; 137 G4double r1 = r + step;
130 kinEnergyFinal = ComputeEnergy(r1,part,m << 138 kinEnergyFinal = ComputeEnergy(r1,part);
131 } 139 }
132 } 140 }
133 return kinEnergyFinal; 141 return kinEnergyFinal;
134 } 142 }
135 143
136 //....oooOO0OOooo........oooOO0OOooo........oo 144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
137 145
138 G4double << 146 G4double G4EnergyLossForExtrapolator::TrueStepLength(G4double kinEnergy,
139 G4EnergyLossForExtrapolator::TrueStepLength(G4 << 147 G4double stepLength,
140 G4double stepLength, << 148 const G4Material* mat,
141 const G4Material* mat, << 149 const G4ParticleDefinition* part)
142 const G4ParticleDefinition* part <<
143 { 150 {
144 G4double res = stepLength; 151 G4double res = stepLength;
145 //G4cout << "## G4EnergyLossForExtrapolator: << 152 if(!isInitialised) Initialisation();
146 // << " " << part->GetParticleName() << <<
147 if(SetupKinematics(part, mat, kinEnergy)) { 153 if(SetupKinematics(part, mat, kinEnergy)) {
148 if(part == electron || part == positron) { 154 if(part == electron || part == positron) {
149 const G4double x = stepLength* << 155 G4double x = stepLength*ComputeValue(kinEnergy, mscElectron);
150 ComputeValue(kinEnergy, GetPhysicsTable(fMsc << 156 if(x < 0.2) res *= (1.0 + 0.5*x + x*x/3.0);
151 //G4cout << " x= " << x << G4endl; << 157 else if(x < 0.9999) res = -std::log(1.0 - x)*stepLength/x;
152 if(x < 0.2) { res *= (1.0 + 0.5* << 158 else res = ComputeRange(kinEnergy,part);
153 else if(x < 0.9999) { res = -G4Log(1.0 - << 159
154 else { res = ComputeRange(kinEnergy, par <<
155 } else { 160 } else {
156 res = ComputeTrueStep(mat,part,kinEnergy 161 res = ComputeTrueStep(mat,part,kinEnergy,stepLength);
157 } 162 }
158 } 163 }
159 return res; 164 return res;
160 } 165 }
161 166
162 //....oooOO0OOooo........oooOO0OOooo........oo 167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
163 168
164 G4bool << 169 G4bool G4EnergyLossForExtrapolator::SetupKinematics(const G4ParticleDefinition* part,
165 G4EnergyLossForExtrapolator::SetupKinematics(c << 170 const G4Material* mat,
166 const G4Material* mat, << 171 G4double kinEnergy)
167 G4double kinEnergy) << 172 {
168 { << 173 if(!part || !mat || kinEnergy < keV) return false;
169 if(mat->GetNumberOfMaterials() != nmat) { In << 174 if(!isInitialised) Initialisation();
170 if(nullptr == part || nullptr == mat || kinE << 175 G4bool flag = false;
171 { return false; } <<
172 if(part != currentParticle) { 176 if(part != currentParticle) {
>> 177 flag = true;
173 currentParticle = part; 178 currentParticle = part;
>> 179 mass = part->GetPDGMass();
174 G4double q = part->GetPDGCharge()/eplus; 180 G4double q = part->GetPDGCharge()/eplus;
175 charge2 = q*q; 181 charge2 = q*q;
176 } 182 }
177 if(mat != currentMaterial) { 183 if(mat != currentMaterial) {
178 size_t i = mat->GetIndex(); << 184 G4int i = mat->GetIndex();
179 if(i >= nmat) { 185 if(i >= nmat) {
180 G4cout << "### G4EnergyLossForExtrapolat << 186 G4cout << "### G4EnergyLossForExtrapolator WARNING:index i= "
181 << i << " above number of materials " < << 187 << i << " is out of table - NO extrapolation" << G4endl;
182 return false; <<
183 } else { 188 } else {
>> 189 flag = true;
184 currentMaterial = mat; 190 currentMaterial = mat;
185 electronDensity = mat->GetElectronDensit 191 electronDensity = mat->GetElectronDensity();
186 radLength = mat->GetRadlen(); 192 radLength = mat->GetRadlen();
>> 193 index = i;
187 } 194 }
188 } 195 }
189 if(kinEnergy != kineticEnergy) { << 196 if(flag || kinEnergy != kineticEnergy) {
190 kineticEnergy = kinEnergy; 197 kineticEnergy = kinEnergy;
191 G4double mass = part->GetPDGMass(); <<
192 G4double tau = kinEnergy/mass; 198 G4double tau = kinEnergy/mass;
193 199
194 gam = tau + 1.0; 200 gam = tau + 1.0;
195 bg2 = tau * (tau + 2.0); 201 bg2 = tau * (tau + 2.0);
196 beta2 = bg2/(gam*gam); 202 beta2 = bg2/(gam*gam);
197 tmax = kinEnergy; 203 tmax = kinEnergy;
198 if(part == electron) tmax *= 0.5; 204 if(part == electron) tmax *= 0.5;
199 else if(part != positron) { 205 else if(part != positron) {
200 G4double r = CLHEP::electron_mass_c2/mas << 206 G4double r = electron_mass_c2/mass;
201 tmax = 2.0*bg2*CLHEP::electron_mass_c2/( << 207 tmax = 2.0*bg2*electron_mass_c2/(1.0 + 2.0*gam*r + r*r);
202 } 208 }
203 tmax = std::min(tmax, maxEnergyTransfer); << 209 if(tmax > maxEnergyTransfer) tmax = maxEnergyTransfer;
204 } 210 }
205 return true; 211 return true;
206 } 212 }
207 213
208 //....oooOO0OOooo........oooOO0OOooo........oo 214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
209 215
210 const G4ParticleDefinition* << 216 void G4EnergyLossForExtrapolator::Initialisation()
211 G4EnergyLossForExtrapolator::FindParticle(cons << 217 {
>> 218 isInitialised = true;
>> 219 if(verbose>1)
>> 220 G4cout << "### G4EnergyLossForExtrapolator::Initialisation" << G4endl;
>> 221 currentParticle = 0;
>> 222 currentMaterial = 0;
>> 223 kineticEnergy = 0.0;
>> 224 electron = G4Electron::Electron();
>> 225 positron = G4Positron::Positron();
>> 226 proton = G4Proton::Proton();
>> 227 muonPlus = G4MuonPlus::MuonPlus();
>> 228 muonMinus= G4MuonMinus::MuonMinus();
>> 229
>> 230 currentParticleName = "";
>> 231
>> 232 linLossLimit = 0.01;
>> 233 emin = 1.*MeV;
>> 234 emax = 10.*TeV;
>> 235 nbins = 70;
>> 236
>> 237 nmat = G4Material::GetNumberOfMaterials();
>> 238 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
>> 239 cuts = new G4ProductionCuts();
>> 240
>> 241 const G4MaterialCutsCouple* couple;
>> 242 for(G4int i=0; i<nmat; i++) {
>> 243 couple = new G4MaterialCutsCouple((*mtable)[i],cuts);
>> 244 couples.push_back(couple);
>> 245 }
>> 246
>> 247 dedxElectron = PrepareTable();
>> 248 dedxPositron = PrepareTable();
>> 249 dedxMuon = PrepareTable();
>> 250 dedxProton = PrepareTable();
>> 251 rangeElectron = PrepareTable();
>> 252 rangePositron = PrepareTable();
>> 253 rangeMuon = PrepareTable();
>> 254 rangeProton = PrepareTable();
>> 255 invRangeElectron = PrepareTable();
>> 256 invRangePositron = PrepareTable();
>> 257 invRangeMuon = PrepareTable();
>> 258 invRangeProton = PrepareTable();
>> 259 mscElectron = PrepareTable();
>> 260
>> 261 G4LossTableBuilder builder;
>> 262
>> 263 if(verbose>1)
>> 264 G4cout << "### G4EnergyLossForExtrapolator Builds electron tables" << G4endl;
>> 265
>> 266 ComputeElectronDEDX(electron, dedxElectron);
>> 267 builder.BuildRangeTable(dedxElectron,rangeElectron);
>> 268 builder.BuildInverseRangeTable(rangeElectron, invRangeElectron);
>> 269
>> 270 if(verbose>1)
>> 271 G4cout << "### G4EnergyLossForExtrapolator Builds positron tables" << G4endl;
>> 272
>> 273 ComputeElectronDEDX(positron, dedxPositron);
>> 274 builder.BuildRangeTable(dedxPositron, rangePositron);
>> 275 builder.BuildInverseRangeTable(rangePositron, invRangePositron);
>> 276
>> 277 if(verbose>1)
>> 278 G4cout << "### G4EnergyLossForExtrapolator Builds muon tables" << G4endl;
>> 279
>> 280 ComputeMuonDEDX(muonPlus, dedxMuon);
>> 281 builder.BuildRangeTable(dedxMuon, rangeMuon);
>> 282 builder.BuildInverseRangeTable(rangeMuon, invRangeMuon);
>> 283
>> 284 if(verbose>1)
>> 285 G4cout << "### G4EnergyLossForExtrapolator Builds proton tables" << G4endl;
>> 286
>> 287 ComputeProtonDEDX(proton, dedxProton);
>> 288 builder.BuildRangeTable(dedxProton, rangeProton);
>> 289 builder.BuildInverseRangeTable(rangeProton, invRangeProton);
>> 290
>> 291 ComputeTrasportXS(electron, mscElectron);
>> 292 }
>> 293
>> 294 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 295
>> 296 G4PhysicsTable* G4EnergyLossForExtrapolator::PrepareTable()
212 { 297 {
213 currentParticle = G4ParticleTable::GetPartic << 298 G4PhysicsTable* table = new G4PhysicsTable();
214 if(nullptr == currentParticle) { << 299
215 G4cout << "### G4EnergyLossForExtrapolator << 300 for(G4int i=0; i<nmat; i++) {
216 << "FindParticle() fails to find " << nam << 301
>> 302 G4PhysicsVector* v = new G4PhysicsLogVector(emin, emax, nbins);
>> 303 v->SetSpline(G4LossTableManager::Instance()->SplineFlag());
>> 304 table->push_back(v);
217 } 305 }
218 return currentParticle; << 306 return table;
219 } 307 }
220 308
221 //....oooOO0OOooo........oooOO0OOooo........oo 309 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
222 310
223 G4double << 311 const G4ParticleDefinition* G4EnergyLossForExtrapolator::FindParticle(const G4String& name)
224 G4EnergyLossForExtrapolator::ComputeDEDX(G4dou << 312 {
225 const G4ParticleDefinition* part, << 313 const G4ParticleDefinition* p = 0;
226 const << 314 if(name != currentParticleName) {
>> 315 p = G4ParticleTable::GetParticleTable()->FindParticle(name);
>> 316 if(!p) {
>> 317 G4cout << "### G4EnergyLossForExtrapolator WARNING:FindParticle fails to find "
>> 318 << name << G4endl;
>> 319 }
>> 320 } else {
>> 321 p = currentParticle;
>> 322 }
>> 323 return p;
>> 324 }
>> 325
>> 326 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 327
>> 328 G4double G4EnergyLossForExtrapolator::ComputeDEDX(G4double kinEnergy,
>> 329 const G4ParticleDefinition* part)
227 { 330 {
228 if(mat->GetNumberOfMaterials() != nmat) { In <<
229 G4double x = 0.0; 331 G4double x = 0.0;
230 if(part == electron) { << 332 if(part == electron) x = ComputeValue(kinEnergy, dedxElectron);
231 x = ComputeValue(ekin, GetPhysicsTable(fDe << 333 else if(part == positron) x = ComputeValue(kinEnergy, dedxPositron);
232 } else if(part == positron) { << 334 else if(part == muonPlus || part == muonMinus) {
233 x = ComputeValue(ekin, GetPhysicsTable(fDe << 335 x = ComputeValue(kinEnergy, dedxMuon);
234 } else if(part == muonPlus || part == muonMi <<
235 x = ComputeValue(ekin, GetPhysicsTable(fDe <<
236 } else { 336 } else {
237 G4double e = ekin*CLHEP::proton_mass_c2/pa << 337 G4double e = kinEnergy*proton_mass_c2/mass;
238 G4double q = part->GetPDGCharge()/CLHEP::e << 338 x = ComputeValue(e, dedxProton)*charge2;
239 x = ComputeValue(e, GetPhysicsTable(fDedxP <<
240 } 339 }
241 return x; 340 return x;
242 } 341 }
243 342
244 //....oooOO0OOooo........oooOO0OOooo........oo 343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
245 344
246 G4double << 345 G4double G4EnergyLossForExtrapolator::ComputeRange(G4double kinEnergy,
247 G4EnergyLossForExtrapolator::ComputeRange(G4do << 346 const G4ParticleDefinition* part)
248 const G4ParticleDefinition* part, <<
249 const G4Material* mat) <<
250 { 347 {
251 if(mat->GetNumberOfMaterials() != nmat) { In <<
252 G4double x = 0.0; 348 G4double x = 0.0;
253 if(part == electron) { << 349 if(part == electron) x = ComputeValue(kinEnergy, rangeElectron);
254 x = ComputeValue(ekin, GetPhysicsTable(fRa << 350 else if(part == positron) x = ComputeValue(kinEnergy, rangePositron);
255 } else if(part == positron) { << 351 else if(part == muonPlus || part == muonMinus)
256 x = ComputeValue(ekin, GetPhysicsTable(fRa << 352 x = ComputeValue(kinEnergy, rangeMuon);
257 } else if(part == muonPlus || part == muonMi << 353 else {
258 x = ComputeValue(ekin, GetPhysicsTable(fRa << 354 G4double massratio = proton_mass_c2/mass;
259 } else { << 355 G4double e = kinEnergy*massratio;
260 G4double massratio = CLHEP::proton_mass_c2 << 356 x = ComputeValue(e, rangeProton)/(charge2*massratio);
261 G4double e = ekin*massratio; <<
262 G4double q = part->GetPDGCharge()/CLHEP::e <<
263 x = ComputeValue(e, GetPhysicsTable(fRange <<
264 /(q*q*massratio); <<
265 } 357 }
266 return x; 358 return x;
267 } 359 }
268 360
269 //....oooOO0OOooo........oooOO0OOooo........oo 361 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
270 362
271 G4double << 363 G4double G4EnergyLossForExtrapolator::ComputeEnergy(G4double range,
272 G4EnergyLossForExtrapolator::ComputeEnergy(G4d << 364 const G4ParticleDefinition* part)
273 const G4ParticleDefinition* part, <<
274 const G4Material* mat) <<
275 { 365 {
276 if(mat->GetNumberOfMaterials() != nmat) { In <<
277 G4double x = 0.0; 366 G4double x = 0.0;
278 if(part == electron) { << 367 if(part == electron) x = ComputeValue(range, invRangeElectron);
279 x = ComputeValue(range,GetPhysicsTable(fIn << 368 else if(part == positron) x = ComputeValue(range, invRangePositron);
280 } else if(part == positron) { << 369 else if(part == muonPlus || part == muonMinus)
281 x = ComputeValue(range,GetPhysicsTable(fIn << 370 x = ComputeValue(range, invRangeMuon);
282 } else if(part == muonPlus || part == muonMi << 371 else {
283 x = ComputeValue(range, GetPhysicsTable(fI << 372 G4double massratio = proton_mass_c2/mass;
284 } else { << 373 G4double r = range*massratio*charge2;
285 G4double massratio = CLHEP::proton_mass_c2 << 374 x = ComputeValue(r, invRangeProton)/massratio;
286 G4double q = part->GetPDGCharge()/CLHEP::e <<
287 G4double r = range*massratio*q*q; <<
288 x = ComputeValue(r, GetPhysicsTable(fInvRa <<
289 } 375 }
290 return x; 376 return x;
291 } 377 }
292 378
293 //....oooOO0OOooo........oooOO0OOooo........oo 379 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
294 << 380
295 G4double << 381 void G4EnergyLossForExtrapolator::ComputeElectronDEDX(const G4ParticleDefinition* part,
296 G4EnergyLossForExtrapolator::EnergyDispersion( << 382 G4PhysicsTable* table)
297 G4double stepLength, <<
298 const G4Material* mat, <<
299 const G4ParticleDefinition* pa <<
300 { 383 {
301 G4double sig2 = 0.0; << 384 G4DataVector v;
302 if(SetupKinematics(part, mat, kinEnergy)) { << 385 G4MollerBhabhaModel* ioni = new G4MollerBhabhaModel();
303 G4double step = ComputeTrueStep(mat,part,k << 386 G4eBremsstrahlungModel* brem = new G4eBremsstrahlungModel();
304 sig2 = (1.0/beta2 - 0.5) << 387 ioni->Initialise(part, v);
305 *CLHEP::twopi_mc2_rcl2*tmax*step*electro << 388 brem->Initialise(part, v);
>> 389
>> 390 mass = electron_mass_c2;
>> 391 charge2 = 1.0;
>> 392 currentParticle = part;
>> 393
>> 394 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
>> 395 if(0<verbose) {
>> 396 G4cout << "G4EnergyLossForExtrapolator::ComputeElectronDEDX for "
>> 397 << part->GetParticleName()
>> 398 << G4endl;
>> 399 }
>> 400 for(G4int i=0; i<nmat; i++) {
>> 401
>> 402 const G4Material* mat = (*mtable)[i];
>> 403 if(1<verbose)
>> 404 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
>> 405 const G4MaterialCutsCouple* couple = couples[i];
>> 406 G4PhysicsVector* aVector = (*table)[i];
>> 407
>> 408 for(G4int j=0; j<=nbins; j++) {
>> 409
>> 410 G4double e = aVector->Energy(j);
>> 411 G4double dedx = ioni->ComputeDEDX(couple,part,e,e) + brem->ComputeDEDX(couple,part,e,e);
>> 412 if(1<verbose) {
>> 413 G4cout << "j= " << j
>> 414 << " e(MeV)= " << e/MeV
>> 415 << " dedx(Mev/cm)= " << dedx*cm/MeV
>> 416 << " dedx(Mev.cm2/g)= " << dedx/((MeV*mat->GetDensity())/(g/cm2)) << G4endl;
>> 417 }
>> 418 aVector->PutValue(j,dedx);
>> 419 }
>> 420 }
>> 421 delete ioni;
>> 422 delete brem;
>> 423 }
>> 424
>> 425 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 426
>> 427 void G4EnergyLossForExtrapolator::ComputeMuonDEDX(const G4ParticleDefinition* part,
>> 428 G4PhysicsTable* table)
>> 429 {
>> 430 G4DataVector v;
>> 431 G4BetheBlochModel* ioni = new G4BetheBlochModel();
>> 432 G4MuPairProductionModel* pair = new G4MuPairProductionModel();
>> 433 G4MuBremsstrahlungModel* brem = new G4MuBremsstrahlungModel();
>> 434 ioni->Initialise(part, v);
>> 435 pair->Initialise(part, v);
>> 436 brem->Initialise(part, v);
>> 437
>> 438 mass = part->GetPDGMass();
>> 439 charge2 = 1.0;
>> 440 currentParticle = part;
>> 441
>> 442 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
>> 443
>> 444 if(0<verbose) {
>> 445 G4cout << "G4EnergyLossForExtrapolator::ComputeMuonDEDX for " << part->GetParticleName()
>> 446 << G4endl;
306 } 447 }
307 return sig2; << 448
>> 449 for(G4int i=0; i<nmat; i++) {
>> 450
>> 451 const G4Material* mat = (*mtable)[i];
>> 452 if(1<verbose)
>> 453 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
>> 454 const G4MaterialCutsCouple* couple = couples[i];
>> 455 G4PhysicsVector* aVector = (*table)[i];
>> 456 for(G4int j=0; j<=nbins; j++) {
>> 457
>> 458 G4double e = aVector->Energy(j);
>> 459 G4double dedx = ioni->ComputeDEDX(couple,part,e,e) +
>> 460 pair->ComputeDEDX(couple,part,e,e) +
>> 461 brem->ComputeDEDX(couple,part,e,e);
>> 462 aVector->PutValue(j,dedx);
>> 463 if(1<verbose) {
>> 464 G4cout << "j= " << j
>> 465 << " e(MeV)= " << e/MeV
>> 466 << " dedx(Mev/cm)= " << dedx*cm/MeV
>> 467 << " dedx(Mev/(g/cm2)= " << dedx/((MeV*mat->GetDensity())/(g/cm2)) << G4endl;
>> 468 }
>> 469 }
>> 470 }
>> 471 delete ioni;
308 } 472 }
309 473
310 //....oooOO0OOooo........oooOO0OOooo........oo 474 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
311 475
312 G4double G4EnergyLossForExtrapolator::AverageS << 476 void G4EnergyLossForExtrapolator::ComputeProtonDEDX(const G4ParticleDefinition* part,
313 G4double kinEnergy, << 477 G4PhysicsTable* table)
314 G4double stepLength, <<
315 const G4Material* mat, <<
316 const G4ParticleDefinition* part) <<
317 { 478 {
318 G4double theta = 0.0; << 479 G4DataVector v;
319 if(SetupKinematics(part, mat, kinEnergy)) { << 480 G4BetheBlochModel* ioni = new G4BetheBlochModel();
320 G4double t = stepLength/radLength; << 481 ioni->Initialise(part, v);
321 G4double y = std::max(0.001, t); << 482
322 theta = 19.23*CLHEP::MeV*std::sqrt(charge2 << 483 mass = part->GetPDGMass();
323 /(beta2*gam*part->GetPDGMass()); << 484 charge2 = 1.0;
>> 485 currentParticle = part;
>> 486
>> 487 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
>> 488
>> 489 if(0<verbose) {
>> 490 G4cout << "G4EnergyLossForExtrapolator::ComputeProtonDEDX for " << part->GetParticleName()
>> 491 << G4endl;
>> 492 }
>> 493
>> 494 for(G4int i=0; i<nmat; i++) {
>> 495
>> 496 const G4Material* mat = (*mtable)[i];
>> 497 if(1<verbose)
>> 498 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
>> 499 const G4MaterialCutsCouple* couple = couples[i];
>> 500 G4PhysicsVector* aVector = (*table)[i];
>> 501 for(G4int j=0; j<=nbins; j++) {
>> 502
>> 503 G4double e = aVector->Energy(j);
>> 504 G4double dedx = ioni->ComputeDEDX(couple,part,e,e);
>> 505 aVector->PutValue(j,dedx);
>> 506 if(1<verbose) {
>> 507 G4cout << "j= " << j
>> 508 << " e(MeV)= " << e/MeV
>> 509 << " dedx(Mev/cm)= " << dedx*cm/MeV
>> 510 << " dedx(Mev.cm2/g)= " << dedx/((mat->GetDensity())/(g/cm2)) << G4endl;
>> 511 }
>> 512 }
324 } 513 }
325 return theta; << 514 delete ioni;
326 } 515 }
327 516
328 //....oooOO0OOooo........oooOO0OOooo........oo 517 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
329 518
330 void G4EnergyLossForExtrapolator::Initialisati << 519 void G4EnergyLossForExtrapolator::ComputeTrasportXS(const G4ParticleDefinition* part,
>> 520 G4PhysicsTable* table)
331 { 521 {
332 if(verbose>0) { << 522 G4DataVector v;
333 G4cout << "### G4EnergyLossForExtrapolator << 523 G4WentzelVIModel* msc = new G4WentzelVIModel();
334 << tables << G4endl; << 524 msc->SetPolarAngleLimit(CLHEP::pi);
>> 525 msc->Initialise(part, v);
>> 526
>> 527 mass = part->GetPDGMass();
>> 528 charge2 = 1.0;
>> 529 currentParticle = part;
>> 530
>> 531 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
>> 532
>> 533 if(0<verbose) {
>> 534 G4cout << "G4EnergyLossForExtrapolator::ComputeProtonDEDX for " << part->GetParticleName()
>> 535 << G4endl;
335 } 536 }
336 electron = G4Electron::Electron(); << 537
337 positron = G4Positron::Positron(); << 538 for(G4int i=0; i<nmat; i++) {
338 proton = G4Proton::Proton(); <<
339 muonPlus = G4MuonPlus::MuonPlus(); <<
340 muonMinus= G4MuonMinus::MuonMinus(); <<
341 539
342 // initialisation for the 1st run << 540 const G4Material* mat = (*mtable)[i];
343 if(nullptr == tables) { << 541 if(1<verbose)
344 #ifdef G4MULTITHREADED << 542 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
345 G4MUTEXLOCK(&extrMutex); << 543 G4PhysicsVector* aVector = (*table)[i];
346 if(nullptr == tables) { << 544 for(G4int j=0; j<=nbins; j++) {
347 #endif << 545
348 isMaster = true; << 546 G4double e = aVector->Energy(j);
349 tables = new G4TablesForExtrapolator(ver << 547 G4double xs = msc->CrossSectionPerVolume(mat,part,e);
350 tables->Initialisation(); << 548 aVector->PutValue(j,xs);
351 nmat = G4Material::GetNumberOfMaterials( << 549 if(1<verbose) {
352 if(verbose > 0) { << 550 G4cout << "j= " << j << " e(MeV)= " << e/MeV
353 G4cout << "### G4EnergyLossForExtrapol << 551 << " xs(1/mm)= " << xs*mm << G4endl;
354 << nmat << " materials Nbins= " << 552 }
355 << nbins << " Emin(MeV)= " << e <<
356 << G4endl; <<
357 } <<
358 #ifdef G4MULTITHREADED <<
359 } 553 }
360 G4MUTEXUNLOCK(&extrMutex); <<
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 <<
373 } 554 }
374 nmat = G4Material::GetNumberOfMaterials(); << 555 delete msc;
375 } 556 }
376 557
377 //....oooOO0OOooo........oooOO0OOooo........oo 558 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 559
378 560