| 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 // GEANT4 Class file 29 // GEANT4 Class file
30 // 30 //
31 // 31 //
32 // File name: G4BraggModel 32 // File name: G4BraggModel
33 // 33 //
34 // Author: Vladimir Ivanchenko 34 // Author: Vladimir Ivanchenko
35 // 35 //
36 // Creation date: 03.01.2002 36 // Creation date: 03.01.2002
37 // 37 //
38 // Modifications: 38 // Modifications:
39 // 39 //
40 // 04-12-02 Fix problem of G4DynamicParticle c 40 // 04-12-02 Fix problem of G4DynamicParticle constructor (V.Ivanchenko)
41 // 23-12-02 Change interface in order to move 41 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
42 // 27-01-03 Make models region aware (V.Ivanch 42 // 27-01-03 Make models region aware (V.Ivanchenko)
43 // 13-02-03 Add name (V.Ivanchenko) 43 // 13-02-03 Add name (V.Ivanchenko)
44 // 04-06-03 Fix compilation warnings (V.Ivanch 44 // 04-06-03 Fix compilation warnings (V.Ivanchenko)
45 // 12-09-04 Add lowestKinEnergy and change ord 45 // 12-09-04 Add lowestKinEnergy and change order of if in DEDX method (VI)
46 // 11-04-05 Major optimisation of internal int 46 // 11-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
47 // 16-06-05 Fix problem of chemical formula (V 47 // 16-06-05 Fix problem of chemical formula (V.Ivantchenko)
48 // 15-02-06 ComputeCrossSectionPerElectron, Co 48 // 15-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
49 // 25-04-06 Add stopping data from PSTAR (V.Iv 49 // 25-04-06 Add stopping data from PSTAR (V.Ivanchenko)
50 // 12-08-08 Added methods GetParticleCharge, G 50 // 12-08-08 Added methods GetParticleCharge, GetChargeSquareRatio,
51 // CorrectionsAlongStep needed for io 51 // CorrectionsAlongStep needed for ions(V.Ivanchenko)
52 52
53 // Class Description: 53 // Class Description:
54 // 54 //
55 // Implementation of energy loss and delta-ele 55 // Implementation of energy loss and delta-electron production by
56 // slow charged heavy particles 56 // slow charged heavy particles
57 57
58 // ------------------------------------------- 58 // -------------------------------------------------------------------
59 // 59 //
60 60
61 61
62 //....oooOO0OOooo........oooOO0OOooo........oo 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
63 //....oooOO0OOooo........oooOO0OOooo........oo 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
64 64
65 #include "G4BraggModel.hh" 65 #include "G4BraggModel.hh"
66 #include "G4PhysicalConstants.hh" 66 #include "G4PhysicalConstants.hh"
67 #include "G4SystemOfUnits.hh" 67 #include "G4SystemOfUnits.hh"
68 #include "Randomize.hh" 68 #include "Randomize.hh"
69 #include "G4Electron.hh" 69 #include "G4Electron.hh"
70 #include "G4ParticleChangeForLoss.hh" 70 #include "G4ParticleChangeForLoss.hh"
71 #include "G4LossTableManager.hh" 71 #include "G4LossTableManager.hh"
72 #include "G4EmCorrections.hh" 72 #include "G4EmCorrections.hh"
73 #include "G4EmParameters.hh" 73 #include "G4EmParameters.hh"
74 #include "G4DeltaAngle.hh" 74 #include "G4DeltaAngle.hh"
75 #include "G4ICRU90StoppingData.hh" 75 #include "G4ICRU90StoppingData.hh"
76 #include "G4NistManager.hh" 76 #include "G4NistManager.hh"
77 #include "G4Log.hh" 77 #include "G4Log.hh"
78 #include "G4Exp.hh" 78 #include "G4Exp.hh"
79 #include "G4AutoLock.hh" 79 #include "G4AutoLock.hh"
80 80
81 //....oooOO0OOooo........oooOO0OOooo........oo 81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82 82
83 G4ICRU90StoppingData* G4BraggModel::fICRU90 = 83 G4ICRU90StoppingData* G4BraggModel::fICRU90 = nullptr;
84 G4PSTARStopping* G4BraggModel::fPSTAR = nullpt 84 G4PSTARStopping* G4BraggModel::fPSTAR = nullptr;
85 85
86 namespace 86 namespace
87 { 87 {
88 G4Mutex ionMutex = G4MUTEX_INITIALIZER; 88 G4Mutex ionMutex = G4MUTEX_INITIALIZER;
89 } 89 }
90 90
91 G4BraggModel::G4BraggModel(const G4ParticleDef 91 G4BraggModel::G4BraggModel(const G4ParticleDefinition* p, const G4String& nam)
92 : G4VEmModel(nam) 92 : G4VEmModel(nam)
93 { 93 {
94 SetHighEnergyLimit(2.0*CLHEP::MeV); 94 SetHighEnergyLimit(2.0*CLHEP::MeV);
95 95
96 lowestKinEnergy = 0.25*CLHEP::keV; 96 lowestKinEnergy = 0.25*CLHEP::keV;
97 theZieglerFactor = CLHEP::eV*CLHEP::cm2*1.0e 97 theZieglerFactor = CLHEP::eV*CLHEP::cm2*1.0e-15;
98 theElectron = G4Electron::Electron(); 98 theElectron = G4Electron::Electron();
99 expStopPower125 = 0.0; 99 expStopPower125 = 0.0;
100 100
101 corr = G4LossTableManager::Instance()->EmCor 101 corr = G4LossTableManager::Instance()->EmCorrections();
102 if(nullptr != p) { SetParticle(p); } 102 if(nullptr != p) { SetParticle(p); }
103 else { SetParticle(theElectron); } 103 else { SetParticle(theElectron); }
104 } 104 }
105 105
106 //....oooOO0OOooo........oooOO0OOooo........oo 106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107 107
108 G4BraggModel::~G4BraggModel() 108 G4BraggModel::~G4BraggModel()
109 { 109 {
110 if(isFirst) { 110 if(isFirst) {
111 delete fPSTAR; 111 delete fPSTAR;
112 fPSTAR = nullptr; 112 fPSTAR = nullptr;
113 } 113 }
114 } 114 }
115 115
116 //....oooOO0OOooo........oooOO0OOooo........oo 116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117 117
118 void G4BraggModel::Initialise(const G4Particle 118 void G4BraggModel::Initialise(const G4ParticleDefinition* p,
119 const G4DataVect 119 const G4DataVector&)
120 { 120 {
121 if(p != particle) { SetParticle(p); } 121 if(p != particle) { SetParticle(p); }
122 122
123 // always false before the run 123 // always false before the run
124 SetDeexcitationFlag(false); 124 SetDeexcitationFlag(false);
125 125
126 // initialise data only once 126 // initialise data only once
127 if(nullptr == fPSTAR) { 127 if(nullptr == fPSTAR) {
128 G4AutoLock l(&ionMutex); 128 G4AutoLock l(&ionMutex);
129 if(nullptr == fPSTAR) { 129 if(nullptr == fPSTAR) {
130 isFirst = true; 130 isFirst = true;
131 fPSTAR = new G4PSTARStopping(); 131 fPSTAR = new G4PSTARStopping();
132 if(G4EmParameters::Instance()->UseICRU90 132 if(G4EmParameters::Instance()->UseICRU90Data()) {
133 fICRU90 = G4NistManager::Instance()->GetICRU 133 fICRU90 = G4NistManager::Instance()->GetICRU90StoppingData();
134 } 134 }
135 } 135 }
136 l.unlock(); 136 l.unlock();
137 } 137 }
138 if(isFirst) { 138 if(isFirst) {
139 if(nullptr != fICRU90) { fICRU90->Initiali 139 if(nullptr != fICRU90) { fICRU90->Initialise(); }
140 fPSTAR->Initialise(); 140 fPSTAR->Initialise();
141 } 141 }
142 142
143 if(nullptr == fParticleChange) { 143 if(nullptr == fParticleChange) {
144 144
145 if(UseAngularGeneratorFlag() && !GetAngula 145 if(UseAngularGeneratorFlag() && !GetAngularDistribution()) {
146 SetAngularDistribution(new G4DeltaAngle( 146 SetAngularDistribution(new G4DeltaAngle());
147 } 147 }
148 G4String pname = particle->GetParticleName 148 G4String pname = particle->GetParticleName();
149 if(particle->GetParticleType() == "nucleus 149 if(particle->GetParticleType() == "nucleus" &&
150 pname != "deuteron" && pname != "triton 150 pname != "deuteron" && pname != "triton" &&
151 pname != "alpha+" && pname != "helium 151 pname != "alpha+" && pname != "helium" &&
152 pname != "hydrogen") { isIon = true; } 152 pname != "hydrogen") { isIon = true; }
153 153
154 fParticleChange = GetParticleChangeForLoss 154 fParticleChange = GetParticleChangeForLoss();
155 } 155 }
156 } 156 }
157 157
158 //....oooOO0OOooo........oooOO0OOooo........oo 158 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
159 159
160 void G4BraggModel::SetParticle(const G4Particl 160 void G4BraggModel::SetParticle(const G4ParticleDefinition* p)
161 { 161 {
162 particle = p; 162 particle = p;
163 mass = particle->GetPDGMass(); 163 mass = particle->GetPDGMass();
164 spin = particle->GetPDGSpin(); 164 spin = particle->GetPDGSpin();
165 G4double q = particle->GetPDGCharge()/CLHEP: 165 G4double q = particle->GetPDGCharge()/CLHEP::eplus;
166 chargeSquare = q*q; 166 chargeSquare = q*q;
167 massRate = mass/CLHEP::proton_mass_c2; 167 massRate = mass/CLHEP::proton_mass_c2;
168 ratio = CLHEP::electron_mass_c2/mass; 168 ratio = CLHEP::electron_mass_c2/mass;
169 } 169 }
170 170
171 //....oooOO0OOooo........oooOO0OOooo........oo 171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
172 172
173 G4double G4BraggModel::GetChargeSquareRatio(co 173 G4double G4BraggModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
174 co 174 const G4Material* mat,
175 G4 175 G4double kinEnergy)
176 { 176 {
177 // this method is called only for ions 177 // this method is called only for ions
178 chargeSquare = corr->EffectiveChargeSquareRa 178 chargeSquare = corr->EffectiveChargeSquareRatio(p, mat, kinEnergy);
179 return chargeSquare; 179 return chargeSquare;
180 } 180 }
181 181
182 //....oooOO0OOooo........oooOO0OOooo........oo 182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
183 183
184 G4double G4BraggModel::MinEnergyCut(const G4Pa 184 G4double G4BraggModel::MinEnergyCut(const G4ParticleDefinition*,
185 const G4Ma 185 const G4MaterialCutsCouple* couple)
186 { 186 {
187 return couple->GetMaterial()->GetIonisation( 187 return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
188 } 188 }
189 189
190 //....oooOO0OOooo........oooOO0OOooo........oo 190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
191 191
192 G4double G4BraggModel::GetParticleCharge(const 192 G4double G4BraggModel::GetParticleCharge(const G4ParticleDefinition* p,
193 const 193 const G4Material* mat,
194 G4dou 194 G4double kineticEnergy)
195 { 195 {
196 // this method is called only for ions, so n 196 // this method is called only for ions, so no check if it is an ion
197 return corr->GetParticleCharge(p,mat,kinetic 197 return corr->GetParticleCharge(p,mat,kineticEnergy);
198 } 198 }
199 199
200 //....oooOO0OOooo........oooOO0OOooo........oo 200 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
201 201
202 G4double G4BraggModel::ComputeCrossSectionPerE 202 G4double G4BraggModel::ComputeCrossSectionPerElectron(
203 con 203 const G4ParticleDefinition* p,
204 204 G4double kineticEnergy,
205 205 G4double cut,
206 206 G4double maxKinEnergy)
207 { 207 {
208 G4double cross = 0.0; 208 G4double cross = 0.0;
209 const G4double tmax = MaxSecondaryEnerg 209 const G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
210 const G4double maxEnergy = std::min(tmax, ma 210 const G4double maxEnergy = std::min(tmax, maxKinEnergy);
211 const G4double cutEnergy = std::max(cut, low 211 const G4double cutEnergy = std::max(cut, lowestKinEnergy*massRate);
212 if(cutEnergy < maxEnergy) { 212 if(cutEnergy < maxEnergy) {
213 213
214 const G4double energy = kineticEnergy + m 214 const G4double energy = kineticEnergy + mass;
215 const G4double energy2 = energy*energy; 215 const G4double energy2 = energy*energy;
216 const G4double beta2 = kineticEnergy*(ki 216 const G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
217 cross = (maxEnergy - cutEnergy)/(cutEnergy 217 cross = (maxEnergy - cutEnergy)/(cutEnergy*maxEnergy)
218 - beta2*G4Log(maxEnergy/cutEnergy)/tmax; 218 - beta2*G4Log(maxEnergy/cutEnergy)/tmax;
219 219
220 if( 0.0 < spin ) { cross += 0.5*(maxEnergy 220 if( 0.0 < spin ) { cross += 0.5*(maxEnergy - cutEnergy)/energy2; }
221 221
222 cross *= CLHEP::twopi_mc2_rcl2*chargeSquar 222 cross *= CLHEP::twopi_mc2_rcl2*chargeSquare/beta2;
223 } 223 }
224 // G4cout << "BR: e= " << kineticEnergy << 224 // G4cout << "BR: e= " << kineticEnergy << " tmin= " << cutEnergy
225 // << " tmax= " << tmax << " cross= 225 // << " tmax= " << tmax << " cross= " << cross << G4endl;
226 return cross; 226 return cross;
227 } 227 }
228 228
229 //....oooOO0OOooo........oooOO0OOooo........oo 229 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
230 230
231 G4double 231 G4double
232 G4BraggModel::ComputeCrossSectionPerAtom(const 232 G4BraggModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition* p,
233 G4dou 233 G4double kineticEnergy,
234 G4dou 234 G4double Z, G4double,
235 G4dou 235 G4double cutEnergy,
236 G4dou 236 G4double maxEnergy)
237 { 237 {
238 return 238 return
239 Z*ComputeCrossSectionPerElectron(p,kinetic 239 Z*ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
240 } 240 }
241 241
242 //....oooOO0OOooo........oooOO0OOooo........oo 242 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
243 243
244 G4double G4BraggModel::CrossSectionPerVolume(c 244 G4double G4BraggModel::CrossSectionPerVolume(const G4Material* material,
245 c 245 const G4ParticleDefinition* p,
246 G 246 G4double kineticEnergy,
247 G 247 G4double cutEnergy,
248 G 248 G4double maxEnergy)
249 { 249 {
250 return material->GetElectronDensity() 250 return material->GetElectronDensity()
251 *ComputeCrossSectionPerElectron(p,kineticE 251 *ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
252 } 252 }
253 253
254 //....oooOO0OOooo........oooOO0OOooo........oo 254 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
255 255
256 G4double G4BraggModel::ComputeDEDXPerVolume(co 256 G4double G4BraggModel::ComputeDEDXPerVolume(const G4Material* material,
257 co 257 const G4ParticleDefinition* p,
258 G4 258 G4double kinEnergy,
259 G4 259 G4double cut)
260 { 260 {
261 const G4double tmax = MaxSecondaryEnergy(p, 261 const G4double tmax = MaxSecondaryEnergy(p, kinEnergy);
262 const G4double tkin = kinEnergy/massRate; << 262 const G4double tlim = lowestKinEnergy*massRate;
263 const G4double cutEnergy = std::max(cut, low << 263 const G4double tmin = std::max(std::min(cut, tmax), tlim);
264 G4double dedx = 0.0; 264 G4double dedx = 0.0;
265 265
266 // tkin is the scaled energy to proton << 266 if(kinEnergy < tlim) {
267 if (tkin < lowestKinEnergy) { << 267 dedx = DEDX(material, lowestKinEnergy)*std::sqrt(kinEnergy/tlim);
268 dedx = DEDX(material, lowestKinEnergy)*std <<
269 } else { 268 } else {
270 dedx = DEDX(material, tkin); << 269 dedx = DEDX(material, kinEnergy);
271 270
272 // correction for low cut value using Beth << 271 if (tmin < tmax) {
273 if (cutEnergy < tmax) { <<
274 const G4double tau = kinEnergy/mass; 272 const G4double tau = kinEnergy/mass;
275 const G4double x = cutEnergy/tmax; << 273 const G4double x = tmin/tmax;
276 274
277 dedx += (G4Log(x)*(tau + 1.)*(tau + 1.)/ 275 dedx += (G4Log(x)*(tau + 1.)*(tau + 1.)/(tau * (tau + 2.0)) + 1.0 - x) *
278 CLHEP::twopi_mc2_rcl2 * material->GetElectro 276 CLHEP::twopi_mc2_rcl2 * material->GetElectronDensity();
279 } 277 }
280 } 278 }
281 dedx = std::max(dedx, 0.0) * chargeSquare; 279 dedx = std::max(dedx, 0.0) * chargeSquare;
282 280
283 //G4cout << "E(MeV)= " << tkin/MeV << " dedx 281 //G4cout << "E(MeV)= " << tkin/MeV << " dedx= " << dedx
284 // << " " << material->GetName() << 282 // << " " << material->GetName() << G4endl;
285 return dedx; 283 return dedx;
286 } 284 }
287 285
288 //....oooOO0OOooo........oooOO0OOooo........oo 286 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
289 287
290 void G4BraggModel::SampleSecondaries(std::vect 288 void G4BraggModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
291 const G4M 289 const G4MaterialCutsCouple* couple,
292 const G4D 290 const G4DynamicParticle* dp,
293 G4double 291 G4double minEnergy,
294 G4double 292 G4double maxEnergy)
295 { 293 {
296 const G4double tmax = MaxSecondaryKinEnergy( 294 const G4double tmax = MaxSecondaryKinEnergy(dp);
297 const G4double xmax = std::min(tmax, maxEner 295 const G4double xmax = std::min(tmax, maxEnergy);
298 const G4double xmin = std::max(lowestKinEner 296 const G4double xmin = std::max(lowestKinEnergy*massRate, std::min(minEnergy, xmax));
299 if(xmin >= xmax) { return; } 297 if(xmin >= xmax) { return; }
300 298
301 G4double kineticEnergy = dp->GetKineticEnerg 299 G4double kineticEnergy = dp->GetKineticEnergy();
302 const G4double energy = kineticEnergy + mas 300 const G4double energy = kineticEnergy + mass;
303 const G4double energy2 = energy*energy; 301 const G4double energy2 = energy*energy;
304 const G4double beta2 = kineticEnergy*(kine 302 const G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
305 const G4double grej = 1.0; 303 const G4double grej = 1.0;
306 G4double deltaKinEnergy, f; 304 G4double deltaKinEnergy, f;
307 305
308 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra 306 CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
309 G4double rndm[2]; 307 G4double rndm[2];
310 308
311 // sampling follows ... 309 // sampling follows ...
312 do { 310 do {
313 rndmEngineMod->flatArray(2, rndm); 311 rndmEngineMod->flatArray(2, rndm);
314 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - rn 312 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - rndm[0]) + xmax*rndm[0]);
315 313
316 f = 1.0 - beta2*deltaKinEnergy/tmax; 314 f = 1.0 - beta2*deltaKinEnergy/tmax;
317 315
318 if(f > grej) { 316 if(f > grej) {
319 G4cout << "G4BraggModel::SampleSeconda 317 G4cout << "G4BraggModel::SampleSecondary Warning! "
320 << "Majorant " << grej << " < " 318 << "Majorant " << grej << " < "
321 << f << " for e= " << deltaKinE 319 << f << " for e= " << deltaKinEnergy
322 << G4endl; 320 << G4endl;
323 } 321 }
324 322
325 // Loop checking, 03-Aug-2015, Vladimir Iv 323 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
326 } while( grej*rndm[1] >= f ); 324 } while( grej*rndm[1] >= f );
327 325
328 G4ThreeVector deltaDirection; 326 G4ThreeVector deltaDirection;
329 327
330 if(UseAngularGeneratorFlag()) { 328 if(UseAngularGeneratorFlag()) {
331 const G4Material* mat = couple->GetMateri 329 const G4Material* mat = couple->GetMaterial();
332 G4int Z = SelectRandomAtomNumber(mat); 330 G4int Z = SelectRandomAtomNumber(mat);
333 331
334 deltaDirection = 332 deltaDirection =
335 GetAngularDistribution()->SampleDirectio 333 GetAngularDistribution()->SampleDirection(dp, deltaKinEnergy, Z, mat);
336 334
337 } else { 335 } else {
338 336
339 G4double deltaMomentum = 337 G4double deltaMomentum =
340 std::sqrt(deltaKinEnergy * (deltaKinEner 338 std::sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
341 G4double cost = deltaKinEnergy * (energy + 339 G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
342 (deltaMomentum * dp->GetTotalMomentum()) 340 (deltaMomentum * dp->GetTotalMomentum());
343 if(cost > 1.0) { cost = 1.0; } 341 if(cost > 1.0) { cost = 1.0; }
344 G4double sint = std::sqrt((1.0 - cost)*(1. 342 G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost));
345 343
346 G4double phi = twopi*rndmEngineMod->flat() 344 G4double phi = twopi*rndmEngineMod->flat();
347 345
348 deltaDirection.set(sint*std::cos(phi),sint 346 deltaDirection.set(sint*std::cos(phi),sint*std::sin(phi), cost) ;
349 deltaDirection.rotateUz(dp->GetMomentumDir 347 deltaDirection.rotateUz(dp->GetMomentumDirection());
350 } 348 }
351 349
352 // create G4DynamicParticle object for delta 350 // create G4DynamicParticle object for delta ray
353 auto delta = new G4DynamicParticle(theElectr 351 auto delta = new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
354 352
355 // Change kinematics of primary particle 353 // Change kinematics of primary particle
356 kineticEnergy -= deltaKinEnergy; 354 kineticEnergy -= deltaKinEnergy;
357 G4ThreeVector finalP = dp->GetMomentum() - d 355 G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum();
358 finalP = finalP.unit(); 356 finalP = finalP.unit();
359 357
360 fParticleChange->SetProposedKineticEnergy(ki 358 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
361 fParticleChange->SetProposedMomentumDirectio 359 fParticleChange->SetProposedMomentumDirection(finalP);
362 360
363 vdp->push_back(delta); 361 vdp->push_back(delta);
364 } 362 }
365 363
366 //....oooOO0OOooo........oooOO0OOooo........oo 364 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
367 365
368 G4double G4BraggModel::MaxSecondaryEnergy(cons 366 G4double G4BraggModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
369 G4do 367 G4double kinEnergy)
370 { 368 {
371 if(pd != particle) { SetParticle(pd); } 369 if(pd != particle) { SetParticle(pd); }
372 G4double tau = kinEnergy/mass; 370 G4double tau = kinEnergy/mass;
373 G4double tmax = 2.0*electron_mass_c2*tau*(ta 371 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
374 (1. + 2.0*(tau + 1.)*ratio + 372 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
375 return tmax; 373 return tmax;
376 } 374 }
377 375
378 //....oooOO0OOooo........oooOO0OOooo........oo 376 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
379 377
380 void G4BraggModel::HasMaterial(const G4Materia 378 void G4BraggModel::HasMaterial(const G4Material* mat)
381 { 379 {
382 const G4String& chFormula = mat->GetChemical 380 const G4String& chFormula = mat->GetChemicalFormula();
383 if(chFormula.empty()) { return; } 381 if(chFormula.empty()) { return; }
384 382
385 // ICRU Report N49, 1993. Power's model for 383 // ICRU Report N49, 1993. Power's model for H
386 static const G4int numberOfMolecula = 11; 384 static const G4int numberOfMolecula = 11;
387 static const G4String molName[numberOfMolecu 385 static const G4String molName[numberOfMolecula] = {
388 "Al_2O_3", "CO_2", 386 "Al_2O_3", "CO_2", "CH_4",
389 "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Pol 387 "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Polypropylene", "(C_8H_8)_N",
390 "C_3H_8", "SiO_2", 388 "C_3H_8", "SiO_2", "H_2O",
391 "H_2O-Gas", "Graphite" } ; 389 "H_2O-Gas", "Graphite" } ;
392 390
393 // Search for the material in the table 391 // Search for the material in the table
394 for (G4int i=0; i<numberOfMolecula; ++i) { 392 for (G4int i=0; i<numberOfMolecula; ++i) {
395 if (chFormula == molName[i]) { 393 if (chFormula == molName[i]) {
396 iMolecula = i; 394 iMolecula = i;
397 return; 395 return;
398 } 396 }
399 } 397 }
400 return; 398 return;
401 } 399 }
402 400
403 //....oooOO0OOooo........oooOO0OOooo........oo 401 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
404 402
405 G4double G4BraggModel::StoppingPower(const G4M 403 G4double G4BraggModel::StoppingPower(const G4Material* material,
406 G4double 404 G4double kineticEnergy)
407 { 405 {
408 G4double ionloss = 0.0 ; 406 G4double ionloss = 0.0 ;
409 407
410 if (iMolecula >= 0) { 408 if (iMolecula >= 0) {
411 409
412 // The data and the fit from: 410 // The data and the fit from:
413 // ICRU Report N49, 1993. Ziegler's model 411 // ICRU Report N49, 1993. Ziegler's model for protons.
414 // Proton kinetic energy for parametrisati 412 // Proton kinetic energy for parametrisation (keV/amu)
415 413
416 G4double T = kineticEnergy/(keV*protonMass 414 G4double T = kineticEnergy/(keV*protonMassAMU) ;
417 415
418 static const G4float a[11][5] = { 416 static const G4float a[11][5] = {
419 {1.187E+1f, 1.343E+1f, 1.069E+4f, 7.723E+2f 417 {1.187E+1f, 1.343E+1f, 1.069E+4f, 7.723E+2f, 2.153E-2f},
420 {7.802E+0f, 8.814E+0f, 8.303E+3f, 7.446E+2f 418 {7.802E+0f, 8.814E+0f, 8.303E+3f, 7.446E+2f, 7.966E-3f},
421 {7.294E+0f, 8.284E+0f, 5.010E+3f, 4.544E+2f 419 {7.294E+0f, 8.284E+0f, 5.010E+3f, 4.544E+2f, 8.153E-3f},
422 {8.646E+0f, 9.800E+0f, 7.066E+3f, 4.581E+2f 420 {8.646E+0f, 9.800E+0f, 7.066E+3f, 4.581E+2f, 9.383E-3f},
423 {1.286E+1f, 1.462E+1f, 5.625E+3f, 2.621E+3f 421 {1.286E+1f, 1.462E+1f, 5.625E+3f, 2.621E+3f, 3.512E-2f},
424 {3.229E+1f, 3.696E+1f, 8.918E+3f, 3.244E+3f 422 {3.229E+1f, 3.696E+1f, 8.918E+3f, 3.244E+3f, 1.273E-1f},
425 {1.604E+1f, 1.825E+1f, 6.967E+3f, 2.307E+3f 423 {1.604E+1f, 1.825E+1f, 6.967E+3f, 2.307E+3f, 3.775E-2f},
426 {8.049E+0f, 9.099E+0f, 9.257E+3f, 3.846E+2f 424 {8.049E+0f, 9.099E+0f, 9.257E+3f, 3.846E+2f, 1.007E-2f},
427 {4.015E+0f, 4.542E+0f, 3.955E+3f, 4.847E+2f 425 {4.015E+0f, 4.542E+0f, 3.955E+3f, 4.847E+2f, 7.904E-3f},
428 {4.571E+0f, 5.173E+0f, 4.346E+3f, 4.779E+2f 426 {4.571E+0f, 5.173E+0f, 4.346E+3f, 4.779E+2f, 8.572E-3f},
429 {2.631E+0f, 2.601E+0f, 1.701E+3f, 1.279E+3f 427 {2.631E+0f, 2.601E+0f, 1.701E+3f, 1.279E+3f, 1.638E-2f} };
430 428
431 static const G4float atomicWeight[11] = { 429 static const G4float atomicWeight[11] = {
432 101.96128f, 44.0098f, 16.0426f, 28.0536f, 430 101.96128f, 44.0098f, 16.0426f, 28.0536f, 42.0804f,
433 104.1512f, 44.665f, 60.0843f, 18.0152f, 18 431 104.1512f, 44.665f, 60.0843f, 18.0152f, 18.0152f, 12.0f};
434 432
435 if ( T < 10.0 ) { 433 if ( T < 10.0 ) {
436 ionloss = ((G4double)(a[iMolecula][0])) 434 ionloss = ((G4double)(a[iMolecula][0])) * std::sqrt(T) ;
437 435
438 } else if ( T < 10000.0 ) { 436 } else if ( T < 10000.0 ) {
439 G4double x1 = (G4double)(a[iMolecula][1] 437 G4double x1 = (G4double)(a[iMolecula][1]);
440 G4double x2 = (G4double)(a[iMolecula][2] 438 G4double x2 = (G4double)(a[iMolecula][2]);
441 G4double x3 = (G4double)(a[iMolecula][3] 439 G4double x3 = (G4double)(a[iMolecula][3]);
442 G4double x4 = (G4double)(a[iMolecula][4] 440 G4double x4 = (G4double)(a[iMolecula][4]);
443 G4double slow = x1 * G4Exp(G4Log(T)* 0. 441 G4double slow = x1 * G4Exp(G4Log(T)* 0.45);
444 G4double shigh = G4Log( 1.0 + x3/T + x4 442 G4double shigh = G4Log( 1.0 + x3/T + x4*T ) * x2/T;
445 ionloss = slow*shigh / (slow + shigh) ; 443 ionloss = slow*shigh / (slow + shigh) ;
446 } 444 }
447 445
448 ionloss = std::max(ionloss, 0.0); 446 ionloss = std::max(ionloss, 0.0);
449 if ( 10 == iMolecula ) { 447 if ( 10 == iMolecula ) {
450 static const G4double invLog10 = 1.0/G4L 448 static const G4double invLog10 = 1.0/G4Log(10.);
451 449
452 if (T < 100.0) { 450 if (T < 100.0) {
453 ionloss *= (1.0+0.023+0.0066*G4Log(T)* 451 ionloss *= (1.0+0.023+0.0066*G4Log(T)*invLog10);
454 } 452 }
455 else if (T < 700.0) { 453 else if (T < 700.0) {
456 ionloss *=(1.0+0.089-0.0248*G4Log(T-99 454 ionloss *=(1.0+0.089-0.0248*G4Log(T-99.)*invLog10);
457 } 455 }
458 else if (T < 10000.0) { 456 else if (T < 10000.0) {
459 ionloss *=(1.0+0.089-0.0248*G4Log(700. 457 ionloss *=(1.0+0.089-0.0248*G4Log(700.-99.)*invLog10);
460 } 458 }
461 } 459 }
462 ionloss /= (G4double)atomicWeight[iMolecul 460 ionloss /= (G4double)atomicWeight[iMolecula];
463 461
464 // pure material (normally not the case for 462 // pure material (normally not the case for this function)
465 } else if(1 == (material->GetNumberOfElement 463 } else if(1 == (material->GetNumberOfElements())) {
466 G4double z = material->GetZ() ; 464 G4double z = material->GetZ() ;
467 ionloss = ElectronicStoppingPower( z, kine 465 ionloss = ElectronicStoppingPower( z, kineticEnergy ) ;
468 } 466 }
469 467
470 return ionloss; 468 return ionloss;
471 } 469 }
472 470
473 //....oooOO0OOooo........oooOO0OOooo........oo 471 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
474 472
475 G4double G4BraggModel::ElectronicStoppingPower 473 G4double G4BraggModel::ElectronicStoppingPower(G4double z,
476 474 G4double kineticEnergy) const
477 { 475 {
478 G4double ionloss ; 476 G4double ionloss ;
479 G4int i = std::min(std::max(G4lrint(z)-1,0), 477 G4int i = std::min(std::max(G4lrint(z)-1,0),91); // index of atom
480 478
481 // The data and the fit from: 479 // The data and the fit from:
482 // ICRU Report 49, 1993. Ziegler's type of p 480 // ICRU Report 49, 1993. Ziegler's type of parametrisations.
483 // Proton kinetic energy for parametrisation 481 // Proton kinetic energy for parametrisation (keV/amu)
484 482
485 G4double T = kineticEnergy/(keV*protonMassAM 483 G4double T = kineticEnergy/(keV*protonMassAMU) ;
486 484
487 static const G4float a[92][5] = { 485 static const G4float a[92][5] = {
488 {1.254E+0f, 1.440E+0f, 2.426E+2f, 1.200E+4f 486 {1.254E+0f, 1.440E+0f, 2.426E+2f, 1.200E+4f, 1.159E-1f},
489 {1.229E+0f, 1.397E+0f, 4.845E+2f, 5.873E+3f 487 {1.229E+0f, 1.397E+0f, 4.845E+2f, 5.873E+3f, 5.225E-2f},
490 {1.411E+0f, 1.600E+0f, 7.256E+2f, 3.013E+3f 488 {1.411E+0f, 1.600E+0f, 7.256E+2f, 3.013E+3f, 4.578E-2f},
491 {2.248E+0f, 2.590E+0f, 9.660E+2f, 1.538E+2f 489 {2.248E+0f, 2.590E+0f, 9.660E+2f, 1.538E+2f, 3.475E-2f},
492 {2.474E+0f, 2.815E+0f, 1.206E+3f, 1.060E+3f 490 {2.474E+0f, 2.815E+0f, 1.206E+3f, 1.060E+3f, 2.855E-2f},
493 {2.631E+0f, 2.601E+0f, 1.701E+3f, 1.279E+3f 491 {2.631E+0f, 2.601E+0f, 1.701E+3f, 1.279E+3f, 1.638E-2f},
494 {2.954E+0f, 3.350E+0f, 1.683E+3f, 1.900E+3f 492 {2.954E+0f, 3.350E+0f, 1.683E+3f, 1.900E+3f, 2.513E-2f},
495 {2.652E+0f, 3.000E+0f, 1.920E+3f, 2.000E+3f 493 {2.652E+0f, 3.000E+0f, 1.920E+3f, 2.000E+3f, 2.230E-2f},
496 {2.085E+0f, 2.352E+0f, 2.157E+3f, 2.634E+3f 494 {2.085E+0f, 2.352E+0f, 2.157E+3f, 2.634E+3f, 1.816E-2f},
497 {1.951E+0f, 2.199E+0f, 2.393E+3f, 2.699E+3f 495 {1.951E+0f, 2.199E+0f, 2.393E+3f, 2.699E+3f, 1.568E-2f},
498 // Z= 11-20 496 // Z= 11-20
499 {2.542E+0f, 2.869E+0f, 2.628E+3f, 1.854E+3f 497 {2.542E+0f, 2.869E+0f, 2.628E+3f, 1.854E+3f, 1.472E-2f},
500 {3.791E+0f, 4.293E+0f, 2.862E+3f, 1.009E+3f 498 {3.791E+0f, 4.293E+0f, 2.862E+3f, 1.009E+3f, 1.397E-2f},
501 {4.154E+0f, 4.739E+0f, 2.766E+3f, 1.645E+2f 499 {4.154E+0f, 4.739E+0f, 2.766E+3f, 1.645E+2f, 2.023E-2f},
502 {4.914E+0f, 5.598E+0f, 3.193E+3f, 2.327E+2f 500 {4.914E+0f, 5.598E+0f, 3.193E+3f, 2.327E+2f, 1.419E-2f},
503 {3.232E+0f, 3.647E+0f, 3.561E+3f, 1.560E+3f 501 {3.232E+0f, 3.647E+0f, 3.561E+3f, 1.560E+3f, 1.267E-2f},
504 {3.447E+0f, 3.891E+0f, 3.792E+3f, 1.219E+3f 502 {3.447E+0f, 3.891E+0f, 3.792E+3f, 1.219E+3f, 1.211E-2f},
505 {5.301E+0f, 6.008E+0f, 3.969E+3f, 6.451E+2f 503 {5.301E+0f, 6.008E+0f, 3.969E+3f, 6.451E+2f, 1.183E-2f},
506 {5.731E+0f, 6.500E+0f, 4.253E+3f, 5.300E+2f 504 {5.731E+0f, 6.500E+0f, 4.253E+3f, 5.300E+2f, 1.123E-2f},
507 {5.152E+0f, 5.833E+0f, 4.482E+3f, 5.457E+2f 505 {5.152E+0f, 5.833E+0f, 4.482E+3f, 5.457E+2f, 1.129E-2f},
508 {5.521E+0f, 6.252E+0f, 4.710E+3f, 5.533E+2f 506 {5.521E+0f, 6.252E+0f, 4.710E+3f, 5.533E+2f, 1.112E-2f},
509 // Z= 21-30 507 // Z= 21-30
510 {5.201E+0f, 5.884E+0f, 4.938E+3f, 5.609E+2f 508 {5.201E+0f, 5.884E+0f, 4.938E+3f, 5.609E+2f, 9.995E-3f},
511 {4.858E+0f, 5.489E+0f, 5.260E+3f, 6.511E+2f 509 {4.858E+0f, 5.489E+0f, 5.260E+3f, 6.511E+2f, 8.930E-3f},
512 {4.479E+0f, 5.055E+0f, 5.391E+3f, 9.523E+2f 510 {4.479E+0f, 5.055E+0f, 5.391E+3f, 9.523E+2f, 9.117E-3f},
513 {3.983E+0f, 4.489E+0f, 5.616E+3f, 1.336E+3f 511 {3.983E+0f, 4.489E+0f, 5.616E+3f, 1.336E+3f, 8.413E-3f},
514 {3.469E+0f, 3.907E+0f, 5.725E+3f, 1.461E+3f 512 {3.469E+0f, 3.907E+0f, 5.725E+3f, 1.461E+3f, 8.829E-3f},
515 {3.519E+0f, 3.963E+0f, 6.065E+3f, 1.243E+3f 513 {3.519E+0f, 3.963E+0f, 6.065E+3f, 1.243E+3f, 7.782E-3f},
516 {3.140E+0f, 3.535E+0f, 6.288E+3f, 1.372E+3f 514 {3.140E+0f, 3.535E+0f, 6.288E+3f, 1.372E+3f, 7.361E-3f},
517 {3.553E+0f, 4.004E+0f, 6.205E+3f, 5.551E+2f 515 {3.553E+0f, 4.004E+0f, 6.205E+3f, 5.551E+2f, 8.763E-3f},
518 {3.696E+0f, 4.194E+0f, 4.649E+3f, 8.113E+1f 516 {3.696E+0f, 4.194E+0f, 4.649E+3f, 8.113E+1f, 2.242E-2f},
519 {4.210E+0f, 4.750E+0f, 6.953E+3f, 2.952E+2f 517 {4.210E+0f, 4.750E+0f, 6.953E+3f, 2.952E+2f, 6.809E-3f},
520 // Z= 31-40 518 // Z= 31-40
521 {5.041E+0f, 5.697E+0f, 7.173E+3f, 2.026E+2f 519 {5.041E+0f, 5.697E+0f, 7.173E+3f, 2.026E+2f, 6.725E-3f},
522 {5.554E+0f, 6.300E+0f, 6.496E+3f, 1.100E+2f 520 {5.554E+0f, 6.300E+0f, 6.496E+3f, 1.100E+2f, 9.689E-3f},
523 {5.323E+0f, 6.012E+0f, 7.611E+3f, 2.925E+2f 521 {5.323E+0f, 6.012E+0f, 7.611E+3f, 2.925E+2f, 6.447E-3f},
524 {5.874E+0f, 6.656E+0f, 7.395E+3f, 1.175E+2f 522 {5.874E+0f, 6.656E+0f, 7.395E+3f, 1.175E+2f, 7.684E-3f},
525 {6.658E+0f, 7.536E+0f, 7.694E+3f, 2.223E+2f 523 {6.658E+0f, 7.536E+0f, 7.694E+3f, 2.223E+2f, 6.509E-3f},
526 {6.413E+0f, 7.240E+0f, 1.185E+4f, 1.537E+2f 524 {6.413E+0f, 7.240E+0f, 1.185E+4f, 1.537E+2f, 2.880E-3f},
527 {5.694E+0f, 6.429E+0f, 8.478E+3f, 2.929E+2f 525 {5.694E+0f, 6.429E+0f, 8.478E+3f, 2.929E+2f, 6.087E-3f},
528 {6.339E+0f, 7.159E+0f, 8.693E+3f, 3.303E+2f 526 {6.339E+0f, 7.159E+0f, 8.693E+3f, 3.303E+2f, 6.003E-3f},
529 {6.407E+0f, 7.234E+0f, 8.907E+3f, 3.678E+2f 527 {6.407E+0f, 7.234E+0f, 8.907E+3f, 3.678E+2f, 5.889E-3f},
530 {6.734E+0f, 7.603E+0f, 9.120E+3f, 4.052E+2f 528 {6.734E+0f, 7.603E+0f, 9.120E+3f, 4.052E+2f, 5.765E-3f},
531 // Z= 41-50 529 // Z= 41-50
532 {6.901E+0f, 7.791E+0f, 9.333E+3f, 4.427E+2f 530 {6.901E+0f, 7.791E+0f, 9.333E+3f, 4.427E+2f, 5.587E-3f},
533 {6.424E+0f, 7.248E+0f, 9.545E+3f, 4.802E+2f 531 {6.424E+0f, 7.248E+0f, 9.545E+3f, 4.802E+2f, 5.376E-3f},
534 {6.799E+0f, 7.671E+0f, 9.756E+3f, 5.176E+2f 532 {6.799E+0f, 7.671E+0f, 9.756E+3f, 5.176E+2f, 5.315E-3f},
535 {6.109E+0f, 6.887E+0f, 9.966E+3f, 5.551E+2f 533 {6.109E+0f, 6.887E+0f, 9.966E+3f, 5.551E+2f, 5.151E-3f},
536 {5.924E+0f, 6.677E+0f, 1.018E+4f, 5.925E+2f 534 {5.924E+0f, 6.677E+0f, 1.018E+4f, 5.925E+2f, 4.919E-3f},
537 {5.238E+0f, 5.900E+0f, 1.038E+4f, 6.300E+2f 535 {5.238E+0f, 5.900E+0f, 1.038E+4f, 6.300E+2f, 4.758E-3f},
538 // {5.623f, 6.354f, 7160.0f, 337.6f 536 // {5.623f, 6.354f, 7160.0f, 337.6f, 0.013940f}, // Ag Ziegler77
539 {5.345E+0f, 6.038E+0f, 6.790E+3f, 3.978E+2f 537 {5.345E+0f, 6.038E+0f, 6.790E+3f, 3.978E+2f, 1.676E-2f}, // Ag ICRU49
540 {5.814E+0f, 6.554E+0f, 1.080E+4f, 3.555E+2f 538 {5.814E+0f, 6.554E+0f, 1.080E+4f, 3.555E+2f, 4.626E-3f},
541 {6.229E+0f, 7.024E+0f, 1.101E+4f, 3.709E+2f 539 {6.229E+0f, 7.024E+0f, 1.101E+4f, 3.709E+2f, 4.540E-3f},
542 {6.409E+0f, 7.227E+0f, 1.121E+4f, 3.864E+2f 540 {6.409E+0f, 7.227E+0f, 1.121E+4f, 3.864E+2f, 4.474E-3f},
543 // Z= 51-60 541 // Z= 51-60
544 {7.500E+0f, 8.480E+0f, 8.608E+3f, 3.480E+2f 542 {7.500E+0f, 8.480E+0f, 8.608E+3f, 3.480E+2f, 9.074E-3f},
545 {6.979E+0f, 7.871E+0f, 1.162E+4f, 3.924E+2f 543 {6.979E+0f, 7.871E+0f, 1.162E+4f, 3.924E+2f, 4.402E-3f},
546 {7.725E+0f, 8.716E+0f, 1.183E+4f, 3.948E+2f 544 {7.725E+0f, 8.716E+0f, 1.183E+4f, 3.948E+2f, 4.376E-3f},
547 {8.337E+0f, 9.425E+0f, 1.051E+4f, 2.696E+2f 545 {8.337E+0f, 9.425E+0f, 1.051E+4f, 2.696E+2f, 6.206E-3f},
548 {7.287E+0f, 8.218E+0f, 1.223E+4f, 3.997E+2f 546 {7.287E+0f, 8.218E+0f, 1.223E+4f, 3.997E+2f, 4.447E-3f},
549 {7.899E+0f, 8.911E+0f, 1.243E+4f, 4.021E+2f 547 {7.899E+0f, 8.911E+0f, 1.243E+4f, 4.021E+2f, 4.511E-3f},
550 {8.041E+0f, 9.071E+0f, 1.263E+4f, 4.045E+2f 548 {8.041E+0f, 9.071E+0f, 1.263E+4f, 4.045E+2f, 4.540E-3f},
551 {7.488E+0f, 8.444E+0f, 1.283E+4f, 4.069E+2f 549 {7.488E+0f, 8.444E+0f, 1.283E+4f, 4.069E+2f, 4.420E-3f},
552 {7.291E+0f, 8.219E+0f, 1.303E+4f, 4.093E+2f 550 {7.291E+0f, 8.219E+0f, 1.303E+4f, 4.093E+2f, 4.298E-3f},
553 {7.098E+0f, 8.000E+0f, 1.323E+4f, 4.118E+2f 551 {7.098E+0f, 8.000E+0f, 1.323E+4f, 4.118E+2f, 4.182E-3f},
554 // Z= 61-70 552 // Z= 61-70
555 {6.909E+0f, 7.786E+0f, 1.343E+4f, 4.142E+2f 553 {6.909E+0f, 7.786E+0f, 1.343E+4f, 4.142E+2f, 4.058E-3f},
556 {6.728E+0f, 7.580E+0f, 1.362E+4f, 4.166E+2f 554 {6.728E+0f, 7.580E+0f, 1.362E+4f, 4.166E+2f, 3.976E-3f},
557 {6.551E+0f, 7.380E+0f, 1.382E+4f, 4.190E+2f 555 {6.551E+0f, 7.380E+0f, 1.382E+4f, 4.190E+2f, 3.877E-3f},
558 {6.739E+0f, 7.592E+0f, 1.402E+4f, 4.214E+2f 556 {6.739E+0f, 7.592E+0f, 1.402E+4f, 4.214E+2f, 3.863E-3f},
559 {6.212E+0f, 6.996E+0f, 1.421E+4f, 4.239E+2f 557 {6.212E+0f, 6.996E+0f, 1.421E+4f, 4.239E+2f, 3.725E-3f},
560 {5.517E+0f, 6.210E+0f, 1.440E+4f, 4.263E+2f 558 {5.517E+0f, 6.210E+0f, 1.440E+4f, 4.263E+2f, 3.632E-3f},
561 {5.220E+0f, 5.874E+0f, 1.460E+4f, 4.287E+2f 559 {5.220E+0f, 5.874E+0f, 1.460E+4f, 4.287E+2f, 3.498E-3f},
562 {5.071E+0f, 5.706E+0f, 1.479E+4f, 4.330E+2f 560 {5.071E+0f, 5.706E+0f, 1.479E+4f, 4.330E+2f, 3.405E-3f},
563 {4.926E+0f, 5.542E+0f, 1.498E+4f, 4.335E+2f 561 {4.926E+0f, 5.542E+0f, 1.498E+4f, 4.335E+2f, 3.342E-3f},
564 {4.788E+0f, 5.386E+0f, 1.517E+4f, 4.359E+2f 562 {4.788E+0f, 5.386E+0f, 1.517E+4f, 4.359E+2f, 3.292E-3f},
565 // Z= 71-80 563 // Z= 71-80
566 {4.893E+0f, 5.505E+0f, 1.536E+4f, 4.384E+2f 564 {4.893E+0f, 5.505E+0f, 1.536E+4f, 4.384E+2f, 3.243E-3f},
567 {5.028E+0f, 5.657E+0f, 1.555E+4f, 4.408E+2f 565 {5.028E+0f, 5.657E+0f, 1.555E+4f, 4.408E+2f, 3.195E-3f},
568 {4.738E+0f, 5.329E+0f, 1.574E+4f, 4.432E+2f 566 {4.738E+0f, 5.329E+0f, 1.574E+4f, 4.432E+2f, 3.186E-3f},
569 {4.587E+0f, 5.160E+0f, 1.541E+4f, 4.153E+2f 567 {4.587E+0f, 5.160E+0f, 1.541E+4f, 4.153E+2f, 3.406E-3f},
570 {5.201E+0f, 5.851E+0f, 1.612E+4f, 4.416E+2f 568 {5.201E+0f, 5.851E+0f, 1.612E+4f, 4.416E+2f, 3.122E-3f},
571 {5.071E+0f, 5.704E+0f, 1.630E+4f, 4.409E+2f 569 {5.071E+0f, 5.704E+0f, 1.630E+4f, 4.409E+2f, 3.082E-3f},
572 {4.946E+0f, 5.563E+0f, 1.649E+4f, 4.401E+2f 570 {4.946E+0f, 5.563E+0f, 1.649E+4f, 4.401E+2f, 2.965E-3f},
573 {4.477E+0f, 5.034E+0f, 1.667E+4f, 4.393E+2f 571 {4.477E+0f, 5.034E+0f, 1.667E+4f, 4.393E+2f, 2.871E-3f},
574 // {4.856f, 5.460f, 18320.0f, 438.5 572 // {4.856f, 5.460f, 18320.0f, 438.5f, 0.002542f}, //Ziegler77
575 {4.844E+0f, 5.458E+0f, 7.852E+3f, 9.758E+2f 573 {4.844E+0f, 5.458E+0f, 7.852E+3f, 9.758E+2f, 2.077E-2f}, //ICRU49
576 {4.307E+0f, 4.843E+0f, 1.704E+4f, 4.878E+2f 574 {4.307E+0f, 4.843E+0f, 1.704E+4f, 4.878E+2f, 2.882E-3f},
577 // Z= 81-90 575 // Z= 81-90
578 {4.723E+0f, 5.311E+0f, 1.722E+4f, 5.370E+2f 576 {4.723E+0f, 5.311E+0f, 1.722E+4f, 5.370E+2f, 2.913E-3f},
579 {5.319E+0f, 5.982E+0f, 1.740E+4f, 5.863E+2f 577 {5.319E+0f, 5.982E+0f, 1.740E+4f, 5.863E+2f, 2.871E-3f},
580 {5.956E+0f, 6.700E+0f, 1.780E+4f, 6.770E+2f 578 {5.956E+0f, 6.700E+0f, 1.780E+4f, 6.770E+2f, 2.660E-3f},
581 {6.158E+0f, 6.928E+0f, 1.777E+4f, 5.863E+2f 579 {6.158E+0f, 6.928E+0f, 1.777E+4f, 5.863E+2f, 2.812E-3f},
582 {6.203E+0f, 6.979E+0f, 1.795E+4f, 5.863E+2f 580 {6.203E+0f, 6.979E+0f, 1.795E+4f, 5.863E+2f, 2.776E-3f},
583 {6.181E+0f, 6.954E+0f, 1.812E+4f, 5.863E+2f 581 {6.181E+0f, 6.954E+0f, 1.812E+4f, 5.863E+2f, 2.748E-3f},
584 {6.949E+0f, 7.820E+0f, 1.830E+4f, 5.863E+2f 582 {6.949E+0f, 7.820E+0f, 1.830E+4f, 5.863E+2f, 2.737E-3f},
585 {7.506E+0f, 8.448E+0f, 1.848E+4f, 5.863E+2f 583 {7.506E+0f, 8.448E+0f, 1.848E+4f, 5.863E+2f, 2.727E-3f},
586 {7.648E+0f, 8.609E+0f, 1.866E+4f, 5.863E+2f 584 {7.648E+0f, 8.609E+0f, 1.866E+4f, 5.863E+2f, 2.697E-3f},
587 {7.711E+0f, 8.679E+0f, 1.883E+4f, 5.863E+2f 585 {7.711E+0f, 8.679E+0f, 1.883E+4f, 5.863E+2f, 2.641E-3f},
588 // Z= 91-92 586 // Z= 91-92
589 {7.407E+0f, 8.336E+0f, 1.901E+4f, 5.863E+2f 587 {7.407E+0f, 8.336E+0f, 1.901E+4f, 5.863E+2f, 2.603E-3f},
590 {7.290E+0f, 8.204E+0f, 1.918E+4f, 5.863E+2f 588 {7.290E+0f, 8.204E+0f, 1.918E+4f, 5.863E+2f, 2.673E-3f}
591 }; 589 };
592 590
593 G4double fac = 1.0 ; 591 G4double fac = 1.0 ;
594 592
595 // Carbon specific case for E < 40 keV 593 // Carbon specific case for E < 40 keV
596 if ( T < 40.0 && 5 == i) { 594 if ( T < 40.0 && 5 == i) {
597 fac = std::sqrt(T*0.025); 595 fac = std::sqrt(T*0.025);
598 T = 40.0; 596 T = 40.0;
599 597
600 // Free electron gas model 598 // Free electron gas model
601 } else if ( T < 10.0 ) { 599 } else if ( T < 10.0 ) {
602 fac = std::sqrt(T*0.1) ; 600 fac = std::sqrt(T*0.1) ;
603 T = 10.0; 601 T = 10.0;
604 } 602 }
605 603
606 // Main parametrisation 604 // Main parametrisation
607 G4double x1 = (G4double)(a[i][1]); 605 G4double x1 = (G4double)(a[i][1]);
608 G4double x2 = (G4double)(a[i][2]); 606 G4double x2 = (G4double)(a[i][2]);
609 G4double x3 = (G4double)(a[i][3]); 607 G4double x3 = (G4double)(a[i][3]);
610 G4double x4 = (G4double)(a[i][4]); 608 G4double x4 = (G4double)(a[i][4]);
611 G4double slow = x1 * G4Exp(G4Log(T) * 0.45) 609 G4double slow = x1 * G4Exp(G4Log(T) * 0.45);
612 G4double shigh = G4Log( 1.0 + x3/T + x4*T ) 610 G4double shigh = G4Log( 1.0 + x3/T + x4*T ) * x2/T;
613 ionloss = slow*shigh*fac / (slow + shigh); 611 ionloss = slow*shigh*fac / (slow + shigh);
614 612
615 ionloss = std::max(ionloss, 0.0); 613 ionloss = std::max(ionloss, 0.0);
616 614
617 return ionloss; 615 return ionloss;
618 } 616 }
619 617
620 //....oooOO0OOooo........oooOO0OOooo........oo 618 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
621 619
622 G4double G4BraggModel::DEDX(const G4Material* 620 G4double G4BraggModel::DEDX(const G4Material* material, G4double kineticEnergy)
623 { 621 {
624 G4double eloss = 0.0; 622 G4double eloss = 0.0;
625 623
626 // check DB 624 // check DB
627 if(material != currentMaterial) { 625 if(material != currentMaterial) {
628 currentMaterial = material; 626 currentMaterial = material;
629 baseMaterial = material->GetBaseMaterial() 627 baseMaterial = material->GetBaseMaterial()
630 ? material->GetBaseMaterial() : material 628 ? material->GetBaseMaterial() : material;
631 iPSTAR = -1; 629 iPSTAR = -1;
632 iMolecula = -1; 630 iMolecula = -1;
633 iICRU90 = fICRU90 ? fICRU90->GetIndex(base 631 iICRU90 = fICRU90 ? fICRU90->GetIndex(baseMaterial) : -1;
634 632
635 if(iICRU90 < 0) { 633 if(iICRU90 < 0) {
636 iPSTAR = fPSTAR->GetIndex(baseMaterial); 634 iPSTAR = fPSTAR->GetIndex(baseMaterial);
637 if(iPSTAR < 0) { HasMaterial(baseMateria 635 if(iPSTAR < 0) { HasMaterial(baseMaterial); }
638 } 636 }
639 //G4cout << "%%% " <<material->GetName() < 637 //G4cout << "%%% " <<material->GetName() << " iMolecula= "
640 // << iMolecula << " iPSTAR= " << i 638 // << iMolecula << " iPSTAR= " << iPSTAR
641 // << " iICRU90= " << iICRU90<< G4e 639 // << " iICRU90= " << iICRU90<< G4endl;
642 } 640 }
643 641
644 // ICRU90 parameterisation 642 // ICRU90 parameterisation
645 if(iICRU90 >= 0) { 643 if(iICRU90 >= 0) {
646 return fICRU90->GetElectronicDEDXforProton 644 return fICRU90->GetElectronicDEDXforProton(iICRU90, kineticEnergy)
647 *material->GetDensity(); 645 *material->GetDensity();
648 } 646 }
649 // PSTAR parameterisation 647 // PSTAR parameterisation
650 if( iPSTAR >= 0 ) { 648 if( iPSTAR >= 0 ) {
651 return fPSTAR->GetElectronicDEDX(iPSTAR, k 649 return fPSTAR->GetElectronicDEDX(iPSTAR, kineticEnergy)
652 *material->GetDensity(); 650 *material->GetDensity();
653 651
654 } 652 }
655 const std::size_t numberOfElements = materia 653 const std::size_t numberOfElements = material->GetNumberOfElements();
656 const G4double* theAtomicNumDensityVector = 654 const G4double* theAtomicNumDensityVector =
657 material->Get 655 material->GetAtomicNumDensityVector();
658 656
659 657
660 if(iMolecula >= 0) { 658 if(iMolecula >= 0) {
661 eloss = StoppingPower(baseMaterial, kineti 659 eloss = StoppingPower(baseMaterial, kineticEnergy)*
662 material->GetDensity 660 material->GetDensity()/amu;
663 661
664 // Pure material ICRU49 paralmeterisation 662 // Pure material ICRU49 paralmeterisation
665 } else if(1 == numberOfElements) { 663 } else if(1 == numberOfElements) {
666 664
667 G4double z = material->GetZ(); 665 G4double z = material->GetZ();
668 eloss = ElectronicStoppingPower(z, kinetic 666 eloss = ElectronicStoppingPower(z, kineticEnergy)
669 * (material->Ge 667 * (material->GetTotNbOfAtomsPerVolume());
670 668
671 669
672 // Experimental data exist only for kinetic 670 // Experimental data exist only for kinetic energy 125 keV
673 } else if( MolecIsInZiegler1988(material) ) 671 } else if( MolecIsInZiegler1988(material) ) {
674 672
675 // Loop over elements - calculation based 673 // Loop over elements - calculation based on Bragg's rule
676 G4double eloss125 = 0.0 ; 674 G4double eloss125 = 0.0 ;
677 const G4ElementVector* theElementVector = 675 const G4ElementVector* theElementVector =
678 material->GetElemen 676 material->GetElementVector();
679 677
680 // Loop for the elements in the material 678 // Loop for the elements in the material
681 for (std::size_t i=0; i<numberOfElements; 679 for (std::size_t i=0; i<numberOfElements; ++i) {
682 const G4Element* element = (*theElementV 680 const G4Element* element = (*theElementVector)[i] ;
683 G4double z = element->GetZ() ; 681 G4double z = element->GetZ() ;
684 eloss += ElectronicStoppingPower(z,ki 682 eloss += ElectronicStoppingPower(z,kineticEnergy)
685 * theAtomi 683 * theAtomicNumDensityVector[i] ;
686 eloss125 += ElectronicStoppingPower(z,12 684 eloss125 += ElectronicStoppingPower(z,125.0*keV)
687 * theAtomi 685 * theAtomicNumDensityVector[i] ;
688 } 686 }
689 687
690 // Chemical factor is taken into account 688 // Chemical factor is taken into account
691 if (eloss125 > 0.0) { << 689 eloss *= ChemicalFactor(kineticEnergy, eloss125) ;
692 eloss *= ChemicalFactor(kineticEnergy, e <<
693 } <<
694 690
695 // Brugg's rule calculation 691 // Brugg's rule calculation
696 } else { 692 } else {
697 const G4ElementVector* theElementVector = 693 const G4ElementVector* theElementVector =
698 material->GetElemen 694 material->GetElementVector() ;
699 695
700 // loop for the elements in the material 696 // loop for the elements in the material
701 for (std::size_t i=0; i<numberOfElements; 697 for (std::size_t i=0; i<numberOfElements; ++i)
702 { 698 {
703 const G4Element* element = (*theElementV 699 const G4Element* element = (*theElementVector)[i] ;
704 eloss += ElectronicStoppingPower(eleme 700 eloss += ElectronicStoppingPower(element->GetZ(), kineticEnergy)
705 * theAtomic 701 * theAtomicNumDensityVector[i];
706 } 702 }
707 } 703 }
708 return eloss*theZieglerFactor; 704 return eloss*theZieglerFactor;
709 } 705 }
710 706
711 //....oooOO0OOooo........oooOO0OOooo........oo 707 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
712 708
713 G4bool G4BraggModel::MolecIsInZiegler1988(cons 709 G4bool G4BraggModel::MolecIsInZiegler1988(const G4Material* material)
714 { 710 {
715 // The list of molecules from 711 // The list of molecules from
716 // J.F.Ziegler and J.M.Manoyan, The stopping 712 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
717 // Nucl. Inst. & Meth. in Phys. Res. B35 (19 713 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
718 714
719 G4String myFormula = G4String(" ") ; 715 G4String myFormula = G4String(" ") ;
720 const G4String chFormula = material->GetChem 716 const G4String chFormula = material->GetChemicalFormula() ;
721 if (myFormula == chFormula ) { return false; 717 if (myFormula == chFormula ) { return false; }
722 718
723 // There are no evidence for difference of 719 // There are no evidence for difference of stopping power depended on
724 // phase of the compound except for water. 720 // phase of the compound except for water. The stopping power of the
725 // water in gas phase can be predicted usin 721 // water in gas phase can be predicted using Bragg's rule.
726 // 722 //
727 // No chemical factor for water-gas 723 // No chemical factor for water-gas
728 724
729 myFormula = G4String("H_2O") ; 725 myFormula = G4String("H_2O") ;
730 const G4State theState = material->GetState( 726 const G4State theState = material->GetState() ;
731 if( theState == kStateGas && myFormula == ch 727 if( theState == kStateGas && myFormula == chFormula) return false ;
732 728
733 729
734 // The coffecient from Table.4 of Ziegler & 730 // The coffecient from Table.4 of Ziegler & Manoyan
735 static const G4float HeEff = 2.8735f; 731 static const G4float HeEff = 2.8735f;
736 732
737 static const std::size_t numberOfMolecula = 733 static const std::size_t numberOfMolecula = 53;
738 static const G4String nameOfMol[numberOfMole 734 static const G4String nameOfMol[numberOfMolecula] = {
739 "H_2O", "C_2H_4O", "C_3H_6O", "C_ 735 "H_2O", "C_2H_4O", "C_3H_6O", "C_2H_2", "C_H_3OH",
740 "C_2H_5OH", "C_3H_7OH", "C_3H_4", "NH 736 "C_2H_5OH", "C_3H_7OH", "C_3H_4", "NH_3", "C_14H_10",
741 "C_6H_6", "C_4H_10", "C_4H_6", "C_ 737 "C_6H_6", "C_4H_10", "C_4H_6", "C_4H_8O", "CCl_4",
742 "CF_4", "C_6H_8", "C_6H_12", "C_ 738 "CF_4", "C_6H_8", "C_6H_12", "C_6H_10O", "C_6H_10",
743 "C_8H_16", "C_5H_10", "C_5H_8", "C_ 739 "C_8H_16", "C_5H_10", "C_5H_8", "C_3H_6-Cyclopropane","C_2H_4F_2",
744 "C_2H_2F_2", "C_4H_8O_2", "C_2H_6", "C_ 740 "C_2H_2F_2", "C_4H_8O_2", "C_2H_6", "C_2F_6", "C_2H_6O",
745 "C_3H_6O", "C_4H_10O", "C_2H_4", "C_ 741 "C_3H_6O", "C_4H_10O", "C_2H_4", "C_2H_4O", "C_2H_4S",
746 "SH_2", "CH_4", "CCLF_3", "CC 742 "SH_2", "CH_4", "CCLF_3", "CCl_2F_2", "CHCl_2F",
747 "(CH_3)_2S", "N_2O", "C_5H_10O", "C_ 743 "(CH_3)_2S", "N_2O", "C_5H_10O", "C_8H_6", "(CH_2)_N",
748 "(C_3H_6)_N","(C_8H_8)_N", "C_3H_8", "C_ 744 "(C_3H_6)_N","(C_8H_8)_N", "C_3H_8", "C_3H_6-Propylene", "C_3H_6O",
749 "C_3H_6S", "C_4H_4S", "C_7H_8" 745 "C_3H_6S", "C_4H_4S", "C_7H_8"
750 }; 746 };
751 747
752 static const G4float expStopping[numberOfMol 748 static const G4float expStopping[numberOfMolecula] = {
753 66.1f, 190.4f, 258.7f, 42.2f, 141.5f, 749 66.1f, 190.4f, 258.7f, 42.2f, 141.5f,
754 210.9f, 279.6f, 198.8f, 31.0f, 267.5f, 750 210.9f, 279.6f, 198.8f, 31.0f, 267.5f,
755 122.8f, 311.4f, 260.3f, 328.9f, 391.3f, 751 122.8f, 311.4f, 260.3f, 328.9f, 391.3f,
756 206.6f, 374.0f, 422.0f, 432.0f, 398.0f, 752 206.6f, 374.0f, 422.0f, 432.0f, 398.0f,
757 554.0f, 353.0f, 326.0f, 74.6f, 220.5f, 753 554.0f, 353.0f, 326.0f, 74.6f, 220.5f,
758 197.4f, 362.0f, 170.0f, 330.5f, 211.3f, 754 197.4f, 362.0f, 170.0f, 330.5f, 211.3f,
759 262.3f, 349.6f, 51.3f, 187.0f, 236.9f, 755 262.3f, 349.6f, 51.3f, 187.0f, 236.9f,
760 121.9f, 35.8f, 247.0f, 292.6f, 268.0f, 756 121.9f, 35.8f, 247.0f, 292.6f, 268.0f,
761 262.3f, 49.0f, 398.9f, 444.0f, 22.91f, 757 262.3f, 49.0f, 398.9f, 444.0f, 22.91f,
762 68.0f, 155.0f, 84.0f, 74.2f, 254.7f, 758 68.0f, 155.0f, 84.0f, 74.2f, 254.7f,
763 306.8f, 324.4f, 420.0f 759 306.8f, 324.4f, 420.0f
764 } ; 760 } ;
765 761
766 static const G4float expCharge[numberOfMolec 762 static const G4float expCharge[numberOfMolecula] = {
767 HeEff, HeEff, HeEff, 1.0f, HeEff, 763 HeEff, HeEff, HeEff, 1.0f, HeEff,
768 HeEff, HeEff, HeEff, 1.0f, 1.0f, 764 HeEff, HeEff, HeEff, 1.0f, 1.0f,
769 1.0f, HeEff, HeEff, HeEff, HeEff, 765 1.0f, HeEff, HeEff, HeEff, HeEff,
770 HeEff, HeEff, HeEff, HeEff, HeEff, 766 HeEff, HeEff, HeEff, HeEff, HeEff,
771 HeEff, HeEff, HeEff, 1.0f, HeEff, 767 HeEff, HeEff, HeEff, 1.0f, HeEff,
772 HeEff, HeEff, HeEff, HeEff, HeEff, 768 HeEff, HeEff, HeEff, HeEff, HeEff,
773 HeEff, HeEff, 1.0f, HeEff, HeEff, 769 HeEff, HeEff, 1.0f, HeEff, HeEff,
774 HeEff, 1.0f, HeEff, HeEff, HeEff, 770 HeEff, 1.0f, HeEff, HeEff, HeEff,
775 HeEff, 1.0f, HeEff, HeEff, 1.0f, 771 HeEff, 1.0f, HeEff, HeEff, 1.0f,
776 1.0f, 1.0f, 1.0f, 1.0f, HeEff, 772 1.0f, 1.0f, 1.0f, 1.0f, HeEff,
777 HeEff, HeEff, HeEff 773 HeEff, HeEff, HeEff
778 } ; 774 } ;
779 775
780 static const G4int numberOfAtomsPerMolecula[ 776 static const G4int numberOfAtomsPerMolecula[numberOfMolecula] = {
781 3, 7, 10, 4, 6, 777 3, 7, 10, 4, 6,
782 9, 12, 7, 4, 24, 778 9, 12, 7, 4, 24,
783 12,14, 10, 13, 5, 779 12,14, 10, 13, 5,
784 5, 14, 18, 17, 17, 780 5, 14, 18, 17, 17,
785 24,15, 13, 9, 8, 781 24,15, 13, 9, 8,
786 6, 14, 8, 8, 9, 782 6, 14, 8, 8, 9,
787 10,15, 6, 7, 7, 783 10,15, 6, 7, 7,
788 3, 5, 5, 5, 5, 784 3, 5, 5, 5, 5,
789 9, 3, 16, 14, 3, 785 9, 3, 16, 14, 3,
790 9, 16, 11, 9, 10, 786 9, 16, 11, 9, 10,
791 10, 9, 15}; 787 10, 9, 15};
792 788
793 // Search for the compaund in the table 789 // Search for the compaund in the table
794 for (std::size_t i=0; i<numberOfMolecula; ++ 790 for (std::size_t i=0; i<numberOfMolecula; ++i) {
795 if(chFormula == nameOfMol[i]) { 791 if(chFormula == nameOfMol[i]) {
796 expStopPower125 = ((G4double)expStopping 792 expStopPower125 = ((G4double)expStopping[i])
797 * (material->GetTotNbOfAtomsPerVolume( 793 * (material->GetTotNbOfAtomsPerVolume()) /
798 ((G4double)(expCharge[i] * numberOfAto 794 ((G4double)(expCharge[i] * numberOfAtomsPerMolecula[i]));
799 return true; 795 return true;
800 } 796 }
801 } 797 }
802 return false; 798 return false;
803 } 799 }
804 800
805 //....oooOO0OOooo........oooOO0OOooo........oo 801 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
806 802
807 G4double G4BraggModel::ChemicalFactor(G4double 803 G4double G4BraggModel::ChemicalFactor(G4double kineticEnergy,
808 G4double 804 G4double eloss125) const
809 { 805 {
810 // Approximation of Chemical Factor accordin 806 // Approximation of Chemical Factor according to
811 // J.F.Ziegler and J.M.Manoyan, The stopping 807 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
812 // Nucl. Inst. & Meth. in Phys. Res. B35 (19 808 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
813 809
814 static const G4double gamma25 = 1.0 + 25.0* 810 static const G4double gamma25 = 1.0 + 25.0*keV /proton_mass_c2;
815 static const G4double gamma125 = 1.0 + 125.0 811 static const G4double gamma125 = 1.0 + 125.0*keV/proton_mass_c2;
816 static const G4double beta25 = std::sqrt(1 812 static const G4double beta25 = std::sqrt(1.0 - 1.0/(gamma25*gamma25));
817 static const G4double beta125 = std::sqrt(1 813 static const G4double beta125 = std::sqrt(1.0 - 1.0/(gamma125*gamma125));
818 static const G4double f12525 = 1.0 + G4Exp 814 static const G4double f12525 = 1.0 + G4Exp( 1.48*(beta125/beta25 - 7.0) );
819 815
820 G4double gamma = 1.0 + kineticEnergy/proton_ 816 G4double gamma = 1.0 + kineticEnergy/proton_mass_c2;
821 G4double beta = std::sqrt(1.0 - 1.0/(gamma* 817 G4double beta = std::sqrt(1.0 - 1.0/(gamma*gamma));
822 818
823 G4double factor = 1.0 + (expStopPower125/elo 819 G4double factor = 1.0 + (expStopPower125/eloss125 - 1.0) * f12525/
824 (1.0 + G4Exp( 1.48 * ( beta/beta25 - 7. 820 (1.0 + G4Exp( 1.48 * ( beta/beta25 - 7.0 ) ) );
825 821
826 return factor ; 822 return factor ;
827 } 823 }
828 824
829 //....oooOO0OOooo........oooOO0OOooo........oo 825 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
830 826
831 827