| 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" <<
74 #include "G4DeltaAngle.hh" 73 #include "G4DeltaAngle.hh"
75 #include "G4ICRU90StoppingData.hh" 74 #include "G4ICRU90StoppingData.hh"
76 #include "G4NistManager.hh" 75 #include "G4NistManager.hh"
77 #include "G4Log.hh" 76 #include "G4Log.hh"
78 #include "G4Exp.hh" 77 #include "G4Exp.hh"
79 #include "G4AutoLock.hh" <<
80 78
81 //....oooOO0OOooo........oooOO0OOooo........oo 79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82 80
83 G4ICRU90StoppingData* G4BraggModel::fICRU90 = << 81 using namespace std;
84 G4PSTARStopping* G4BraggModel::fPSTAR = nullpt <<
85 82
86 namespace << 83 G4PSTARStopping* G4BraggModel::fPSTAR = nullptr;
87 { <<
88 G4Mutex ionMutex = G4MUTEX_INITIALIZER; <<
89 } <<
90 84
91 G4BraggModel::G4BraggModel(const G4ParticleDef 85 G4BraggModel::G4BraggModel(const G4ParticleDefinition* p, const G4String& nam)
92 : G4VEmModel(nam) << 86 : G4VEmModel(nam),
>> 87 particle(nullptr),
>> 88 fICRU90(nullptr),
>> 89 currentMaterial(nullptr),
>> 90 baseMaterial(nullptr),
>> 91 protonMassAMU(1.007276),
>> 92 iMolecula(-1),
>> 93 iPSTAR(-1),
>> 94 iICRU90(-1),
>> 95 isIon(false)
93 { 96 {
94 SetHighEnergyLimit(2.0*CLHEP::MeV); << 97 fParticleChange = nullptr;
>> 98 SetHighEnergyLimit(2.0*MeV);
95 99
96 lowestKinEnergy = 0.25*CLHEP::keV; << 100 lowestKinEnergy = 1.0*keV;
97 theZieglerFactor = CLHEP::eV*CLHEP::cm2*1.0e << 101 theZieglerFactor = eV*cm2*1.0e-15;
98 theElectron = G4Electron::Electron(); 102 theElectron = G4Electron::Electron();
99 expStopPower125 = 0.0; 103 expStopPower125 = 0.0;
100 104
101 corr = G4LossTableManager::Instance()->EmCor 105 corr = G4LossTableManager::Instance()->EmCorrections();
102 if(nullptr != p) { SetParticle(p); } << 106 if(p) { SetParticle(p); }
103 else { SetParticle(theElectron); } << 107 else { SetParticle(theElectron); }
104 } 108 }
105 109
106 //....oooOO0OOooo........oooOO0OOooo........oo 110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107 111
108 G4BraggModel::~G4BraggModel() 112 G4BraggModel::~G4BraggModel()
109 { 113 {
110 if(isFirst) { << 114 if(IsMaster()) { delete fPSTAR; fPSTAR = nullptr; }
111 delete fPSTAR; <<
112 fPSTAR = nullptr; <<
113 } <<
114 } 115 }
115 116
116 //....oooOO0OOooo........oooOO0OOooo........oo 117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117 118
118 void G4BraggModel::Initialise(const G4Particle 119 void G4BraggModel::Initialise(const G4ParticleDefinition* p,
119 const G4DataVect 120 const G4DataVector&)
120 { 121 {
121 if(p != particle) { SetParticle(p); } 122 if(p != particle) { SetParticle(p); }
122 123
123 // always false before the run 124 // always false before the run
124 SetDeexcitationFlag(false); 125 SetDeexcitationFlag(false);
125 126
126 // initialise data only once << 127 if(IsMaster()) {
127 if(nullptr == fPSTAR) { << 128 if(nullptr == fPSTAR) { fPSTAR = new G4PSTARStopping(); }
128 G4AutoLock l(&ionMutex); << 129 if(particle->GetPDGMass() < GeV) { fPSTAR->Initialise(); }
129 if(nullptr == fPSTAR) { << 130 if(G4EmParameters::Instance()->UseICRU90Data()) {
130 isFirst = true; << 131 if(!fICRU90) {
131 fPSTAR = new G4PSTARStopping(); <<
132 if(G4EmParameters::Instance()->UseICRU90 <<
133 fICRU90 = G4NistManager::Instance()->GetICRU 132 fICRU90 = G4NistManager::Instance()->GetICRU90StoppingData();
134 } << 133 } else if(particle->GetPDGMass() < GeV) { fICRU90->Initialise(); }
135 } << 134 }
136 l.unlock(); <<
137 } <<
138 if(isFirst) { <<
139 if(nullptr != fICRU90) { fICRU90->Initiali <<
140 fPSTAR->Initialise(); <<
141 } 135 }
142 136
143 if(nullptr == fParticleChange) { 137 if(nullptr == fParticleChange) {
144 138
145 if(UseAngularGeneratorFlag() && !GetAngula 139 if(UseAngularGeneratorFlag() && !GetAngularDistribution()) {
146 SetAngularDistribution(new G4DeltaAngle( 140 SetAngularDistribution(new G4DeltaAngle());
147 } 141 }
148 G4String pname = particle->GetParticleName 142 G4String pname = particle->GetParticleName();
149 if(particle->GetParticleType() == "nucleus 143 if(particle->GetParticleType() == "nucleus" &&
150 pname != "deuteron" && pname != "triton 144 pname != "deuteron" && pname != "triton" &&
151 pname != "alpha+" && pname != "helium 145 pname != "alpha+" && pname != "helium" &&
152 pname != "hydrogen") { isIon = true; } 146 pname != "hydrogen") { isIon = true; }
153 147
154 fParticleChange = GetParticleChangeForLoss 148 fParticleChange = GetParticleChangeForLoss();
155 } 149 }
156 } 150 }
157 151
158 //....oooOO0OOooo........oooOO0OOooo........oo 152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
159 153
160 void G4BraggModel::SetParticle(const G4Particl <<
161 { <<
162 particle = p; <<
163 mass = particle->GetPDGMass(); <<
164 spin = particle->GetPDGSpin(); <<
165 G4double q = particle->GetPDGCharge()/CLHEP: <<
166 chargeSquare = q*q; <<
167 massRate = mass/CLHEP::proton_mass_c2; <<
168 ratio = CLHEP::electron_mass_c2/mass; <<
169 } <<
170 <<
171 //....oooOO0OOooo........oooOO0OOooo........oo <<
172 <<
173 G4double G4BraggModel::GetChargeSquareRatio(co 154 G4double G4BraggModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
174 co 155 const G4Material* mat,
175 G4 << 156 G4double kineticEnergy)
176 { 157 {
177 // this method is called only for ions 158 // this method is called only for ions
178 chargeSquare = corr->EffectiveChargeSquareRa << 159 G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,kineticEnergy);
179 return chargeSquare; << 160 GetModelOfFluctuations()->SetParticleAndCharge(p, q2);
180 } << 161 return q2*corr->EffectiveChargeCorrection(p,mat,kineticEnergy);
181 <<
182 //....oooOO0OOooo........oooOO0OOooo........oo <<
183 <<
184 G4double G4BraggModel::MinEnergyCut(const G4Pa <<
185 const G4Ma <<
186 { <<
187 return couple->GetMaterial()->GetIonisation( <<
188 } 162 }
189 163
190 //....oooOO0OOooo........oooOO0OOooo........oo 164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
191 165
192 G4double G4BraggModel::GetParticleCharge(const 166 G4double G4BraggModel::GetParticleCharge(const G4ParticleDefinition* p,
193 const 167 const G4Material* mat,
194 G4dou 168 G4double kineticEnergy)
195 { 169 {
196 // this method is called only for ions, so n 170 // this method is called only for ions, so no check if it is an ion
197 return corr->GetParticleCharge(p,mat,kinetic 171 return corr->GetParticleCharge(p,mat,kineticEnergy);
198 } 172 }
199 173
200 //....oooOO0OOooo........oooOO0OOooo........oo 174 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
201 175
202 G4double G4BraggModel::ComputeCrossSectionPerE 176 G4double G4BraggModel::ComputeCrossSectionPerElectron(
203 con 177 const G4ParticleDefinition* p,
204 178 G4double kineticEnergy,
205 << 179 G4double cutEnergy,
206 180 G4double maxKinEnergy)
207 { 181 {
208 G4double cross = 0.0; << 182 G4double cross = 0.0;
209 const G4double tmax = MaxSecondaryEnerg << 183 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
210 const G4double maxEnergy = std::min(tmax, ma << 184 G4double maxEnergy = std::min(tmax,maxKinEnergy);
211 const G4double cutEnergy = std::max(cut, low <<
212 if(cutEnergy < maxEnergy) { 185 if(cutEnergy < maxEnergy) {
213 186
214 const G4double energy = kineticEnergy + m << 187 G4double energy = kineticEnergy + mass;
215 const G4double energy2 = energy*energy; << 188 G4double energy2 = energy*energy;
216 const G4double beta2 = kineticEnergy*(ki << 189 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
217 cross = (maxEnergy - cutEnergy)/(cutEnergy 190 cross = (maxEnergy - cutEnergy)/(cutEnergy*maxEnergy)
218 - beta2*G4Log(maxEnergy/cutEnergy)/tmax; 191 - beta2*G4Log(maxEnergy/cutEnergy)/tmax;
219 192
220 if( 0.0 < spin ) { cross += 0.5*(maxEnergy 193 if( 0.0 < spin ) { cross += 0.5*(maxEnergy - cutEnergy)/energy2; }
221 194
222 cross *= CLHEP::twopi_mc2_rcl2*chargeSquar << 195 cross *= twopi_mc2_rcl2*chargeSquare/beta2;
223 } 196 }
224 // G4cout << "BR: e= " << kineticEnergy << << 197 // G4cout << "BR: e= " << kineticEnergy << " tmin= " << cutEnergy
225 // << " tmax= " << tmax << " cross= << 198 // << " tmax= " << tmax << " cross= " << cross << G4endl;
>> 199
226 return cross; 200 return cross;
227 } 201 }
228 202
229 //....oooOO0OOooo........oooOO0OOooo........oo 203 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
230 204
231 G4double 205 G4double
232 G4BraggModel::ComputeCrossSectionPerAtom(const 206 G4BraggModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition* p,
233 G4dou 207 G4double kineticEnergy,
234 G4dou 208 G4double Z, G4double,
235 G4dou 209 G4double cutEnergy,
236 G4dou 210 G4double maxEnergy)
237 { 211 {
238 return 212 return
239 Z*ComputeCrossSectionPerElectron(p,kinetic 213 Z*ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
240 } 214 }
241 215
242 //....oooOO0OOooo........oooOO0OOooo........oo 216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
243 217
244 G4double G4BraggModel::CrossSectionPerVolume(c 218 G4double G4BraggModel::CrossSectionPerVolume(const G4Material* material,
245 c 219 const G4ParticleDefinition* p,
246 G 220 G4double kineticEnergy,
247 G 221 G4double cutEnergy,
248 G 222 G4double maxEnergy)
249 { 223 {
250 return material->GetElectronDensity() 224 return material->GetElectronDensity()
251 *ComputeCrossSectionPerElectron(p,kineticE 225 *ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
252 } 226 }
253 227
254 //....oooOO0OOooo........oooOO0OOooo........oo 228 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
255 229
256 G4double G4BraggModel::ComputeDEDXPerVolume(co 230 G4double G4BraggModel::ComputeDEDXPerVolume(const G4Material* material,
257 co 231 const G4ParticleDefinition* p,
258 G4 << 232 G4double kineticEnergy,
259 G4 << 233 G4double cutEnergy)
260 { 234 {
261 const G4double tmax = MaxSecondaryEnergy(p, << 235 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
262 const G4double tkin = kinEnergy/massRate; << 236 G4double tkin = kineticEnergy/massRate;
263 const G4double cutEnergy = std::max(cut, low <<
264 G4double dedx = 0.0; 237 G4double dedx = 0.0;
265 238
266 // tkin is the scaled energy to proton << 239 if(tkin < lowestKinEnergy) {
267 if (tkin < lowestKinEnergy) { << 240 dedx = DEDX(material, lowestKinEnergy)*sqrt(tkin/lowestKinEnergy);
268 dedx = DEDX(material, lowestKinEnergy)*std <<
269 } else { 241 } else {
270 dedx = DEDX(material, tkin); << 242 dedx = DEDX(material, tkin);
>> 243 }
271 244
272 // correction for low cut value using Beth << 245 if (cutEnergy < tmax) {
273 if (cutEnergy < tmax) { <<
274 const G4double tau = kinEnergy/mass; <<
275 const G4double x = cutEnergy/tmax; <<
276 246
277 dedx += (G4Log(x)*(tau + 1.)*(tau + 1.)/ << 247 G4double tau = kineticEnergy/mass;
278 CLHEP::twopi_mc2_rcl2 * material->GetElectro << 248 G4double gam = tau + 1.0;
279 } << 249 G4double bg2 = tau * (tau+2.0);
>> 250 G4double beta2 = bg2/(gam*gam);
>> 251 G4double x = cutEnergy/tmax;
>> 252
>> 253 dedx += (G4Log(x) + (1.0 - x)*beta2) * twopi_mc2_rcl2
>> 254 * (material->GetElectronDensity())/beta2;
280 } 255 }
>> 256
>> 257 // now compute the total ionization loss
>> 258
281 dedx = std::max(dedx, 0.0) * chargeSquare; 259 dedx = std::max(dedx, 0.0) * chargeSquare;
282 260
283 //G4cout << "E(MeV)= " << tkin/MeV << " dedx 261 //G4cout << "E(MeV)= " << tkin/MeV << " dedx= " << dedx
284 // << " " << material->GetName() << 262 // << " " << material->GetName() << G4endl;
>> 263
285 return dedx; 264 return dedx;
286 } 265 }
287 266
288 //....oooOO0OOooo........oooOO0OOooo........oo 267 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
289 268
290 void G4BraggModel::SampleSecondaries(std::vect << 269 void G4BraggModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
291 const G4M 270 const G4MaterialCutsCouple* couple,
292 const G4D 271 const G4DynamicParticle* dp,
293 G4double << 272 G4double xmin,
294 G4double 273 G4double maxEnergy)
295 { 274 {
296 const G4double tmax = MaxSecondaryKinEnergy( << 275 G4double tmax = MaxSecondaryKinEnergy(dp);
297 const G4double xmax = std::min(tmax, maxEner << 276 G4double xmax = std::min(tmax, maxEnergy);
298 const G4double xmin = std::max(lowestKinEner <<
299 if(xmin >= xmax) { return; } 277 if(xmin >= xmax) { return; }
300 278
301 G4double kineticEnergy = dp->GetKineticEnerg 279 G4double kineticEnergy = dp->GetKineticEnergy();
302 const G4double energy = kineticEnergy + mas << 280 G4double energy = kineticEnergy + mass;
303 const G4double energy2 = energy*energy; << 281 G4double energy2 = energy*energy;
304 const G4double beta2 = kineticEnergy*(kine << 282 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
305 const G4double grej = 1.0; << 283 G4double grej = 1.0;
306 G4double deltaKinEnergy, f; 284 G4double deltaKinEnergy, f;
307 285
308 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra 286 CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
309 G4double rndm[2]; 287 G4double rndm[2];
310 288
311 // sampling follows ... 289 // sampling follows ...
312 do { 290 do {
313 rndmEngineMod->flatArray(2, rndm); 291 rndmEngineMod->flatArray(2, rndm);
314 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - rn 292 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - rndm[0]) + xmax*rndm[0]);
315 293
316 f = 1.0 - beta2*deltaKinEnergy/tmax; 294 f = 1.0 - beta2*deltaKinEnergy/tmax;
317 295
318 if(f > grej) { 296 if(f > grej) {
319 G4cout << "G4BraggModel::SampleSeconda 297 G4cout << "G4BraggModel::SampleSecondary Warning! "
320 << "Majorant " << grej << " < " 298 << "Majorant " << grej << " < "
321 << f << " for e= " << deltaKinE 299 << f << " for e= " << deltaKinEnergy
322 << G4endl; 300 << G4endl;
323 } 301 }
324 302
325 // Loop checking, 03-Aug-2015, Vladimir Iv 303 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
326 } while( grej*rndm[1] >= f ); 304 } while( grej*rndm[1] >= f );
327 305
328 G4ThreeVector deltaDirection; 306 G4ThreeVector deltaDirection;
329 307
330 if(UseAngularGeneratorFlag()) { 308 if(UseAngularGeneratorFlag()) {
331 const G4Material* mat = couple->GetMateri 309 const G4Material* mat = couple->GetMaterial();
332 G4int Z = SelectRandomAtomNumber(mat); 310 G4int Z = SelectRandomAtomNumber(mat);
333 311
334 deltaDirection = 312 deltaDirection =
335 GetAngularDistribution()->SampleDirectio 313 GetAngularDistribution()->SampleDirection(dp, deltaKinEnergy, Z, mat);
336 314
337 } else { 315 } else {
338 316
339 G4double deltaMomentum = 317 G4double deltaMomentum =
340 std::sqrt(deltaKinEnergy * (deltaKinEner << 318 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
341 G4double cost = deltaKinEnergy * (energy + 319 G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
342 (deltaMomentum * dp->GetTotalMomentum()) 320 (deltaMomentum * dp->GetTotalMomentum());
343 if(cost > 1.0) { cost = 1.0; } 321 if(cost > 1.0) { cost = 1.0; }
344 G4double sint = std::sqrt((1.0 - cost)*(1. << 322 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
345 323
346 G4double phi = twopi*rndmEngineMod->flat() 324 G4double phi = twopi*rndmEngineMod->flat();
347 325
348 deltaDirection.set(sint*std::cos(phi),sint << 326 deltaDirection.set(sint*cos(phi),sint*sin(phi), cost) ;
349 deltaDirection.rotateUz(dp->GetMomentumDir 327 deltaDirection.rotateUz(dp->GetMomentumDirection());
350 } 328 }
351 329
352 // create G4DynamicParticle object for delta 330 // create G4DynamicParticle object for delta ray
353 auto delta = new G4DynamicParticle(theElectr << 331 G4DynamicParticle* delta =
>> 332 new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
354 333
355 // Change kinematics of primary particle 334 // Change kinematics of primary particle
356 kineticEnergy -= deltaKinEnergy; 335 kineticEnergy -= deltaKinEnergy;
357 G4ThreeVector finalP = dp->GetMomentum() - d 336 G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum();
358 finalP = finalP.unit(); << 337 finalP = finalP.unit();
359 338
360 fParticleChange->SetProposedKineticEnergy(ki 339 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
361 fParticleChange->SetProposedMomentumDirectio 340 fParticleChange->SetProposedMomentumDirection(finalP);
362 341
363 vdp->push_back(delta); 342 vdp->push_back(delta);
364 } 343 }
365 344
366 //....oooOO0OOooo........oooOO0OOooo........oo 345 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
367 346
368 G4double G4BraggModel::MaxSecondaryEnergy(cons 347 G4double G4BraggModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
369 G4do 348 G4double kinEnergy)
370 { 349 {
371 if(pd != particle) { SetParticle(pd); } 350 if(pd != particle) { SetParticle(pd); }
372 G4double tau = kinEnergy/mass; 351 G4double tau = kinEnergy/mass;
373 G4double tmax = 2.0*electron_mass_c2*tau*(ta 352 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
374 (1. + 2.0*(tau + 1.)*ratio + 353 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
375 return tmax; 354 return tmax;
376 } 355 }
377 356
378 //....oooOO0OOooo........oooOO0OOooo........oo 357 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
379 358
380 void G4BraggModel::HasMaterial(const G4Materia 359 void G4BraggModel::HasMaterial(const G4Material* mat)
381 { 360 {
382 const G4String& chFormula = mat->GetChemical 361 const G4String& chFormula = mat->GetChemicalFormula();
383 if(chFormula.empty()) { return; } 362 if(chFormula.empty()) { return; }
384 363
385 // ICRU Report N49, 1993. Power's model for 364 // ICRU Report N49, 1993. Power's model for H
386 static const G4int numberOfMolecula = 11; << 365 static const size_t numberOfMolecula = 11;
387 static const G4String molName[numberOfMolecu 366 static const G4String molName[numberOfMolecula] = {
388 "Al_2O_3", "CO_2", 367 "Al_2O_3", "CO_2", "CH_4",
389 "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Pol 368 "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Polypropylene", "(C_8H_8)_N",
390 "C_3H_8", "SiO_2", 369 "C_3H_8", "SiO_2", "H_2O",
391 "H_2O-Gas", "Graphite" } ; 370 "H_2O-Gas", "Graphite" } ;
392 371
393 // Search for the material in the table 372 // Search for the material in the table
394 for (G4int i=0; i<numberOfMolecula; ++i) { << 373 for (size_t i=0; i<numberOfMolecula; ++i) {
395 if (chFormula == molName[i]) { 374 if (chFormula == molName[i]) {
396 iMolecula = i; 375 iMolecula = i;
397 return; 376 return;
398 } 377 }
399 } 378 }
400 return; 379 return;
401 } 380 }
402 381
403 //....oooOO0OOooo........oooOO0OOooo........oo 382 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
404 383
405 G4double G4BraggModel::StoppingPower(const G4M 384 G4double G4BraggModel::StoppingPower(const G4Material* material,
406 G4double 385 G4double kineticEnergy)
407 { 386 {
408 G4double ionloss = 0.0 ; 387 G4double ionloss = 0.0 ;
409 388
410 if (iMolecula >= 0) { 389 if (iMolecula >= 0) {
411 390
412 // The data and the fit from: 391 // The data and the fit from:
413 // ICRU Report N49, 1993. Ziegler's model 392 // ICRU Report N49, 1993. Ziegler's model for protons.
414 // Proton kinetic energy for parametrisati 393 // Proton kinetic energy for parametrisation (keV/amu)
415 394
416 G4double T = kineticEnergy/(keV*protonMass 395 G4double T = kineticEnergy/(keV*protonMassAMU) ;
417 396
418 static const G4float a[11][5] = { 397 static const G4float a[11][5] = {
419 {1.187E+1f, 1.343E+1f, 1.069E+4f, 7.723E+2f 398 {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 399 {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 400 {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 401 {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 402 {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 403 {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 404 {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 405 {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 406 {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 407 {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 408 {2.631E+0f, 2.601E+0f, 1.701E+3f, 1.279E+3f, 1.638E-2f} };
430 409
431 static const G4float atomicWeight[11] = { 410 static const G4float atomicWeight[11] = {
432 101.96128f, 44.0098f, 16.0426f, 28.0536f, 411 101.96128f, 44.0098f, 16.0426f, 28.0536f, 42.0804f,
433 104.1512f, 44.665f, 60.0843f, 18.0152f, 18 412 104.1512f, 44.665f, 60.0843f, 18.0152f, 18.0152f, 12.0f};
434 413
435 if ( T < 10.0 ) { 414 if ( T < 10.0 ) {
436 ionloss = ((G4double)(a[iMolecula][0])) << 415 ionloss = ((G4double)(a[iMolecula][0])) * sqrt(T) ;
437 416
438 } else if ( T < 10000.0 ) { 417 } else if ( T < 10000.0 ) {
439 G4double x1 = (G4double)(a[iMolecula][1] 418 G4double x1 = (G4double)(a[iMolecula][1]);
440 G4double x2 = (G4double)(a[iMolecula][2] 419 G4double x2 = (G4double)(a[iMolecula][2]);
441 G4double x3 = (G4double)(a[iMolecula][3] 420 G4double x3 = (G4double)(a[iMolecula][3]);
442 G4double x4 = (G4double)(a[iMolecula][4] 421 G4double x4 = (G4double)(a[iMolecula][4]);
443 G4double slow = x1 * G4Exp(G4Log(T)* 0. 422 G4double slow = x1 * G4Exp(G4Log(T)* 0.45);
444 G4double shigh = G4Log( 1.0 + x3/T + x4 423 G4double shigh = G4Log( 1.0 + x3/T + x4*T ) * x2/T;
445 ionloss = slow*shigh / (slow + shigh) ; 424 ionloss = slow*shigh / (slow + shigh) ;
446 } 425 }
447 426
448 ionloss = std::max(ionloss, 0.0); 427 ionloss = std::max(ionloss, 0.0);
449 if ( 10 == iMolecula ) { 428 if ( 10 == iMolecula ) {
450 static const G4double invLog10 = 1.0/G4L 429 static const G4double invLog10 = 1.0/G4Log(10.);
451 430
452 if (T < 100.0) { 431 if (T < 100.0) {
453 ionloss *= (1.0+0.023+0.0066*G4Log(T)* 432 ionloss *= (1.0+0.023+0.0066*G4Log(T)*invLog10);
454 } 433 }
455 else if (T < 700.0) { 434 else if (T < 700.0) {
456 ionloss *=(1.0+0.089-0.0248*G4Log(T-99 435 ionloss *=(1.0+0.089-0.0248*G4Log(T-99.)*invLog10);
457 } 436 }
458 else if (T < 10000.0) { 437 else if (T < 10000.0) {
459 ionloss *=(1.0+0.089-0.0248*G4Log(700. 438 ionloss *=(1.0+0.089-0.0248*G4Log(700.-99.)*invLog10);
460 } 439 }
461 } 440 }
462 ionloss /= (G4double)atomicWeight[iMolecul 441 ionloss /= (G4double)atomicWeight[iMolecula];
463 442
464 // pure material (normally not the case for 443 // pure material (normally not the case for this function)
465 } else if(1 == (material->GetNumberOfElement 444 } else if(1 == (material->GetNumberOfElements())) {
466 G4double z = material->GetZ() ; 445 G4double z = material->GetZ() ;
467 ionloss = ElectronicStoppingPower( z, kine 446 ionloss = ElectronicStoppingPower( z, kineticEnergy ) ;
468 } 447 }
469 448
470 return ionloss; 449 return ionloss;
471 } 450 }
472 451
473 //....oooOO0OOooo........oooOO0OOooo........oo 452 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
474 453
475 G4double G4BraggModel::ElectronicStoppingPower 454 G4double G4BraggModel::ElectronicStoppingPower(G4double z,
476 455 G4double kineticEnergy) const
477 { 456 {
478 G4double ionloss ; 457 G4double ionloss ;
479 G4int i = std::min(std::max(G4lrint(z)-1,0), 458 G4int i = std::min(std::max(G4lrint(z)-1,0),91); // index of atom
480 459
481 // The data and the fit from: 460 // The data and the fit from:
482 // ICRU Report 49, 1993. Ziegler's type of p 461 // ICRU Report 49, 1993. Ziegler's type of parametrisations.
483 // Proton kinetic energy for parametrisation 462 // Proton kinetic energy for parametrisation (keV/amu)
484 463
485 G4double T = kineticEnergy/(keV*protonMassAM 464 G4double T = kineticEnergy/(keV*protonMassAMU) ;
486 465
487 static const G4float a[92][5] = { 466 static const G4float a[92][5] = {
488 {1.254E+0f, 1.440E+0f, 2.426E+2f, 1.200E+4f 467 {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 468 {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 469 {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 470 {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 471 {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 472 {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 473 {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 474 {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 475 {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 476 {1.951E+0f, 2.199E+0f, 2.393E+3f, 2.699E+3f, 1.568E-2f},
498 // Z= 11-20 477 // Z= 11-20
499 {2.542E+0f, 2.869E+0f, 2.628E+3f, 1.854E+3f 478 {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 479 {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 480 {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 481 {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 482 {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 483 {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 484 {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 485 {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 486 {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 487 {5.521E+0f, 6.252E+0f, 4.710E+3f, 5.533E+2f, 1.112E-2f},
509 // Z= 21-30 488 // Z= 21-30
510 {5.201E+0f, 5.884E+0f, 4.938E+3f, 5.609E+2f 489 {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 490 {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 491 {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 492 {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 493 {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 494 {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 495 {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 496 {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 497 {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 498 {4.210E+0f, 4.750E+0f, 6.953E+3f, 2.952E+2f, 6.809E-3f},
520 // Z= 31-40 499 // Z= 31-40
521 {5.041E+0f, 5.697E+0f, 7.173E+3f, 2.026E+2f 500 {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 501 {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 502 {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 503 {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 504 {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 505 {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 506 {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 507 {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 508 {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 509 {6.734E+0f, 7.603E+0f, 9.120E+3f, 4.052E+2f, 5.765E-3f},
531 // Z= 41-50 510 // Z= 41-50
532 {6.901E+0f, 7.791E+0f, 9.333E+3f, 4.427E+2f 511 {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 512 {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 513 {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 514 {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 515 {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 516 {5.238E+0f, 5.900E+0f, 1.038E+4f, 6.300E+2f, 4.758E-3f},
538 // {5.623f, 6.354f, 7160.0f, 337.6f 517 // {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 518 {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 519 {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 520 {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 521 {6.409E+0f, 7.227E+0f, 1.121E+4f, 3.864E+2f, 4.474E-3f},
543 // Z= 51-60 522 // Z= 51-60
544 {7.500E+0f, 8.480E+0f, 8.608E+3f, 3.480E+2f 523 {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 524 {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 525 {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 526 {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 527 {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 528 {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 529 {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 530 {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 531 {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 532 {7.098E+0f, 8.000E+0f, 1.323E+4f, 4.118E+2f, 4.182E-3f},
554 // Z= 61-70 533 // Z= 61-70
555 {6.909E+0f, 7.786E+0f, 1.343E+4f, 4.142E+2f 534 {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 535 {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 536 {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 537 {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 538 {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 539 {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 540 {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 541 {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 542 {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 543 {4.788E+0f, 5.386E+0f, 1.517E+4f, 4.359E+2f, 3.292E-3f},
565 // Z= 71-80 544 // Z= 71-80
566 {4.893E+0f, 5.505E+0f, 1.536E+4f, 4.384E+2f 545 {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 546 {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 547 {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 548 {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 549 {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 550 {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 551 {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 552 {4.477E+0f, 5.034E+0f, 1.667E+4f, 4.393E+2f, 2.871E-3f},
574 // {4.856f, 5.460f, 18320.0f, 438.5 553 // {4.856f, 5.460f, 18320.0f, 438.5f, 0.002542f}, //Ziegler77
575 {4.844E+0f, 5.458E+0f, 7.852E+3f, 9.758E+2f 554 {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 555 {4.307E+0f, 4.843E+0f, 1.704E+4f, 4.878E+2f, 2.882E-3f},
577 // Z= 81-90 556 // Z= 81-90
578 {4.723E+0f, 5.311E+0f, 1.722E+4f, 5.370E+2f 557 {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 558 {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 559 {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 560 {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 561 {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 562 {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 563 {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 564 {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 565 {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 566 {7.711E+0f, 8.679E+0f, 1.883E+4f, 5.863E+2f, 2.641E-3f},
588 // Z= 91-92 567 // Z= 91-92
589 {7.407E+0f, 8.336E+0f, 1.901E+4f, 5.863E+2f 568 {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 569 {7.290E+0f, 8.204E+0f, 1.918E+4f, 5.863E+2f, 2.673E-3f}
591 }; 570 };
592 571
593 G4double fac = 1.0 ; 572 G4double fac = 1.0 ;
594 573
595 // Carbon specific case for E < 40 keV 574 // Carbon specific case for E < 40 keV
596 if ( T < 40.0 && 5 == i) { 575 if ( T < 40.0 && 5 == i) {
597 fac = std::sqrt(T*0.025); << 576 fac = sqrt(T*0.025);
598 T = 40.0; 577 T = 40.0;
599 578
600 // Free electron gas model 579 // Free electron gas model
601 } else if ( T < 10.0 ) { 580 } else if ( T < 10.0 ) {
602 fac = std::sqrt(T*0.1) ; << 581 fac = sqrt(T*0.1) ;
603 T = 10.0; 582 T = 10.0;
604 } 583 }
605 584
606 // Main parametrisation 585 // Main parametrisation
607 G4double x1 = (G4double)(a[i][1]); 586 G4double x1 = (G4double)(a[i][1]);
608 G4double x2 = (G4double)(a[i][2]); 587 G4double x2 = (G4double)(a[i][2]);
609 G4double x3 = (G4double)(a[i][3]); 588 G4double x3 = (G4double)(a[i][3]);
610 G4double x4 = (G4double)(a[i][4]); 589 G4double x4 = (G4double)(a[i][4]);
611 G4double slow = x1 * G4Exp(G4Log(T) * 0.45) 590 G4double slow = x1 * G4Exp(G4Log(T) * 0.45);
612 G4double shigh = G4Log( 1.0 + x3/T + x4*T ) 591 G4double shigh = G4Log( 1.0 + x3/T + x4*T ) * x2/T;
613 ionloss = slow*shigh*fac / (slow + shigh); 592 ionloss = slow*shigh*fac / (slow + shigh);
614 593
615 ionloss = std::max(ionloss, 0.0); 594 ionloss = std::max(ionloss, 0.0);
616 595
617 return ionloss; 596 return ionloss;
618 } 597 }
619 598
620 //....oooOO0OOooo........oooOO0OOooo........oo 599 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
621 600
622 G4double G4BraggModel::DEDX(const G4Material* 601 G4double G4BraggModel::DEDX(const G4Material* material, G4double kineticEnergy)
623 { 602 {
624 G4double eloss = 0.0; 603 G4double eloss = 0.0;
625 604
626 // check DB 605 // check DB
627 if(material != currentMaterial) { 606 if(material != currentMaterial) {
628 currentMaterial = material; 607 currentMaterial = material;
629 baseMaterial = material->GetBaseMaterial() 608 baseMaterial = material->GetBaseMaterial()
630 ? material->GetBaseMaterial() : material 609 ? material->GetBaseMaterial() : material;
631 iPSTAR = -1; 610 iPSTAR = -1;
632 iMolecula = -1; 611 iMolecula = -1;
633 iICRU90 = fICRU90 ? fICRU90->GetIndex(base 612 iICRU90 = fICRU90 ? fICRU90->GetIndex(baseMaterial) : -1;
634 613
635 if(iICRU90 < 0) { 614 if(iICRU90 < 0) {
636 iPSTAR = fPSTAR->GetIndex(baseMaterial); 615 iPSTAR = fPSTAR->GetIndex(baseMaterial);
637 if(iPSTAR < 0) { HasMaterial(baseMateria 616 if(iPSTAR < 0) { HasMaterial(baseMaterial); }
638 } 617 }
639 //G4cout << "%%% " <<material->GetName() < 618 //G4cout << "%%% " <<material->GetName() << " iMolecula= "
640 // << iMolecula << " iPSTAR= " << i 619 // << iMolecula << " iPSTAR= " << iPSTAR
641 // << " iICRU90= " << iICRU90<< G4e 620 // << " iICRU90= " << iICRU90<< G4endl;
642 } 621 }
643 622
644 // ICRU90 parameterisation 623 // ICRU90 parameterisation
645 if(iICRU90 >= 0) { 624 if(iICRU90 >= 0) {
646 return fICRU90->GetElectronicDEDXforProton 625 return fICRU90->GetElectronicDEDXforProton(iICRU90, kineticEnergy)
647 *material->GetDensity(); 626 *material->GetDensity();
648 } 627 }
649 // PSTAR parameterisation 628 // PSTAR parameterisation
650 if( iPSTAR >= 0 ) { 629 if( iPSTAR >= 0 ) {
651 return fPSTAR->GetElectronicDEDX(iPSTAR, k 630 return fPSTAR->GetElectronicDEDX(iPSTAR, kineticEnergy)
652 *material->GetDensity(); 631 *material->GetDensity();
653 632
654 } 633 }
655 const std::size_t numberOfElements = materia << 634 const G4int numberOfElements = material->GetNumberOfElements();
656 const G4double* theAtomicNumDensityVector = 635 const G4double* theAtomicNumDensityVector =
657 material->Get 636 material->GetAtomicNumDensityVector();
658 637
659 638
660 if(iMolecula >= 0) { 639 if(iMolecula >= 0) {
661 eloss = StoppingPower(baseMaterial, kineti 640 eloss = StoppingPower(baseMaterial, kineticEnergy)*
662 material->GetDensity 641 material->GetDensity()/amu;
663 642
664 // Pure material ICRU49 paralmeterisation 643 // Pure material ICRU49 paralmeterisation
665 } else if(1 == numberOfElements) { 644 } else if(1 == numberOfElements) {
666 645
667 G4double z = material->GetZ(); 646 G4double z = material->GetZ();
668 eloss = ElectronicStoppingPower(z, kinetic 647 eloss = ElectronicStoppingPower(z, kineticEnergy)
669 * (material->Ge 648 * (material->GetTotNbOfAtomsPerVolume());
670 649
671 650
672 // Experimental data exist only for kinetic 651 // Experimental data exist only for kinetic energy 125 keV
673 } else if( MolecIsInZiegler1988(material) ) 652 } else if( MolecIsInZiegler1988(material) ) {
674 653
675 // Loop over elements - calculation based 654 // Loop over elements - calculation based on Bragg's rule
676 G4double eloss125 = 0.0 ; 655 G4double eloss125 = 0.0 ;
677 const G4ElementVector* theElementVector = 656 const G4ElementVector* theElementVector =
678 material->GetElemen 657 material->GetElementVector();
679 658
680 // Loop for the elements in the material 659 // Loop for the elements in the material
681 for (std::size_t i=0; i<numberOfElements; << 660 for (G4int i=0; i<numberOfElements; ++i) {
682 const G4Element* element = (*theElementV 661 const G4Element* element = (*theElementVector)[i] ;
683 G4double z = element->GetZ() ; 662 G4double z = element->GetZ() ;
684 eloss += ElectronicStoppingPower(z,ki 663 eloss += ElectronicStoppingPower(z,kineticEnergy)
685 * theAtomi 664 * theAtomicNumDensityVector[i] ;
686 eloss125 += ElectronicStoppingPower(z,12 665 eloss125 += ElectronicStoppingPower(z,125.0*keV)
687 * theAtomi 666 * theAtomicNumDensityVector[i] ;
688 } 667 }
689 668
690 // Chemical factor is taken into account 669 // Chemical factor is taken into account
691 if (eloss125 > 0.0) { << 670 eloss *= ChemicalFactor(kineticEnergy, eloss125) ;
692 eloss *= ChemicalFactor(kineticEnergy, e <<
693 } <<
694 671
695 // Brugg's rule calculation 672 // Brugg's rule calculation
696 } else { 673 } else {
697 const G4ElementVector* theElementVector = 674 const G4ElementVector* theElementVector =
698 material->GetElemen 675 material->GetElementVector() ;
699 676
700 // loop for the elements in the material 677 // loop for the elements in the material
701 for (std::size_t i=0; i<numberOfElements; << 678 for (G4int i=0; i<numberOfElements; ++i)
702 { 679 {
703 const G4Element* element = (*theElementV 680 const G4Element* element = (*theElementVector)[i] ;
704 eloss += ElectronicStoppingPower(eleme 681 eloss += ElectronicStoppingPower(element->GetZ(), kineticEnergy)
705 * theAtomic 682 * theAtomicNumDensityVector[i];
706 } 683 }
707 } 684 }
708 return eloss*theZieglerFactor; 685 return eloss*theZieglerFactor;
709 } 686 }
710 687
711 //....oooOO0OOooo........oooOO0OOooo........oo 688 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
712 689
713 G4bool G4BraggModel::MolecIsInZiegler1988(cons 690 G4bool G4BraggModel::MolecIsInZiegler1988(const G4Material* material)
714 { 691 {
715 // The list of molecules from 692 // The list of molecules from
716 // J.F.Ziegler and J.M.Manoyan, The stopping 693 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
717 // Nucl. Inst. & Meth. in Phys. Res. B35 (19 694 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
718 695
719 G4String myFormula = G4String(" ") ; 696 G4String myFormula = G4String(" ") ;
720 const G4String chFormula = material->GetChem 697 const G4String chFormula = material->GetChemicalFormula() ;
721 if (myFormula == chFormula ) { return false; 698 if (myFormula == chFormula ) { return false; }
722 699
723 // There are no evidence for difference of 700 // There are no evidence for difference of stopping power depended on
724 // phase of the compound except for water. 701 // phase of the compound except for water. The stopping power of the
725 // water in gas phase can be predicted usin 702 // water in gas phase can be predicted using Bragg's rule.
726 // 703 //
727 // No chemical factor for water-gas 704 // No chemical factor for water-gas
728 705
729 myFormula = G4String("H_2O") ; 706 myFormula = G4String("H_2O") ;
730 const G4State theState = material->GetState( 707 const G4State theState = material->GetState() ;
731 if( theState == kStateGas && myFormula == ch 708 if( theState == kStateGas && myFormula == chFormula) return false ;
732 709
733 710
734 // The coffecient from Table.4 of Ziegler & 711 // The coffecient from Table.4 of Ziegler & Manoyan
735 static const G4float HeEff = 2.8735f; 712 static const G4float HeEff = 2.8735f;
736 713
737 static const std::size_t numberOfMolecula = << 714 static const size_t numberOfMolecula = 53;
738 static const G4String nameOfMol[numberOfMole 715 static const G4String nameOfMol[numberOfMolecula] = {
739 "H_2O", "C_2H_4O", "C_3H_6O", "C_ 716 "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 717 "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_ 718 "C_6H_6", "C_4H_10", "C_4H_6", "C_4H_8O", "CCl_4",
742 "CF_4", "C_6H_8", "C_6H_12", "C_ 719 "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_ 720 "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_ 721 "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_ 722 "C_3H_6O", "C_4H_10O", "C_2H_4", "C_2H_4O", "C_2H_4S",
746 "SH_2", "CH_4", "CCLF_3", "CC 723 "SH_2", "CH_4", "CCLF_3", "CCl_2F_2", "CHCl_2F",
747 "(CH_3)_2S", "N_2O", "C_5H_10O", "C_ 724 "(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_ 725 "(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" 726 "C_3H_6S", "C_4H_4S", "C_7H_8"
750 }; 727 };
751 728
752 static const G4float expStopping[numberOfMol 729 static const G4float expStopping[numberOfMolecula] = {
753 66.1f, 190.4f, 258.7f, 42.2f, 141.5f, 730 66.1f, 190.4f, 258.7f, 42.2f, 141.5f,
754 210.9f, 279.6f, 198.8f, 31.0f, 267.5f, 731 210.9f, 279.6f, 198.8f, 31.0f, 267.5f,
755 122.8f, 311.4f, 260.3f, 328.9f, 391.3f, 732 122.8f, 311.4f, 260.3f, 328.9f, 391.3f,
756 206.6f, 374.0f, 422.0f, 432.0f, 398.0f, 733 206.6f, 374.0f, 422.0f, 432.0f, 398.0f,
757 554.0f, 353.0f, 326.0f, 74.6f, 220.5f, 734 554.0f, 353.0f, 326.0f, 74.6f, 220.5f,
758 197.4f, 362.0f, 170.0f, 330.5f, 211.3f, 735 197.4f, 362.0f, 170.0f, 330.5f, 211.3f,
759 262.3f, 349.6f, 51.3f, 187.0f, 236.9f, 736 262.3f, 349.6f, 51.3f, 187.0f, 236.9f,
760 121.9f, 35.8f, 247.0f, 292.6f, 268.0f, 737 121.9f, 35.8f, 247.0f, 292.6f, 268.0f,
761 262.3f, 49.0f, 398.9f, 444.0f, 22.91f, 738 262.3f, 49.0f, 398.9f, 444.0f, 22.91f,
762 68.0f, 155.0f, 84.0f, 74.2f, 254.7f, 739 68.0f, 155.0f, 84.0f, 74.2f, 254.7f,
763 306.8f, 324.4f, 420.0f 740 306.8f, 324.4f, 420.0f
764 } ; 741 } ;
765 742
766 static const G4float expCharge[numberOfMolec << 743 static const G4float expCharge[53] = {
767 HeEff, HeEff, HeEff, 1.0f, HeEff, 744 HeEff, HeEff, HeEff, 1.0f, HeEff,
768 HeEff, HeEff, HeEff, 1.0f, 1.0f, 745 HeEff, HeEff, HeEff, 1.0f, 1.0f,
769 1.0f, HeEff, HeEff, HeEff, HeEff, 746 1.0f, HeEff, HeEff, HeEff, HeEff,
770 HeEff, HeEff, HeEff, HeEff, HeEff, 747 HeEff, HeEff, HeEff, HeEff, HeEff,
771 HeEff, HeEff, HeEff, 1.0f, HeEff, 748 HeEff, HeEff, HeEff, 1.0f, HeEff,
772 HeEff, HeEff, HeEff, HeEff, HeEff, 749 HeEff, HeEff, HeEff, HeEff, HeEff,
773 HeEff, HeEff, 1.0f, HeEff, HeEff, 750 HeEff, HeEff, 1.0f, HeEff, HeEff,
774 HeEff, 1.0f, HeEff, HeEff, HeEff, 751 HeEff, 1.0f, HeEff, HeEff, HeEff,
775 HeEff, 1.0f, HeEff, HeEff, 1.0f, 752 HeEff, 1.0f, HeEff, HeEff, 1.0f,
776 1.0f, 1.0f, 1.0f, 1.0f, HeEff, 753 1.0f, 1.0f, 1.0f, 1.0f, HeEff,
777 HeEff, HeEff, HeEff 754 HeEff, HeEff, HeEff
778 } ; 755 } ;
779 756
780 static const G4int numberOfAtomsPerMolecula[ << 757 static const G4int numberOfAtomsPerMolecula[53] = {
781 3, 7, 10, 4, 6, 758 3, 7, 10, 4, 6,
782 9, 12, 7, 4, 24, 759 9, 12, 7, 4, 24,
783 12,14, 10, 13, 5, 760 12,14, 10, 13, 5,
784 5, 14, 18, 17, 17, 761 5, 14, 18, 17, 17,
785 24,15, 13, 9, 8, 762 24,15, 13, 9, 8,
786 6, 14, 8, 8, 9, 763 6, 14, 8, 8, 9,
787 10,15, 6, 7, 7, 764 10,15, 6, 7, 7,
788 3, 5, 5, 5, 5, 765 3, 5, 5, 5, 5,
789 9, 3, 16, 14, 3, 766 9, 3, 16, 14, 3,
790 9, 16, 11, 9, 10, 767 9, 16, 11, 9, 10,
791 10, 9, 15}; 768 10, 9, 15};
792 769
793 // Search for the compaund in the table 770 // Search for the compaund in the table
794 for (std::size_t i=0; i<numberOfMolecula; ++ << 771 for (size_t i=0; i<numberOfMolecula; ++i) {
795 if(chFormula == nameOfMol[i]) { 772 if(chFormula == nameOfMol[i]) {
796 expStopPower125 = ((G4double)expStopping 773 expStopPower125 = ((G4double)expStopping[i])
797 * (material->GetTotNbOfAtomsPerVolume( 774 * (material->GetTotNbOfAtomsPerVolume()) /
798 ((G4double)(expCharge[i] * numberOfAto 775 ((G4double)(expCharge[i] * numberOfAtomsPerMolecula[i]));
799 return true; 776 return true;
800 } 777 }
801 } 778 }
802 return false; 779 return false;
803 } 780 }
804 781
805 //....oooOO0OOooo........oooOO0OOooo........oo 782 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
806 783
807 G4double G4BraggModel::ChemicalFactor(G4double 784 G4double G4BraggModel::ChemicalFactor(G4double kineticEnergy,
808 G4double 785 G4double eloss125) const
809 { 786 {
810 // Approximation of Chemical Factor accordin 787 // Approximation of Chemical Factor according to
811 // J.F.Ziegler and J.M.Manoyan, The stopping 788 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
812 // Nucl. Inst. & Meth. in Phys. Res. B35 (19 789 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
813 790
814 static const G4double gamma25 = 1.0 + 25.0* 791 static const G4double gamma25 = 1.0 + 25.0*keV /proton_mass_c2;
815 static const G4double gamma125 = 1.0 + 125.0 792 static const G4double gamma125 = 1.0 + 125.0*keV/proton_mass_c2;
816 static const G4double beta25 = std::sqrt(1 << 793 static const G4double beta25 = sqrt(1.0 - 1.0/(gamma25*gamma25));
817 static const G4double beta125 = std::sqrt(1 << 794 static const G4double beta125 = sqrt(1.0 - 1.0/(gamma125*gamma125));
818 static const G4double f12525 = 1.0 + G4Exp 795 static const G4double f12525 = 1.0 + G4Exp( 1.48*(beta125/beta25 - 7.0) );
819 796
820 G4double gamma = 1.0 + kineticEnergy/proton_ 797 G4double gamma = 1.0 + kineticEnergy/proton_mass_c2;
821 G4double beta = std::sqrt(1.0 - 1.0/(gamma* << 798 G4double beta = sqrt(1.0 - 1.0/(gamma*gamma));
822 799
823 G4double factor = 1.0 + (expStopPower125/elo 800 G4double factor = 1.0 + (expStopPower125/eloss125 - 1.0) * f12525/
824 (1.0 + G4Exp( 1.48 * ( beta/beta25 - 7. 801 (1.0 + G4Exp( 1.48 * ( beta/beta25 - 7.0 ) ) );
825 802
826 return factor ; 803 return factor ;
827 } 804 }
828 805
829 //....oooOO0OOooo........oooOO0OOooo........oo 806 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
830 807
831 808