| 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$
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 G4PSTARStopping* G4BraggModel::fPSTAR = nullptr;
87 { <<
88 G4Mutex ionMutex = G4MUTEX_INITIALIZER; <<
89 } <<
90 83
91 G4BraggModel::G4BraggModel(const G4ParticleDef 84 G4BraggModel::G4BraggModel(const G4ParticleDefinition* p, const G4String& nam)
92 : G4VEmModel(nam) << 85 : G4VEmModel(nam),
>> 86 particle(0),
>> 87 currentMaterial(0),
>> 88 protonMassAMU(1.007276),
>> 89 iMolecula(-1),
>> 90 iPSTAR(-1),
>> 91 isIon(false)
93 { 92 {
94 SetHighEnergyLimit(2.0*CLHEP::MeV); << 93 fParticleChange = nullptr;
>> 94 SetHighEnergyLimit(2.0*MeV);
95 95
96 lowestKinEnergy = 0.25*CLHEP::keV; << 96 lowestKinEnergy = 1.0*keV;
97 theZieglerFactor = CLHEP::eV*CLHEP::cm2*1.0e << 97 theZieglerFactor = eV*cm2*1.0e-15;
98 theElectron = G4Electron::Electron(); 98 theElectron = G4Electron::Electron();
99 expStopPower125 = 0.0; 99 expStopPower125 = 0.0;
100 100
101 corr = G4LossTableManager::Instance()->EmCor 101 corr = G4LossTableManager::Instance()->EmCorrections();
102 if(nullptr != p) { SetParticle(p); } << 102 if(p) { SetParticle(p); }
103 else { SetParticle(theElectron); } << 103 else { SetParticle(theElectron); }
104 } 104 }
105 105
106 //....oooOO0OOooo........oooOO0OOooo........oo 106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107 107
108 G4BraggModel::~G4BraggModel() 108 G4BraggModel::~G4BraggModel()
109 { 109 {
110 if(isFirst) { << 110 if(IsMaster()) { delete fPSTAR; fPSTAR = nullptr; }
111 delete fPSTAR; <<
112 fPSTAR = nullptr; <<
113 } <<
114 } 111 }
115 112
116 //....oooOO0OOooo........oooOO0OOooo........oo 113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117 114
118 void G4BraggModel::Initialise(const G4Particle 115 void G4BraggModel::Initialise(const G4ParticleDefinition* p,
119 const G4DataVect 116 const G4DataVector&)
120 { 117 {
121 if(p != particle) { SetParticle(p); } 118 if(p != particle) { SetParticle(p); }
122 119
123 // always false before the run 120 // always false before the run
124 SetDeexcitationFlag(false); 121 SetDeexcitationFlag(false);
125 122
126 // initialise data only once << 123 if(IsMaster()) {
127 if(nullptr == fPSTAR) { << 124 if(nullptr == fPSTAR) { fPSTAR = new G4PSTARStopping(); }
128 G4AutoLock l(&ionMutex); << 125 if(particle->GetPDGMass() < GeV) { fPSTAR->Initialise(); }
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 } 126 }
142 127
143 if(nullptr == fParticleChange) { 128 if(nullptr == fParticleChange) {
144 129
145 if(UseAngularGeneratorFlag() && !GetAngula 130 if(UseAngularGeneratorFlag() && !GetAngularDistribution()) {
146 SetAngularDistribution(new G4DeltaAngle( 131 SetAngularDistribution(new G4DeltaAngle());
147 } 132 }
148 G4String pname = particle->GetParticleName 133 G4String pname = particle->GetParticleName();
149 if(particle->GetParticleType() == "nucleus 134 if(particle->GetParticleType() == "nucleus" &&
150 pname != "deuteron" && pname != "triton 135 pname != "deuteron" && pname != "triton" &&
151 pname != "alpha+" && pname != "helium 136 pname != "alpha+" && pname != "helium" &&
152 pname != "hydrogen") { isIon = true; } 137 pname != "hydrogen") { isIon = true; }
153 138
154 fParticleChange = GetParticleChangeForLoss 139 fParticleChange = GetParticleChangeForLoss();
155 } 140 }
156 } 141 }
157 142
158 //....oooOO0OOooo........oooOO0OOooo........oo 143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
159 144
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 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 }
>> 251
>> 252 // now compute the total ionization loss
>> 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;
>> 258
285 return dedx; 259 return dedx;
286 } 260 }
287 261
288 //....oooOO0OOooo........oooOO0OOooo........oo 262 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
289 263
290 void G4BraggModel::SampleSecondaries(std::vect << 264 void G4BraggModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
291 const G4M 265 const G4MaterialCutsCouple* couple,
292 const G4D 266 const G4DynamicParticle* dp,
293 G4double << 267 G4double xmin,
294 G4double 268 G4double maxEnergy)
295 { 269 {
296 const G4double tmax = MaxSecondaryKinEnergy( << 270 G4double tmax = MaxSecondaryKinEnergy(dp);
297 const G4double xmax = std::min(tmax, maxEner << 271 G4double xmax = std::min(tmax, maxEnergy);
298 const G4double xmin = std::max(lowestKinEner <<
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 G4double energy = kineticEnergy + mass;
303 const G4double energy2 = energy*energy; << 276 G4double energy2 = energy*energy;
304 const G4double beta2 = kineticEnergy*(kine << 277 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
305 const G4double grej = 1.0; << 278 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 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 = 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*cos(phi),sint*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 G4DynamicParticle* delta =
>> 327 new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
354 328
355 // Change kinematics of primary particle 329 // Change kinematics of primary particle
356 kineticEnergy -= deltaKinEnergy; 330 kineticEnergy -= deltaKinEnergy;
357 G4ThreeVector finalP = dp->GetMomentum() - d 331 G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum();
358 finalP = finalP.unit(); << 332 finalP = finalP.unit();
359 333
360 fParticleChange->SetProposedKineticEnergy(ki 334 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
361 fParticleChange->SetProposedMomentumDirectio 335 fParticleChange->SetProposedMomentumDirection(finalP);
362 336
363 vdp->push_back(delta); 337 vdp->push_back(delta);
364 } 338 }
365 339
366 //....oooOO0OOooo........oooOO0OOooo........oo 340 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
367 341
368 G4double G4BraggModel::MaxSecondaryEnergy(cons 342 G4double G4BraggModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
369 G4do 343 G4double kinEnergy)
370 { 344 {
371 if(pd != particle) { SetParticle(pd); } 345 if(pd != particle) { SetParticle(pd); }
372 G4double tau = kinEnergy/mass; 346 G4double tau = kinEnergy/mass;
373 G4double tmax = 2.0*electron_mass_c2*tau*(ta 347 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
374 (1. + 2.0*(tau + 1.)*ratio + 348 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
375 return tmax; 349 return tmax;
376 } 350 }
377 351
378 //....oooOO0OOooo........oooOO0OOooo........oo 352 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
379 353
380 void G4BraggModel::HasMaterial(const G4Materia << 354 G4bool G4BraggModel::HasMaterial(const G4Material*)
381 { 355 {
382 const G4String& chFormula = mat->GetChemical << 356 return false;
383 if(chFormula.empty()) { return; } << 357 /*
>> 358 G4String chFormula = material->GetChemicalFormula();
>> 359 if("" == chFormula) { return false; }
384 360
385 // ICRU Report N49, 1993. Power's model for 361 // ICRU Report N49, 1993. Power's model for H
386 static const G4int numberOfMolecula = 11; << 362 static const size_t numberOfMolecula = 11;
387 static const G4String molName[numberOfMolecu 363 static const G4String molName[numberOfMolecula] = {
388 "Al_2O_3", "CO_2", 364 "Al_2O_3", "CO_2", "CH_4",
389 "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Pol 365 "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Polypropylene", "(C_8H_8)_N",
390 "C_3H_8", "SiO_2", 366 "C_3H_8", "SiO_2", "H_2O",
391 "H_2O-Gas", "Graphite" } ; 367 "H_2O-Gas", "Graphite" } ;
392 368
393 // Search for the material in the table 369 // Search for the material in the table
394 for (G4int i=0; i<numberOfMolecula; ++i) { << 370 for (size_t i=0; i<numberOfMolecula; ++i) {
395 if (chFormula == molName[i]) { 371 if (chFormula == molName[i]) {
396 iMolecula = i; << 372 iPSTAR = fPSTAR->GetIndex(matName[i]);
397 return; << 373 break;
398 } 374 }
399 } 375 }
400 return; << 376 return (iPSTAR >= 0);
>> 377 */
401 } 378 }
402 379
403 //....oooOO0OOooo........oooOO0OOooo........oo 380 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
404 381
405 G4double G4BraggModel::StoppingPower(const G4M 382 G4double G4BraggModel::StoppingPower(const G4Material* material,
406 G4double << 383 G4double kineticEnergy)
407 { 384 {
408 G4double ionloss = 0.0 ; 385 G4double ionloss = 0.0 ;
409 386
410 if (iMolecula >= 0) { 387 if (iMolecula >= 0) {
411 388
412 // The data and the fit from: 389 // The data and the fit from:
413 // ICRU Report N49, 1993. Ziegler's model 390 // ICRU Report N49, 1993. Ziegler's model for protons.
414 // Proton kinetic energy for parametrisati 391 // Proton kinetic energy for parametrisation (keV/amu)
415 392
416 G4double T = kineticEnergy/(keV*protonMass 393 G4double T = kineticEnergy/(keV*protonMassAMU) ;
417 394
418 static const G4float a[11][5] = { 395 static const G4float a[11][5] = {
419 {1.187E+1f, 1.343E+1f, 1.069E+4f, 7.723E+2f 396 {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 397 {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 398 {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 399 {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 400 {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 401 {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 402 {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 403 {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 404 {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 405 {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 406 {2.631E+0f, 2.601E+0f, 1.701E+3f, 1.279E+3f, 1.638E-2f} };
430 407
431 static const G4float atomicWeight[11] = { 408 static const G4float atomicWeight[11] = {
432 101.96128f, 44.0098f, 16.0426f, 28.0536f, 409 101.96128f, 44.0098f, 16.0426f, 28.0536f, 42.0804f,
433 104.1512f, 44.665f, 60.0843f, 18.0152f, 18 410 104.1512f, 44.665f, 60.0843f, 18.0152f, 18.0152f, 12.0f};
434 411
435 if ( T < 10.0 ) { 412 if ( T < 10.0 ) {
436 ionloss = ((G4double)(a[iMolecula][0])) << 413 ionloss = ((G4double)(a[iMolecula][0])) * sqrt(T) ;
437 414
438 } else if ( T < 10000.0 ) { 415 } else if ( T < 10000.0 ) {
439 G4double x1 = (G4double)(a[iMolecula][1] 416 G4double x1 = (G4double)(a[iMolecula][1]);
440 G4double x2 = (G4double)(a[iMolecula][2] 417 G4double x2 = (G4double)(a[iMolecula][2]);
441 G4double x3 = (G4double)(a[iMolecula][3] 418 G4double x3 = (G4double)(a[iMolecula][3]);
442 G4double x4 = (G4double)(a[iMolecula][4] 419 G4double x4 = (G4double)(a[iMolecula][4]);
443 G4double slow = x1 * G4Exp(G4Log(T)* 0. 420 G4double slow = x1 * G4Exp(G4Log(T)* 0.45);
444 G4double shigh = G4Log( 1.0 + x3/T + x4 421 G4double shigh = G4Log( 1.0 + x3/T + x4*T ) * x2/T;
445 ionloss = slow*shigh / (slow + shigh) ; 422 ionloss = slow*shigh / (slow + shigh) ;
446 } 423 }
447 424
448 ionloss = std::max(ionloss, 0.0); 425 ionloss = std::max(ionloss, 0.0);
449 if ( 10 == iMolecula ) { 426 if ( 10 == iMolecula ) {
450 static const G4double invLog10 = 1.0/G4L 427 static const G4double invLog10 = 1.0/G4Log(10.);
451 428
452 if (T < 100.0) { 429 if (T < 100.0) {
453 ionloss *= (1.0+0.023+0.0066*G4Log(T)* 430 ionloss *= (1.0+0.023+0.0066*G4Log(T)*invLog10);
454 } 431 }
455 else if (T < 700.0) { 432 else if (T < 700.0) {
456 ionloss *=(1.0+0.089-0.0248*G4Log(T-99 433 ionloss *=(1.0+0.089-0.0248*G4Log(T-99.)*invLog10);
457 } 434 }
458 else if (T < 10000.0) { 435 else if (T < 10000.0) {
459 ionloss *=(1.0+0.089-0.0248*G4Log(700. 436 ionloss *=(1.0+0.089-0.0248*G4Log(700.-99.)*invLog10);
460 } 437 }
461 } 438 }
462 ionloss /= (G4double)atomicWeight[iMolecul 439 ionloss /= (G4double)atomicWeight[iMolecula];
463 440
464 // pure material (normally not the case for 441 // pure material (normally not the case for this function)
465 } else if(1 == (material->GetNumberOfElement 442 } else if(1 == (material->GetNumberOfElements())) {
466 G4double z = material->GetZ() ; 443 G4double z = material->GetZ() ;
467 ionloss = ElectronicStoppingPower( z, kine 444 ionloss = ElectronicStoppingPower( z, kineticEnergy ) ;
468 } 445 }
469 446
470 return ionloss; 447 return ionloss;
471 } 448 }
472 449
473 //....oooOO0OOooo........oooOO0OOooo........oo 450 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
474 451
475 G4double G4BraggModel::ElectronicStoppingPower 452 G4double G4BraggModel::ElectronicStoppingPower(G4double z,
476 453 G4double kineticEnergy) const
477 { 454 {
478 G4double ionloss ; 455 G4double ionloss ;
479 G4int i = std::min(std::max(G4lrint(z)-1,0), << 456 G4int i = G4lrint(z)-1 ; // index of atom
480 << 457 if(i < 0) i = 0 ;
>> 458 if(i > 91) i = 91 ;
>> 459
481 // The data and the fit from: 460 // The data and the fit from:
482 // ICRU Report 49, 1993. Ziegler's type of p 461 // ICRU Report 49, 1993. Ziegler's type of parametrisations.
483 // Proton kinetic energy for parametrisation 462 // Proton kinetic energy for parametrisation (keV/amu)
484 463
485 G4double T = kineticEnergy/(keV*protonMassAM 464 G4double T = kineticEnergy/(keV*protonMassAMU) ;
486 465
487 static const G4float a[92][5] = { 466 static const G4float a[92][5] = {
488 {1.254E+0f, 1.440E+0f, 2.426E+2f, 1.200E+4f 467 {1.254E+0f, 1.440E+0f, 2.426E+2f, 1.200E+4f, 1.159E-1f},
489 {1.229E+0f, 1.397E+0f, 4.845E+2f, 5.873E+3f 468 {1.229E+0f, 1.397E+0f, 4.845E+2f, 5.873E+3f, 5.225E-2f},
490 {1.411E+0f, 1.600E+0f, 7.256E+2f, 3.013E+3f 469 {1.411E+0f, 1.600E+0f, 7.256E+2f, 3.013E+3f, 4.578E-2f},
491 {2.248E+0f, 2.590E+0f, 9.660E+2f, 1.538E+2f 470 {2.248E+0f, 2.590E+0f, 9.660E+2f, 1.538E+2f, 3.475E-2f},
492 {2.474E+0f, 2.815E+0f, 1.206E+3f, 1.060E+3f 471 {2.474E+0f, 2.815E+0f, 1.206E+3f, 1.060E+3f, 2.855E-2f},
493 {2.631E+0f, 2.601E+0f, 1.701E+3f, 1.279E+3f 472 {2.631E+0f, 2.601E+0f, 1.701E+3f, 1.279E+3f, 1.638E-2f},
494 {2.954E+0f, 3.350E+0f, 1.683E+3f, 1.900E+3f 473 {2.954E+0f, 3.350E+0f, 1.683E+3f, 1.900E+3f, 2.513E-2f},
495 {2.652E+0f, 3.000E+0f, 1.920E+3f, 2.000E+3f 474 {2.652E+0f, 3.000E+0f, 1.920E+3f, 2.000E+3f, 2.230E-2f},
496 {2.085E+0f, 2.352E+0f, 2.157E+3f, 2.634E+3f 475 {2.085E+0f, 2.352E+0f, 2.157E+3f, 2.634E+3f, 1.816E-2f},
497 {1.951E+0f, 2.199E+0f, 2.393E+3f, 2.699E+3f 476 {1.951E+0f, 2.199E+0f, 2.393E+3f, 2.699E+3f, 1.568E-2f},
498 // Z= 11-20 477 // Z= 11-20
499 {2.542E+0f, 2.869E+0f, 2.628E+3f, 1.854E+3f 478 {2.542E+0f, 2.869E+0f, 2.628E+3f, 1.854E+3f, 1.472E-2f},
500 {3.791E+0f, 4.293E+0f, 2.862E+3f, 1.009E+3f 479 {3.791E+0f, 4.293E+0f, 2.862E+3f, 1.009E+3f, 1.397E-2f},
501 {4.154E+0f, 4.739E+0f, 2.766E+3f, 1.645E+2f 480 {4.154E+0f, 4.739E+0f, 2.766E+3f, 1.645E+2f, 2.023E-2f},
502 {4.914E+0f, 5.598E+0f, 3.193E+3f, 2.327E+2f 481 {4.914E+0f, 5.598E+0f, 3.193E+3f, 2.327E+2f, 1.419E-2f},
503 {3.232E+0f, 3.647E+0f, 3.561E+3f, 1.560E+3f 482 {3.232E+0f, 3.647E+0f, 3.561E+3f, 1.560E+3f, 1.267E-2f},
504 {3.447E+0f, 3.891E+0f, 3.792E+3f, 1.219E+3f 483 {3.447E+0f, 3.891E+0f, 3.792E+3f, 1.219E+3f, 1.211E-2f},
505 {5.301E+0f, 6.008E+0f, 3.969E+3f, 6.451E+2f 484 {5.301E+0f, 6.008E+0f, 3.969E+3f, 6.451E+2f, 1.183E-2f},
506 {5.731E+0f, 6.500E+0f, 4.253E+3f, 5.300E+2f 485 {5.731E+0f, 6.500E+0f, 4.253E+3f, 5.300E+2f, 1.123E-2f},
507 {5.152E+0f, 5.833E+0f, 4.482E+3f, 5.457E+2f 486 {5.152E+0f, 5.833E+0f, 4.482E+3f, 5.457E+2f, 1.129E-2f},
508 {5.521E+0f, 6.252E+0f, 4.710E+3f, 5.533E+2f 487 {5.521E+0f, 6.252E+0f, 4.710E+3f, 5.533E+2f, 1.112E-2f},
509 // Z= 21-30 488 // Z= 21-30
510 {5.201E+0f, 5.884E+0f, 4.938E+3f, 5.609E+2f 489 {5.201E+0f, 5.884E+0f, 4.938E+3f, 5.609E+2f, 9.995E-3f},
511 {4.858E+0f, 5.489E+0f, 5.260E+3f, 6.511E+2f 490 {4.858E+0f, 5.489E+0f, 5.260E+3f, 6.511E+2f, 8.930E-3f},
512 {4.479E+0f, 5.055E+0f, 5.391E+3f, 9.523E+2f 491 {4.479E+0f, 5.055E+0f, 5.391E+3f, 9.523E+2f, 9.117E-3f},
513 {3.983E+0f, 4.489E+0f, 5.616E+3f, 1.336E+3f 492 {3.983E+0f, 4.489E+0f, 5.616E+3f, 1.336E+3f, 8.413E-3f},
514 {3.469E+0f, 3.907E+0f, 5.725E+3f, 1.461E+3f 493 {3.469E+0f, 3.907E+0f, 5.725E+3f, 1.461E+3f, 8.829E-3f},
515 {3.519E+0f, 3.963E+0f, 6.065E+3f, 1.243E+3f 494 {3.519E+0f, 3.963E+0f, 6.065E+3f, 1.243E+3f, 7.782E-3f},
516 {3.140E+0f, 3.535E+0f, 6.288E+3f, 1.372E+3f 495 {3.140E+0f, 3.535E+0f, 6.288E+3f, 1.372E+3f, 7.361E-3f},
517 {3.553E+0f, 4.004E+0f, 6.205E+3f, 5.551E+2f 496 {3.553E+0f, 4.004E+0f, 6.205E+3f, 5.551E+2f, 8.763E-3f},
518 {3.696E+0f, 4.194E+0f, 4.649E+3f, 8.113E+1f 497 {3.696E+0f, 4.194E+0f, 4.649E+3f, 8.113E+1f, 2.242E-2f},
519 {4.210E+0f, 4.750E+0f, 6.953E+3f, 2.952E+2f 498 {4.210E+0f, 4.750E+0f, 6.953E+3f, 2.952E+2f, 6.809E-3f},
520 // Z= 31-40 499 // Z= 31-40
521 {5.041E+0f, 5.697E+0f, 7.173E+3f, 2.026E+2f 500 {5.041E+0f, 5.697E+0f, 7.173E+3f, 2.026E+2f, 6.725E-3f},
522 {5.554E+0f, 6.300E+0f, 6.496E+3f, 1.100E+2f 501 {5.554E+0f, 6.300E+0f, 6.496E+3f, 1.100E+2f, 9.689E-3f},
523 {5.323E+0f, 6.012E+0f, 7.611E+3f, 2.925E+2f 502 {5.323E+0f, 6.012E+0f, 7.611E+3f, 2.925E+2f, 6.447E-3f},
524 {5.874E+0f, 6.656E+0f, 7.395E+3f, 1.175E+2f 503 {5.874E+0f, 6.656E+0f, 7.395E+3f, 1.175E+2f, 7.684E-3f},
525 {6.658E+0f, 7.536E+0f, 7.694E+3f, 2.223E+2f 504 {6.658E+0f, 7.536E+0f, 7.694E+3f, 2.223E+2f, 6.509E-3f},
526 {6.413E+0f, 7.240E+0f, 1.185E+4f, 1.537E+2f 505 {6.413E+0f, 7.240E+0f, 1.185E+4f, 1.537E+2f, 2.880E-3f},
527 {5.694E+0f, 6.429E+0f, 8.478E+3f, 2.929E+2f 506 {5.694E+0f, 6.429E+0f, 8.478E+3f, 2.929E+2f, 6.087E-3f},
528 {6.339E+0f, 7.159E+0f, 8.693E+3f, 3.303E+2f 507 {6.339E+0f, 7.159E+0f, 8.693E+3f, 3.303E+2f, 6.003E-3f},
529 {6.407E+0f, 7.234E+0f, 8.907E+3f, 3.678E+2f 508 {6.407E+0f, 7.234E+0f, 8.907E+3f, 3.678E+2f, 5.889E-3f},
530 {6.734E+0f, 7.603E+0f, 9.120E+3f, 4.052E+2f 509 {6.734E+0f, 7.603E+0f, 9.120E+3f, 4.052E+2f, 5.765E-3f},
531 // Z= 41-50 510 // Z= 41-50
532 {6.901E+0f, 7.791E+0f, 9.333E+3f, 4.427E+2f 511 {6.901E+0f, 7.791E+0f, 9.333E+3f, 4.427E+2f, 5.587E-3f},
533 {6.424E+0f, 7.248E+0f, 9.545E+3f, 4.802E+2f 512 {6.424E+0f, 7.248E+0f, 9.545E+3f, 4.802E+2f, 5.376E-3f},
534 {6.799E+0f, 7.671E+0f, 9.756E+3f, 5.176E+2f 513 {6.799E+0f, 7.671E+0f, 9.756E+3f, 5.176E+2f, 5.315E-3f},
535 {6.109E+0f, 6.887E+0f, 9.966E+3f, 5.551E+2f 514 {6.109E+0f, 6.887E+0f, 9.966E+3f, 5.551E+2f, 5.151E-3f},
536 {5.924E+0f, 6.677E+0f, 1.018E+4f, 5.925E+2f 515 {5.924E+0f, 6.677E+0f, 1.018E+4f, 5.925E+2f, 4.919E-3f},
537 {5.238E+0f, 5.900E+0f, 1.038E+4f, 6.300E+2f 516 {5.238E+0f, 5.900E+0f, 1.038E+4f, 6.300E+2f, 4.758E-3f},
538 // {5.623f, 6.354f, 7160.0f, 337.6f 517 // {5.623f, 6.354f, 7160.0f, 337.6f, 0.013940f}, // Ag Ziegler77
539 {5.345E+0f, 6.038E+0f, 6.790E+3f, 3.978E+2f 518 {5.345E+0f, 6.038E+0f, 6.790E+3f, 3.978E+2f, 1.676E-2f}, // Ag ICRU49
540 {5.814E+0f, 6.554E+0f, 1.080E+4f, 3.555E+2f 519 {5.814E+0f, 6.554E+0f, 1.080E+4f, 3.555E+2f, 4.626E-3f},
541 {6.229E+0f, 7.024E+0f, 1.101E+4f, 3.709E+2f 520 {6.229E+0f, 7.024E+0f, 1.101E+4f, 3.709E+2f, 4.540E-3f},
542 {6.409E+0f, 7.227E+0f, 1.121E+4f, 3.864E+2f 521 {6.409E+0f, 7.227E+0f, 1.121E+4f, 3.864E+2f, 4.474E-3f},
543 // Z= 51-60 522 // Z= 51-60
544 {7.500E+0f, 8.480E+0f, 8.608E+3f, 3.480E+2f 523 {7.500E+0f, 8.480E+0f, 8.608E+3f, 3.480E+2f, 9.074E-3f},
545 {6.979E+0f, 7.871E+0f, 1.162E+4f, 3.924E+2f 524 {6.979E+0f, 7.871E+0f, 1.162E+4f, 3.924E+2f, 4.402E-3f},
546 {7.725E+0f, 8.716E+0f, 1.183E+4f, 3.948E+2f 525 {7.725E+0f, 8.716E+0f, 1.183E+4f, 3.948E+2f, 4.376E-3f},
547 {8.337E+0f, 9.425E+0f, 1.051E+4f, 2.696E+2f 526 {8.337E+0f, 9.425E+0f, 1.051E+4f, 2.696E+2f, 6.206E-3f},
548 {7.287E+0f, 8.218E+0f, 1.223E+4f, 3.997E+2f 527 {7.287E+0f, 8.218E+0f, 1.223E+4f, 3.997E+2f, 4.447E-3f},
549 {7.899E+0f, 8.911E+0f, 1.243E+4f, 4.021E+2f 528 {7.899E+0f, 8.911E+0f, 1.243E+4f, 4.021E+2f, 4.511E-3f},
550 {8.041E+0f, 9.071E+0f, 1.263E+4f, 4.045E+2f 529 {8.041E+0f, 9.071E+0f, 1.263E+4f, 4.045E+2f, 4.540E-3f},
551 {7.488E+0f, 8.444E+0f, 1.283E+4f, 4.069E+2f 530 {7.488E+0f, 8.444E+0f, 1.283E+4f, 4.069E+2f, 4.420E-3f},
552 {7.291E+0f, 8.219E+0f, 1.303E+4f, 4.093E+2f 531 {7.291E+0f, 8.219E+0f, 1.303E+4f, 4.093E+2f, 4.298E-3f},
553 {7.098E+0f, 8.000E+0f, 1.323E+4f, 4.118E+2f 532 {7.098E+0f, 8.000E+0f, 1.323E+4f, 4.118E+2f, 4.182E-3f},
554 // Z= 61-70 533 // Z= 61-70
555 {6.909E+0f, 7.786E+0f, 1.343E+4f, 4.142E+2f 534 {6.909E+0f, 7.786E+0f, 1.343E+4f, 4.142E+2f, 4.058E-3f},
556 {6.728E+0f, 7.580E+0f, 1.362E+4f, 4.166E+2f 535 {6.728E+0f, 7.580E+0f, 1.362E+4f, 4.166E+2f, 3.976E-3f},
557 {6.551E+0f, 7.380E+0f, 1.382E+4f, 4.190E+2f 536 {6.551E+0f, 7.380E+0f, 1.382E+4f, 4.190E+2f, 3.877E-3f},
558 {6.739E+0f, 7.592E+0f, 1.402E+4f, 4.214E+2f 537 {6.739E+0f, 7.592E+0f, 1.402E+4f, 4.214E+2f, 3.863E-3f},
559 {6.212E+0f, 6.996E+0f, 1.421E+4f, 4.239E+2f 538 {6.212E+0f, 6.996E+0f, 1.421E+4f, 4.239E+2f, 3.725E-3f},
560 {5.517E+0f, 6.210E+0f, 1.440E+4f, 4.263E+2f 539 {5.517E+0f, 6.210E+0f, 1.440E+4f, 4.263E+2f, 3.632E-3f},
561 {5.220E+0f, 5.874E+0f, 1.460E+4f, 4.287E+2f 540 {5.220E+0f, 5.874E+0f, 1.460E+4f, 4.287E+2f, 3.498E-3f},
562 {5.071E+0f, 5.706E+0f, 1.479E+4f, 4.330E+2f 541 {5.071E+0f, 5.706E+0f, 1.479E+4f, 4.330E+2f, 3.405E-3f},
563 {4.926E+0f, 5.542E+0f, 1.498E+4f, 4.335E+2f 542 {4.926E+0f, 5.542E+0f, 1.498E+4f, 4.335E+2f, 3.342E-3f},
564 {4.788E+0f, 5.386E+0f, 1.517E+4f, 4.359E+2f 543 {4.788E+0f, 5.386E+0f, 1.517E+4f, 4.359E+2f, 3.292E-3f},
565 // Z= 71-80 544 // Z= 71-80
566 {4.893E+0f, 5.505E+0f, 1.536E+4f, 4.384E+2f 545 {4.893E+0f, 5.505E+0f, 1.536E+4f, 4.384E+2f, 3.243E-3f},
567 {5.028E+0f, 5.657E+0f, 1.555E+4f, 4.408E+2f 546 {5.028E+0f, 5.657E+0f, 1.555E+4f, 4.408E+2f, 3.195E-3f},
568 {4.738E+0f, 5.329E+0f, 1.574E+4f, 4.432E+2f 547 {4.738E+0f, 5.329E+0f, 1.574E+4f, 4.432E+2f, 3.186E-3f},
569 {4.587E+0f, 5.160E+0f, 1.541E+4f, 4.153E+2f 548 {4.587E+0f, 5.160E+0f, 1.541E+4f, 4.153E+2f, 3.406E-3f},
570 {5.201E+0f, 5.851E+0f, 1.612E+4f, 4.416E+2f 549 {5.201E+0f, 5.851E+0f, 1.612E+4f, 4.416E+2f, 3.122E-3f},
571 {5.071E+0f, 5.704E+0f, 1.630E+4f, 4.409E+2f 550 {5.071E+0f, 5.704E+0f, 1.630E+4f, 4.409E+2f, 3.082E-3f},
572 {4.946E+0f, 5.563E+0f, 1.649E+4f, 4.401E+2f 551 {4.946E+0f, 5.563E+0f, 1.649E+4f, 4.401E+2f, 2.965E-3f},
573 {4.477E+0f, 5.034E+0f, 1.667E+4f, 4.393E+2f 552 {4.477E+0f, 5.034E+0f, 1.667E+4f, 4.393E+2f, 2.871E-3f},
574 // {4.856f, 5.460f, 18320.0f, 438.5 553 // {4.856f, 5.460f, 18320.0f, 438.5f, 0.002542f}, //Ziegler77
575 {4.844E+0f, 5.458E+0f, 7.852E+3f, 9.758E+2f 554 {4.844E+0f, 5.458E+0f, 7.852E+3f, 9.758E+2f, 2.077E-2f}, //ICRU49
576 {4.307E+0f, 4.843E+0f, 1.704E+4f, 4.878E+2f 555 {4.307E+0f, 4.843E+0f, 1.704E+4f, 4.878E+2f, 2.882E-3f},
577 // Z= 81-90 556 // Z= 81-90
578 {4.723E+0f, 5.311E+0f, 1.722E+4f, 5.370E+2f 557 {4.723E+0f, 5.311E+0f, 1.722E+4f, 5.370E+2f, 2.913E-3f},
579 {5.319E+0f, 5.982E+0f, 1.740E+4f, 5.863E+2f 558 {5.319E+0f, 5.982E+0f, 1.740E+4f, 5.863E+2f, 2.871E-3f},
580 {5.956E+0f, 6.700E+0f, 1.780E+4f, 6.770E+2f 559 {5.956E+0f, 6.700E+0f, 1.780E+4f, 6.770E+2f, 2.660E-3f},
581 {6.158E+0f, 6.928E+0f, 1.777E+4f, 5.863E+2f 560 {6.158E+0f, 6.928E+0f, 1.777E+4f, 5.863E+2f, 2.812E-3f},
582 {6.203E+0f, 6.979E+0f, 1.795E+4f, 5.863E+2f 561 {6.203E+0f, 6.979E+0f, 1.795E+4f, 5.863E+2f, 2.776E-3f},
583 {6.181E+0f, 6.954E+0f, 1.812E+4f, 5.863E+2f 562 {6.181E+0f, 6.954E+0f, 1.812E+4f, 5.863E+2f, 2.748E-3f},
584 {6.949E+0f, 7.820E+0f, 1.830E+4f, 5.863E+2f 563 {6.949E+0f, 7.820E+0f, 1.830E+4f, 5.863E+2f, 2.737E-3f},
585 {7.506E+0f, 8.448E+0f, 1.848E+4f, 5.863E+2f 564 {7.506E+0f, 8.448E+0f, 1.848E+4f, 5.863E+2f, 2.727E-3f},
586 {7.648E+0f, 8.609E+0f, 1.866E+4f, 5.863E+2f 565 {7.648E+0f, 8.609E+0f, 1.866E+4f, 5.863E+2f, 2.697E-3f},
587 {7.711E+0f, 8.679E+0f, 1.883E+4f, 5.863E+2f 566 {7.711E+0f, 8.679E+0f, 1.883E+4f, 5.863E+2f, 2.641E-3f},
588 // Z= 91-92 567 // Z= 91-92
589 {7.407E+0f, 8.336E+0f, 1.901E+4f, 5.863E+2f 568 {7.407E+0f, 8.336E+0f, 1.901E+4f, 5.863E+2f, 2.603E-3f},
590 {7.290E+0f, 8.204E+0f, 1.918E+4f, 5.863E+2f 569 {7.290E+0f, 8.204E+0f, 1.918E+4f, 5.863E+2f, 2.673E-3f}
591 }; 570 };
592 571
593 G4double fac = 1.0 ; 572 G4double fac = 1.0 ;
594 573
595 // Carbon specific case for E < 40 keV 574 // Carbon specific case for E < 40 keV
596 if ( T < 40.0 && 5 == i) { 575 if ( T < 40.0 && 5 == i) {
597 fac = std::sqrt(T*0.025); << 576 fac = sqrt(T*0.025);
598 T = 40.0; 577 T = 40.0;
599 578
600 // Free electron gas model 579 // Free electron gas model
601 } else if ( T < 10.0 ) { 580 } else if ( T < 10.0 ) {
602 fac = std::sqrt(T*0.1) ; << 581 fac = sqrt(T*0.1) ;
603 T = 10.0; 582 T = 10.0;
604 } 583 }
605 584
606 // Main parametrisation 585 // Main parametrisation
607 G4double x1 = (G4double)(a[i][1]); 586 G4double x1 = (G4double)(a[i][1]);
608 G4double x2 = (G4double)(a[i][2]); 587 G4double x2 = (G4double)(a[i][2]);
609 G4double x3 = (G4double)(a[i][3]); 588 G4double x3 = (G4double)(a[i][3]);
610 G4double x4 = (G4double)(a[i][4]); 589 G4double x4 = (G4double)(a[i][4]);
611 G4double slow = x1 * G4Exp(G4Log(T) * 0.45) 590 G4double slow = x1 * G4Exp(G4Log(T) * 0.45);
612 G4double shigh = G4Log( 1.0 + x3/T + x4*T ) 591 G4double shigh = G4Log( 1.0 + x3/T + x4*T ) * x2/T;
613 ionloss = slow*shigh*fac / (slow + shigh); 592 ionloss = slow*shigh*fac / (slow + shigh);
614 593
615 ionloss = std::max(ionloss, 0.0); 594 ionloss = std::max(ionloss, 0.0);
616 595
617 return ionloss; 596 return ionloss;
618 } 597 }
619 598
620 //....oooOO0OOooo........oooOO0OOooo........oo 599 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
621 600
622 G4double G4BraggModel::DEDX(const G4Material* 601 G4double G4BraggModel::DEDX(const G4Material* material, G4double kineticEnergy)
623 { 602 {
624 G4double eloss = 0.0; 603 G4double eloss = 0.0;
625 604
626 // check DB 605 // check DB
627 if(material != currentMaterial) { 606 if(material != currentMaterial) {
628 currentMaterial = material; 607 currentMaterial = material;
629 baseMaterial = material->GetBaseMaterial() <<
630 ? material->GetBaseMaterial() : material <<
631 iPSTAR = -1; 608 iPSTAR = -1;
632 iMolecula = -1; 609 iMolecula = -1;
633 iICRU90 = fICRU90 ? fICRU90->GetIndex(base << 610 if( !HasMaterial(material) ) { iPSTAR = fPSTAR->GetIndex(material); }
634 << 611
635 if(iICRU90 < 0) { <<
636 iPSTAR = fPSTAR->GetIndex(baseMaterial); <<
637 if(iPSTAR < 0) { HasMaterial(baseMateria <<
638 } <<
639 //G4cout << "%%% " <<material->GetName() < 612 //G4cout << "%%% " <<material->GetName() << " iMolecula= "
640 // << iMolecula << " iPSTAR= " << i << 613 // << iMolecula << " iPSTAR= " << iPSTAR << G4endl;
641 // << " iICRU90= " << iICRU90<< G4e <<
642 } <<
643 614
644 // ICRU90 parameterisation <<
645 if(iICRU90 >= 0) { <<
646 return fICRU90->GetElectronicDEDXforProton <<
647 *material->GetDensity(); <<
648 } 615 }
649 // PSTAR parameterisation <<
650 if( iPSTAR >= 0 ) { <<
651 return fPSTAR->GetElectronicDEDX(iPSTAR, k <<
652 *material->GetDensity(); <<
653 616
654 } << 617 const G4int numberOfElements = material->GetNumberOfElements();
655 const std::size_t numberOfElements = materia <<
656 const G4double* theAtomicNumDensityVector = 618 const G4double* theAtomicNumDensityVector =
657 material->Get 619 material->GetAtomicNumDensityVector();
658 620
>> 621 if( iPSTAR >= 0 ) {
>> 622 return
>> 623 fPSTAR->GetElectronicDEDX(iPSTAR, kineticEnergy)*material->GetDensity();
659 624
660 if(iMolecula >= 0) { << 625 } else if(iMolecula >= 0) {
661 eloss = StoppingPower(baseMaterial, kineti << 626
>> 627 eloss = StoppingPower(material, kineticEnergy)*
662 material->GetDensity 628 material->GetDensity()/amu;
663 629
664 // Pure material ICRU49 paralmeterisation << 630 // Pure material ICRU49 paralmeterisation
665 } else if(1 == numberOfElements) { 631 } else if(1 == numberOfElements) {
666 632
667 G4double z = material->GetZ(); 633 G4double z = material->GetZ();
668 eloss = ElectronicStoppingPower(z, kinetic 634 eloss = ElectronicStoppingPower(z, kineticEnergy)
669 * (material->Ge 635 * (material->GetTotNbOfAtomsPerVolume());
670 636
671 637
672 // Experimental data exist only for kinetic 638 // Experimental data exist only for kinetic energy 125 keV
673 } else if( MolecIsInZiegler1988(material) ) 639 } else if( MolecIsInZiegler1988(material) ) {
674 640
675 // Loop over elements - calculation based 641 // Loop over elements - calculation based on Bragg's rule
676 G4double eloss125 = 0.0 ; 642 G4double eloss125 = 0.0 ;
677 const G4ElementVector* theElementVector = 643 const G4ElementVector* theElementVector =
678 material->GetElemen 644 material->GetElementVector();
679 645
680 // Loop for the elements in the material 646 // Loop for the elements in the material
681 for (std::size_t i=0; i<numberOfElements; << 647 for (G4int i=0; i<numberOfElements; ++i) {
682 const G4Element* element = (*theElementV 648 const G4Element* element = (*theElementVector)[i] ;
683 G4double z = element->GetZ() ; 649 G4double z = element->GetZ() ;
684 eloss += ElectronicStoppingPower(z,ki 650 eloss += ElectronicStoppingPower(z,kineticEnergy)
685 * theAtomi 651 * theAtomicNumDensityVector[i] ;
686 eloss125 += ElectronicStoppingPower(z,12 652 eloss125 += ElectronicStoppingPower(z,125.0*keV)
687 * theAtomi 653 * theAtomicNumDensityVector[i] ;
688 } 654 }
689 655
690 // Chemical factor is taken into account 656 // Chemical factor is taken into account
691 if (eloss125 > 0.0) { << 657 eloss *= ChemicalFactor(kineticEnergy, eloss125) ;
692 eloss *= ChemicalFactor(kineticEnergy, e <<
693 } <<
694 658
695 // Brugg's rule calculation 659 // Brugg's rule calculation
696 } else { 660 } else {
697 const G4ElementVector* theElementVector = 661 const G4ElementVector* theElementVector =
698 material->GetElemen 662 material->GetElementVector() ;
699 663
700 // loop for the elements in the material 664 // loop for the elements in the material
701 for (std::size_t i=0; i<numberOfElements; << 665 for (G4int i=0; i<numberOfElements; ++i)
702 { 666 {
703 const G4Element* element = (*theElementV 667 const G4Element* element = (*theElementVector)[i] ;
704 eloss += ElectronicStoppingPower(eleme 668 eloss += ElectronicStoppingPower(element->GetZ(), kineticEnergy)
705 * theAtomic 669 * theAtomicNumDensityVector[i];
706 } 670 }
707 } 671 }
708 return eloss*theZieglerFactor; 672 return eloss*theZieglerFactor;
709 } 673 }
710 674
711 //....oooOO0OOooo........oooOO0OOooo........oo 675 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
712 676
713 G4bool G4BraggModel::MolecIsInZiegler1988(cons 677 G4bool G4BraggModel::MolecIsInZiegler1988(const G4Material* material)
714 { 678 {
715 // The list of molecules from 679 // The list of molecules from
716 // J.F.Ziegler and J.M.Manoyan, The stopping 680 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
717 // Nucl. Inst. & Meth. in Phys. Res. B35 (19 681 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
718 682
719 G4String myFormula = G4String(" ") ; 683 G4String myFormula = G4String(" ") ;
720 const G4String chFormula = material->GetChem 684 const G4String chFormula = material->GetChemicalFormula() ;
721 if (myFormula == chFormula ) { return false; 685 if (myFormula == chFormula ) { return false; }
722 686
723 // There are no evidence for difference of 687 // There are no evidence for difference of stopping power depended on
724 // phase of the compound except for water. 688 // phase of the compound except for water. The stopping power of the
725 // water in gas phase can be predicted usin 689 // water in gas phase can be predicted using Bragg's rule.
726 // 690 //
727 // No chemical factor for water-gas 691 // No chemical factor for water-gas
728 692
729 myFormula = G4String("H_2O") ; 693 myFormula = G4String("H_2O") ;
730 const G4State theState = material->GetState( 694 const G4State theState = material->GetState() ;
731 if( theState == kStateGas && myFormula == ch 695 if( theState == kStateGas && myFormula == chFormula) return false ;
732 696
733 697
734 // The coffecient from Table.4 of Ziegler & 698 // The coffecient from Table.4 of Ziegler & Manoyan
735 static const G4float HeEff = 2.8735f; 699 static const G4float HeEff = 2.8735f;
736 700
737 static const std::size_t numberOfMolecula = << 701 static const size_t numberOfMolecula = 53;
738 static const G4String nameOfMol[numberOfMole 702 static const G4String nameOfMol[numberOfMolecula] = {
739 "H_2O", "C_2H_4O", "C_3H_6O", "C_ 703 "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 704 "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_ 705 "C_6H_6", "C_4H_10", "C_4H_6", "C_4H_8O", "CCl_4",
742 "CF_4", "C_6H_8", "C_6H_12", "C_ 706 "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_ 707 "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_ 708 "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_ 709 "C_3H_6O", "C_4H_10O", "C_2H_4", "C_2H_4O", "C_2H_4S",
746 "SH_2", "CH_4", "CCLF_3", "CC 710 "SH_2", "CH_4", "CCLF_3", "CCl_2F_2", "CHCl_2F",
747 "(CH_3)_2S", "N_2O", "C_5H_10O", "C_ 711 "(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_ 712 "(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" 713 "C_3H_6S", "C_4H_4S", "C_7H_8"
750 }; 714 };
751 715
752 static const G4float expStopping[numberOfMol 716 static const G4float expStopping[numberOfMolecula] = {
753 66.1f, 190.4f, 258.7f, 42.2f, 141.5f, 717 66.1f, 190.4f, 258.7f, 42.2f, 141.5f,
754 210.9f, 279.6f, 198.8f, 31.0f, 267.5f, 718 210.9f, 279.6f, 198.8f, 31.0f, 267.5f,
755 122.8f, 311.4f, 260.3f, 328.9f, 391.3f, 719 122.8f, 311.4f, 260.3f, 328.9f, 391.3f,
756 206.6f, 374.0f, 422.0f, 432.0f, 398.0f, 720 206.6f, 374.0f, 422.0f, 432.0f, 398.0f,
757 554.0f, 353.0f, 326.0f, 74.6f, 220.5f, 721 554.0f, 353.0f, 326.0f, 74.6f, 220.5f,
758 197.4f, 362.0f, 170.0f, 330.5f, 211.3f, 722 197.4f, 362.0f, 170.0f, 330.5f, 211.3f,
759 262.3f, 349.6f, 51.3f, 187.0f, 236.9f, 723 262.3f, 349.6f, 51.3f, 187.0f, 236.9f,
760 121.9f, 35.8f, 247.0f, 292.6f, 268.0f, 724 121.9f, 35.8f, 247.0f, 292.6f, 268.0f,
761 262.3f, 49.0f, 398.9f, 444.0f, 22.91f, 725 262.3f, 49.0f, 398.9f, 444.0f, 22.91f,
762 68.0f, 155.0f, 84.0f, 74.2f, 254.7f, 726 68.0f, 155.0f, 84.0f, 74.2f, 254.7f,
763 306.8f, 324.4f, 420.0f 727 306.8f, 324.4f, 420.0f
764 } ; 728 } ;
765 729
766 static const G4float expCharge[numberOfMolec << 730 static const G4float expCharge[53] = {
767 HeEff, HeEff, HeEff, 1.0f, HeEff, 731 HeEff, HeEff, HeEff, 1.0f, HeEff,
768 HeEff, HeEff, HeEff, 1.0f, 1.0f, 732 HeEff, HeEff, HeEff, 1.0f, 1.0f,
769 1.0f, HeEff, HeEff, HeEff, HeEff, 733 1.0f, HeEff, HeEff, HeEff, HeEff,
770 HeEff, HeEff, HeEff, HeEff, HeEff, 734 HeEff, HeEff, HeEff, HeEff, HeEff,
771 HeEff, HeEff, HeEff, 1.0f, HeEff, 735 HeEff, HeEff, HeEff, 1.0f, HeEff,
772 HeEff, HeEff, HeEff, HeEff, HeEff, 736 HeEff, HeEff, HeEff, HeEff, HeEff,
773 HeEff, HeEff, 1.0f, HeEff, HeEff, 737 HeEff, HeEff, 1.0f, HeEff, HeEff,
774 HeEff, 1.0f, HeEff, HeEff, HeEff, 738 HeEff, 1.0f, HeEff, HeEff, HeEff,
775 HeEff, 1.0f, HeEff, HeEff, 1.0f, 739 HeEff, 1.0f, HeEff, HeEff, 1.0f,
776 1.0f, 1.0f, 1.0f, 1.0f, HeEff, 740 1.0f, 1.0f, 1.0f, 1.0f, HeEff,
777 HeEff, HeEff, HeEff 741 HeEff, HeEff, HeEff
778 } ; 742 } ;
779 743
780 static const G4int numberOfAtomsPerMolecula[ << 744 static const G4int numberOfAtomsPerMolecula[53] = {
781 3, 7, 10, 4, 6, 745 3, 7, 10, 4, 6,
782 9, 12, 7, 4, 24, 746 9, 12, 7, 4, 24,
783 12,14, 10, 13, 5, 747 12,14, 10, 13, 5,
784 5, 14, 18, 17, 17, 748 5, 14, 18, 17, 17,
785 24,15, 13, 9, 8, 749 24,15, 13, 9, 8,
786 6, 14, 8, 8, 9, 750 6, 14, 8, 8, 9,
787 10,15, 6, 7, 7, 751 10,15, 6, 7, 7,
788 3, 5, 5, 5, 5, 752 3, 5, 5, 5, 5,
789 9, 3, 16, 14, 3, 753 9, 3, 16, 14, 3,
790 9, 16, 11, 9, 10, 754 9, 16, 11, 9, 10,
791 10, 9, 15}; 755 10, 9, 15};
792 756
793 // Search for the compaund in the table 757 // Search for the compaund in the table
794 for (std::size_t i=0; i<numberOfMolecula; ++ << 758 for (size_t i=0; i<numberOfMolecula; ++i) {
795 if(chFormula == nameOfMol[i]) { 759 if(chFormula == nameOfMol[i]) {
796 expStopPower125 = ((G4double)expStopping 760 expStopPower125 = ((G4double)expStopping[i])
797 * (material->GetTotNbOfAtomsPerVolume( 761 * (material->GetTotNbOfAtomsPerVolume()) /
798 ((G4double)(expCharge[i] * numberOfAto 762 ((G4double)(expCharge[i] * numberOfAtomsPerMolecula[i]));
799 return true; 763 return true;
800 } 764 }
801 } 765 }
802 return false; 766 return false;
803 } 767 }
804 768
805 //....oooOO0OOooo........oooOO0OOooo........oo 769 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
806 770
807 G4double G4BraggModel::ChemicalFactor(G4double 771 G4double G4BraggModel::ChemicalFactor(G4double kineticEnergy,
808 G4double 772 G4double eloss125) const
809 { 773 {
810 // Approximation of Chemical Factor accordin 774 // Approximation of Chemical Factor according to
811 // J.F.Ziegler and J.M.Manoyan, The stopping 775 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
812 // Nucl. Inst. & Meth. in Phys. Res. B35 (19 776 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
813 777
814 static const G4double gamma25 = 1.0 + 25.0* 778 static const G4double gamma25 = 1.0 + 25.0*keV /proton_mass_c2;
815 static const G4double gamma125 = 1.0 + 125.0 779 static const G4double gamma125 = 1.0 + 125.0*keV/proton_mass_c2;
816 static const G4double beta25 = std::sqrt(1 << 780 static const G4double beta25 = sqrt(1.0 - 1.0/(gamma25*gamma25));
817 static const G4double beta125 = std::sqrt(1 << 781 static const G4double beta125 = sqrt(1.0 - 1.0/(gamma125*gamma125));
818 static const G4double f12525 = 1.0 + G4Exp 782 static const G4double f12525 = 1.0 + G4Exp( 1.48*(beta125/beta25 - 7.0) );
819 783
820 G4double gamma = 1.0 + kineticEnergy/proton_ 784 G4double gamma = 1.0 + kineticEnergy/proton_mass_c2;
821 G4double beta = std::sqrt(1.0 - 1.0/(gamma* << 785 G4double beta = sqrt(1.0 - 1.0/(gamma*gamma));
822 786
823 G4double factor = 1.0 + (expStopPower125/elo 787 G4double factor = 1.0 + (expStopPower125/eloss125 - 1.0) * f12525/
824 (1.0 + G4Exp( 1.48 * ( beta/beta25 - 7. 788 (1.0 + G4Exp( 1.48 * ( beta/beta25 - 7.0 ) ) );
825 789
826 return factor ; 790 return factor ;
827 } 791 }
828 792
829 //....oooOO0OOooo........oooOO0OOooo........oo 793 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
830 794
831 795