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