| 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 // Modified by Z. Francis, S. Incerti to handl 27 // Modified by Z. Francis, S. Incerti to handle HZE
28 // && inverse rudd function sampling 26-10-201 28 // && inverse rudd function sampling 26-10-2010
29 // 29 //
30 // Rewitten by V.Ivanchenko 21.05.2023 30 // Rewitten by V.Ivanchenko 21.05.2023
31 // 31 //
32 32
33 #include "G4EmCorrections.hh" 33 #include "G4EmCorrections.hh"
34 #include "G4DNARuddIonisationExtendedModel.hh" 34 #include "G4DNARuddIonisationExtendedModel.hh"
35 #include "G4PhysicalConstants.hh" 35 #include "G4PhysicalConstants.hh"
36 #include "G4SystemOfUnits.hh" 36 #include "G4SystemOfUnits.hh"
37 #include "G4UAtomicDeexcitation.hh" 37 #include "G4UAtomicDeexcitation.hh"
38 #include "G4LossTableManager.hh" 38 #include "G4LossTableManager.hh"
39 #include "G4NistManager.hh" 39 #include "G4NistManager.hh"
40 #include "G4DNAChemistryManager.hh" 40 #include "G4DNAChemistryManager.hh"
41 #include "G4DNAMolecularMaterial.hh" 41 #include "G4DNAMolecularMaterial.hh"
42 42
43 #include "G4IonTable.hh" 43 #include "G4IonTable.hh"
44 #include "G4DNARuddAngle.hh" 44 #include "G4DNARuddAngle.hh"
45 #include "G4DeltaAngle.hh" 45 #include "G4DeltaAngle.hh"
46 #include "G4Exp.hh" 46 #include "G4Exp.hh"
47 #include "G4Log.hh" 47 #include "G4Log.hh"
48 #include "G4Pow.hh" 48 #include "G4Pow.hh"
49 #include "G4Alpha.hh" 49 #include "G4Alpha.hh"
50 #include "G4Proton.hh" 50 #include "G4Proton.hh"
>> 51 #include "G4AutoLock.hh"
51 52
52 //....oooOO0OOooo........oooOO0OOooo........oo 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53 54
54 G4DNACrossSectionDataSet* G4DNARuddIonisationE 55 G4DNACrossSectionDataSet* G4DNARuddIonisationExtendedModel::xsdata[] = {nullptr};
55 G4DNACrossSectionDataSet* G4DNARuddIonisationE 56 G4DNACrossSectionDataSet* G4DNARuddIonisationExtendedModel::xshelium = nullptr;
56 G4DNACrossSectionDataSet* G4DNARuddIonisationE 57 G4DNACrossSectionDataSet* G4DNARuddIonisationExtendedModel::xsalphaplus = nullptr;
57 const std::vector<G4double>* G4DNARuddIonisati 58 const std::vector<G4double>* G4DNARuddIonisationExtendedModel::fpWaterDensity = nullptr;
58 59
59 namespace 60 namespace
60 { 61 {
>> 62 G4Mutex ionDNAMutex = G4MUTEX_INITIALIZER;
61 const G4double scaleFactor = CLHEP::m*CLHEP: 63 const G4double scaleFactor = CLHEP::m*CLHEP::m;
62 const G4double tolerance = 1*CLHEP::eV; 64 const G4double tolerance = 1*CLHEP::eV;
63 const G4double Ry = 13.6*CLHEP::eV; 65 const G4double Ry = 13.6*CLHEP::eV;
>> 66 const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
64 67
65 // Following values provided by M. Dingfelde 68 // Following values provided by M. Dingfelder (priv. comm)
66 const G4double Bj[5] = {12.60*CLHEP::eV, 14. 69 const G4double Bj[5] = {12.60*CLHEP::eV, 14.70*CLHEP::eV, 18.40*CLHEP::eV,
67 32.20*CLHEP::eV, 539 70 32.20*CLHEP::eV, 539*CLHEP::eV};
68 } 71 }
69 72
70 //....oooOO0OOooo........oooOO0OOooo........oo 73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71 74
72 G4DNARuddIonisationExtendedModel::G4DNARuddIon 75 G4DNARuddIonisationExtendedModel::G4DNARuddIonisationExtendedModel(const G4ParticleDefinition*,
73 76 const G4String& nam)
74 : G4VEmModel(nam) 77 : G4VEmModel(nam)
75 { 78 {
76 fEmCorrections = G4LossTableManager::Instanc 79 fEmCorrections = G4LossTableManager::Instance()->EmCorrections();
77 fGpow = G4Pow::GetInstance(); 80 fGpow = G4Pow::GetInstance();
78 fLowestEnergy = 100*CLHEP::eV; 81 fLowestEnergy = 100*CLHEP::eV;
79 fLimitEnergy = 1*CLHEP::keV; 82 fLimitEnergy = 1*CLHEP::keV;
80 83
81 // Mark this model as "applicable" for atomi 84 // Mark this model as "applicable" for atomic deexcitation
82 SetDeexcitationFlag(true); 85 SetDeexcitationFlag(true);
83 86
84 // Define default angular generator 87 // Define default angular generator
85 SetAngularDistribution(new G4DNARuddAngle()) 88 SetAngularDistribution(new G4DNARuddAngle());
86 <<
87 if (nullptr == xshelium) { LoadData(); } <<
88 } 89 }
89 90
90 //....oooOO0OOooo........oooOO0OOooo........oo 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91 92
92 G4DNARuddIonisationExtendedModel::~G4DNARuddIo 93 G4DNARuddIonisationExtendedModel::~G4DNARuddIonisationExtendedModel()
93 { 94 {
94 if(isFirst) { 95 if(isFirst) {
95 for(auto & i : xsdata) { delete i; } 96 for(auto & i : xsdata) { delete i; }
96 } 97 }
97 } 98 }
98 99
99 //....oooOO0OOooo........oooOO0OOooo........oo 100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100 101
101 void G4DNARuddIonisationExtendedModel::LoadDat <<
102 { <<
103 // initialisation of static data once <<
104 isFirst = true; <<
105 G4String filename("dna/sigma_ionisation_h_ru <<
106 xsdata[0] = new G4DNACrossSectionDataSet(new <<
107 xsdata[0]->LoadData(filename); <<
108 <<
109 filename = "dna/sigma_ionisation_p_rudd"; <<
110 xsdata[1] = new G4DNACrossSectionDataSet(new <<
111 xsdata[1]->LoadData(filename); <<
112 <<
113 filename = "dna/sigma_ionisation_alphapluspl <<
114 xsdata[2] = new G4DNACrossSectionDataSet(new <<
115 xsdata[2]->LoadData(filename); <<
116 <<
117 filename = "dna/sigma_ionisation_li_rudd"; <<
118 xsdata[3] = new G4DNACrossSectionDataSet(new <<
119 xsdata[3]->LoadData(filename); <<
120 <<
121 filename = "dna/sigma_ionisation_be_rudd"; <<
122 xsdata[4] = new G4DNACrossSectionDataSet(new <<
123 xsdata[4]->LoadData(filename); <<
124 <<
125 filename = "dna/sigma_ionisation_b_rudd"; <<
126 xsdata[5] = new G4DNACrossSectionDataSet(new <<
127 xsdata[5]->LoadData(filename); <<
128 <<
129 filename = "dna/sigma_ionisation_c_rudd"; <<
130 xsdata[6] = new G4DNACrossSectionDataSet(new <<
131 xsdata[6]->LoadData(filename); <<
132 <<
133 filename = "dna/sigma_ionisation_n_rudd"; <<
134 xsdata[7] = new G4DNACrossSectionDataSet(new <<
135 xsdata[7]->LoadData(filename); <<
136 <<
137 filename = "dna/sigma_ionisation_o_rudd"; <<
138 xsdata[8] = new G4DNACrossSectionDataSet(new <<
139 xsdata[8]->LoadData(filename); <<
140 <<
141 filename = "dna/sigma_ionisation_si_rudd"; <<
142 xsdata[14] = new G4DNACrossSectionDataSet(ne <<
143 xsdata[14]->LoadData(filename); <<
144 <<
145 filename = "dna/sigma_ionisation_fe_rudd"; <<
146 xsdata[26] = new G4DNACrossSectionDataSet(ne <<
147 xsdata[26]->LoadData(filename); <<
148 <<
149 filename = "dna/sigma_ionisation_alphaplus_r <<
150 xsalphaplus = new G4DNACrossSectionDataSet(n <<
151 xsalphaplus->LoadData(filename); <<
152 <<
153 filename = "dna/sigma_ionisation_he_rudd"; <<
154 xshelium = new G4DNACrossSectionDataSet(new <<
155 xshelium->LoadData(filename); <<
156 <<
157 // to avoid possible threading problem fill <<
158 auto water = G4NistManager::Instance()->Find <<
159 fpWaterDensity = <<
160 G4DNAMolecularMaterial::Instance()->GetNum <<
161 } <<
162 <<
163 //....oooOO0OOooo........oooOO0OOooo........oo <<
164 <<
165 void G4DNARuddIonisationExtendedModel::Initial 102 void G4DNARuddIonisationExtendedModel::Initialise(const G4ParticleDefinition* p,
166 103 const G4DataVector&)
167 { 104 {
168 if (p != fParticle) { SetParticle(p); } << 105 if(p != fParticle) { SetParticle(p); }
169 106
170 // particle change object may be externally << 107 // initialisation of static data once
171 if (nullptr == fParticleChangeForGamma) { << 108 if(nullptr == xsdata[0]) {
172 fParticleChangeForGamma = GetParticleChang << 109 G4AutoLock l(&ionDNAMutex);
>> 110 if(nullptr == xsdata[0]) {
>> 111 isFirst = true;
>> 112 G4String filename("dna/sigma_ionisation_h_rudd");
>> 113 xsdata[0] = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, CLHEP::eV, scaleFactor);
>> 114 xsdata[0]->LoadData(filename);
>> 115
>> 116 filename = "dna/sigma_ionisation_p_rudd";
>> 117 xsdata[1] = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, CLHEP::eV, scaleFactor);
>> 118 xsdata[1]->LoadData(filename);
>> 119
>> 120 filename = "dna/sigma_ionisation_alphaplusplus_rudd";
>> 121 xsdata[2] = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, CLHEP::eV, scaleFactor);
>> 122 xsdata[2]->LoadData(filename);
>> 123
>> 124 filename = "dna/sigma_ionisation_li_rudd";
>> 125 xsdata[3] = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, CLHEP::eV, scaleFactor);
>> 126 xsdata[3]->LoadData(filename);
>> 127
>> 128 filename = "dna/sigma_ionisation_be_rudd";
>> 129 xsdata[4] = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, CLHEP::eV, scaleFactor);
>> 130 xsdata[4]->LoadData(filename);
>> 131
>> 132 filename = "dna/sigma_ionisation_b_rudd";
>> 133 xsdata[5] = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, CLHEP::eV, scaleFactor);
>> 134 xsdata[5]->LoadData(filename);
>> 135
>> 136 filename = "dna/sigma_ionisation_c_rudd";
>> 137 xsdata[6] = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, CLHEP::eV, scaleFactor);
>> 138 xsdata[6]->LoadData(filename);
>> 139
>> 140 filename = "dna/sigma_ionisation_n_rudd";
>> 141 xsdata[7] = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, CLHEP::eV, scaleFactor);
>> 142 xsdata[7]->LoadData(filename);
>> 143
>> 144 filename = "dna/sigma_ionisation_o_rudd";
>> 145 xsdata[8] = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, CLHEP::eV, scaleFactor);
>> 146 xsdata[8]->LoadData(filename);
>> 147
>> 148 filename = "dna/sigma_ionisation_si_rudd";
>> 149 xsdata[14] = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, CLHEP::eV, scaleFactor);
>> 150 xsdata[14]->LoadData(filename);
>> 151
>> 152 filename = "dna/sigma_ionisation_fe_rudd";
>> 153 xsdata[26] = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, CLHEP::eV, scaleFactor);
>> 154 xsdata[26]->LoadData(filename);
>> 155 filename = "dna/sigma_ionisation_alphaplus_rudd";
>> 156 xsalphaplus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, CLHEP::eV, scaleFactor);
>> 157 xsalphaplus->LoadData(filename);
>> 158
>> 159 filename = "dna/sigma_ionisation_he_rudd";
>> 160 xshelium = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, CLHEP::eV, scaleFactor);
>> 161 xshelium->LoadData(filename);
>> 162 }
>> 163 // to avoid possible threading problem fill this vector only once
>> 164 auto water = G4NistManager::Instance()->FindMaterial("G4_WATER");
>> 165 fpWaterDensity =
>> 166 G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(water);
>> 167
>> 168 l.unlock();
173 } 169 }
174 170
175 // initialisation once in each thread 171 // initialisation once in each thread
176 if (!isInitialised) { << 172 if(nullptr == fParticleChangeForGamma) {
177 isInitialised = true; << 173 fParticleChangeForGamma = GetParticleChangeForGamma();
178 const G4String& pname = fParticle->GetPart 174 const G4String& pname = fParticle->GetParticleName();
179 if (pname == "proton") { << 175 if(pname == "proton") {
180 idx = 1; 176 idx = 1;
181 xscurrent = xsdata[1]; 177 xscurrent = xsdata[1];
182 fElow = fLowestEnergy; 178 fElow = fLowestEnergy;
183 } else if (pname == "hydrogen") { << 179 } else if(pname == "hydrogen") {
184 idx = 0; 180 idx = 0;
185 xscurrent = xsdata[0]; 181 xscurrent = xsdata[0];
186 fElow = fLowestEnergy; 182 fElow = fLowestEnergy;
187 } else if (pname == "alpha") { << 183 } else if(pname == "alpha") {
188 idx = 1; 184 idx = 1;
189 xscurrent = xsdata[2]; 185 xscurrent = xsdata[2];
190 isHelium = true; 186 isHelium = true;
191 fElow = fLimitEnergy; 187 fElow = fLimitEnergy;
192 } else if (pname == "alpha+") { << 188 } else if(pname == "alpha+") {
193 idx = 1; 189 idx = 1;
194 isHelium = true; 190 isHelium = true;
195 xscurrent = xsalphaplus; 191 xscurrent = xsalphaplus;
196 fElow = fLimitEnergy; 192 fElow = fLimitEnergy;
197 // The following values are provided by 193 // The following values are provided by M. Dingfelder (priv. comm)
198 slaterEffectiveCharge[0]=2.0; 194 slaterEffectiveCharge[0]=2.0;
199 slaterEffectiveCharge[1]=2.0; 195 slaterEffectiveCharge[1]=2.0;
200 slaterEffectiveCharge[2]=2.0; 196 slaterEffectiveCharge[2]=2.0;
201 sCoefficient[0]=0.7; 197 sCoefficient[0]=0.7;
202 sCoefficient[1]=0.15; 198 sCoefficient[1]=0.15;
203 sCoefficient[2]=0.15; 199 sCoefficient[2]=0.15;
204 } else if (pname == "helium") { << 200 } else if(pname == "helium") {
205 idx = 0; 201 idx = 0;
206 isHelium = true; 202 isHelium = true;
207 fElow = fLimitEnergy; 203 fElow = fLimitEnergy;
208 xscurrent = xshelium; 204 xscurrent = xshelium;
209 slaterEffectiveCharge[0]=1.7; 205 slaterEffectiveCharge[0]=1.7;
210 slaterEffectiveCharge[1]=1.15; 206 slaterEffectiveCharge[1]=1.15;
211 slaterEffectiveCharge[2]=1.15; 207 slaterEffectiveCharge[2]=1.15;
212 sCoefficient[0]=0.5; 208 sCoefficient[0]=0.5;
213 sCoefficient[1]=0.25; 209 sCoefficient[1]=0.25;
214 sCoefficient[2]=0.25; 210 sCoefficient[2]=0.25;
215 } else { 211 } else {
216 isIon = true; 212 isIon = true;
217 idx = -1; <<
218 xscurrent = xsdata[1]; <<
219 fElow = fLowestEnergy; <<
220 } 213 }
221 // defined stationary mode 214 // defined stationary mode
222 statCode = G4EmParameters::Instance()->DNA 215 statCode = G4EmParameters::Instance()->DNAStationary();
223 216
224 // initialise atomic de-excitation 217 // initialise atomic de-excitation
225 fAtomDeexcitation = G4LossTableManager::In 218 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
226 219
227 if (verbose > 0) { 220 if (verbose > 0) {
228 G4cout << "### G4DNARuddIonisationExtend 221 G4cout << "### G4DNARuddIonisationExtendedModel::Initialise(..) " << pname
229 << "/n idx=" << idx << " isIon=" << << 222 << "/n idx=" << idx << " Amass=" << fAmass
230 << " isHelium=" << isHelium << G4endl; << 223 << " isIon=" << isIon << " isHelium=" << isHelium << G4endl;
231 } 224 }
232 } 225 }
233 } 226 }
234 227
235 //....oooOO0OOooo........oooOO0OOooo........oo 228 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
236 229
237 void G4DNARuddIonisationExtendedModel::SetPart 230 void G4DNARuddIonisationExtendedModel::SetParticle(const G4ParticleDefinition* p)
238 { 231 {
239 fParticle = p; 232 fParticle = p;
240 fMass = p->GetPDGMass(); 233 fMass = p->GetPDGMass();
241 fMassRate = (isIon) ? CLHEP::proton_mass_c2/ << 234 fAmass = p->GetAtomicMass();
>> 235
>> 236 // for generic ions idx is dynamic, -1 means that data for the ion does not exist
>> 237 if(isIon) {
>> 238 G4int i = p->GetAtomicNumber();
>> 239 idx = -1;
>> 240 if (i < RUDDZMAX && nullptr != xsdata[i]) {
>> 241 idx = i;
>> 242 fElow = fAmass*fLowestEnergy;
>> 243 }
>> 244 }
242 } 245 }
243 246
244 //....oooOO0OOooo........oooOO0OOooo........oo 247 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
245 248
246 G4double 249 G4double
247 G4DNARuddIonisationExtendedModel::CrossSection 250 G4DNARuddIonisationExtendedModel::CrossSectionPerVolume(const G4Material* material,
248 251 const G4ParticleDefinition* part,
249 252 G4double kinE,
250 253 G4double, G4double)
251 { 254 {
252 // check if model is applicable for given ma 255 // check if model is applicable for given material
253 G4double density = (material->GetIndex() < f 256 G4double density = (material->GetIndex() < fpWaterDensity->size())
254 ? (*fpWaterDensity)[material->GetIndex()] 257 ? (*fpWaterDensity)[material->GetIndex()] : 0.0;
255 if (0.0 == density) { return 0.0; } 258 if (0.0 == density) { return 0.0; }
256 259
257 // ion may be different 260 // ion may be different
258 if (fParticle != part) { SetParticle(part); 261 if (fParticle != part) { SetParticle(part); }
259 262
>> 263 // initilise mass rate
>> 264 fMassRate = 1.0;
>> 265
260 // ion shoud be stopped - check on kinetic e 266 // ion shoud be stopped - check on kinetic energy and not scaled energy
261 if (kinE < fLowestEnergy) { return DBL_MAX; 267 if (kinE < fLowestEnergy) { return DBL_MAX; }
262 268
263 G4double e = kinE*fMassRate; << 269 G4double sigma = 0.;
264 << 270
265 G4double sigma = (e > fElow) ? xscurrent->Fi << 271 // use ion table if available for given energy
266 : xscurrent->FindValue(fElow) * e / fElow; << 272 // for proton, hydrogen, alpha, alpha+, and helium no scaling to proton x-section
>> 273 if (idx == 0 || idx == 1) {
>> 274 sigma = (kinE > fElow) ? xscurrent->FindValue(kinE)
>> 275 : xscurrent->FindValue(fElow)*kinE/fElow;
>> 276
>> 277 // for ions with data above limit energy
>> 278 } else if (idx > 1) {
>> 279 sigma = (kinE > fElow) ? xsdata[idx]->FindValue(kinE)
>> 280 : xsdata[idx]->FindValue(fElow)*kinE/fElow;
267 281
268 if (idx == -1) { << 282 // scaling from proton
>> 283 } else {
>> 284 fMassRate = CLHEP::proton_mass_c2/fMass;
>> 285 G4double e = kinE*fMassRate;
>> 286 sigma = (e > fLowestEnergy) ? xsdata[1]->FindValue(e)
>> 287 : xsdata[1]->FindValue(fLowestEnergy)*e/fLowestEnergy;
269 sigma *= fEmCorrections->EffectiveChargeSq 288 sigma *= fEmCorrections->EffectiveChargeSquareRatio(part, material, kinE);
270 } 289 }
271 <<
272 sigma *= density; 290 sigma *= density;
273 291
274 if (verbose > 1) { 292 if (verbose > 1) {
275 G4cout << "G4DNARuddIonisationExtendedMode 293 G4cout << "G4DNARuddIonisationExtendedModel for " << part->GetParticleName()
276 << " Ekin(keV)=" << kinE/CLHEP::keV 294 << " Ekin(keV)=" << kinE/CLHEP::keV
277 << " sigma(cm^2)=" << sigma/CLHEP:: 295 << " sigma(cm^2)=" << sigma/CLHEP::cm2 << G4endl;
278 } 296 }
279 return sigma; 297 return sigma;
280 } 298 }
281 299
282 //....oooOO0OOooo........oooOO0OOooo........oo 300 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
283 301
284 void 302 void
285 G4DNARuddIonisationExtendedModel::SampleSecond 303 G4DNARuddIonisationExtendedModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
286 304 const G4MaterialCutsCouple* couple,
287 305 const G4DynamicParticle* dpart,
288 306 G4double, G4double)
289 { 307 {
290 const G4ParticleDefinition* pd = dpart->GetD 308 const G4ParticleDefinition* pd = dpart->GetDefinition();
291 if (fParticle != pd) { SetParticle(pd); } 309 if (fParticle != pd) { SetParticle(pd); }
292 310
293 // stop ion with energy below low energy lim 311 // stop ion with energy below low energy limit
294 G4double kinE = dpart->GetKineticEnergy(); 312 G4double kinE = dpart->GetKineticEnergy();
295 // ion shoud be stopped - check on kinetic e 313 // ion shoud be stopped - check on kinetic energy and not scaled energy
296 if (kinE <= fLowestEnergy) { 314 if (kinE <= fLowestEnergy) {
297 fParticleChangeForGamma->SetProposedKineti 315 fParticleChangeForGamma->SetProposedKineticEnergy(0.);
298 fParticleChangeForGamma->ProposeTrackStatu 316 fParticleChangeForGamma->ProposeTrackStatus(fStopButAlive);
299 fParticleChangeForGamma->ProposeLocalEnerg 317 fParticleChangeForGamma->ProposeLocalEnergyDeposit(kinE);
300 return; 318 return;
301 } 319 }
302 320
303 G4int shell = SelectShell(kinE*fMassRate); << 321 G4int shell = SelectShell(kinE);
304 G4double bindingEnergy = (useDNAWaterStructu 322 G4double bindingEnergy = (useDNAWaterStructure)
305 ? waterStructure.IonisationEnergy(shell) : 323 ? waterStructure.IonisationEnergy(shell) : Bj[shell];
306 324
307 //Si: additional protection if tcs interpola 325 //Si: additional protection if tcs interpolation method is modified
308 if (kinE < bindingEnergy) { return; } << 326 if (kinE < bindingEnergy) return;
309 << 327
310 G4double esec = SampleElectronEnergy(kinE, s << 328 G4double esec = SampleElectronEnergy(kinE, bindingEnergy, shell);
311 G4double esum = 0.0; 329 G4double esum = 0.0;
312 330
313 // sample deexcitation 331 // sample deexcitation
314 // here we assume that H2O electronic levels 332 // here we assume that H2O electronic levels are the same as Oxygen.
315 // this can be considered true with a rough 333 // this can be considered true with a rough 10% error in energy on K-shell,
316 G4int Z = 8; 334 G4int Z = 8;
317 G4ThreeVector deltaDir = 335 G4ThreeVector deltaDir =
318 GetAngularDistribution()->SampleDirectionF 336 GetAngularDistribution()->SampleDirectionForShell(dpart, esec, Z, shell, couple->GetMaterial());
319 337
320 // SI: only atomic deexcitation from K shell 338 // SI: only atomic deexcitation from K shell is considered
321 if(fAtomDeexcitation != nullptr && shell == 339 if(fAtomDeexcitation != nullptr && shell == 4) {
322 auto as = G4AtomicShellEnumerator(0); 340 auto as = G4AtomicShellEnumerator(0);
323 auto ashell = fAtomDeexcitation->GetAtomic 341 auto ashell = fAtomDeexcitation->GetAtomicShell(Z, as);
324 fAtomDeexcitation->GenerateParticles(fvect 342 fAtomDeexcitation->GenerateParticles(fvect, ashell, Z, 0, 0);
325 343
326 // compute energy sum from de-excitation 344 // compute energy sum from de-excitation
327 for (auto const & ptr : *fvect) { << 345 std::size_t nn = fvect->size();
328 esum += ptr->GetKineticEnergy(); << 346 for (std::size_t i=0; i<nn; ++i) {
>> 347 esum += (*fvect)[i]->GetKineticEnergy();
329 } 348 }
330 } 349 }
331 // check energy balance 350 // check energy balance
332 // remaining excitation energy of water mole 351 // remaining excitation energy of water molecule
333 G4double exc = bindingEnergy - esum; 352 G4double exc = bindingEnergy - esum;
334 353
335 // remaining projectile energy 354 // remaining projectile energy
336 G4double scatteredEnergy = kinE - bindingEne 355 G4double scatteredEnergy = kinE - bindingEnergy - esec;
337 if(scatteredEnergy < -tolerance || exc < -to 356 if(scatteredEnergy < -tolerance || exc < -tolerance) {
338 G4cout << "G4DNARuddIonisationExtendedMode 357 G4cout << "G4DNARuddIonisationExtendedModel::SampleSecondaries: "
339 << "negative final E(keV)=" << scat 358 << "negative final E(keV)=" << scatteredEnergy/CLHEP::keV << " Ein(keV)="
340 << kinE/CLHEP::keV << " " << pd->G 359 << kinE/CLHEP::keV << " " << pd->GetParticleName()
341 << " Edelta(keV)=" << esec/CLHEP::k 360 << " Edelta(keV)=" << esec/CLHEP::keV << " MeV, Exc(keV)=" << exc/CLHEP::keV
342 << G4endl; 361 << G4endl;
343 } 362 }
344 363
345 // projectile 364 // projectile
346 if (!statCode) { 365 if (!statCode) {
347 fParticleChangeForGamma->SetProposedKineti 366 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
348 fParticleChangeForGamma->ProposeLocalEnerg 367 fParticleChangeForGamma->ProposeLocalEnergyDeposit(exc);
349 } else { 368 } else {
350 fParticleChangeForGamma->SetProposedKineti 369 fParticleChangeForGamma->SetProposedKineticEnergy(kinE);
351 fParticleChangeForGamma->ProposeLocalEnerg 370 fParticleChangeForGamma->ProposeLocalEnergyDeposit(kinE - scatteredEnergy);
352 } 371 }
353 372
354 // delta-electron 373 // delta-electron
355 auto dp = new G4DynamicParticle(G4Electron: 374 auto dp = new G4DynamicParticle(G4Electron::Electron(), deltaDir, esec);
356 fvect->push_back(dp); 375 fvect->push_back(dp);
357 376
358 // create radical 377 // create radical
359 const G4Track* theIncomingTrack = fParticleC 378 const G4Track* theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
360 G4DNAChemistryManager::Instance()->CreateWat 379 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule, shell,
361 theIncomingTrack); 380 theIncomingTrack);
362 } 381 }
363 382
364 //....oooOO0OOooo........oooOO0OOooo........oo 383 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
365 384
366 G4int G4DNARuddIonisationExtendedModel::Select 385 G4int G4DNARuddIonisationExtendedModel::SelectShell(G4double e)
367 { 386 {
368 G4double sum = 0.0; 387 G4double sum = 0.0;
369 G4double xs; 388 G4double xs;
370 for (G4int i=0; i<5; ++i) { << 389 for(G4int i=0; i<5; ++i) {
371 auto ptr = xscurrent->GetComponent(i); << 390 if (idx == 0 || idx == 1) {
372 xs = (e > fElow) ? ptr->FindValue(e) : ptr << 391 auto ptr = xscurrent->GetComponent(i);
>> 392 xs = (e > fElow) ? ptr->FindValue(e) : ptr->FindValue(fElow)*e/fElow;
>> 393
>> 394 } else if (idx > 1) {
>> 395 auto ptr = xsdata[idx]->GetComponent(i);
>> 396 xs = (e > fElow) ? ptr->FindValue(e) : ptr->FindValue(fElow)*e/fElow;
>> 397
>> 398 } else {
>> 399 // use scaling from proton
>> 400 auto ptr = xsdata[1]->GetComponent(i);
>> 401 G4double x = e*fMassRate;
>> 402 xs = (x >= fLowestEnergy) ? ptr->FindValue(x)
>> 403 : ptr->FindValue(fLowestEnergy)*x/fLowestEnergy;
>> 404 }
373 sum += xs; 405 sum += xs;
374 fTemp[i] = sum; 406 fTemp[i] = sum;
375 } 407 }
376 sum *= G4UniformRand(); 408 sum *= G4UniformRand();
377 for (G4int i=0; i<5; ++i) { << 409 for(G4int i=0; i<5; ++i) {
378 if (sum <= fTemp[i]) { return i; } << 410 if(sum <= fTemp[i]) { return i; }
379 } 411 }
380 return 0; 412 return 0;
381 } 413 }
382 414
383 //....oooOO0OOooo........oooOO0OOooo........oo 415 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
384 416
385 G4double G4DNARuddIonisationExtendedModel::Max << 417 G4double G4DNARuddIonisationExtendedModel::SampleElectronEnergy(G4double kine,
>> 418 G4double eexc,
>> 419 G4int shell)
386 { 420 {
387 // kinematic limit 421 // kinematic limit
388 G4double tau = kine/fMass; 422 G4double tau = kine/fMass;
389 G4double gam = 1.0 + tau; <<
390 G4double emax = 2.0*CLHEP::electron_mass_c2* 423 G4double emax = 2.0*CLHEP::electron_mass_c2*tau*(tau + 2.0);
391 <<
392 // Initialisation of sampling <<
393 G4double A1, B1, C1, D1, E1, A2, B2, C2, D2; <<
394 if (shell == 4) { <<
395 //Data For Liquid Water K SHELL from Dingf <<
396 A1 = 1.25; <<
397 B1 = 0.5; <<
398 C1 = 1.00; <<
399 D1 = 1.00; <<
400 E1 = 3.00; <<
401 A2 = 1.10; <<
402 B2 = 1.30; <<
403 C2 = 1.00; <<
404 D2 = 0.00; <<
405 alphaConst = 0.66; <<
406 } else { <<
407 //Data For Liquid Water from Dingfelder (P <<
408 A1 = 1.02; <<
409 B1 = 82.0; <<
410 C1 = 0.45; <<
411 D1 = -0.80; <<
412 E1 = 0.38; <<
413 A2 = 1.07; <<
414 // Value provided by M. Dingfelder (priv. <<
415 B2 = 11.6; <<
416 C2 = 0.60; <<
417 D2 = 0.04; <<
418 alphaConst = 0.64; <<
419 } <<
420 bEnergy = Bj[shell]; <<
421 G4double v2 = 0.25*emax/(bEnergy*gam*gam); <<
422 v = std::sqrt(v2); <<
423 u = Ry/bEnergy; <<
424 wc = 4.*v2 - 2.*v - 0.25*u; <<
425 <<
426 G4double L1 = (C1 * fGpow->powA(v, D1)) / (1 <<
427 G4double L2 = C2 * fGpow->powA(v, D2); <<
428 G4double H1 = (A1 * G4Log(1. + v2)) / (v2 + <<
429 G4double H2 = (A2 / v2) + (B2 / (v2 * v2)); <<
430 <<
431 F1 = L1 + H1; <<
432 F2 = (L2 * H2) / (L2 + H2); <<
433 return emax; <<
434 } <<
435 <<
436 //....oooOO0OOooo........oooOO0OOooo........oo <<
437 <<
438 G4double G4DNARuddIonisationExtendedModel::Sam <<
439 <<
440 { <<
441 G4double emax = MaxEnergy(kine, shell); <<
442 // compute cumulative probability function 424 // compute cumulative probability function
443 G4double step = 1*CLHEP::eV; 425 G4double step = 1*CLHEP::eV;
444 auto nn = (G4int)(emax/step); << 426 auto nn = (G4int)(emax/step);
445 nn = std::min(std::max(nn, 10), 100); << 427 nn = std::max(nn, 10);
446 step = emax/(G4double)nn; 428 step = emax/(G4double)nn;
447 429
448 // find max probability 430 // find max probability
449 G4double pmax = ProbabilityFunction(kine, 0. << 431 G4double pmax = ProbabilityFunction(kine, 0.0, eexc, shell);
450 //G4cout << "## E(keV)=" << kine/keV << " em << 432 //G4cout << "E(keV)=" << kine/keV << " emax=" << emax/keV
451 // << " pmax(0)=" << pmax << " shell=" 433 // << " pmax(0)=" << pmax << " shell=" << shell << " nn=" << nn << G4endl;
452 434
>> 435 G4double e2 = 0.0; // backup energy
453 G4double e0 = 0.0; // energy with max probab 436 G4double e0 = 0.0; // energy with max probability
454 // 2 areas after point with max probability <<
455 G4double e1 = emax; <<
456 G4double e2 = emax; <<
457 G4double p1 = 0.0; <<
458 G4double p2 = 0.0; <<
459 const G4double f = 0.25; <<
460 <<
461 // find max probability <<
462 G4double e = 0.0; 437 G4double e = 0.0;
463 G4double p = 0.0; <<
464 for (G4int i=0; i<nn; ++i) { 438 for (G4int i=0; i<nn; ++i) {
465 e += step; 439 e += step;
466 p = ProbabilityFunction(kine, e, shell); << 440 G4double prob = ProbabilityFunction(kine, e, eexc, shell);
467 if (p > pmax) { << 441 if (prob < pmax) {
468 pmax = p; << 442 e2 = 2*e;
469 e0 = e; <<
470 } else { <<
471 break; 443 break;
472 } 444 }
>> 445 pmax = prob;
>> 446 e0 = e;
473 } 447 }
474 // increase step to be more effective << 448 //G4cout << " E0(keV)=" << e0/keV << " pmax=" << pmax << G4endl;
475 step *= 2.0; <<
476 // 2-nd area <<
477 for (G4int i=0; i<nn; ++i) { <<
478 e += step; <<
479 if (std::abs(e - emax) < step) { <<
480 e1 = emax; <<
481 break; <<
482 } <<
483 p = ProbabilityFunction(kine, e, shell); <<
484 if (p < f*pmax) { <<
485 p1 = p; <<
486 e1 = e; <<
487 break; <<
488 } <<
489 } <<
490 // 3-d area <<
491 if (e < emax) { <<
492 for (G4int i=0; i<nn; ++i) { <<
493 e += step; <<
494 if (std::abs(e - emax) < step) { <<
495 e2 = emax; <<
496 break; <<
497 } <<
498 p = ProbabilityFunction(kine, e, shell); <<
499 if (p < f*p1) { <<
500 p2 = p; <<
501 e2 = e; <<
502 break; <<
503 } <<
504 } <<
505 } <<
506 pmax *= 1.05; 449 pmax *= 1.05;
507 // regression method with 3 regions << 450 // regression method with two regions
508 G4double s0 = pmax*e1; << 451 G4double e1 = emax;
509 G4double s1 = s0 + p1 * (e2 - e1); << 452 G4double p1 = 0.0;
510 G4double s2 = s1 + p2 * (emax - e2); << 453 if (2*e0 < emax) {
511 s0 = (s0 == s1) ? 1.0 : s0 / s2; << 454 e1 = e0 + 0.25*(emax - e0);
512 s1 = (s1 == s2) ? 1.0 : s1 / s2; << 455 p1 = ProbabilityFunction(kine, e1, eexc, shell);
513 << 456 }
514 //G4cout << "pmax=" << pmax << " e1(keV)=" < << 457 G4double s2 = p1*(emax - e1);
515 // << " p2=" << p2 << " s0=" << s0 << " s1 << 458 s2 /= (s2 + e1*pmax);
>> 459 G4double s1 = 1.0 - s2;
516 460
517 // sampling 461 // sampling
518 G4int count = 0; 462 G4int count = 0;
519 G4double ymax, y, deltae; 463 G4double ymax, y, deltae;
520 for (G4int i = 0; i<100000; ++i) { 464 for (G4int i = 0; i<100000; ++i) {
521 G4double q = G4UniformRand(); 465 G4double q = G4UniformRand();
522 if (q <= s0) { << 466 if (q <= s1) {
523 ymax = pmax; 467 ymax = pmax;
524 deltae = e1 * q / s0; << 468 deltae = e1 * q / s1;
525 } else if (q <= s1) { <<
526 ymax = p1; <<
527 deltae = e1 + (e2 - e1) * (q - s0) / (s1 <<
528 } else { 469 } else {
529 ymax = p2; << 470 ymax = p1;
530 deltae = e2 + (emax - e2) * (q - s1) / ( << 471 deltae = e1 + (emax - e1)* (q - s1) / s2;
531 } 472 }
532 y = ProbabilityFunction(kine, deltae, shel << 473 y = ProbabilityFunction(kine, deltae, eexc, shell);
533 //G4cout << " " << i << ". deltae=" << 474 //G4cout << " " << i << ". deltae=" << deltae/CLHEP::keV
534 // << " y=" << y << " ymax=" << ymax << 475 // << " y=" << y << " ymax=" << ymax << G4endl;
535 if (y > ymax && count < 10) { 476 if (y > ymax && count < 10) {
536 ++count; 477 ++count;
537 G4cout << "G4DNARuddIonisationExtendedMo 478 G4cout << "G4DNARuddIonisationExtendedModel::SampleElectronEnergy warning: "
538 << fParticle->GetParticleName() << " E( 479 << fParticle->GetParticleName() << " E(keV)=" << kine/CLHEP::keV
539 << " Edelta(keV)=" << deltae/CLHEP::keV 480 << " Edelta(keV)=" << deltae/CLHEP::keV
540 << " y=" << y << " ymax=" << ymax << " 481 << " y=" << y << " ymax=" << ymax << " n=" << i << G4endl;
541 } 482 }
542 if (ymax * G4UniformRand() <= y) { << 483 if (ymax * G4UniformRand() < y) {
543 return deltae; 484 return deltae;
544 } 485 }
545 } 486 }
546 deltae = std::min(e0 + step, 0.5*emax); << 487 deltae = e2;
547 return deltae; 488 return deltae;
548 } 489 }
549 490
550 //....oooOO0OOooo........oooOO0OOooo........oo 491 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
551 492
552 G4double G4DNARuddIonisationExtendedModel::Pro 493 G4double G4DNARuddIonisationExtendedModel::ProbabilityFunction(G4double kine,
553 494 G4double deltae,
>> 495 G4double bindingEnergy,
554 496 G4int shell)
555 { 497 {
556 // Shells ids are 0 1 2 3 4 (4 is k shell) 498 // Shells ids are 0 1 2 3 4 (4 is k shell)
557 // !!Attention, "energyTransfer" here is the 499 // !!Attention, "energyTransfer" here is the energy transfered to the electron which means
558 // that the secondary kinetic en 500 // that the secondary kinetic energy is w = energyTransfer - bindingEnergy
559 // 501 //
560 // ds S F1(nu) + 502 // ds S F1(nu) + w * F2(nu)
561 // ---- = G(k) * ---- ----------------- 503 // ---- = G(k) * ---- -------------------------------------------
562 // dw Bj (1+w)^3 * [1 + e 504 // dw Bj (1+w)^3 * [1 + exp{alpha * (w - wc) / nu}]
563 // 505 //
564 // w is the secondary electron kinetic Energ 506 // w is the secondary electron kinetic Energy in eV
565 // 507 //
566 // All the other parameters can be found in 508 // All the other parameters can be found in Rudd's Papers
567 // 509 //
568 // M.Eugene Rudd, 1988, User-Friendly model 510 // M.Eugene Rudd, 1988, User-Friendly model for the energy distribution of
569 // electrons from protons or electron collis 511 // electrons from protons or electron collisions. Nucl. Tracks Rad. Meas.Vol 16 N0 2/3 pp 219-218
570 // 512 //
>> 513 G4double A1, B1, C1, D1, E1, A2, B2, C2, D2, alphaConst;
>> 514 if (shell == 4) {
>> 515 //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
>> 516 A1 = 1.25;
>> 517 B1 = 0.5;
>> 518 C1 = 1.00;
>> 519 D1 = 1.00;
>> 520 E1 = 3.00;
>> 521 A2 = 1.10;
>> 522 B2 = 1.30;
>> 523 C2 = 1.00;
>> 524 D2 = 0.00;
>> 525 alphaConst = 0.66;
>> 526 } else {
>> 527 //Data For Liquid Water from Dingfelder (Protons in Water)
>> 528 A1 = 1.02;
>> 529 B1 = 82.0;
>> 530 C1 = 0.45;
>> 531 D1 = -0.80;
>> 532 E1 = 0.38;
>> 533 A2 = 1.07;
>> 534 // Value provided by M. Dingfelder (priv. comm)
>> 535 B2 = 11.6;
>> 536 C2 = 0.60;
>> 537 D2 = 0.04;
>> 538 alphaConst = 0.64;
>> 539 }
>> 540 G4double bEnergy = Bj[shell];
571 G4double w = deltae/bEnergy; 541 G4double w = deltae/bEnergy;
>> 542 G4double u = Ry/bEnergy;
>> 543 G4double tau = kine/fMass;
>> 544 G4double gam = 1.0 + tau;;
>> 545
>> 546 G4double v2 = 0.5*CLHEP::electron_mass_c2*tau*(tau + 2.0)/(bEnergy*gam*gam);
>> 547 G4double v = std::sqrt(v2);
>> 548 G4double wc = 4.*v2 - 2.*v - 0.25*u;
>> 549
572 G4double x = alphaConst*(w - wc)/v; 550 G4double x = alphaConst*(w - wc)/v;
573 G4double y = (x > -15.) ? 1.0 + G4Exp(x) : 1 551 G4double y = (x > -15.) ? 1.0 + G4Exp(x) : 1.0;
574 552
575 G4double res = CorrectionFactor(kine, shell) << 553 G4double L1 = (C1 * fGpow->powA(v, D1)) / (1. + E1 * fGpow->powA(v, (D1 + 4.)));
>> 554 G4double L2 = C2 * fGpow->powA(v, D2);
>> 555 G4double H1 = (A1 * G4Log(1. + v2)) / (v2 + (B1 / v2));
>> 556 G4double H2 = (A2 / v2) + (B2 / (v2 * v2));
>> 557
>> 558 G4double F1 = L1 + H1;
>> 559 G4double F2 = (L2 * H2) / (L2 + H2);
>> 560
>> 561 G4double res = CorrectionFactor(kine, shell) * (F1 + w*F2) * Gj[shell] /
576 (fGpow->powN((1. + w)/u, 3) * y); 562 (fGpow->powN((1. + w)/u, 3) * y);
577 563
578 if (isHelium) { << 564 if(isHelium) {
579 G4double energyTransfer = deltae + bEnergy << 565 G4double energyTransfer = deltae + bindingEnergy;
580 G4double Zeff = 2.0 - 566 G4double Zeff = 2.0 -
581 (sCoefficient[0] * S_1s(kine, energyTran 567 (sCoefficient[0] * S_1s(kine, energyTransfer, slaterEffectiveCharge[0], 1.) +
582 sCoefficient[1] * S_2s(kine, energyTran 568 sCoefficient[1] * S_2s(kine, energyTransfer, slaterEffectiveCharge[1], 2.) +
583 sCoefficient[2] * S_2p(kine, energyTran 569 sCoefficient[2] * S_2p(kine, energyTransfer, slaterEffectiveCharge[2], 2.) );
584 570
585 res *= Zeff * Zeff; 571 res *= Zeff * Zeff;
586 } 572 }
587 return std::max(res, 0.0); 573 return std::max(res, 0.0);
588 } 574 }
589 575
590 //....oooOO0OOooo........oooOO0OOooo........oo 576 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
591 577
592 G4double G4DNARuddIonisationExtendedModel::Com 578 G4double G4DNARuddIonisationExtendedModel::ComputeProbabilityFunction(
593 const G4ParticleDefinition* p, G4doub 579 const G4ParticleDefinition* p, G4double e, G4double deltae, G4int shell)
594 { 580 {
595 if (fParticle != p) { SetParticle(p); } 581 if (fParticle != p) { SetParticle(p); }
596 MaxEnergy(e, shell); << 582 G4double bEnergy = (useDNAWaterStructure)
597 return ProbabilityFunction(e, deltae, shell) << 583 ? waterStructure.IonisationEnergy(shell) : Bj[shell];
>> 584 return ProbabilityFunction(e, deltae, bEnergy, shell);
598 } 585 }
599 586
600 //....oooOO0OOooo........oooOO0OOooo........oo 587 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
601 588
602 G4double G4DNARuddIonisationExtendedModel::S_1 589 G4double G4DNARuddIonisationExtendedModel::S_1s(G4double kine,
603 590 G4double energyTransfer,
604 591 G4double slaterEffCharge,
605 592 G4double shellNumber)
606 { 593 {
607 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2) 594 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
608 // Dingfelder, in Chattanooga 2005 proceedin 595 // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
609 596
610 G4double r = Rh(kine, energyTransfer, slater 597 G4double r = Rh(kine, energyTransfer, slaterEffCharge, shellNumber);
611 G4double value = 1. - G4Exp(-2 * r) * ( ( 2. 598 G4double value = 1. - G4Exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
612 return value; 599 return value;
613 } 600 }
614 601
615 //....oooOO0OOooo........oooOO0OOooo........oo 602 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
616 603
617 G4double G4DNARuddIonisationExtendedModel::S_2 604 G4double G4DNARuddIonisationExtendedModel::S_2s(G4double kine,
618 605 G4double energyTransfer,
619 606 G4double slaterEffCharge,
620 607 G4double shellNumber)
621 { 608 {
622 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4) 609 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
623 // Dingfelder, in Chattanooga 2005 proceedin 610 // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
624 611
625 G4double r = Rh(kine, energyTransfer, slater 612 G4double r = Rh(kine, energyTransfer, slaterEffCharge, shellNumber);
626 G4double value = 613 G4double value =
627 1. - G4Exp(-2 * r) * (((2. * r * r + 2.) * 614 1. - G4Exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
628 615
629 return value; 616 return value;
630 } 617 }
631 618
632 //....oooOO0OOooo........oooOO0OOooo........oo 619 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
633 620
634 G4double G4DNARuddIonisationExtendedModel::S_2 621 G4double G4DNARuddIonisationExtendedModel::S_2p(G4double kine,
635 622 G4double energyTransfer,
636 623 G4double slaterEffCharge,
637 624 G4double shellNumber)
638 { 625 {
639 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^ 626 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
640 // Dingfelder, in Chattanooga 2005 proceedin 627 // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
641 628
642 G4double r = Rh(kine, energyTransfer, slater 629 G4double r = Rh(kine, energyTransfer, slaterEffCharge, shellNumber);
643 G4double value = 630 G4double value =
644 1. - G4Exp(-2 * r) * (((( 2./3. * r + 4./3 631 1. - G4Exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
645 632
646 return value; 633 return value;
647 } 634 }
648 635
649 //....oooOO0OOooo........oooOO0OOooo........oo 636 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
650 637
651 G4double G4DNARuddIonisationExtendedModel::Rh( 638 G4double G4DNARuddIonisationExtendedModel::Rh(G4double ekin, G4double etrans,
652 639 G4double q, G4double shell)
653 { 640 {
654 // The following values are provided by M. D 641 // The following values are provided by M. Dingfelder (priv. comm)
655 // Dingfelder, in Chattanooga 2005 proceedin 642 // Dingfelder, in Chattanooga 2005 proceedings, p 4
656 643
657 G4double escaled = CLHEP::electron_mass_c2/f 644 G4double escaled = CLHEP::electron_mass_c2/fMass * ekin;
658 const G4double H = 13.60569172 * CLHEP::eV; 645 const G4double H = 13.60569172 * CLHEP::eV;
659 G4double value = 2.0*std::sqrt(escaled / H)* 646 G4double value = 2.0*std::sqrt(escaled / H)*q*H /(etrans*shell);
660 647
661 return value; 648 return value;
662 } 649 }
663 650
664 //....oooOO0OOooo........oooOO0OOooo........oo 651 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
665 652
666 G4double 653 G4double
667 G4DNARuddIonisationExtendedModel::CorrectionFa 654 G4DNARuddIonisationExtendedModel::CorrectionFactor(G4double kine, G4int shell)
668 { 655 {
669 // ZF Shortened 656 // ZF Shortened
670 G4double res = 1.0; 657 G4double res = 1.0;
671 if (shell < 4 && 0 != idx) { 658 if (shell < 4 && 0 != idx) {
672 const G4double ln10 = fGpow->logZ(10); 659 const G4double ln10 = fGpow->logZ(10);
673 G4double x = 2.0*((G4Log(kine/CLHEP::eV)/l 660 G4double x = 2.0*((G4Log(kine/CLHEP::eV)/ln10) - 4.2);
674 // The following values are provided by M. 661 // The following values are provided by M. Dingfelder (priv. comm)
675 res = 0.6/(1.0 + G4Exp(x)) + 0.9; 662 res = 0.6/(1.0 + G4Exp(x)) + 0.9;
676 } 663 }
677 return res; 664 return res;
678 } 665 }
679 666
680 //....oooOO0OOooo........oooOO0OOooo........oo 667 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
681 668