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