| 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: G4DNABornIonisationModel2.cc 87631 2014-12-14 12:42:05Z matkara $
26 // 27 //
27 28
28 #include "G4DNABornIonisationModel2.hh" 29 #include "G4DNABornIonisationModel2.hh"
29 #include "G4PhysicalConstants.hh" 30 #include "G4PhysicalConstants.hh"
30 #include "G4SystemOfUnits.hh" 31 #include "G4SystemOfUnits.hh"
31 #include "G4UAtomicDeexcitation.hh" 32 #include "G4UAtomicDeexcitation.hh"
32 #include "G4LossTableManager.hh" 33 #include "G4LossTableManager.hh"
33 #include "G4DNAChemistryManager.hh" 34 #include "G4DNAChemistryManager.hh"
34 #include "G4DNAMolecularMaterial.hh" 35 #include "G4DNAMolecularMaterial.hh"
35 #include "G4DNABornAngle.hh" 36 #include "G4DNABornAngle.hh"
36 #include "G4DeltaAngle.hh" 37 #include "G4DeltaAngle.hh"
37 #include "G4Exp.hh" 38 #include "G4Exp.hh"
38 39
39 //....oooOO0OOooo........oooOO0OOooo........oo 40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
40 41
41 using namespace std; 42 using namespace std;
42 43
43 //....oooOO0OOooo........oooOO0OOooo........oo 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44 45
45 G4DNABornIonisationModel2::G4DNABornIonisation 46 G4DNABornIonisationModel2::G4DNABornIonisationModel2(const G4ParticleDefinition*,
46 47 const G4String& nam) :
47 G4VEmModel(nam) << 48 G4VEmModel(nam), isInitialised(false)
48 { 49 {
49 verboseLevel = 0; 50 verboseLevel = 0;
50 // Verbosity scale: 51 // Verbosity scale:
51 // 0 = nothing 52 // 0 = nothing
52 // 1 = warning for energy non-conservation 53 // 1 = warning for energy non-conservation
53 // 2 = details of energy budget 54 // 2 = details of energy budget
54 // 3 = calculation of cross sections, file o 55 // 3 = calculation of cross sections, file openings, sampling of atoms
55 // 4 = entering in methods 56 // 4 = entering in methods
56 57
57 if (verboseLevel > 0) 58 if (verboseLevel > 0)
58 { 59 {
59 G4cout << "Born ionisation model is constr 60 G4cout << "Born ionisation model is constructed " << G4endl;
60 } 61 }
61 62
62 // Mark this model as "applicable" for atomi 63 // Mark this model as "applicable" for atomic deexcitation
63 64
64 SetDeexcitationFlag(true); 65 SetDeexcitationFlag(true);
65 fAtomDeexcitation = nullptr; << 66 fAtomDeexcitation = 0;
66 fParticleChangeForGamma = nullptr; << 67 fParticleChangeForGamma = 0;
67 fpMolWaterDensity = nullptr; << 68 fpMolWaterDensity = 0;
68 fTableData = nullptr; << 69 fTableData = 0;
69 fLowEnergyLimit = 0; 70 fLowEnergyLimit = 0;
70 fHighEnergyLimit = 0; 71 fHighEnergyLimit = 0;
71 fParticleDef = nullptr; << 72 fParticleDef = 0;
72 73
73 // Define default angular generator 74 // Define default angular generator
74 75
75 SetAngularDistribution(new G4DNABornAngle()) 76 SetAngularDistribution(new G4DNABornAngle());
76 77
77 // Selection of computation method 78 // Selection of computation method
78 79
79 fasterCode = false; 80 fasterCode = false;
80 81
81 // Selection of stationary mode 82 // Selection of stationary mode
82 83
83 statCode = false; 84 statCode = false;
84 85
85 // Selection of SP scaling 86 // Selection of SP scaling
86 87
87 spScaling = true; 88 spScaling = true;
88 } 89 }
89 90
90 //....oooOO0OOooo........oooOO0OOooo........oo 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91 92
92 G4DNABornIonisationModel2::~G4DNABornIonisatio 93 G4DNABornIonisationModel2::~G4DNABornIonisationModel2()
93 { 94 {
94 // Cross section 95 // Cross section
95 96
96 << 97 if (fTableData)
97 delete fTableData; 98 delete fTableData;
98 99
99 // Final state 100 // Final state
100 101
101 fVecm.clear(); 102 fVecm.clear();
102 } 103 }
103 104
104 //....oooOO0OOooo........oooOO0OOooo........oo 105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
105 106
106 void G4DNABornIonisationModel2::Initialise(con 107 void G4DNABornIonisationModel2::Initialise(const G4ParticleDefinition* particle,
107 con 108 const G4DataVector& /*cuts*/)
108 { 109 {
109 110
110 if (verboseLevel > 3) 111 if (verboseLevel > 3)
111 { 112 {
112 G4cout << "Calling G4DNABornIonisationMode 113 G4cout << "Calling G4DNABornIonisationModel2::Initialise()" << G4endl;
113 } 114 }
114 115
115 if(fParticleDef != nullptr && particle != fP << 116 if(fParticleDef != 0 && particle != fParticleDef)
116 { 117 {
117 G4ExceptionDescription description; 118 G4ExceptionDescription description;
118 description << "You are trying to initiali 119 description << "You are trying to initialized G4DNABornIonisationModel2 "
119 "for particle " 120 "for particle "
120 << particle->GetParticleName() 121 << particle->GetParticleName()
121 << G4endl; 122 << G4endl;
122 description << "G4DNABornIonisationModel2 << 123 description << "G4DNABornIonisationModel2 was already initiliased "
123 "for particle:" << fParticleDef->GetPartic 124 "for particle:" << fParticleDef->GetParticleName() << G4endl;
124 G4Exception("G4DNABornIonisationModel2::In 125 G4Exception("G4DNABornIonisationModel2::Initialise","bornIonInit",
125 FatalException,description); 126 FatalException,description);
126 } 127 }
127 128
128 fParticleDef = particle; 129 fParticleDef = particle;
129 130
130 // Energy limits 131 // Energy limits
131 const char* path = G4FindDataDir("G4LEDATA") << 132 char *path = getenv("G4LEDATA");
132 133
133 // *** 134 // ***
134 135
135 G4String particleName = particle->GetParticl 136 G4String particleName = particle->GetParticleName();
136 std::ostringstream fullFileName; 137 std::ostringstream fullFileName;
137 fullFileName << path; 138 fullFileName << path;
138 139
139 if(particleName == "e-") 140 if(particleName == "e-")
140 { 141 {
141 fTableFile = "dna/sigma_ionisation_e_born" 142 fTableFile = "dna/sigma_ionisation_e_born";
142 fLowEnergyLimit = 11.*eV; 143 fLowEnergyLimit = 11.*eV;
143 fHighEnergyLimit = 1.*MeV; 144 fHighEnergyLimit = 1.*MeV;
144 145
145 if (fasterCode) 146 if (fasterCode)
146 { 147 {
147 fullFileName << "/dna/sigmadiff_cumulate 148 fullFileName << "/dna/sigmadiff_cumulated_ionisation_e_born_hp.dat";
148 } 149 }
149 else 150 else
150 { 151 {
151 fullFileName << "/dna/sigmadiff_ionisati 152 fullFileName << "/dna/sigmadiff_ionisation_e_born.dat";
152 } 153 }
153 } 154 }
154 else if(particleName == "proton") 155 else if(particleName == "proton")
155 { 156 {
156 fTableFile = "dna/sigma_ionisation_p_born" 157 fTableFile = "dna/sigma_ionisation_p_born";
157 fLowEnergyLimit = 500. * keV; 158 fLowEnergyLimit = 500. * keV;
158 fHighEnergyLimit = 100. * MeV; 159 fHighEnergyLimit = 100. * MeV;
159 160
160 if (fasterCode) 161 if (fasterCode)
161 { 162 {
162 fullFileName << "/dna/sigmadiff_cumulate 163 fullFileName << "/dna/sigmadiff_cumulated_ionisation_p_born_hp.dat";
163 } 164 }
164 else 165 else
165 { 166 {
166 fullFileName << "/dna/sigmadiff_ionisati 167 fullFileName << "/dna/sigmadiff_ionisation_p_born.dat";
167 } 168 }
168 } 169 }
169 170
170 // Cross section 171 // Cross section
171 172
172 G4double scaleFactor = (1.e-22 / 3.343) * m* 173 G4double scaleFactor = (1.e-22 / 3.343) * m*m;
173 fTableData = new G4DNACrossSectionDataSet(ne 174 fTableData = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
174 fTableData->LoadData(fTableFile); 175 fTableData->LoadData(fTableFile);
175 176
176 // Final state 177 // Final state
177 178
178 std::ifstream diffCrossSection(fullFileName. 179 std::ifstream diffCrossSection(fullFileName.str().c_str());
179 180
180 if (!diffCrossSection) 181 if (!diffCrossSection)
181 { 182 {
182 G4ExceptionDescription description; 183 G4ExceptionDescription description;
183 description << "Missing data file:" << G4e 184 description << "Missing data file:" << G4endl << fullFileName.str() << G4endl;
184 G4Exception("G4DNABornIonisationModel2::In 185 G4Exception("G4DNABornIonisationModel2::Initialise","em0003",
185 FatalException,description); 186 FatalException,description);
186 } 187 }
187 188
188 // Clear the arrays for re-initialization ca 189 // Clear the arrays for re-initialization case (MT mode)
189 // March 25th, 2014 - Vaclav Stepan, Sebasti 190 // March 25th, 2014 - Vaclav Stepan, Sebastien Incerti
190 191
191 fTdummyVec.clear(); 192 fTdummyVec.clear();
192 fVecm.clear(); 193 fVecm.clear();
193 194
194 for (int j=0; j<5; j++) 195 for (int j=0; j<5; j++)
195 { 196 {
196 fProbaShellMap[j].clear(); 197 fProbaShellMap[j].clear();
197 fDiffCrossSectionData[j].clear(); 198 fDiffCrossSectionData[j].clear();
198 fNrjTransfData[j].clear(); 199 fNrjTransfData[j].clear();
199 } 200 }
200 201
201 // 202 //
202 203
203 fTdummyVec.push_back(0.); 204 fTdummyVec.push_back(0.);
204 while(!diffCrossSection.eof()) 205 while(!diffCrossSection.eof())
205 { 206 {
206 G4double tDummy; 207 G4double tDummy;
207 G4double eDummy; 208 G4double eDummy;
208 diffCrossSection>>tDummy>>eDummy; 209 diffCrossSection>>tDummy>>eDummy;
209 if (tDummy != fTdummyVec.back()) fTdummyVe 210 if (tDummy != fTdummyVec.back()) fTdummyVec.push_back(tDummy);
210 211
211 G4double tmp; 212 G4double tmp;
212 for (int j=0; j<5; j++) 213 for (int j=0; j<5; j++)
213 { 214 {
214 diffCrossSection>> tmp; 215 diffCrossSection>> tmp;
215 216
216 fDiffCrossSectionData[j][tDummy][eDummy] 217 fDiffCrossSectionData[j][tDummy][eDummy] = tmp;
217 218
218 if (fasterCode) 219 if (fasterCode)
219 { 220 {
220 fNrjTransfData[j][tDummy][fDiffCrossSe 221 fNrjTransfData[j][tDummy][fDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
221 fProbaShellMap[j][tDummy].push_back(fD 222 fProbaShellMap[j][tDummy].push_back(fDiffCrossSectionData[j][tDummy][eDummy]);
222 } 223 }
223 224
224 // SI - only if eof is not reached 225 // SI - only if eof is not reached
225 if (!diffCrossSection.eof() && !fasterCo 226 if (!diffCrossSection.eof() && !fasterCode) fDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
226 227
227 if (!fasterCode) fVecm[tDummy].push_back 228 if (!fasterCode) fVecm[tDummy].push_back(eDummy);
228 229
229 } 230 }
230 } 231 }
231 232
232 // 233 //
233 SetLowEnergyLimit(fLowEnergyLimit); 234 SetLowEnergyLimit(fLowEnergyLimit);
234 SetHighEnergyLimit(fHighEnergyLimit); 235 SetHighEnergyLimit(fHighEnergyLimit);
235 236
236 if( verboseLevel>0 ) 237 if( verboseLevel>0 )
237 { 238 {
238 G4cout << "Born ionisation model is initia 239 G4cout << "Born ionisation model is initialized " << G4endl
239 << "Energy range: " 240 << "Energy range: "
240 << LowEnergyLimit() / eV << " eV - " 241 << LowEnergyLimit() / eV << " eV - "
241 << HighEnergyLimit() / keV << " keV for " 242 << HighEnergyLimit() / keV << " keV for "
242 << particle->GetParticleName() 243 << particle->GetParticleName()
243 << G4endl; 244 << G4endl;
244 } 245 }
245 246
246 // Initialize water density pointer 247 // Initialize water density pointer
247 248
248 fpMolWaterDensity = G4DNAMolecularMaterial:: 249 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->
249 GetNumMolPerVolTableFor(G4Material::GetMater 250 GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
250 251
251 // AD 252 // AD
252 253
253 fAtomDeexcitation = G4LossTableManager::Inst 254 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
254 255
255 if (isInitialised) 256 if (isInitialised)
256 { return;} 257 { return;}
257 fParticleChangeForGamma = GetParticleChangeF 258 fParticleChangeForGamma = GetParticleChangeForGamma();
258 isInitialised = true; 259 isInitialised = true;
259 } 260 }
260 261
261 //....oooOO0OOooo........oooOO0OOooo........oo 262 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
262 263
263 G4double G4DNABornIonisationModel2::CrossSecti 264 G4double G4DNABornIonisationModel2::CrossSectionPerVolume(const G4Material* material,
264 265 const G4ParticleDefinition* particleDefinition,
265 266 G4double ekin,
266 267 G4double,
267 268 G4double)
268 { 269 {
269 if (verboseLevel > 3) 270 if (verboseLevel > 3)
270 { 271 {
271 G4cout << "Calling CrossSectionPerVolume() 272 G4cout << "Calling CrossSectionPerVolume() of G4DNABornIonisationModel2"
272 << G4endl; 273 << G4endl;
273 } 274 }
274 275
275 if (particleDefinition != fParticleDef) retu 276 if (particleDefinition != fParticleDef) return 0;
276 277
277 // Calculate total cross section for model 278 // Calculate total cross section for model
278 279
279 G4double sigma=0; 280 G4double sigma=0;
280 281
281 G4double waterDensity = (*fpMolWaterDensity) 282 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
282 283
283 if (ekin >= fLowEnergyLimit && ekin <= fHigh << 284 if(waterDensity!= 0.0)
284 { 285 {
285 sigma = fTableData->FindValue(ekin); << 286 if (ekin >= fLowEnergyLimit && ekin < fHighEnergyLimit)
>> 287 {
>> 288 sigma = fTableData->FindValue(ekin);
286 289
287 // ICRU49 electronic SP scaling - ZF, SI << 290 // ICRU49 electronic SP scaling - ZF, SI
288 291
289 if (particleDefinition == G4Proton::Proton << 292 if (particleDefinition == G4Proton::ProtonDefinition() && ekin < 70*MeV && spScaling)
290 { << 293 {
291 G4double A = 1.39241700556072800000E-009 << 294 G4double A = 1.39241700556072800000E-009 ;
292 G4double B = -8.52610412942622630000E-00 << 295 G4double B = -8.52610412942622630000E-002 ;
293 sigma = sigma * G4Exp(A*(ekin/eV)+B); << 296 sigma = sigma * G4Exp(A*(ekin/eV)+B);
>> 297 }
>> 298 //
294 } 299 }
295 // <<
296 } <<
297 300
298 if (verboseLevel > 2) << 301 if (verboseLevel > 2)
299 { << 302 {
300 G4cout << "_______________________________ << 303 G4cout << "__________________________________" << G4endl;
301 G4cout << "G4DNABornIonisationModel2 - XS << 304 G4cout << "G4DNABornIonisationModel2 - XS INFO START" << G4endl;
302 G4cout << "Kinetic energy(eV)=" << ekin/eV << 305 G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
303 G4cout << "Cross section per water molecul << 306 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
304 G4cout << "Cross section per water molecul << 307 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
305 G4cout << "G4DNABornIonisationModel2 - XS << 308 G4cout << "G4DNABornIonisationModel2 - XS INFO END" << G4endl;
306 } << 309 }
>> 310 } // if (waterMaterial)
307 311
308 return sigma*waterDensity; 312 return sigma*waterDensity;
309 } 313 }
310 314
311 //....oooOO0OOooo........oooOO0OOooo........oo 315 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
312 316
313 void G4DNABornIonisationModel2::SampleSecondar 317 void G4DNABornIonisationModel2::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
314 318 const G4MaterialCutsCouple* couple,
315 319 const G4DynamicParticle* particle,
316 320 G4double,
317 321 G4double)
318 { 322 {
319 323
320 if (verboseLevel > 3) 324 if (verboseLevel > 3)
321 { 325 {
322 G4cout << "Calling SampleSecondaries() of 326 G4cout << "Calling SampleSecondaries() of G4DNABornIonisationModel2"
323 << G4endl; 327 << G4endl;
324 } 328 }
325 329
326 G4double k = particle->GetKineticEnergy(); 330 G4double k = particle->GetKineticEnergy();
327 331
328 if (k >= fLowEnergyLimit && k <= fHighEnergy << 332 if (k >= fLowEnergyLimit && k < fHighEnergyLimit)
329 { 333 {
330 G4ParticleMomentum primaryDirection = part 334 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
331 G4double particleMass = particle->GetDefin 335 G4double particleMass = particle->GetDefinition()->GetPDGMass();
332 G4double totalEnergy = k + particleMass; 336 G4double totalEnergy = k + particleMass;
333 G4double pSquare = k * (totalEnergy + part 337 G4double pSquare = k * (totalEnergy + particleMass);
334 G4double totalMomentum = std::sqrt(pSquare 338 G4double totalMomentum = std::sqrt(pSquare);
335 339
336 G4int ionizationShell = 0; 340 G4int ionizationShell = 0;
337 341
338 if (!fasterCode) ionizationShell = RandomS 342 if (!fasterCode) ionizationShell = RandomSelect(k);
339 343
340 // SI: The following protection is necessa 344 // SI: The following protection is necessary to avoid infinite loops :
341 // sigmadiff_ionisation_e_born.dat has no 345 // sigmadiff_ionisation_e_born.dat has non zero partial xs at 18 eV for shell 3 (ionizationShell ==2)
342 // sigmadiff_cumulated_ionisation_e_born. 346 // sigmadiff_cumulated_ionisation_e_born.dat has zero cumulated partial xs at 18 eV for shell 3 (ionizationShell ==2)
343 // this is due to the fact that the max a 347 // this is due to the fact that the max allowed transfered energy is (18+10.79)/2=17.025 eV and only transfered energies
344 // strictly above this value have non zer 348 // strictly above this value have non zero partial xs in sigmadiff_ionisation_e_born.dat (starting at trans = 17.12 eV)
345 349
346 if (fasterCode) 350 if (fasterCode)
347 do 351 do
348 { 352 {
349 ionizationShell = RandomSelect(k); 353 ionizationShell = RandomSelect(k);
350 } while (k<19*eV && ionizationShell==2 && << 354 }while (k<19*eV && ionizationShell==2 && particle->GetDefinition()==G4Electron::ElectronDefinition());
351 355
352 G4double secondaryKinetic=-1000*eV; 356 G4double secondaryKinetic=-1000*eV;
353 357
354 if (!fasterCode) << 358 if (fasterCode == false)
355 { 359 {
356 secondaryKinetic = RandomizeEjectedElect 360 secondaryKinetic = RandomizeEjectedElectronEnergy(particle->GetDefinition(),k,ionizationShell);
357 } 361 }
358 else 362 else
359 { 363 {
360 secondaryKinetic = RandomizeEjectedElect 364 secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(particle->GetDefinition(),k,ionizationShell);
361 } 365 }
362 366
363 G4int Z = 8; 367 G4int Z = 8;
364 368
365 G4ThreeVector deltaDirection = 369 G4ThreeVector deltaDirection =
366 GetAngularDistribution()->SampleDirectionF 370 GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
367 Z, ionizationShell, 371 Z, ionizationShell,
368 couple->GetMaterial()); 372 couple->GetMaterial());
369 373
370 if (secondaryKinetic>0) 374 if (secondaryKinetic>0)
371 { 375 {
372 auto dp = new G4DynamicParticle (G4Elec << 376 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic);
373 fvect->push_back(dp); 377 fvect->push_back(dp);
374 } 378 }
375 379
376 if (particle->GetDefinition() == G4Electro 380 if (particle->GetDefinition() == G4Electron::ElectronDefinition())
377 { 381 {
378 G4double deltaTotalMomentum = std::sqrt( 382 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
379 383
380 G4double finalPx = totalMomentum*primary 384 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
381 G4double finalPy = totalMomentum*primary 385 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
382 G4double finalPz = totalMomentum*primary 386 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
383 G4double finalMomentum = std::sqrt(final 387 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
384 finalPx /= finalMomentum; 388 finalPx /= finalMomentum;
385 finalPy /= finalMomentum; 389 finalPy /= finalMomentum;
386 finalPz /= finalMomentum; 390 finalPz /= finalMomentum;
387 391
388 G4ThreeVector direction; 392 G4ThreeVector direction;
389 direction.set(finalPx,finalPy,finalPz); 393 direction.set(finalPx,finalPy,finalPz);
390 394
391 fParticleChangeForGamma->ProposeMomentum 395 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit());
392 } 396 }
393 397
394 else fParticleChangeForGamma->ProposeMomen 398 else fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
395 399
396 // AM: sample deexcitation 400 // AM: sample deexcitation
397 // here we assume that H_{2}O electronic l 401 // here we assume that H_{2}O electronic levels are the same as Oxygen.
398 // this can be considered true with a roug 402 // this can be considered true with a rough 10% error in energy on K-shell,
399 403
400 std::size_t secNumberInit = 0; << 404 G4int secNumberInit = 0;// need to know at a certain point the energy of secondaries
401 std::size_t secNumberFinal = 0; << 405 G4int secNumberFinal = 0;// So I'll make the diference and then sum the energies
402 406
403 G4double bindingEnergy = 0; 407 G4double bindingEnergy = 0;
404 bindingEnergy = waterStructure.IonisationE 408 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
405 409
406 // SI: additional protection if tcs interp 410 // SI: additional protection if tcs interpolation method is modified
407 if (k<bindingEnergy) return; 411 if (k<bindingEnergy) return;
408 // 412 //
409 413
410 G4double scatteredEnergy = k-bindingEnergy << 414 if(fAtomDeexcitation)
411 <<
412 // SI: only atomic deexcitation from K she <<
413 if((fAtomDeexcitation != nullptr) && ioniz <<
414 { 415 {
415 const G4AtomicShell* shell = << 416 G4AtomicShellEnumerator as = fKShell;
416 fAtomDeexcitation->GetAtomicShell(Z, G << 417
>> 418 if (ionizationShell <5 && ionizationShell >1)
>> 419 {
>> 420 as = G4AtomicShellEnumerator(4-ionizationShell);
>> 421 }
>> 422 else if (ionizationShell <2)
>> 423 {
>> 424 as = G4AtomicShellEnumerator(3);
>> 425 }
>> 426
>> 427 // FOR DEBUG ONLY
>> 428 // if (ionizationShell == 4) {
>> 429 //
>> 430 // G4cout << "Z: " << Z << " as: " << as
>> 431 // << " ionizationShell: " << ionizationShell << " bindingEnergy: "<< bindingEnergy/eV << G4endl;
>> 432 // G4cout << "Press <Enter> key to continue..." << G4endl;
>> 433 // G4cin.ignore();
>> 434 // }
>> 435
>> 436 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
417 secNumberInit = fvect->size(); 437 secNumberInit = fvect->size();
418 fAtomDeexcitation->GenerateParticles(fve 438 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
419 secNumberFinal = fvect->size(); 439 secNumberFinal = fvect->size();
>> 440 }
420 441
421 if(secNumberFinal > secNumberInit) << 442 // Note that secondaryKinetic is the energy of the delta ray, not of all secondaries.
422 { <<
423 for (std::size_t i=secNumberInit; i<secNumbe <<
424 { <<
425 //Check if there is enough residual <<
426 if (bindingEnergy >= ((*fvect)[i])-> <<
427 { <<
428 //Ok, this is a valid secondary: <<
429 bindingEnergy -= ((*fvect)[i])->GetKine <<
430 } <<
431 else <<
432 { <<
433 //Invalid secondary: not enough energy <<
434 //Keep its energy in the local deposit <<
435 delete (*fvect)[i]; <<
436 (*fvect)[i]=nullptr; <<
437 } <<
438 } <<
439 } <<
440 443
>> 444 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
>> 445 G4double deexSecEnergy = 0;
>> 446 for (G4int j=secNumberInit; j < secNumberFinal; j++)
>> 447 {
>> 448 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
441 } 449 }
442 450
443 //This should never happen <<
444 if(bindingEnergy < 0.0) <<
445 G4Exception("G4DNAEmfietzoglouIonisatioMo <<
446 "em2050",FatalException,"Nega <<
447 <<
448 //bindingEnergy has been decreased <<
449 //by the amount of energy taken away by de <<
450 if (!statCode) 451 if (!statCode)
451 { 452 {
452 fParticleChangeForGamma->SetProposedKine 453 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
453 fParticleChangeForGamma->ProposeLocalEne << 454 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy-secondaryKinetic-deexSecEnergy);
454 } 455 }
455 else 456 else
456 { 457 {
457 fParticleChangeForGamma->SetProposedKine 458 fParticleChangeForGamma->SetProposedKineticEnergy(k);
458 fParticleChangeForGamma->ProposeLocalEne 459 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy);
459 } 460 }
460 << 461
461 // TEST ////////////////////////// 462 // TEST //////////////////////////
462 // if (secondaryKinetic<0) abort(); 463 // if (secondaryKinetic<0) abort();
463 // if (scatteredEnergy<0) abort(); 464 // if (scatteredEnergy<0) abort();
464 // if (k-scatteredEnergy-secondaryKinetic- 465 // if (k-scatteredEnergy-secondaryKinetic-deexSecEnergy<0) abort();
465 // if (k-scatteredEnergy<0) abort(); 466 // if (k-scatteredEnergy<0) abort();
466 ///////////////////////////////// 467 /////////////////////////////////
467 468
468 const G4Track * theIncomingTrack = fPartic 469 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
469 G4DNAChemistryManager::Instance()->CreateW 470 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
470 ionizationShell, 471 ionizationShell,
471 theIncomingTrack); 472 theIncomingTrack);
472 } 473 }
473 } 474 }
474 475
475 //....oooOO0OOooo........oooOO0OOooo........oo 476 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
476 477
477 G4double G4DNABornIonisationModel2::RandomizeE 478 G4double G4DNABornIonisationModel2::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
478 479 G4double k,
479 480 G4int shell)
480 { 481 {
481 // G4cout << "*** SLOW computation for " << 482 // G4cout << "*** SLOW computation for " << " " << particleDefinition->GetParticleName() << G4endl;
482 483
483 if (particleDefinition == G4Electron::Electr 484 if (particleDefinition == G4Electron::ElectronDefinition())
484 { 485 {
485 G4double maximumEnergyTransfer = 0.; 486 G4double maximumEnergyTransfer = 0.;
486 if ((k + waterStructure.IonisationEnergy(s 487 if ((k + waterStructure.IonisationEnergy(shell)) / 2. > k)
487 maximumEnergyTransfer = k; 488 maximumEnergyTransfer = k;
488 else 489 else
489 maximumEnergyTransfer = (k + waterStruct 490 maximumEnergyTransfer = (k + waterStructure.IonisationEnergy(shell)) / 2.;
490 491
491 // SI : original method 492 // SI : original method
492 /* 493 /*
493 G4double crossSectionMaximum = 0.; 494 G4double crossSectionMaximum = 0.;
494 for(G4double value=waterStructure.Ionisat 495 for(G4double value=waterStructure.IonisationEnergy(shell); value<=maximumEnergyTransfer; value+=0.1*eV)
495 { 496 {
496 G4double differentialCrossSection = Diffe 497 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
497 if(differentialCrossSection >= crossSecti 498 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
498 } 499 }
499 */ 500 */
500 501
501 // SI : alternative method 502 // SI : alternative method
502 G4double crossSectionMaximum = 0.; 503 G4double crossSectionMaximum = 0.;
503 504
504 G4double minEnergy = waterStructure.Ionisa 505 G4double minEnergy = waterStructure.IonisationEnergy(shell);
505 G4double maxEnergy = maximumEnergyTransfer 506 G4double maxEnergy = maximumEnergyTransfer;
506 G4int nEnergySteps = 50; 507 G4int nEnergySteps = 50;
507 508
508 G4double value(minEnergy); 509 G4double value(minEnergy);
509 G4double stpEnergy(std::pow(maxEnergy / va 510 G4double stpEnergy(std::pow(maxEnergy / value,
510 1. / static_ca 511 1. / static_cast<G4double>(nEnergySteps - 1)));
511 G4int step(nEnergySteps); 512 G4int step(nEnergySteps);
512 while (step > 0) 513 while (step > 0)
513 { 514 {
514 step--; 515 step--;
515 G4double differentialCrossSection = 516 G4double differentialCrossSection =
516 DifferentialCrossSection(particleDef 517 DifferentialCrossSection(particleDefinition,
517 k / eV, 518 k / eV,
518 value / eV, 519 value / eV,
519 shell); 520 shell);
520 if (differentialCrossSection >= crossSec 521 if (differentialCrossSection >= crossSectionMaximum)
521 crossSectionMaximum = differentialCros 522 crossSectionMaximum = differentialCrossSection;
522 value *= stpEnergy; 523 value *= stpEnergy;
523 } 524 }
524 // 525 //
525 526
526 G4double secondaryElectronKineticEnergy = 527 G4double secondaryElectronKineticEnergy = 0.;
527 do 528 do
528 { 529 {
529 secondaryElectronKineticEnergy = G4Unifo 530 secondaryElectronKineticEnergy = G4UniformRand()* (maximumEnergyTransfer-waterStructure.IonisationEnergy(shell));
530 } while(G4UniformRand()*crossSectionMaximu << 531 }while(G4UniformRand()*crossSectionMaximum >
531 DifferentialCrossSection(particleDefin 532 DifferentialCrossSection(particleDefinition, k/eV,
532 (secondaryElectronKineticEnergy+wa 533 (secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
533 534
534 return secondaryElectronKineticEnergy; 535 return secondaryElectronKineticEnergy;
535 536
536 } 537 }
537 538
538 if (particleDefinition == G4Proton::ProtonDe << 539 else if (particleDefinition == G4Proton::ProtonDefinition())
539 { 540 {
540 G4double maximumKineticEnergyTransfer = 4. 541 G4double maximumKineticEnergyTransfer = 4.
541 * (electron_mass_c2 / proton_mass_c2) 542 * (electron_mass_c2 / proton_mass_c2) * k;
542 543
543 G4double crossSectionMaximum = 0.; 544 G4double crossSectionMaximum = 0.;
544 for (G4double value = waterStructure.Ionis 545 for (G4double value = waterStructure.IonisationEnergy(shell);
545 value <= 4. * waterStructure.Ionisatio 546 value <= 4. * waterStructure.IonisationEnergy(shell); value += 0.1 * eV)
546 { 547 {
547 G4double differentialCrossSection = 548 G4double differentialCrossSection =
548 DifferentialCrossSection(particleDef 549 DifferentialCrossSection(particleDefinition,
549 k / eV, 550 k / eV,
550 value / eV, 551 value / eV,
551 shell); 552 shell);
552 if (differentialCrossSection >= crossSec 553 if (differentialCrossSection >= crossSectionMaximum)
553 crossSectionMaximum = differentialCros 554 crossSectionMaximum = differentialCrossSection;
554 } 555 }
555 556
556 G4double secondaryElectronKineticEnergy = 557 G4double secondaryElectronKineticEnergy = 0.;
557 do 558 do
558 { 559 {
559 secondaryElectronKineticEnergy = G4Unifo 560 secondaryElectronKineticEnergy = G4UniformRand()* maximumKineticEnergyTransfer;
560 } while(G4UniformRand()*crossSectionMaximu << 561 }while(G4UniformRand()*crossSectionMaximum >=
561 DifferentialCrossSection(particleDefin 562 DifferentialCrossSection(particleDefinition, k/eV,
562 (secondaryElectronKineticEnergy+wa 563 (secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
563 564
564 return secondaryElectronKineticEnergy; 565 return secondaryElectronKineticEnergy;
565 } 566 }
566 567
567 return 0; 568 return 0;
568 } 569 }
569 570
570 //....oooOO0OOooo........oooOO0OOooo........oo 571 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
571 572
572 // The following section is not used anymore b 573 // The following section is not used anymore but is kept for memory
573 // GetAngularDistribution()->SampleDirectionFo 574 // GetAngularDistribution()->SampleDirectionForShell is used instead
574 575
575 /* 576 /*
576 void G4DNABornIonisationModel2::RandomizeEjec 577 void G4DNABornIonisationModel2::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
577 G4double k, 578 G4double k,
578 G4double secKinetic, 579 G4double secKinetic,
579 G4double & cosTheta, 580 G4double & cosTheta,
580 G4double & phi ) 581 G4double & phi )
581 { 582 {
582 if (particleDefinition == G4Electron::Electro 583 if (particleDefinition == G4Electron::ElectronDefinition())
583 { 584 {
584 phi = twopi * G4UniformRand(); 585 phi = twopi * G4UniformRand();
585 if (secKinetic < 50.*eV) cosTheta = (2.*G4Uni 586 if (secKinetic < 50.*eV) cosTheta = (2.*G4UniformRand())-1.;
586 else if (secKinetic <= 200.*eV) 587 else if (secKinetic <= 200.*eV)
587 { 588 {
588 if (G4UniformRand() <= 0.1) cosTheta = (2.*G4 589 if (G4UniformRand() <= 0.1) cosTheta = (2.*G4UniformRand())-1.;
589 else cosTheta = G4UniformRand()*(std::sqrt(2. 590 else cosTheta = G4UniformRand()*(std::sqrt(2.)/2);
590 } 591 }
591 else 592 else
592 { 593 {
593 G4double sin2O = (1.-secKinetic/k) / (1.+secK 594 G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
594 cosTheta = std::sqrt(1.-sin2O); 595 cosTheta = std::sqrt(1.-sin2O);
595 } 596 }
596 } 597 }
597 598
598 else if (particleDefinition == G4Proton::Prot 599 else if (particleDefinition == G4Proton::ProtonDefinition())
599 { 600 {
600 G4double maxSecKinetic = 4.* (electron_mass_c 601 G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
601 phi = twopi * G4UniformRand(); 602 phi = twopi * G4UniformRand();
602 603
603 // cosTheta = std::sqrt(secKinetic / maxSecKi 604 // cosTheta = std::sqrt(secKinetic / maxSecKinetic);
604 605
605 // Restriction below 100 eV from Emfietzoglou 606 // Restriction below 100 eV from Emfietzoglou (2000)
606 607
607 if (secKinetic>100*eV) cosTheta = std::sqrt(s 608 if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
608 else cosTheta = (2.*G4UniformRand())-1.; 609 else cosTheta = (2.*G4UniformRand())-1.;
609 610
610 } 611 }
611 } 612 }
612 */ 613 */
613 614
614 //....oooOO0OOooo........oooOO0OOooo........oo 615 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
615 G4double G4DNABornIonisationModel2::Differenti 616 G4double G4DNABornIonisationModel2::DifferentialCrossSection(G4ParticleDefinition * /*particleDefinition*/,
616 617 G4double k,
617 618 G4double energyTransfer,
618 619 G4int ionizationLevelIndex)
619 { 620 {
620 G4double sigma = 0.; 621 G4double sigma = 0.;
621 622
622 if (energyTransfer >= waterStructure.Ionisat 623 if (energyTransfer >= waterStructure.IonisationEnergy(ionizationLevelIndex)/eV)
623 { 624 {
624 G4double valueT1 = 0; 625 G4double valueT1 = 0;
625 G4double valueT2 = 0; 626 G4double valueT2 = 0;
626 G4double valueE21 = 0; 627 G4double valueE21 = 0;
627 G4double valueE22 = 0; 628 G4double valueE22 = 0;
628 G4double valueE12 = 0; 629 G4double valueE12 = 0;
629 G4double valueE11 = 0; 630 G4double valueE11 = 0;
630 631
631 G4double xs11 = 0; 632 G4double xs11 = 0;
632 G4double xs12 = 0; 633 G4double xs12 = 0;
633 G4double xs21 = 0; 634 G4double xs21 = 0;
634 G4double xs22 = 0; 635 G4double xs22 = 0;
635 636
636 // Protection against out of boundary acce <<
637 if (k==fTdummyVec.back()) k=k*(1.-1e-12); <<
638 // <<
639 <<
640 // k should be in eV and energy transfer e 637 // k should be in eV and energy transfer eV also
641 638
642 auto t2 = std::upper_bound(fTdummyVec.begi << 639 std::vector<G4double>::iterator t2 = std::upper_bound(fTdummyVec.begin(),
643 640 fTdummyVec.end(),
644 641 k);
645 642
646 auto t1 = t2 - 1; << 643 std::vector<G4double>::iterator t1 = t2 - 1;
647 644
648 // SI : the following condition avoids sit 645 // SI : the following condition avoids situations where energyTransfer >last vector element
649 646
650 if (energyTransfer <= fVecm[(*t1)].back() 647 if (energyTransfer <= fVecm[(*t1)].back()
651 && energyTransfer <= fVecm[(*t2)].back 648 && energyTransfer <= fVecm[(*t2)].back())
652 { 649 {
653 auto e12 = std::upper_bound(fVecm[(*t1)] << 650 std::vector<G4double>::iterator e12 = std::upper_bound(fVecm[(*t1)].begin(),
654 651 fVecm[(*t1)].end(),
655 652 energyTransfer);
656 auto e11 = e12 - 1; << 653 std::vector<G4double>::iterator e11 = e12 - 1;
657 654
658 auto e22 = std::upper_bound(fVecm[(*t2)] << 655 std::vector<G4double>::iterator e22 = std::upper_bound(fVecm[(*t2)].begin(),
659 656 fVecm[(*t2)].end(),
660 657 energyTransfer);
661 auto e21 = e22 - 1; << 658 std::vector<G4double>::iterator e21 = e22 - 1;
662 659
663 valueT1 = *t1; 660 valueT1 = *t1;
664 valueT2 = *t2; 661 valueT2 = *t2;
665 valueE21 = *e21; 662 valueE21 = *e21;
666 valueE22 = *e22; 663 valueE22 = *e22;
667 valueE12 = *e12; 664 valueE12 = *e12;
668 valueE11 = *e11; 665 valueE11 = *e11;
669 666
670 xs11 = fDiffCrossSectionData[ionizationL 667 xs11 = fDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
671 xs12 = fDiffCrossSectionData[ionizationL 668 xs12 = fDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
672 xs21 = fDiffCrossSectionData[ionizationL 669 xs21 = fDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
673 xs22 = fDiffCrossSectionData[ionizationL 670 xs22 = fDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
674 671
675 } 672 }
676 673
677 G4double xsProduct = xs11 * xs12 * xs21 * 674 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
678 if (xsProduct != 0.) 675 if (xsProduct != 0.)
679 { 676 {
680 sigma = QuadInterpolator(valueE11, 677 sigma = QuadInterpolator(valueE11,
681 valueE12, 678 valueE12,
682 valueE21, 679 valueE21,
683 valueE22, 680 valueE22,
684 xs11, 681 xs11,
685 xs12, 682 xs12,
686 xs21, 683 xs21,
687 xs22, 684 xs22,
688 valueT1, 685 valueT1,
689 valueT2, 686 valueT2,
690 k, 687 k,
691 energyTransfer) 688 energyTransfer);
692 } 689 }
693 } 690 }
694 691
695 return sigma; 692 return sigma;
696 } 693 }
697 694
698 //....oooOO0OOooo........oooOO0OOooo........oo 695 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
699 696
700 G4double G4DNABornIonisationModel2::Interpolat 697 G4double G4DNABornIonisationModel2::Interpolate(G4double e1,
701 698 G4double e2,
702 699 G4double e,
703 700 G4double xs1,
704 701 G4double xs2)
705 { 702 {
706 G4double value = 0.; 703 G4double value = 0.;
707 704
708 // Log-log interpolation by default 705 // Log-log interpolation by default
709 706
710 if (e1 != 0 && e2 != 0 && (std::log10(e2) - 707 if (e1 != 0 && e2 != 0 && (std::log10(e2) - std::log10(e1)) != 0
711 && !fasterCode) 708 && !fasterCode)
712 { 709 {
713 G4double a = (std::log10(xs2) - std::log10 710 G4double a = (std::log10(xs2) - std::log10(xs1))
714 / (std::log10(e2) - std::log10(e1)); 711 / (std::log10(e2) - std::log10(e1));
715 G4double b = std::log10(xs2) - a * std::lo 712 G4double b = std::log10(xs2) - a * std::log10(e2);
716 G4double sigma = a * std::log10(e) + b; 713 G4double sigma = a * std::log10(e) + b;
717 value = (std::pow(10., sigma)); 714 value = (std::pow(10., sigma));
718 } 715 }
719 716
720 // Switch to lin-lin interpolation 717 // Switch to lin-lin interpolation
721 /* 718 /*
722 if ((e2-e1)!=0) 719 if ((e2-e1)!=0)
723 { 720 {
724 G4double d1 = xs1; 721 G4double d1 = xs1;
725 G4double d2 = xs2; 722 G4double d2 = xs2;
726 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1) 723 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
727 } 724 }
728 */ 725 */
729 726
730 // Switch to log-lin interpolation for faste 727 // Switch to log-lin interpolation for faster code
731 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 & 728 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 && fasterCode)
732 { 729 {
733 G4double d1 = std::log10(xs1); 730 G4double d1 = std::log10(xs1);
734 G4double d2 = std::log10(xs2); 731 G4double d2 = std::log10(xs2);
735 value = std::pow(10., (d1 + (d2 - d1) * (e 732 value = std::pow(10., (d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
736 } 733 }
737 734
738 // Switch to lin-lin interpolation for faste 735 // Switch to lin-lin interpolation for faster code
739 // in case one of xs1 or xs2 (=cum proba) va 736 // in case one of xs1 or xs2 (=cum proba) value is zero
740 737
741 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) 738 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) && fasterCode)
742 { 739 {
743 G4double d1 = xs1; 740 G4double d1 = xs1;
744 G4double d2 = xs2; 741 G4double d2 = xs2;
745 value = (d1 + (d2 - d1) * (e - e1) / (e2 - 742 value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
746 } 743 }
747 744
748 /* 745 /*
749 G4cout 746 G4cout
750 << e1 << " " 747 << e1 << " "
751 << e2 << " " 748 << e2 << " "
752 << e << " " 749 << e << " "
753 << xs1 << " " 750 << xs1 << " "
754 << xs2 << " " 751 << xs2 << " "
755 << value 752 << value
756 << G4endl; 753 << G4endl;
757 */ 754 */
758 755
759 return value; 756 return value;
760 } 757 }
761 758
762 //....oooOO0OOooo........oooOO0OOooo........oo 759 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
763 760
764 G4double G4DNABornIonisationModel2::QuadInterp 761 G4double G4DNABornIonisationModel2::QuadInterpolator(G4double e11,
765 762 G4double e12,
766 763 G4double e21,
767 764 G4double e22,
768 765 G4double xs11,
769 766 G4double xs12,
770 767 G4double xs21,
771 768 G4double xs22,
772 769 G4double t1,
773 770 G4double t2,
774 771 G4double t,
775 772 G4double e)
776 { 773 {
777 G4double interpolatedvalue1 = Interpolate(e1 774 G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
778 G4double interpolatedvalue2 = Interpolate(e2 775 G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
779 G4double value = Interpolate(t1, 776 G4double value = Interpolate(t1,
780 t2, 777 t2,
781 t, 778 t,
782 interpolatedval 779 interpolatedvalue1,
783 interpolatedval 780 interpolatedvalue2);
784 781
785 return value; 782 return value;
786 } 783 }
787 784
788 //....oooOO0OOooo........oooOO0OOooo........oo 785 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
789 786
790 G4double G4DNABornIonisationModel2::GetPartial 787 G4double G4DNABornIonisationModel2::GetPartialCrossSection(const G4Material*,
791 788 G4int level,
792 789 const G4ParticleDefinition* /*particle*/,
793 790 G4double kineticEnergy)
794 { 791 {
795 return fTableData->GetComponent(level)->Find 792 return fTableData->GetComponent(level)->FindValue(kineticEnergy);
796 } 793 }
797 794
798 //....oooOO0OOooo........oooOO0OOooo........oo 795 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
799 796
800 G4int G4DNABornIonisationModel2::RandomSelect( 797 G4int G4DNABornIonisationModel2::RandomSelect(G4double k)
801 { 798 {
802 G4int level = 0; 799 G4int level = 0;
803 800
804 auto valuesBuffer = new G4double[fTableData << 801 G4double* valuesBuffer = new G4double[fTableData->NumberOfComponents()];
805 const auto n = (G4int)fTableData->NumberOfC << 802 const size_t n(fTableData->NumberOfComponents());
806 G4int i(n); << 803 size_t i(n);
807 G4double value = 0.; 804 G4double value = 0.;
808 805
809 while (i > 0) 806 while (i > 0)
810 { 807 {
811 i--; 808 i--;
812 valuesBuffer[i] = fTableData->GetComponent 809 valuesBuffer[i] = fTableData->GetComponent(i)->FindValue(k);
813 value += valuesBuffer[i]; 810 value += valuesBuffer[i];
814 } 811 }
815 812
816 value *= G4UniformRand(); 813 value *= G4UniformRand();
817 814
818 i = n; 815 i = n;
819 816
820 while (i > 0) 817 while (i > 0)
821 { 818 {
822 i--; 819 i--;
823 820
824 if (valuesBuffer[i] > value) 821 if (valuesBuffer[i] > value)
825 { 822 {
826 delete[] valuesBuffer; 823 delete[] valuesBuffer;
827 return i; 824 return i;
828 } 825 }
829 value -= valuesBuffer[i]; 826 value -= valuesBuffer[i];
830 } 827 }
831 828
832 << 829 if (valuesBuffer)
833 delete[] valuesBuffer; 830 delete[] valuesBuffer;
834 831
835 return level; 832 return level;
836 } 833 }
837 834
838 //....oooOO0OOooo........oooOO0OOooo........oo 835 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
839 836
840 G4double G4DNABornIonisationModel2::RandomizeE 837 G4double G4DNABornIonisationModel2::RandomizeEjectedElectronEnergyFromCumulatedDcs(G4ParticleDefinition* particleDefinition,
841 838 G4double k,
842 839 G4int shell)
843 { 840 {
844 // G4cout << "*** FAST computation for " << 841 // G4cout << "*** FAST computation for " << " " << particleDefinition->GetParticleName() << G4endl;
845 842
846 G4double secondaryElectronKineticEnergy = 0. 843 G4double secondaryElectronKineticEnergy = 0.;
847 844
848 G4double random = G4UniformRand(); 845 G4double random = G4UniformRand();
849 846
850 secondaryElectronKineticEnergy = TransferedE 847 secondaryElectronKineticEnergy = TransferedEnergy(particleDefinition,
851 848 k / eV,
852 849 shell,
853 850 random) * eV
854 - waterStructure.IonisationEnergy(shell) 851 - waterStructure.IonisationEnergy(shell);
855 852
856 // G4cout << TransferedEnergy(particleDefini 853 // G4cout << TransferedEnergy(particleDefinition, k/eV, shell, random) << G4endl;
857 if (secondaryElectronKineticEnergy < 0.) 854 if (secondaryElectronKineticEnergy < 0.)
858 return 0.; 855 return 0.;
859 // 856 //
860 857
861 return secondaryElectronKineticEnergy; 858 return secondaryElectronKineticEnergy;
862 } 859 }
863 860
864 //....oooOO0OOooo........oooOO0OOooo........oo 861 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
865 862
866 G4double G4DNABornIonisationModel2::Transfered 863 G4double G4DNABornIonisationModel2::TransferedEnergy(G4ParticleDefinition* /*particleDefinition*/,
867 864 G4double k,
868 865 G4int ionizationLevelIndex,
869 866 G4double random)
870 { 867 {
871 868
872 G4double nrj = 0.; 869 G4double nrj = 0.;
873 870
874 G4double valueK1 = 0; 871 G4double valueK1 = 0;
875 G4double valueK2 = 0; 872 G4double valueK2 = 0;
876 G4double valuePROB21 = 0; 873 G4double valuePROB21 = 0;
877 G4double valuePROB22 = 0; 874 G4double valuePROB22 = 0;
878 G4double valuePROB12 = 0; 875 G4double valuePROB12 = 0;
879 G4double valuePROB11 = 0; 876 G4double valuePROB11 = 0;
880 877
881 G4double nrjTransf11 = 0; 878 G4double nrjTransf11 = 0;
882 G4double nrjTransf12 = 0; 879 G4double nrjTransf12 = 0;
883 G4double nrjTransf21 = 0; 880 G4double nrjTransf21 = 0;
884 G4double nrjTransf22 = 0; 881 G4double nrjTransf22 = 0;
885 882
886 // Protection against out of boundary access <<
887 if (k==fTdummyVec.back()) k=k*(1.-1e-12); <<
888 // <<
889 <<
890 // k should be in eV 883 // k should be in eV
891 auto k2 = std::upper_bound(fTdummyVec.begin( << 884 std::vector<G4double>::iterator k2 = std::upper_bound(fTdummyVec.begin(),
892 885 fTdummyVec.end(),
893 886 k);
894 auto k1 = k2 - 1; << 887 std::vector<G4double>::iterator k1 = k2 - 1;
895 888
896 /* 889 /*
897 G4cout << "----> k=" << k 890 G4cout << "----> k=" << k
898 << " " << *k1 891 << " " << *k1
899 << " " << *k2 892 << " " << *k2
900 << " " << random 893 << " " << random
901 << " " << ionizationLevelIndex 894 << " " << ionizationLevelIndex
902 << " " << eProbaShellMap[ionizationLevelInd 895 << " " << eProbaShellMap[ionizationLevelIndex][(*k1)].back()
903 << " " << eProbaShellMap[ionizationLevelInd 896 << " " << eProbaShellMap[ionizationLevelIndex][(*k2)].back()
904 << G4endl; 897 << G4endl;
905 */ 898 */
906 899
907 // SI : the following condition avoids situa 900 // SI : the following condition avoids situations where random >last vector element
908 if (random <= fProbaShellMap[ionizationLevel 901 if (random <= fProbaShellMap[ionizationLevelIndex][(*k1)].back()
909 && random <= fProbaShellMap[ionizationLe 902 && random <= fProbaShellMap[ionizationLevelIndex][(*k2)].back())
910 { 903 {
911 auto prob12 = << 904 std::vector<G4double>::iterator prob12 =
912 std::upper_bound(fProbaShellMap[ioniza 905 std::upper_bound(fProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
913 fProbaShellMap[ioniza 906 fProbaShellMap[ionizationLevelIndex][(*k1)].end(),
914 random); 907 random);
915 908
916 auto prob11 = prob12 - 1; << 909 std::vector<G4double>::iterator prob11 = prob12 - 1;
917 910
918 auto prob22 = << 911 std::vector<G4double>::iterator prob22 =
919 std::upper_bound(fProbaShellMap[ioniza 912 std::upper_bound(fProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
920 fProbaShellMap[ioniza 913 fProbaShellMap[ionizationLevelIndex][(*k2)].end(),
921 random); 914 random);
922 915
923 auto prob21 = prob22 - 1; << 916 std::vector<G4double>::iterator prob21 = prob22 - 1;
924 917
925 valueK1 = *k1; 918 valueK1 = *k1;
926 valueK2 = *k2; 919 valueK2 = *k2;
927 valuePROB21 = *prob21; 920 valuePROB21 = *prob21;
928 valuePROB22 = *prob22; 921 valuePROB22 = *prob22;
929 valuePROB12 = *prob12; 922 valuePROB12 = *prob12;
930 valuePROB11 = *prob11; 923 valuePROB11 = *prob11;
931 924
932 /* 925 /*
933 G4cout << " " << random << " " << 926 G4cout << " " << random << " " << valuePROB11 << " "
934 << valuePROB12 << " " << valuePROB21 << " 927 << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
935 */ 928 */
936 929
937 nrjTransf11 = fNrjTransfData[ionizationLev 930 nrjTransf11 = fNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
938 nrjTransf12 = fNrjTransfData[ionizationLev 931 nrjTransf12 = fNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
939 nrjTransf21 = fNrjTransfData[ionizationLev 932 nrjTransf21 = fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
940 nrjTransf22 = fNrjTransfData[ionizationLev 933 nrjTransf22 = fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
941 934
942 /* 935 /*
943 G4cout << " " << ionizationLevelIn 936 G4cout << " " << ionizationLevelIndex << " "
944 << random << " " <<valueK1 << " " << valu 937 << random << " " <<valueK1 << " " << valueK2 << G4endl;
945 938
946 G4cout << " " << random << " " << 939 G4cout << " " << random << " " << nrjTransf11 << " "
947 << nrjTransf12 << " " << nrjTransf21 << " 940 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
948 */ 941 */
949 942
950 } 943 }
951 // Avoids cases where cum xs is zero for k1 944 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
952 if (random > fProbaShellMap[ionizationLevelI 945 if (random > fProbaShellMap[ionizationLevelIndex][(*k1)].back())
953 { 946 {
954 auto prob22 = << 947 std::vector<G4double>::iterator prob22 =
955 std::upper_bound(fProbaShellMap[ioniza 948 std::upper_bound(fProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
956 fProbaShellMap[ioniza 949 fProbaShellMap[ionizationLevelIndex][(*k2)].end(),
957 random); 950 random);
958 951
959 auto prob21 = prob22 - 1; << 952 std::vector<G4double>::iterator prob21 = prob22 - 1;
960 953
961 valueK1 = *k1; 954 valueK1 = *k1;
962 valueK2 = *k2; 955 valueK2 = *k2;
963 valuePROB21 = *prob21; 956 valuePROB21 = *prob21;
964 valuePROB22 = *prob22; 957 valuePROB22 = *prob22;
965 958
966 // G4cout << " " << random << " " < 959 // G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
967 960
968 nrjTransf21 = fNrjTransfData[ionizationLev 961 nrjTransf21 = fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
969 nrjTransf22 = fNrjTransfData[ionizationLev 962 nrjTransf22 = fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
970 963
971 G4double interpolatedvalue2 = Interpolate( 964 G4double interpolatedvalue2 = Interpolate(valuePROB21,
972 965 valuePROB22,
973 966 random,
974 967 nrjTransf21,
975 968 nrjTransf22);
976 969
977 // zeros are explicitly set 970 // zeros are explicitly set
978 971
979 G4double value = Interpolate(valueK1, valu 972 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
980 973
981 /* 974 /*
982 G4cout << " " << ionizationLevelIn 975 G4cout << " " << ionizationLevelIndex << " "
983 << random << " " <<valueK1 << " " << valu 976 << random << " " <<valueK1 << " " << valueK2 << G4endl;
984 977
985 G4cout << " " << random << " " << 978 G4cout << " " << random << " " << nrjTransf11 << " "
986 << nrjTransf12 << " " << nrjTransf21 << " 979 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
987 980
988 G4cout << "ici" << " " << value << G4endl 981 G4cout << "ici" << " " << value << G4endl;
989 */ 982 */
990 983
991 return value; 984 return value;
992 } 985 }
993 986
994 // End electron and proton cases 987 // End electron and proton cases
995 988
996 G4double nrjTransfProduct = nrjTransf11 * nr 989 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21
997 * nrjTransf22; 990 * nrjTransf22;
998 //G4cout << "nrjTransfProduct=" << nrjTransf 991 //G4cout << "nrjTransfProduct=" << nrjTransfProduct << G4endl;
999 992
1000 if (nrjTransfProduct != 0.) 993 if (nrjTransfProduct != 0.)
1001 { 994 {
1002 nrj = QuadInterpolator(valuePROB11, 995 nrj = QuadInterpolator(valuePROB11,
1003 valuePROB12, 996 valuePROB12,
1004 valuePROB21, 997 valuePROB21,
1005 valuePROB22, 998 valuePROB22,
1006 nrjTransf11, 999 nrjTransf11,
1007 nrjTransf12, 1000 nrjTransf12,
1008 nrjTransf21, 1001 nrjTransf21,
1009 nrjTransf22, 1002 nrjTransf22,
1010 valueK1, 1003 valueK1,
1011 valueK2, 1004 valueK2,
1012 k, 1005 k,
1013 random); 1006 random);
1014 } 1007 }
1015 // G4cout << nrj << endl; 1008 // G4cout << nrj << endl;
1016 1009
1017 return nrj; 1010 return nrj;
1018 } 1011 }
1019 1012