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