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