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