| 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: G4BraggIonModel.cc 93567 2015-10-26 14:51:41Z gcosmo $
>> 27 //
26 // ------------------------------------------- 28 // -------------------------------------------------------------------
27 // 29 //
28 // GEANT4 Class file 30 // GEANT4 Class file
29 // 31 //
30 // 32 //
31 // File name: G4BraggIonModel 33 // File name: G4BraggIonModel
32 // 34 //
33 // Author: Vladimir Ivanchenko 35 // Author: Vladimir Ivanchenko
34 // 36 //
35 // Creation date: 13.10.2004 37 // Creation date: 13.10.2004
36 // 38 //
37 // Modifications: 39 // Modifications:
38 // 11-05-05 Major optimisation of internal int 40 // 11-05-05 Major optimisation of internal interfaces (V.Ivantchenko)
39 // 29-11-05 Do not use G4Alpha class (V.Ivantc 41 // 29-11-05 Do not use G4Alpha class (V.Ivantchenko)
40 // 15-02-06 ComputeCrossSectionPerElectron, Co 42 // 15-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
41 // 25-04-06 Add stopping data from ASTAR (V.Iv 43 // 25-04-06 Add stopping data from ASTAR (V.Ivanchenko)
42 // 23-10-06 Reduce lowestKinEnergy to 0.25 keV 44 // 23-10-06 Reduce lowestKinEnergy to 0.25 keV (V.Ivanchenko)
43 // 12-08-08 Added methods GetParticleCharge, G 45 // 12-08-08 Added methods GetParticleCharge, GetChargeSquareRatio,
44 // CorrectionsAlongStep needed for io 46 // CorrectionsAlongStep needed for ions(V.Ivanchenko)
45 // 47 //
46 48
47 // Class Description: 49 // Class Description:
48 // 50 //
49 // Implementation of energy loss and delta-ele 51 // Implementation of energy loss and delta-electron production by
50 // slow charged heavy particles 52 // slow charged heavy particles
51 53
52 // ------------------------------------------- 54 // -------------------------------------------------------------------
53 // 55 //
54 56
55 //....oooOO0OOooo........oooOO0OOooo........oo 57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 //....oooOO0OOooo........oooOO0OOooo........oo 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57 59
58 #include "G4BraggIonModel.hh" 60 #include "G4BraggIonModel.hh"
59 #include "G4PhysicalConstants.hh" 61 #include "G4PhysicalConstants.hh"
60 #include "G4SystemOfUnits.hh" 62 #include "G4SystemOfUnits.hh"
61 #include "Randomize.hh" 63 #include "Randomize.hh"
62 #include "G4Electron.hh" 64 #include "G4Electron.hh"
63 #include "G4ParticleChangeForLoss.hh" 65 #include "G4ParticleChangeForLoss.hh"
>> 66 #include "G4LossTableManager.hh"
64 #include "G4EmCorrections.hh" 67 #include "G4EmCorrections.hh"
65 #include "G4DeltaAngle.hh" 68 #include "G4DeltaAngle.hh"
66 #include "G4ICRU90StoppingData.hh" <<
67 #include "G4ASTARStopping.hh" <<
68 #include "G4PSTARStopping.hh" <<
69 #include "G4NistManager.hh" <<
70 #include "G4Log.hh" 69 #include "G4Log.hh"
71 #include "G4Exp.hh" 70 #include "G4Exp.hh"
72 #include "G4AutoLock.hh" <<
73 71
74 //....oooOO0OOooo........oooOO0OOooo........oo 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
75 73
76 G4ASTARStopping* G4BraggIonModel::fASTAR = nul << 74 using namespace std;
77 75
78 namespace << 76 G4ASTARStopping* G4BraggIonModel::fASTAR = nullptr;
79 { <<
80 G4Mutex alphaMutex = G4MUTEX_INITIALIZER; <<
81 } <<
82 77
83 G4BraggIonModel::G4BraggIonModel(const G4Parti 78 G4BraggIonModel::G4BraggIonModel(const G4ParticleDefinition* p,
84 const G4Strin 79 const G4String& nam)
85 : G4BraggModel(p, nam) << 80 : G4VEmModel(nam),
86 { << 81 corr(nullptr),
87 HeMass = 3.727417*CLHEP::GeV; << 82 particle(nullptr),
88 massFactor = 1000.*CLHEP::amu_c2/HeMass; << 83 fParticleChange(nullptr),
>> 84 currentMaterial(nullptr),
>> 85 iMolecula(-1),
>> 86 iASTAR(-1),
>> 87 isIon(false),
>> 88 isInitialised(false)
>> 89 {
>> 90 SetHighEnergyLimit(2.0*MeV);
>> 91
>> 92 HeMass = 3.727417*GeV;
>> 93 rateMassHe2p = HeMass/proton_mass_c2;
>> 94 lowestKinEnergy = 1.0*keV/rateMassHe2p;
>> 95 massFactor = 1000.*amu_c2/HeMass;
>> 96 theZieglerFactor = eV*cm2*1.0e-15;
>> 97 theElectron = G4Electron::Electron();
>> 98 corrFactor = 1.0;
>> 99 if(p) { SetParticle(p); }
>> 100 else { SetParticle(theElectron); }
89 } 101 }
90 102
91 //....oooOO0OOooo........oooOO0OOooo........oo 103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
92 104
93 G4BraggIonModel::~G4BraggIonModel() 105 G4BraggIonModel::~G4BraggIonModel()
94 { 106 {
95 if(isFirstAlpha) { << 107 if(IsMaster()) { delete fASTAR; fASTAR = nullptr; }
96 delete fASTAR; <<
97 fASTAR = nullptr; <<
98 } <<
99 } 108 }
100 109
101 //....oooOO0OOooo........oooOO0OOooo........oo 110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
102 111
103 void G4BraggIonModel::Initialise(const G4Parti 112 void G4BraggIonModel::Initialise(const G4ParticleDefinition* p,
104 const G4DataV << 113 const G4DataVector&)
105 { 114 {
106 G4BraggModel::Initialise(p, ref); << 115 if(p != particle) { SetParticle(p); }
107 const G4String& pname = particle->GetParticl << 116
108 if(pname == "alpha") { isAlpha = true; } << 117 corrFactor = chargeSquare;
109 if(isAlpha && fASTAR == nullptr) { << 118
110 G4AutoLock l(&alphaMutex); << 119 // always false before the run
111 if(fASTAR == nullptr) { << 120 SetDeexcitationFlag(false);
112 isFirstAlpha = true; << 121
113 fASTAR = new G4ASTARStopping(); << 122 if(!isInitialised) {
>> 123 isInitialised = true;
>> 124
>> 125 if(UseAngularGeneratorFlag() && !GetAngularDistribution()) {
>> 126 SetAngularDistribution(new G4DeltaAngle());
114 } 127 }
115 l.unlock(); << 128 G4String pname = particle->GetParticleName();
116 } << 129 if(particle->GetParticleType() == "nucleus" &&
117 if(isFirstAlpha) { << 130 pname != "deuteron" && pname != "triton" &&
118 fASTAR->Initialise(); << 131 pname != "alpha+" && pname != "helium" &&
>> 132 pname != "hydrogen") { isIon = true; }
>> 133
>> 134 corr = G4LossTableManager::Instance()->EmCorrections();
>> 135
>> 136 fParticleChange = GetParticleChangeForLoss();
>> 137 if(!fASTAR) { fASTAR = new G4ASTARStopping(); }
119 } 138 }
>> 139 if(IsMaster() && particle->GetPDGMass() < GeV) { fASTAR->Initialise(); }
120 } 140 }
121 141
>> 142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 143
>> 144 G4double G4BraggIonModel::MinEnergyCut(const G4ParticleDefinition*,
>> 145 const G4MaterialCutsCouple* couple)
>> 146 {
>> 147 return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
>> 148 }
122 149
123 //....oooOO0OOooo........oooOO0OOooo........oo 150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
124 151
125 G4double G4BraggIonModel::GetChargeSquareRatio 152 G4double G4BraggIonModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
126 153 const G4Material* mat,
127 << 154 G4double kineticEnergy)
128 { 155 {
129 // this method is called only for ions, so n << 156 //G4cout<<"G4BraggIonModel::GetChargeSquareRatio e= "<<kineticEnergy<<G4endl;
130 if(isAlpha) { return 1.0; } << 157 // this method is called only for ions
131 return G4BraggModel::GetChargeSquareRatio(p, << 158 G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,kineticEnergy);
>> 159 corrFactor = q2*corr->EffectiveChargeCorrection(p,mat,kineticEnergy);
>> 160 return corrFactor;
>> 161 }
>> 162
>> 163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 164
>> 165 G4double G4BraggIonModel::GetParticleCharge(const G4ParticleDefinition* p,
>> 166 const G4Material* mat,
>> 167 G4double kineticEnergy)
>> 168 {
>> 169 //G4cout<<"G4BraggIonModel::GetParticleCharge e= "<<kineticEnergy <<
>> 170 // " q= " << corr->GetParticleCharge(p,mat,kineticEnergy) <<G4endl;
>> 171 // this method is called only for ions
>> 172 return corr->GetParticleCharge(p,mat,kineticEnergy);
>> 173 }
>> 174
>> 175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 176
>> 177 G4double G4BraggIonModel::ComputeCrossSectionPerElectron(
>> 178 const G4ParticleDefinition* p,
>> 179 G4double kineticEnergy,
>> 180 G4double cutEnergy,
>> 181 G4double maxKinEnergy)
>> 182 {
>> 183 G4double cross = 0.0;
>> 184 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
>> 185 G4double maxEnergy = std::min(tmax,maxKinEnergy);
>> 186 if(cutEnergy < tmax) {
>> 187
>> 188 G4double energy = kineticEnergy + mass;
>> 189 G4double energy2 = energy*energy;
>> 190 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
>> 191 cross = (maxEnergy - cutEnergy)/(cutEnergy*maxEnergy)
>> 192 - beta2*G4Log(maxEnergy/cutEnergy)/tmax;
>> 193
>> 194 if( 0.0 < spin ) { cross += 0.5*(maxEnergy - cutEnergy)/energy2; }
>> 195
>> 196 cross *= twopi_mc2_rcl2*chargeSquare/beta2;
>> 197 }
>> 198 // G4cout << "BR: e= " << kineticEnergy << " tmin= " << cutEnergy
>> 199 // << " tmax= " << tmax << " cross= " << cross << G4endl;
>> 200
>> 201 return cross;
132 } 202 }
133 203
134 //....oooOO0OOooo........oooOO0OOooo........oo 204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
135 205
136 G4double G4BraggIonModel::ComputeCrossSectionP 206 G4double G4BraggIonModel::ComputeCrossSectionPerAtom(
137 con 207 const G4ParticleDefinition* p,
138 << 208 G4double kineticEnergy,
139 209 G4double Z, G4double,
140 210 G4double cutEnergy,
141 211 G4double maxEnergy)
142 { 212 {
143 G4double sigma = << 213 G4double cross = Z*ComputeCrossSectionPerElectron
144 Z*ComputeCrossSectionPerElectron(p,kinEner << 214 (p,kineticEnergy,cutEnergy,maxEnergy);
145 if(isAlpha) { << 215 return cross;
146 sigma *= (HeEffChargeSquare(Z, kinEnergy/C <<
147 } <<
148 return sigma; <<
149 } 216 }
150 217
151 //....oooOO0OOooo........oooOO0OOooo........oo 218 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
152 219
153 G4double G4BraggIonModel::CrossSectionPerVolum 220 G4double G4BraggIonModel::CrossSectionPerVolume(
154 con 221 const G4Material* material,
155 con 222 const G4ParticleDefinition* p,
156 << 223 G4double kineticEnergy,
157 224 G4double cutEnergy,
158 225 G4double maxEnergy)
159 { 226 {
160 G4double sigma = material->GetElectronDensit << 227 G4double eDensity = material->GetElectronDensity();
161 ComputeCrossSectionPerElectron(p,kinEnergy << 228 G4double cross = eDensity*ComputeCrossSectionPerElectron
162 if(isAlpha) { << 229 (p,kineticEnergy,cutEnergy,maxEnergy);
163 const G4double zeff = material->GetTotNbOf << 230 return cross;
164 material->GetTotNbOfAtomsPerVolume(); <<
165 sigma *= (HeEffChargeSquare(zeff, kinEnerg <<
166 } <<
167 return sigma; <<
168 } 231 }
169 232
170 //....oooOO0OOooo........oooOO0OOooo........oo 233 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
171 234
172 G4double G4BraggIonModel::ComputeDEDXPerVolume 235 G4double G4BraggIonModel::ComputeDEDXPerVolume(const G4Material* material,
173 236 const G4ParticleDefinition* p,
174 237 G4double kineticEnergy,
175 << 238 G4double cutEnergy)
176 { 239 {
177 const G4double tmax = MaxSecondaryEnergy(p, << 240 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
178 const G4double tlim = lowestKinEnergy*massRa << 241 G4double tmin = min(cutEnergy, tmax);
179 const G4double tmin = std::max(std::min(cut, << 242 G4double tkin = kineticEnergy/massRate;
180 G4double dedx = 0.0; << 243 G4double dedx = 0.0;
181 244
182 if(kineticEnergy < tlim) { << 245 if(tkin < lowestKinEnergy) {
183 dedx = HeDEDX(material, tlim)*std::sqrt(ki << 246 dedx = DEDX(material, lowestKinEnergy)*sqrt(tkin/lowestKinEnergy);
184 } else { 247 } else {
185 dedx = HeDEDX(material, kineticEnergy); << 248 dedx = DEDX(material, tkin);
>> 249 }
186 250
187 if (tmin < tmax) { << 251 if (cutEnergy < tmax) {
188 const G4double tau = kineticEnergy/mass; << 252
189 const G4double x = tmin/tmax; << 253 G4double tau = kineticEnergy/mass;
190 << 254 G4double gam = tau + 1.0;
191 G4double del = << 255 G4double bg2 = tau * (tau+2.0);
192 (G4Log(x)*(tau + 1.)*(tau + 1.)/(tau * << 256 G4double beta2 = bg2/(gam*gam);
193 CLHEP::twopi_mc2_rcl2*material->GetElectronD << 257 G4double x = tmin/tmax;
194 if(isAlpha) { << 258
195 const G4double zeff = material->GetTotNbOfEl << 259 dedx += (G4Log(x) + (1.0 - x)*beta2) * twopi_mc2_rcl2
196 material->GetTotNbOfAtomsPerVolume(); << 260 * (material->GetElectronDensity())/beta2;
197 heChargeSquare = HeEffChargeSquare(zeff, kin <<
198 del *= heChargeSquare; <<
199 } <<
200 dedx += del; <<
201 } <<
202 } 261 }
203 dedx = std::max(dedx, 0.0); << 262
204 /* << 263 // now compute the total ionization loss
205 G4cout << "BraggIon: " << material->GetNam << 264
206 << " E(MeV)=" << kineticEnergy/MeV << 265 if (dedx < 0.0) dedx = 0.0 ;
207 << " Tmin(MeV)=" << tmin << " dedx( << 266
208 << dedx*gram/(MeV*cm2*material->Get << 267 dedx *= chargeSquare;
209 << " q2=" << chargeSquare << G4endl << 268 /*
>> 269 G4cout << "Bragg: tkin(MeV) = " << tkin/MeV << " dedx(MeVxcm^2/g) = "
>> 270 << dedx*gram/(MeV*cm2*material->GetDensity())
>> 271 << " q2 = " << chargeSquare << G4endl;
210 */ 272 */
211 return dedx; 273 return dedx;
212 } 274 }
213 275
214 //....oooOO0OOooo........oooOO0OOooo........oo 276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
215 277
216 void G4BraggIonModel::CorrectionsAlongStep(con 278 void G4BraggIonModel::CorrectionsAlongStep(const G4MaterialCutsCouple* couple,
217 con 279 const G4DynamicParticle* dp,
218 con << 280 G4double& eloss,
219 G4d << 281 G4double&,
>> 282 G4double /*length*/)
220 { 283 {
221 // no correction for alpha << 284 // this method is called only for ions
222 if(isAlpha) { return; } << 285 const G4ParticleDefinition* p = dp->GetDefinition();
>> 286 const G4Material* mat = couple->GetMaterial();
>> 287 G4double preKinEnergy = dp->GetKineticEnergy();
>> 288 G4double e = preKinEnergy - eloss*0.5;
>> 289 if(e < 0.0) { e = preKinEnergy*0.5; }
223 290
224 // no correction at a small step at the last << 291 G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,e);
225 const G4double preKinEnergy = dp->GetKinetic << 292 GetModelOfFluctuations()->SetParticleAndCharge(p, q2);
226 if(eloss >= preKinEnergy || eloss < preKinEn << 293 G4double qfactor = q2*corr->EffectiveChargeCorrection(p,mat,e)/corrFactor;
>> 294 eloss *= qfactor;
227 295
228 // corrections only for ions << 296 //G4cout << "G4BraggIonModel::CorrectionsAlongStep e= " << e
229 const G4ParticleDefinition* p = dp->GetDefin << 297 // << " qfactor= " << qfactor << " " << p->GetParticleName() <<G4endl;
230 if(p != particle) { SetParticle(p); } << 298 }
231 299
232 // effective energy and charge at a step << 300 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
233 const G4Material* mat = couple->GetMaterial( << 301
234 const G4double e = std::max(preKinEnergy - e << 302 void G4BraggIonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
235 const G4double q20 = corr->EffectiveChargeSq << 303 const G4MaterialCutsCouple* couple,
236 const G4double q2 = corr->EffectiveChargeSqu << 304 const G4DynamicParticle* dp,
237 const G4double qfactor = q2/q20; << 305 G4double xmin,
238 /* << 306 G4double maxEnergy)
239 G4cout << "G4BraggIonModel::CorrectionsAlo << 307 {
240 << preKinEnergy << " Eeff(MeV)=" << e << 308 G4double tmax = MaxSecondaryKinEnergy(dp);
241 << " eloss=" << eloss << " elossnew=" << e << 309 G4double xmax = std::min(tmax, maxEnergy);
242 << " qfactor=" << qfactor << " Qpre=" << q << 310 if(xmin >= xmax) { return; }
243 << p->GetParticleName() <<G4endl; << 311
244 */ << 312 G4double kineticEnergy = dp->GetKineticEnergy();
245 eloss *= qfactor; << 313 G4double energy = kineticEnergy + mass;
>> 314 G4double energy2 = energy*energy;
>> 315 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
>> 316 G4double grej = 1.0;
>> 317 G4double deltaKinEnergy, f;
>> 318
>> 319 CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
>> 320 G4double rndm[2];
>> 321
>> 322 // sampling follows ...
>> 323 do {
>> 324 rndmEngineMod->flatArray(2, rndm);
>> 325 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - rndm[0]) + xmax*rndm[0]);
>> 326
>> 327 f = 1.0 - beta2*deltaKinEnergy/tmax;
>> 328
>> 329 if(f > grej) {
>> 330 G4cout << "G4BraggIonModel::SampleSecondary Warning! "
>> 331 << "Majorant " << grej << " < "
>> 332 << f << " for e= " << deltaKinEnergy
>> 333 << G4endl;
>> 334 }
>> 335
>> 336 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
>> 337 } while( grej*rndm[1] >= f );
>> 338
>> 339 G4ThreeVector deltaDirection;
>> 340
>> 341 if(UseAngularGeneratorFlag()) {
>> 342 const G4Material* mat = couple->GetMaterial();
>> 343 G4int Z = SelectRandomAtomNumber(mat);
>> 344
>> 345 deltaDirection =
>> 346 GetAngularDistribution()->SampleDirection(dp, deltaKinEnergy, Z, mat);
>> 347
>> 348 } else {
>> 349
>> 350 G4double deltaMomentum =
>> 351 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
>> 352 G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
>> 353 (deltaMomentum * dp->GetTotalMomentum());
>> 354 if(cost > 1.0) { cost = 1.0; }
>> 355 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
>> 356
>> 357 G4double phi = twopi*rndmEngineMod->flat();
>> 358
>> 359 deltaDirection.set(sint*cos(phi),sint*sin(phi), cost) ;
>> 360 deltaDirection.rotateUz(dp->GetMomentumDirection());
>> 361 }
>> 362
>> 363 // create G4DynamicParticle object for delta ray
>> 364 G4DynamicParticle* delta =
>> 365 new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
>> 366
>> 367 vdp->push_back(delta);
>> 368
>> 369 // Change kinematics of primary particle
>> 370 kineticEnergy -= deltaKinEnergy;
>> 371 G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum();
>> 372 finalP = finalP.unit();
>> 373
>> 374 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
>> 375 fParticleChange->SetProposedMomentumDirection(finalP);
246 } 376 }
247 377
248 //....oooOO0OOooo........oooOO0OOooo........oo 378 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
249 379
250 G4int G4BraggIonModel::HasMaterialForHe(const << 380 G4double G4BraggIonModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
>> 381 G4double kinEnergy)
251 { 382 {
252 const G4String& chFormula = mat->GetChemical << 383 if(pd != particle) { SetParticle(pd); }
253 if(chFormula.empty()) { return -1; } << 384 G4double tau = kinEnergy/mass;
>> 385 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
>> 386 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
>> 387 return tmax;
>> 388 }
>> 389
>> 390 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 391
>> 392 G4bool G4BraggIonModel::HasMaterial(const G4Material*)
>> 393 {
>> 394 return false;
>> 395 /*
>> 396 G4String chFormula = material->GetChemicalFormula();
>> 397 if("" == chFormula) { return false; }
254 398
255 // ICRU Report N49, 1993. Ziegler model for 399 // ICRU Report N49, 1993. Ziegler model for He.
256 400
257 static const G4int numberOfMolecula = 11; << 401 static const size_t numberOfMolecula = 11;
258 static const G4String molName[numberOfMolecu 402 static const G4String molName[numberOfMolecula] = {
259 "CaF_2", "Cellulose_Nitrate", "LiF", "Po 403 "CaF_2", "Cellulose_Nitrate", "LiF", "Policarbonate",
260 "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Po 404 "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Polymethly_Methacralate",
261 "Polysterene", "SiO_2", "NaI", "H_2O", 405 "Polysterene", "SiO_2", "NaI", "H_2O",
262 "Graphite" }; 406 "Graphite" };
263 407
264 // Search for the material in the table 408 // Search for the material in the table
265 for (G4int i=0; i<numberOfMolecula; ++i) { << 409 for (size_t i=0; i<numberOfMolecula; ++i) {
266 if (chFormula == molName[i]) { << 410 if (chFormula == molName[i]) {
267 return i; << 411 iASTAR = fASTAR->GetIndex(matName[i]);
>> 412 break;
268 } 413 }
269 } 414 }
270 return -1; << 415 return (iASTAR >= 0);
>> 416 */
271 } 417 }
272 418
273 //....oooOO0OOooo........oooOO0OOooo........oo 419 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
274 420
275 G4double G4BraggIonModel::HeStoppingPower(cons << 421 G4double G4BraggIonModel::StoppingPower(const G4Material* material,
>> 422 G4double kineticEnergy)
276 { 423 {
277 G4double ionloss = 0.0; << 424 G4double ionloss = 0.0 ;
>> 425
278 if (iMolecula >= 0) { 426 if (iMolecula >= 0) {
279 427
280 // The data and the fit from: 428 // The data and the fit from:
281 // ICRU Report N49, 1993. Ziegler's model 429 // ICRU Report N49, 1993. Ziegler's model for alpha
282 // He energy in internal units of parametr 430 // He energy in internal units of parametrisation formula (MeV)
283 // Input scaled energy of a proton or Gene <<
284 G4double T = kineticEnergy/(massRate*CLHEP <<
285 <<
286 static const G4float a[11][5] = { <<
287 {9.43672f, 0.54398f, 84.341f, 1.3705f, <<
288 {67.1503f, 0.41409f, 404.512f, 148.97f, <<
289 {5.11203f, 0.453f, 36.718f, 50.6f, <<
290 {61.793f, 0.48445f, 361.537f, 57.889f, <<
291 {7.83464f, 0.49804f, 160.452f, 3.192f, <<
292 {19.729f, 0.52153f, 162.341f, 58.35f, <<
293 {26.4648f, 0.50112f, 188.913f, 30.079f, <<
294 {7.8655f, 0.5205f, 63.96f, 51.32f, <<
295 {8.8965f, 0.5148f, 339.36f, 1.7205f, <<
296 {2.959f, 0.53255f, 34.247f, 60.655f, <<
297 {3.80133f, 0.41590f, 12.9966f, 117.83f, <<
298 431
299 static const G4double atomicWeight[11] = { << 432 G4double T = kineticEnergy*rateMassHe2p/MeV ;
300 101.96128f, 44.0098f, 16.0426f, 28.0536 <<
301 104.1512f, 44.665f, 60.0843f, 18.0152 <<
302 433
303 const G4int i = iMolecula; << 434 static const G4double a[11][5] = {
>> 435 {9.43672, 0.54398, 84.341, 1.3705, 57.422},
>> 436 {67.1503, 0.41409, 404.512, 148.97, 20.99},
>> 437 {5.11203, 0.453, 36.718, 50.6, 28.058},
>> 438 {61.793, 0.48445, 361.537, 57.889, 50.674},
>> 439 {7.83464, 0.49804, 160.452, 3.192, 0.71922},
>> 440 {19.729, 0.52153, 162.341, 58.35, 25.668},
>> 441 {26.4648, 0.50112, 188.913, 30.079, 16.509},
>> 442 {7.8655, 0.5205, 63.96, 51.32, 67.775},
>> 443 {8.8965, 0.5148, 339.36, 1.7205, 0.70423},
>> 444 {2.959, 0.53255, 34.247, 60.655, 15.153},
>> 445 {3.80133, 0.41590, 12.9966, 117.83, 242.28} };
304 446
305 G4double slow = (G4double)(a[i][0]); << 447 static const G4double atomicWeight[11] = {
>> 448 101.96128, 44.0098, 16.0426, 28.0536, 42.0804,
>> 449 104.1512, 44.665, 60.0843, 18.0152, 18.0152, 12.0};
306 450
307 G4double x1 = (G4double)(a[i][1]); << 451 G4int i = iMolecula;
308 G4double x2 = (G4double)(a[i][2]); <<
309 G4double x3 = (G4double)(a[i][3]); <<
310 G4double x4 = (G4double)(a[i][4]); <<
311 452
312 // Free electron gas model 453 // Free electron gas model
313 if ( T < 0.001 ) { 454 if ( T < 0.001 ) {
314 G4double shigh = G4Log( 1.0 + x3*1000.0 << 455 G4double slow = a[i][0] ;
>> 456 G4double shigh = G4Log( 1.0 + a[i][3]*1000.0 + a[i][4]*0.001 )
>> 457 * a[i][2]*1000.0 ;
315 ionloss = slow*shigh / (slow + shigh) ; 458 ionloss = slow*shigh / (slow + shigh) ;
316 ionloss *= std::sqrt(T*1000.0) ; << 459 ionloss *= sqrt(T*1000.0) ;
317 460
318 // Main parametrisation 461 // Main parametrisation
319 } else { 462 } else {
320 slow *= G4Exp(G4Log(T*1000.0)*x1) ; << 463 G4double slow = a[i][0] * G4Exp(G4Log(T*1000.0)*a[i][1]) ;
321 G4double shigh = G4Log( 1.0 + x3/T + x4* << 464 G4double shigh = G4Log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ;
322 ionloss = slow*shigh / (slow + shigh) ; 465 ionloss = slow*shigh / (slow + shigh) ;
323 /* 466 /*
324 G4cout << "## " << i << ". T= " << T 467 G4cout << "## " << i << ". T= " << T << " slow= " << slow
325 << " a0= " << a[i][0] << " a1= " << a 468 << " a0= " << a[i][0] << " a1= " << a[i][1]
326 << " shigh= " << shigh 469 << " shigh= " << shigh
327 << " dedx= " << ionloss << " q^2= " < 470 << " dedx= " << ionloss << " q^2= " << HeEffChargeSquare(z, T*MeV)
328 << G4endl; 471 << G4endl;
329 */ 472 */
330 } 473 }
331 ionloss = std::max(ionloss, 0.0) * atomicW << 474 if ( ionloss < 0.0) ionloss = 0.0 ;
>> 475
>> 476 // He effective charge
>> 477 G4double aa = atomicWeight[iMolecula];
>> 478 ionloss /= (HeEffChargeSquare(0.5*aa, T)*aa);
>> 479
>> 480 // pure material (normally not the case for this function)
>> 481 } else if(1 == (material->GetNumberOfElements())) {
>> 482 G4double z = material->GetZ() ;
>> 483 ionloss = ElectronicStoppingPower( z, kineticEnergy ) ;
332 } 484 }
>> 485
333 return ionloss; 486 return ionloss;
334 } 487 }
335 488
336 //....oooOO0OOooo........oooOO0OOooo........oo 489 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
337 490
338 G4double G4BraggIonModel::HeElectronicStopping << 491 G4double G4BraggIonModel::ElectronicStoppingPower(G4double z,
339 const G4double kinet << 492 G4double kineticEnergy) const
340 { 493 {
341 G4double ionloss ; 494 G4double ionloss ;
342 G4int i = std::min(z-1, 91); // index of at << 495 G4int i = G4lrint(z)-1 ; // index of atom
343 //G4cout << "ElectronicStoppingPower z=" << << 496 if(i < 0) i = 0 ;
344 // << " E=" << kineticEnergy << G4endl; << 497 if(i > 91) i = 91 ;
>> 498
345 // The data and the fit from: 499 // The data and the fit from:
346 // ICRU Report 49, 1993. Ziegler's type of p 500 // ICRU Report 49, 1993. Ziegler's type of parametrisations.
347 // Proton kinetic energy for parametrisation 501 // Proton kinetic energy for parametrisation (keV/amu)
348 // He energy in internal units of parametris << 502
349 //G4double T = kineticEnergy*rateMassHe2p/CL << 503 // He energy in internal units of parametrisation formula (MeV)
350 G4double T = kineticEnergy/CLHEP::MeV; << 504 G4double T = kineticEnergy*rateMassHe2p/MeV ;
351 << 505
352 static const G4float a[92][5] = { << 506 static const G4double a[92][5] = {
353 { 0.35485f, 0.6456f, 6.01525f, 20.8933f, << 507 {0.35485, 0.6456, 6.01525, 20.8933, 4.3515
354 },{ 0.58f, 0.59f, 6.3f, 130.0f, << 508 },{ 0.58, 0.59, 6.3, 130.0, 44.07
355 },{ 1.42f, 0.49f, 12.25f, 32.0f, << 509 },{ 1.42, 0.49, 12.25, 32.0, 9.161
356 },{ 2.206f, 0.51f, 15.32f, 0.25f, << 510 },{ 2.206, 0.51, 15.32, 0.25, 8.995 //Be Ziegler77
357 // },{ 2.1895f, 0.47183,7.2362f, 134 << 511 // },{ 2.1895, 0.47183,7.2362, 134.30, 197.96 //Be from ICRU
358 },{ 3.691f, 0.4128f, 18.48f, 50.72f, << 512 },{ 3.691, 0.4128, 18.48, 50.72, 9.0
359 },{ 3.83523f, 0.42993f,12.6125f, 227.41f, << 513 },{ 3.83523, 0.42993,12.6125, 227.41, 188.97
360 // },{ 1.9259f, 0.5550f, 27.15125f, 26 << 514 },{ 1.9259, 0.5550, 27.15125, 26.0665, 6.2768
361 },{ 1.9259f, 0.5550f, 27.1513f, 26.0665f, << 515 },{ 2.81015, 0.4759, 50.0253, 10.556, 1.0382
362 },{ 2.81015f, 0.4759f, 50.0253f, 10.556f, << 516 },{ 1.533, 0.531, 40.44, 18.41, 2.718
363 },{ 1.533f, 0.531f, 40.44f, 18.41f, << 517 },{ 2.303, 0.4861, 37.01, 37.96, 5.092
364 },{ 2.303f, 0.4861f, 37.01f, 37.96f, <<
365 // Z= 11-20 518 // Z= 11-20
366 },{ 9.894f, 0.3081f, 23.65f, 0.384f, << 519 },{ 9.894, 0.3081, 23.65, 0.384, 92.93
367 },{ 4.3f, 0.47f, 34.3f, 3.3f, << 520 },{ 4.3, 0.47, 34.3, 3.3, 12.74
368 },{ 2.5f, 0.625f, 45.7f, 0.1f, << 521 },{ 2.5, 0.625, 45.7, 0.1, 4.359
369 },{ 2.1f, 0.65f, 49.34f, 1.788f, << 522 },{ 2.1, 0.65, 49.34, 1.788, 4.133
370 },{ 1.729f, 0.6562f, 53.41f, 2.405f, << 523 },{ 1.729, 0.6562, 53.41, 2.405, 3.845
371 },{ 1.402f, 0.6791f, 58.98f, 3.528f, << 524 },{ 1.402, 0.6791, 58.98, 3.528, 3.211
372 },{ 1.117f, 0.7044f, 69.69f, 3.705f, << 525 },{ 1.117, 0.7044, 69.69, 3.705, 2.156
373 },{ 2.291f, 0.6284f, 73.88f, 4.478f, << 526 },{ 2.291, 0.6284, 73.88, 4.478, 2.066
374 },{ 8.554f, 0.3817f, 83.61f, 11.84f, << 527 },{ 8.554, 0.3817, 83.61, 11.84, 1.875
375 },{ 6.297f, 0.4622f, 65.39f, 10.14f, << 528 },{ 6.297, 0.4622, 65.39, 10.14, 5.036
376 // Z= 21-30 529 // Z= 21-30
377 },{ 5.307f, 0.4918f, 61.74f, 12.4f, << 530 },{ 5.307, 0.4918, 61.74, 12.4, 6.665
378 },{ 4.71f, 0.5087f, 65.28f, 8.806f, << 531 },{ 4.71, 0.5087, 65.28, 8.806, 5.948
379 },{ 6.151f, 0.4524f, 83.0f, 18.31f, << 532 },{ 6.151, 0.4524, 83.0, 18.31, 2.71
380 },{ 6.57f, 0.4322f, 84.76f, 15.53f, << 533 },{ 6.57, 0.4322, 84.76, 15.53, 2.779
381 },{ 5.738f, 0.4492f, 84.6f, 14.18f, << 534 },{ 5.738, 0.4492, 84.6, 14.18, 3.101
382 },{ 5.013f, 0.4707f, 85.8f, 16.55f, << 535 },{ 5.013, 0.4707, 85.8, 16.55, 3.211
383 },{ 4.32f, 0.4947f, 76.14f, 10.85f, << 536 },{ 4.32, 0.4947, 76.14, 10.85, 5.441
384 },{ 4.652f, 0.4571f, 80.73f, 22.0f, << 537 },{ 4.652, 0.4571, 80.73, 22.0, 4.952
385 },{ 3.114f, 0.5236f, 76.67f, 7.62f, << 538 },{ 3.114, 0.5236, 76.67, 7.62, 6.385
386 },{ 3.114f, 0.5236f, 76.67f, 7.62f, << 539 },{ 3.114, 0.5236, 76.67, 7.62, 7.502
387 // Z= 31-40 540 // Z= 31-40
388 },{ 3.114f, 0.5236f, 76.67f, 7.62f, << 541 },{ 3.114, 0.5236, 76.67, 7.62, 8.514
389 },{ 5.746f, 0.4662f, 79.24f, 1.185f, << 542 },{ 5.746, 0.4662, 79.24, 1.185, 7.993
390 },{ 2.792f, 0.6346f, 106.1f, 0.2986f, << 543 },{ 2.792, 0.6346, 106.1, 0.2986, 2.331
391 },{ 4.667f, 0.5095f, 124.3f, 2.102f, << 544 },{ 4.667, 0.5095, 124.3, 2.102, 1.667
392 },{ 2.44f, 0.6346f, 105.0f, 0.83f, << 545 },{ 2.44, 0.6346, 105.0, 0.83, 2.851
393 },{ 1.413f, 0.7377f, 147.9f, 1.466f, << 546 },{ 1.413, 0.7377, 147.9, 1.466, 1.016
394 },{ 11.72f, 0.3826f, 102.8f, 9.231f, << 547 },{ 11.72, 0.3826, 102.8, 9.231, 4.371
395 },{ 7.126f, 0.4804f, 119.3f, 5.784f, << 548 },{ 7.126, 0.4804, 119.3, 5.784, 2.454
396 },{ 11.61f, 0.3955f, 146.7f, 7.031f, << 549 },{ 11.61, 0.3955, 146.7, 7.031, 1.423
397 },{ 10.99f, 0.41f, 163.9f, 7.1f, << 550 },{ 10.99, 0.41, 163.9, 7.1, 1.052
398 // Z= 41-50 551 // Z= 41-50
399 },{ 9.241f, 0.4275f, 163.1f, 7.954f, << 552 },{ 9.241, 0.4275, 163.1, 7.954, 1.102
400 },{ 9.276f, 0.418f, 157.1f, 8.038f, << 553 },{ 9.276, 0.418, 157.1, 8.038, 1.29
401 },{ 3.999f, 0.6152f, 97.6f, 1.297f, << 554 },{ 3.999, 0.6152, 97.6, 1.297, 5.792
402 },{ 4.306f, 0.5658f, 97.99f, 5.514f, << 555 },{ 4.306, 0.5658, 97.99, 5.514, 5.754
403 },{ 3.615f, 0.6197f, 86.26f, 0.333f, << 556 },{ 3.615, 0.6197, 86.26, 0.333, 8.689
404 },{ 5.8f, 0.49f, 147.2f, 6.903f, << 557 },{ 5.8, 0.49, 147.2, 6.903, 1.289
405 },{ 5.6f, 0.49f, 130.0f, 10.0f, << 558 },{ 5.6, 0.49, 130.0, 10.0, 2.844
406 },{ 3.55f, 0.6068f, 124.7f, 1.112f, << 559 },{ 3.55, 0.6068, 124.7, 1.112, 3.119
407 },{ 3.6f, 0.62f, 105.8f, 0.1692f, << 560 },{ 3.6, 0.62, 105.8, 0.1692, 6.026
408 },{ 5.4f, 0.53f, 103.1f, 3.931f, << 561 },{ 5.4, 0.53, 103.1, 3.931, 7.767
409 // Z= 51-60 562 // Z= 51-60
410 },{ 3.97f, 0.6459f, 131.8f, 0.2233f, << 563 },{ 3.97, 0.6459, 131.8, 0.2233, 2.723
411 },{ 3.65f, 0.64f, 126.8f, 0.6834f, << 564 },{ 3.65, 0.64, 126.8, 0.6834, 3.411
412 },{ 3.118f, 0.6519f, 164.9f, 1.208f, << 565 },{ 3.118, 0.6519, 164.9, 1.208, 1.51
413 },{ 3.949f, 0.6209f, 200.5f, 1.878f, << 566 },{ 3.949, 0.6209, 200.5, 1.878, 0.9126
414 },{ 14.4f, 0.3923f, 152.5f, 8.354f, << 567 },{ 14.4, 0.3923, 152.5, 8.354, 2.597
415 },{ 10.99f, 0.4599f, 138.4f, 4.811f, << 568 },{ 10.99, 0.4599, 138.4, 4.811, 3.726
416 },{ 16.6f, 0.3773f, 224.1f, 6.28f, << 569 },{ 16.6, 0.3773, 224.1, 6.28, 0.9121
417 },{ 10.54f, 0.4533f, 159.3f, 4.832f, << 570 },{ 10.54, 0.4533, 159.3, 4.832, 2.529
418 },{ 10.33f, 0.4502f, 162.0f, 5.132f, << 571 },{ 10.33, 0.4502, 162.0, 5.132, 2.444
419 },{ 10.15f, 0.4471f, 165.6f, 5.378f, << 572 },{ 10.15, 0.4471, 165.6, 5.378, 2.328
420 // Z= 61-70 573 // Z= 61-70
421 },{ 9.976f, 0.4439f, 168.0f, 5.721f, << 574 },{ 9.976, 0.4439, 168.0, 5.721, 2.258
422 },{ 9.804f, 0.4408f, 176.2f, 5.675f, << 575 },{ 9.804, 0.4408, 176.2, 5.675, 1.997
423 },{ 14.22f, 0.363f, 228.4f, 7.024f, << 576 },{ 14.22, 0.363, 228.4, 7.024, 1.016
424 },{ 9.952f, 0.4318f, 233.5f, 5.065f, << 577 },{ 9.952, 0.4318, 233.5, 5.065, 0.9244
425 },{ 9.272f, 0.4345f, 210.0f, 4.911f, << 578 },{ 9.272, 0.4345, 210.0, 4.911, 1.258
426 },{ 10.13f, 0.4146f, 225.7f, 5.525f, << 579 },{ 10.13, 0.4146, 225.7, 5.525, 1.055
427 },{ 8.949f, 0.4304f, 213.3f, 5.071f, << 580 },{ 8.949, 0.4304, 213.3, 5.071, 1.221
428 },{ 11.94f, 0.3783f, 247.2f, 6.655f, << 581 },{ 11.94, 0.3783, 247.2, 6.655, 0.849
429 },{ 8.472f, 0.4405f, 195.5f, 4.051f, << 582 },{ 8.472, 0.4405, 195.5, 4.051, 1.604
430 },{ 8.301f, 0.4399f, 203.7f, 3.667f, << 583 },{ 8.301, 0.4399, 203.7, 3.667, 1.459
431 // Z= 71-80 584 // Z= 71-80
432 },{ 6.567f, 0.4858f, 193.0f, 2.65f, << 585 },{ 6.567, 0.4858, 193.0, 2.65, 1.66
433 },{ 5.951f, 0.5016f, 196.1f, 2.662f, << 586 },{ 5.951, 0.5016, 196.1, 2.662, 1.589
434 },{ 7.495f, 0.4523f, 251.4f, 3.433f, << 587 },{ 7.495, 0.4523, 251.4, 3.433, 0.8619
435 },{ 6.335f, 0.4825f, 255.1f, 2.834f, << 588 },{ 6.335, 0.4825, 255.1, 2.834, 0.8228
436 },{ 4.314f, 0.5558f, 214.8f, 2.354f, << 589 },{ 4.314, 0.5558, 214.8, 2.354, 1.263
437 },{ 4.02f, 0.5681f, 219.9f, 2.402f, << 590 },{ 4.02, 0.5681, 219.9, 2.402, 1.191
438 },{ 3.836f, 0.5765f, 210.2f, 2.742f, << 591 },{ 3.836, 0.5765, 210.2, 2.742, 1.305
439 },{ 4.68f, 0.5247f, 244.7f, 2.749f, << 592 },{ 4.68, 0.5247, 244.7, 2.749, 0.8962
440 },{ 2.892f, 0.6204f, 208.6f, 2.415f, << 593 },{ 2.892, 0.6204, 208.6, 2.415, 1.416 //Au Z77
441 // },{ 3.223f, 0.5883f, 232.7f, 2.9 << 594 // },{ 3.223, 0.5883, 232.7, 2.954, 1.05 //Au ICRU
442 },{ 2.892f, 0.6204f, 208.6f, 2.415f, << 595 },{ 2.892, 0.6204, 208.6, 2.415, 1.416
443 // Z= 81-90 596 // Z= 81-90
444 },{ 4.728f, 0.5522f, 217.0f, 3.091f, << 597 },{ 4.728, 0.5522, 217.0, 3.091, 1.386
445 },{ 6.18f, 0.52f, 170.0f, 4.0f, << 598 },{ 6.18, 0.52, 170.0, 4.0, 3.224
446 },{ 9.0f, 0.47f, 198.0f, 3.8f, << 599 },{ 9.0, 0.47, 198.0, 3.8, 2.032
447 },{ 2.324f, 0.6997f, 216.0f, 1.599f, << 600 },{ 2.324, 0.6997, 216.0, 1.599, 1.399
448 },{ 1.961f, 0.7286f, 223.0f, 1.621f, << 601 },{ 1.961, 0.7286, 223.0, 1.621, 1.296
449 },{ 1.75f, 0.7427f, 350.1f, 0.9789f, << 602 },{ 1.75, 0.7427, 350.1, 0.9789, 0.5507
450 },{ 10.31f, 0.4613f, 261.2f, 4.738f, << 603 },{ 10.31, 0.4613, 261.2, 4.738, 0.9899
451 },{ 7.962f, 0.519f, 235.7f, 4.347f, << 604 },{ 7.962, 0.519, 235.7, 4.347, 1.313
452 },{ 6.227f, 0.5645f, 231.9f, 3.961f, << 605 },{ 6.227, 0.5645, 231.9, 3.961, 1.379
453 },{ 5.246f, 0.5947f, 228.6f, 4.027f, << 606 },{ 5.246, 0.5947, 228.6, 4.027, 1.432
454 // Z= 91-92 607 // Z= 91-92
455 },{ 5.408f, 0.5811f, 235.7f, 3.961f, << 608 },{ 5.408, 0.5811, 235.7, 3.961, 1.358
456 },{ 5.218f, 0.5828f, 245.0f, 3.838f, << 609 },{ 5.218, 0.5828, 245.0, 3.838, 1.25}
457 }; 610 };
458 611
459 G4double slow = (G4double)(a[i][0]); <<
460 <<
461 G4double x1 = (G4double)(a[i][1]); <<
462 G4double x2 = (G4double)(a[i][2]); <<
463 G4double x3 = (G4double)(a[i][3]); <<
464 G4double x4 = (G4double)(a[i][4]); <<
465 <<
466 // Free electron gas model 612 // Free electron gas model
467 if ( T < 0.001 ) { 613 if ( T < 0.001 ) {
468 G4double shigh = G4Log( 1.0 + x3*1000.0 + << 614 G4double slow = a[i][0] ;
469 ionloss = slow*shigh*std::sqrt(T*1000.0) << 615 G4double shigh = G4Log( 1.0 + a[i][3]*1000.0 + a[i][4]*0.001 )
>> 616 * a[i][2]*1000.0 ;
>> 617 ionloss = slow*shigh / (slow + shigh) ;
>> 618 ionloss *= sqrt(T*1000.0) ;
470 619
471 // Main parametrisation 620 // Main parametrisation
472 } else { 621 } else {
473 slow *= G4Exp(G4Log(T*1000.0)*x1); << 622 G4double slow = a[i][0] * G4Exp(G4Log(T*1000.0)*a[i][1]) ;
474 G4double shigh = G4Log( 1.0 + x3/T + x4*T << 623 G4double shigh = G4Log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ;
475 ionloss = slow*shigh / (slow + shigh) ; 624 ionloss = slow*shigh / (slow + shigh) ;
476 /* 625 /*
477 G4cout << "## " << i << ". T= " << T << " 626 G4cout << "## " << i << ". T= " << T << " slow= " << slow
478 << " a0= " << a[i][0] << " a1= " << 627 << " a0= " << a[i][0] << " a1= " << a[i][1]
479 << " shigh= " << shigh 628 << " shigh= " << shigh
480 << " dedx= " << ionloss << " q^2= " << 629 << " dedx= " << ionloss << " q^2= " << HeEffChargeSquare(z, T*MeV)
481 << G4endl; 630 << G4endl;
482 */ 631 */
483 } 632 }
484 ionloss = std::max(ionloss, 0.0); << 633 if ( ionloss < 0.0) { ionloss = 0.0; }
>> 634
>> 635 // He effective charge
>> 636 ionloss /= HeEffChargeSquare(z, T);
>> 637
485 return ionloss; 638 return ionloss;
486 } 639 }
487 640
488 //....oooOO0OOooo........oooOO0OOooo........oo 641 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
489 642
490 G4double G4BraggIonModel::HeDEDX(const G4Mater << 643 G4double G4BraggIonModel::DEDX(const G4Material* material,
491 const G4double << 644 G4double kineticEnergy)
492 { 645 {
493 // aEnergy is energy of alpha <<
494 G4double eloss = 0.0; 646 G4double eloss = 0.0;
495 // check DB 647 // check DB
496 if(material != currentMaterial) { 648 if(material != currentMaterial) {
497 currentMaterial = material; 649 currentMaterial = material;
498 baseMaterial = material->GetBaseMaterial() <<
499 ? material->GetBaseMaterial() : material <<
500 iPSTAR = -1; <<
501 iASTAR = -1; 650 iASTAR = -1;
502 iMolecula = -1; 651 iMolecula = -1;
503 iICRU90 = (nullptr != fICRU90) ? fICRU90-> << 652 if( !HasMaterial(material) ) { iASTAR = fASTAR->GetIndex(material); }
504 <<
505 if(iICRU90 < 0) { <<
506 if(isAlpha) { <<
507 iASTAR = fASTAR->GetIndex(baseMaterial); <<
508 if(iASTAR < 0) { iMolecula = HasMaterialForH <<
509 } else { <<
510 iPSTAR = fPSTAR->GetIndex(baseMaterial); <<
511 } <<
512 } <<
513 /* <<
514 G4cout << "%%% " <<material->GetName() << <<
515 << iMolecula << " iASTAR= " << iAS <<
516 << " iICRU90= " << iICRU90<< G4end <<
517 */ <<
518 } <<
519 // ICRU90 <<
520 if(iICRU90 >= 0) { <<
521 eloss = (isAlpha) <<
522 ? fICRU90->GetElectronicDEDXforAlpha(iIC <<
523 : fICRU90->GetElectronicDEDXforProton(iI <<
524 if(eloss > 0.0) { return eloss*material->G <<
525 } <<
526 // PSTAR parameterisation <<
527 if( iPSTAR >= 0 ) { <<
528 return fPSTAR->GetElectronicDEDX(iPSTAR, a <<
529 *material->GetDensity(); <<
530 } <<
531 // ASTAR <<
532 if( iASTAR >= 0 ) { <<
533 eloss = fASTAR->GetElectronicDEDX(iASTAR, <<
534 /* <<
535 G4cout << "ASTAR: E=" << aEnergy <<
536 << " dedx=" << eloss*material->GetDensity <<
537 << " " << particle->GetParticleName() << <<
538 */ <<
539 if(eloss > 0.0) { return eloss*material->G <<
540 } 653 }
541 654
542 const std::size_t numberOfElements = materia << 655 const G4int numberOfElements = material->GetNumberOfElements();
543 const G4ElementVector* theElmVector = materi <<
544 const G4double* theAtomicNumDensityVector = 656 const G4double* theAtomicNumDensityVector =
545 material->GetAtomicNumDensityVector(); << 657 material->GetAtomicNumDensityVector();
>> 658
>> 659 if( iASTAR >= 0 ) {
>> 660 G4double T = kineticEnergy*rateMassHe2p;
>> 661 G4int zeff = G4lrint(material->GetTotNbOfElectPerVolume()/
>> 662 material->GetTotNbOfAtomsPerVolume());
>> 663 return fASTAR->GetElectronicDEDX(iASTAR, T)*material->GetDensity()/
>> 664 HeEffChargeSquare(zeff, T/MeV);
>> 665
>> 666 } else if(iMolecula >= 0) {
546 667
547 // molecular data use proton stopping power << 668 eloss = StoppingPower(material, kineticEnergy)*
548 // element data from ICRU49 include data for << 669 material->GetDensity()/amu;
549 if(iMolecula >= 0) { <<
550 const G4double zeff = material->GetTotNbOf <<
551 material->GetTotNbOfAtomsPerVolume(); <<
552 heChargeSquare = HeEffChargeSquare(zeff, a <<
553 eloss = HeStoppingPower(aEnergy)*heChargeS <<
554 670
555 // pure material << 671 // pure material
556 } else if(1 == numberOfElements) { 672 } else if(1 == numberOfElements) {
557 673
558 const G4Element* element = (*theElmVector) << 674 G4double z = material->GetZ();
559 eloss = HeElectronicStoppingPower(element- << 675 eloss = ElectronicStoppingPower(z, kineticEnergy)
560 * (material->GetTotNbOfAtomsPerVolume()) << 676 * (material->GetTotNbOfAtomsPerVolume());
561 677
562 // Brugg's rule calculation 678 // Brugg's rule calculation
563 } else { 679 } else {
>> 680 const G4ElementVector* theElementVector =
>> 681 material->GetElementVector() ;
>> 682
564 // loop for the elements in the material 683 // loop for the elements in the material
565 for (std::size_t i=0; i<numberOfElements; << 684 for (G4int i=0; i<numberOfElements; i++)
566 const G4Element* element = (*theElmVecto << 685 {
567 eloss += HeElectronicStoppingPower(eleme << 686 const G4Element* element = (*theElementVector)[i] ;
568 * theAtomicNumDensityVector[i]; << 687 eloss += ElectronicStoppingPower(element->GetZ(), kineticEnergy)
>> 688 * theAtomicNumDensityVector[i];
569 } 689 }
570 } 690 }
571 return eloss*theZieglerFactor; 691 return eloss*theZieglerFactor;
572 } 692 }
573 693
574 //....oooOO0OOooo........oooOO0OOooo........oo 694 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
575 695
576 G4double << 696 G4double G4BraggIonModel::HeEffChargeSquare(G4double z,
577 G4BraggIonModel::HeEffChargeSquare(const G4dou << 697 G4double kinEnergyHeInMeV) const
578 const G4dou <<
579 { 698 {
580 // The aproximation of He effective charge f 699 // The aproximation of He effective charge from:
581 // J.F.Ziegler, J.P. Biersack, U. Littmark 700 // J.F.Ziegler, J.P. Biersack, U. Littmark
582 // The Stopping and Range of Ions in Matter, 701 // The Stopping and Range of Ions in Matter,
583 // Vol.1, Pergamon Press, 1985 702 // Vol.1, Pergamon Press, 1985
584 703
585 static const G4double c[6] = {0.2865, 0.126 704 static const G4double c[6] = {0.2865, 0.1266, -0.001429,
586 0.02402,-0.011 705 0.02402,-0.01135, 0.001475};
587 706
588 G4double e = std::max(0.0, G4Log(kinEnergyHe 707 G4double e = std::max(0.0, G4Log(kinEnergyHeInMeV*massFactor));
589 G4double x = c[0] ; 708 G4double x = c[0] ;
590 G4double y = 1.0 ; 709 G4double y = 1.0 ;
591 for (G4int i=1; i<6; ++i) { 710 for (G4int i=1; i<6; ++i) {
592 y *= e; << 711 y *= e ;
593 x += y * c[i]; << 712 x += y * c[i] ;
594 } 713 }
595 714
596 G4double w = 7.6 - e ; 715 G4double w = 7.6 - e ;
597 w = 1.0 + (0.007 + 0.00005*z) * G4Exp( -w*w 716 w = 1.0 + (0.007 + 0.00005*z) * G4Exp( -w*w ) ;
598 w = 4.0 * (1.0 - G4Exp(-x)) * w * w ; 717 w = 4.0 * (1.0 - G4Exp(-x)) * w * w ;
599 718
600 return w; 719 return w;
601 } 720 }
602 721
603 //....oooOO0OOooo........oooOO0OOooo........oo 722 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
604 723
605 724