| 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: G4DNABornIonisationModel1.cc 87631 2014-12-14 12:42:05Z matkara $
26 // 27 //
27 28
28 #include "G4DNABornIonisationModel1.hh" 29 #include "G4DNABornIonisationModel1.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 G4DNABornIonisationModel1::G4DNABornIonisation 45 G4DNABornIonisationModel1::G4DNABornIonisationModel1(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 SetDeexcitationFlag(true); 63 SetDeexcitationFlag(true);
64 fAtomDeexcitation = nullptr; << 64 fAtomDeexcitation = 0;
65 fParticleChangeForGamma = nullptr; << 65 fParticleChangeForGamma = 0;
66 fpMolWaterDensity = nullptr; << 66 fpMolWaterDensity = 0;
67 67
68 // Define default angular generator << 68 // define default angular generator
69 SetAngularDistribution(new G4DNABornAngle()) 69 SetAngularDistribution(new G4DNABornAngle());
70 70
71 // Selection of computation method 71 // Selection of computation method
72 72
73 fasterCode = false; 73 fasterCode = false;
74 74
75 // Selection of stationary mode 75 // Selection of stationary mode
76 76
77 statCode = false; 77 statCode = false;
78 78
79 // Selection of SP scaling 79 // Selection of SP scaling
80 80
81 spScaling = true; 81 spScaling = true;
82 } 82 }
83 83
84 //....oooOO0OOooo........oooOO0OOooo........oo 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85 85
86 G4DNABornIonisationModel1::~G4DNABornIonisatio 86 G4DNABornIonisationModel1::~G4DNABornIonisationModel1()
87 { 87 {
88 // Cross section 88 // Cross section
89 89
90 std::map<G4String, G4DNACrossSectionDataSet* 90 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
91 for (pos = tableData.begin(); pos != tableDa 91 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
92 { 92 {
93 G4DNACrossSectionDataSet* table = pos->sec 93 G4DNACrossSectionDataSet* table = pos->second;
94 delete table; 94 delete table;
95 } 95 }
96 96
97 // Final state 97 // Final state
98 98
99 eVecm.clear(); 99 eVecm.clear();
100 pVecm.clear(); 100 pVecm.clear();
101 } 101 }
102 102
103 //....oooOO0OOooo........oooOO0OOooo........oo 103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
104 104
105 void G4DNABornIonisationModel1::Initialise(con 105 void G4DNABornIonisationModel1::Initialise(const G4ParticleDefinition* particle,
106 con 106 const G4DataVector& /*cuts*/)
107 { 107 {
108 108
109 if (verboseLevel > 3) 109 if (verboseLevel > 3)
110 { 110 {
111 G4cout << "Calling G4DNABornIonisationMode 111 G4cout << "Calling G4DNABornIonisationModel1::Initialise()" << G4endl;
112 } 112 }
113 113
114 // Energy limits 114 // Energy limits
115 115
116 G4String fileElectron("dna/sigma_ionisation_ 116 G4String fileElectron("dna/sigma_ionisation_e_born");
117 G4String fileProton("dna/sigma_ionisation_p_ 117 G4String fileProton("dna/sigma_ionisation_p_born");
118 118
119 G4ParticleDefinition* electronDef = G4Electr 119 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
120 G4ParticleDefinition* protonDef = G4Proton:: 120 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
121 121
122 G4String electron; 122 G4String electron;
123 G4String proton; 123 G4String proton;
124 124
125 G4double scaleFactor = (1.e-22 / 3.343) * m* 125 G4double scaleFactor = (1.e-22 / 3.343) * m*m;
126 126
127 const char *path = G4FindDataDir("G4LEDATA") << 127 char *path = getenv("G4LEDATA");
128 128
129 // *** ELECTRON 129 // *** ELECTRON
130 130
131 electron = electronDef->GetParticleName(); 131 electron = electronDef->GetParticleName();
132 132
133 tableFile[electron] = fileElectron; 133 tableFile[electron] = fileElectron;
134 134
135 lowEnergyLimit[electron] = 11. * eV; 135 lowEnergyLimit[electron] = 11. * eV;
136 highEnergyLimit[electron] = 1. * MeV; 136 highEnergyLimit[electron] = 1. * MeV;
137 137
138 // Cross section 138 // Cross section
139 139
140 auto tableE = new G4DNACrossSectionDataSet( << 140 G4DNACrossSectionDataSet* tableE = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
141 tableE->LoadData(fileElectron); 141 tableE->LoadData(fileElectron);
142 142
143 tableData[electron] = tableE; 143 tableData[electron] = tableE;
144 144
145 // Final state 145 // Final state
146 146
147 std::ostringstream eFullFileName; 147 std::ostringstream eFullFileName;
148 148
149 if (fasterCode) eFullFileName << path << "/d 149 if (fasterCode) eFullFileName << path << "/dna/sigmadiff_cumulated_ionisation_e_born_hp.dat";
150 if (!fasterCode) eFullFileName << path << "/ 150 if (!fasterCode) eFullFileName << path << "/dna/sigmadiff_ionisation_e_born.dat";
151 151
152 std::ifstream eDiffCrossSection(eFullFileNam 152 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
153 153
154 if (!eDiffCrossSection) 154 if (!eDiffCrossSection)
155 { 155 {
156 if (fasterCode) G4Exception("G4DNABornIoni 156 if (fasterCode) G4Exception("G4DNABornIonisationModel1::Initialise","em0003",
157 FatalException,"Missing data file:/dna 157 FatalException,"Missing data file:/dna/sigmadiff_cumulated_ionisation_e_born_hp.dat");
158 158
159 if (!fasterCode) G4Exception("G4DNABornIon 159 if (!fasterCode) G4Exception("G4DNABornIonisationModel1::Initialise","em0003",
160 FatalException,"Missing data file:/dna 160 FatalException,"Missing data file:/dna/sigmadiff_ionisation_e_born.dat");
161 } 161 }
162 162
>> 163 //
>> 164
163 // Clear the arrays for re-initialization ca 165 // Clear the arrays for re-initialization case (MT mode)
164 // March 25th, 2014 - Vaclav Stepan, Sebasti 166 // March 25th, 2014 - Vaclav Stepan, Sebastien Incerti
165 167
166 eTdummyVec.clear(); 168 eTdummyVec.clear();
167 pTdummyVec.clear(); 169 pTdummyVec.clear();
168 170
169 eVecm.clear(); 171 eVecm.clear();
170 pVecm.clear(); 172 pVecm.clear();
171 173
172 for (G4int j=0; j<5; j++) << 174 for (int j=0; j<5; j++)
173 { 175 {
174 eProbaShellMap[j].clear(); 176 eProbaShellMap[j].clear();
175 pProbaShellMap[j].clear(); 177 pProbaShellMap[j].clear();
176 178
177 eDiffCrossSectionData[j].clear(); 179 eDiffCrossSectionData[j].clear();
178 pDiffCrossSectionData[j].clear(); 180 pDiffCrossSectionData[j].clear();
179 181
180 eNrjTransfData[j].clear(); 182 eNrjTransfData[j].clear();
181 pNrjTransfData[j].clear(); 183 pNrjTransfData[j].clear();
182 } 184 }
183 <<
184 // 185 //
185 186
186 eTdummyVec.push_back(0.); 187 eTdummyVec.push_back(0.);
187 while(!eDiffCrossSection.eof()) 188 while(!eDiffCrossSection.eof())
188 { 189 {
189 G4double tDummy; << 190 double tDummy;
190 G4double eDummy; << 191 double eDummy;
191 eDiffCrossSection>>tDummy>>eDummy; 192 eDiffCrossSection>>tDummy>>eDummy;
192 if (tDummy != eTdummyVec.back()) eTdummyVe 193 if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
193 194
194 G4double tmp; << 195 double tmp;
195 for (G4int j=0; j<5; j++) << 196 for (int j=0; j<5; j++)
196 { 197 {
197 eDiffCrossSection>> tmp; 198 eDiffCrossSection>> tmp;
198 199
199 eDiffCrossSectionData[j][tDummy][eDummy] 200 eDiffCrossSectionData[j][tDummy][eDummy] = tmp;
200 201
201 if (fasterCode) 202 if (fasterCode)
202 { 203 {
203 eNrjTransfData[j][tDummy][eDiffCrossSe 204 eNrjTransfData[j][tDummy][eDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
204 eProbaShellMap[j][tDummy].push_back(eD 205 eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]);
205 } 206 }
206 207
207 // SI - only if eof is not reached 208 // SI - only if eof is not reached
208 if (!eDiffCrossSection.eof() && !fasterC 209 if (!eDiffCrossSection.eof() && !fasterCode) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
209 210
210 if (!fasterCode) eVecm[tDummy].push_back 211 if (!fasterCode) eVecm[tDummy].push_back(eDummy);
211 212
212 } 213 }
213 } 214 }
214 215
215 // *** PROTON 216 // *** PROTON
216 217
217 proton = protonDef->GetParticleName(); 218 proton = protonDef->GetParticleName();
218 219
219 tableFile[proton] = fileProton; 220 tableFile[proton] = fileProton;
220 221
221 lowEnergyLimit[proton] = 500. * keV; 222 lowEnergyLimit[proton] = 500. * keV;
222 highEnergyLimit[proton] = 100. * MeV; 223 highEnergyLimit[proton] = 100. * MeV;
223 224
224 // Cross section 225 // Cross section
225 226
226 auto tableP = new G4DNACrossSectionDataSet( << 227 G4DNACrossSectionDataSet* tableP = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
227 tableP->LoadData(fileProton); 228 tableP->LoadData(fileProton);
228 229
229 tableData[proton] = tableP; 230 tableData[proton] = tableP;
230 231
231 // Final state 232 // Final state
232 233
233 std::ostringstream pFullFileName; 234 std::ostringstream pFullFileName;
234 235
235 if (fasterCode) pFullFileName << path << "/d 236 if (fasterCode) pFullFileName << path << "/dna/sigmadiff_cumulated_ionisation_p_born_hp.dat";
236 237
237 if (!fasterCode) pFullFileName << path << "/ 238 if (!fasterCode) pFullFileName << path << "/dna/sigmadiff_ionisation_p_born.dat";
238 239
239 std::ifstream pDiffCrossSection(pFullFileNam 240 std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
240 241
241 if (!pDiffCrossSection) 242 if (!pDiffCrossSection)
242 { 243 {
243 if (fasterCode) G4Exception("G4DNABornIoni 244 if (fasterCode) G4Exception("G4DNABornIonisationModel1::Initialise","em0003",
244 FatalException,"Missing data file:/dna 245 FatalException,"Missing data file:/dna/sigmadiff_cumulated_ionisation_p_born_hp.dat");
245 246
246 if (!fasterCode) G4Exception("G4DNABornIon 247 if (!fasterCode) G4Exception("G4DNABornIonisationModel1::Initialise","em0003",
247 FatalException,"Missing data file:/dna 248 FatalException,"Missing data file:/dna/sigmadiff_ionisation_p_born.dat");
248 } 249 }
249 250
250 pTdummyVec.push_back(0.); 251 pTdummyVec.push_back(0.);
251 while(!pDiffCrossSection.eof()) 252 while(!pDiffCrossSection.eof())
252 { 253 {
253 G4double tDummy; << 254 double tDummy;
254 G4double eDummy; << 255 double eDummy;
255 pDiffCrossSection>>tDummy>>eDummy; 256 pDiffCrossSection>>tDummy>>eDummy;
256 if (tDummy != pTdummyVec.back()) pTdummyVe 257 if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
257 for (G4int j=0; j<5; j++) << 258 for (int j=0; j<5; j++)
258 { 259 {
259 pDiffCrossSection>>pDiffCrossSectionData 260 pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
260 261
261 if (fasterCode) 262 if (fasterCode)
262 { 263 {
263 pNrjTransfData[j][tDummy][pDiffCrossSe 264 pNrjTransfData[j][tDummy][pDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
264 pProbaShellMap[j][tDummy].push_back(pD 265 pProbaShellMap[j][tDummy].push_back(pDiffCrossSectionData[j][tDummy][eDummy]);
265 } 266 }
266 267
267 // SI - only if eof is not reached ! 268 // SI - only if eof is not reached !
268 if (!pDiffCrossSection.eof() && !fasterC 269 if (!pDiffCrossSection.eof() && !fasterCode) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
269 270
270 if (!fasterCode) pVecm[tDummy].push_back 271 if (!fasterCode) pVecm[tDummy].push_back(eDummy);
271 } 272 }
272 } 273 }
273 274
274 // 275 //
275 276
276 if (particle==electronDef) 277 if (particle==electronDef)
277 { 278 {
278 SetLowEnergyLimit(lowEnergyLimit[electron] 279 SetLowEnergyLimit(lowEnergyLimit[electron]);
279 SetHighEnergyLimit(highEnergyLimit[electro 280 SetHighEnergyLimit(highEnergyLimit[electron]);
280 } 281 }
281 282
282 if (particle==protonDef) 283 if (particle==protonDef)
283 { 284 {
284 SetLowEnergyLimit(lowEnergyLimit[proton]); 285 SetLowEnergyLimit(lowEnergyLimit[proton]);
285 SetHighEnergyLimit(highEnergyLimit[proton] 286 SetHighEnergyLimit(highEnergyLimit[proton]);
286 } 287 }
287 288
288 if( verboseLevel>0 ) 289 if( verboseLevel>0 )
289 { 290 {
290 G4cout << "Born ionisation model is initia 291 G4cout << "Born ionisation model is initialized " << G4endl
291 << "Energy range: " 292 << "Energy range: "
292 << LowEnergyLimit() / eV << " eV - " 293 << LowEnergyLimit() / eV << " eV - "
293 << HighEnergyLimit() / keV << " keV for " 294 << HighEnergyLimit() / keV << " keV for "
294 << particle->GetParticleName() 295 << particle->GetParticleName()
295 << G4endl; 296 << G4endl;
296 } 297 }
297 298
298 // Initialize water density pointer 299 // Initialize water density pointer
299 <<
300 fpMolWaterDensity = G4DNAMolecularMaterial:: 300 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->
301 GetNumMolPerVolTableFor(G4Material::GetMater 301 GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
302 302
303 // AD <<
304 <<
305 fAtomDeexcitation = G4LossTableManager::Inst <<
306 <<
307 // 303 //
>> 304 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
308 305
309 if (isInitialised) 306 if (isInitialised)
310 { return;} 307 { return;}
311 fParticleChangeForGamma = GetParticleChangeF 308 fParticleChangeForGamma = GetParticleChangeForGamma();
312 isInitialised = true; 309 isInitialised = true;
313 } 310 }
314 311
315 //....oooOO0OOooo........oooOO0OOooo........oo 312 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
316 313
317 G4double G4DNABornIonisationModel1::CrossSecti 314 G4double G4DNABornIonisationModel1::CrossSectionPerVolume(const G4Material* material,
318 315 const G4ParticleDefinition* particleDefinition,
319 316 G4double ekin,
320 317 G4double,
321 318 G4double)
322 { 319 {
323 if (verboseLevel > 3) 320 if (verboseLevel > 3)
324 { 321 {
325 G4cout << "Calling CrossSectionPerVolume() 322 G4cout << "Calling CrossSectionPerVolume() of G4DNABornIonisationModel1"
326 << G4endl; 323 << G4endl;
>> 324
327 } 325 }
328 326
329 if ( 327 if (
330 particleDefinition != G4Proton::ProtonDe 328 particleDefinition != G4Proton::ProtonDefinition()
331 && 329 &&
332 particleDefinition != G4Electron::Electr 330 particleDefinition != G4Electron::ElectronDefinition()
333 ) 331 )
334 332
335 return 0; 333 return 0;
336 334
337 // Calculate total cross section for model 335 // Calculate total cross section for model
338 336
339 G4double lowLim = 0; 337 G4double lowLim = 0;
340 G4double highLim = 0; 338 G4double highLim = 0;
341 G4double sigma=0; 339 G4double sigma=0;
342 340
343 G4double waterDensity = (*fpMolWaterDensity) 341 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
344 342
345 const G4String& particleName = particleDefin << 343 if(waterDensity!= 0.0)
346 << 344 // if (material == nistwater || material->GetBaseMaterial() == nistwater)
347 std::map< G4String,G4double,std::less<G4Stri <<
348 pos1 = lowEnergyLimit.find(particleName); <<
349 if (pos1 != lowEnergyLimit.end()) <<
350 { 345 {
351 lowLim = pos1->second; << 346 const G4String& particleName = particleDefinition->GetParticleName();
352 } <<
353 347
354 std::map< G4String,G4double,std::less<G4Stri << 348 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
355 pos2 = highEnergyLimit.find(particleName); << 349 pos1 = lowEnergyLimit.find(particleName);
356 if (pos2 != highEnergyLimit.end()) << 350 if (pos1 != lowEnergyLimit.end())
357 { << 351 {
358 highLim = pos2->second; << 352 lowLim = pos1->second;
359 } << 353 }
360 354
361 if (ekin >= lowLim && ekin <= highLim) << 355 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
362 { << 356 pos2 = highEnergyLimit.find(particleName);
363 std::map< G4String,G4DNACrossSectionDataSe << 357 if (pos2 != highEnergyLimit.end())
364 pos = tableData.find(particleName); << 358 {
>> 359 highLim = pos2->second;
>> 360 }
365 361
366 if (pos != tableData.end()) << 362 if (ekin >= lowLim && ekin < highLim)
367 { 363 {
368 G4DNACrossSectionDataSet* table = pos->s << 364 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
369 if (table != nullptr) << 365 pos = tableData.find(particleName);
>> 366
>> 367 if (pos != tableData.end())
370 { 368 {
371 sigma = table->FindValue(ekin); << 369 G4DNACrossSectionDataSet* table = pos->second;
>> 370 if (table != 0)
>> 371 {
>> 372 sigma = table->FindValue(ekin);
372 373
373 // ICRU49 electronic SP scaling - ZF, << 374 // ICRU49 electronic SP scaling - ZF, SI
374 375
375 if (particleDefinition == G4Proton::Pr << 376 if (particleDefinition == G4Proton::ProtonDefinition() && ekin < 70*MeV && spScaling)
376 { << 377 {
377 G4double A = 1.39241700556072800000E- << 378 G4double A = 1.39241700556072800000E-009 ;
378 G4double B = -8.52610412942622630000E << 379 G4double B = -8.52610412942622630000E-002 ;
379 sigma = sigma * G4Exp(A*(ekin/eV)+B); << 380 sigma = sigma * std::exp(A*(ekin/eV)+B);
380 } << 381 }
381 // << 382 //
382 383
>> 384 }
>> 385 }
>> 386 else
>> 387 {
>> 388 G4Exception("G4DNABornIonisationModel1::CrossSectionPerVolume","em0002",
>> 389 FatalException,"Model not applicable to particle type.");
383 } 390 }
384 } 391 }
385 else << 392
>> 393 if (verboseLevel > 2)
386 { 394 {
387 G4Exception("G4DNABornIonisationModel1:: << 395 G4cout << "__________________________________" << G4endl;
388 FatalException,"Model not applicable << 396 G4cout << "G4DNABornIonisationModel1 - XS INFO START" << G4endl;
>> 397 G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
>> 398 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
>> 399 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
>> 400 G4cout << "G4DNABornIonisationModel1 - XS INFO END" << G4endl;
389 } 401 }
390 } << 402 } // if (waterMaterial)
391 <<
392 if (verboseLevel > 2) <<
393 { <<
394 G4cout << "_______________________________ <<
395 G4cout << "G4DNABornIonisationModel1 - XS <<
396 G4cout << "Kinetic energy(eV)=" << ekin/eV <<
397 G4cout << "Cross section per water molecul <<
398 G4cout << "Cross section per water molecul <<
399 G4cout << "G4DNABornIonisationModel1 - XS <<
400 } <<
401 403
402 return sigma*waterDensity; 404 return sigma*waterDensity;
403 } 405 }
404 406
405 //....oooOO0OOooo........oooOO0OOooo........oo 407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
406 408
407 void G4DNABornIonisationModel1::SampleSecondar 409 void G4DNABornIonisationModel1::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
408 410 const G4MaterialCutsCouple* couple,
409 411 const G4DynamicParticle* particle,
410 412 G4double,
411 413 G4double)
412 { 414 {
413 415
414 if (verboseLevel > 3) 416 if (verboseLevel > 3)
415 { 417 {
416 G4cout << "Calling SampleSecondaries() of 418 G4cout << "Calling SampleSecondaries() of G4DNABornIonisationModel1"
417 << G4endl; 419 << G4endl;
418 } 420 }
419 421
420 G4double lowLim = 0; 422 G4double lowLim = 0;
421 G4double highLim = 0; 423 G4double highLim = 0;
422 424
423 G4double k = particle->GetKineticEnergy(); 425 G4double k = particle->GetKineticEnergy();
424 426
425 const G4String& particleName = particle->Get 427 const G4String& particleName = particle->GetDefinition()->GetParticleName();
426 428
427 std::map< G4String,G4double,std::less<G4Stri 429 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
428 pos1 = lowEnergyLimit.find(particleName); 430 pos1 = lowEnergyLimit.find(particleName);
429 431
430 if (pos1 != lowEnergyLimit.end()) 432 if (pos1 != lowEnergyLimit.end())
431 { 433 {
432 lowLim = pos1->second; 434 lowLim = pos1->second;
433 } 435 }
434 436
435 std::map< G4String,G4double,std::less<G4Stri 437 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
436 pos2 = highEnergyLimit.find(particleName); 438 pos2 = highEnergyLimit.find(particleName);
437 439
438 if (pos2 != highEnergyLimit.end()) 440 if (pos2 != highEnergyLimit.end())
439 { 441 {
440 highLim = pos2->second; 442 highLim = pos2->second;
441 } 443 }
442 444
443 if (k >= lowLim && k <= highLim) << 445 if (k >= lowLim && k < highLim)
444 { 446 {
445 G4ParticleMomentum primaryDirection = part 447 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
446 G4double particleMass = particle->GetDefin 448 G4double particleMass = particle->GetDefinition()->GetPDGMass();
447 G4double totalEnergy = k + particleMass; 449 G4double totalEnergy = k + particleMass;
448 G4double pSquare = k * (totalEnergy + part 450 G4double pSquare = k * (totalEnergy + particleMass);
449 G4double totalMomentum = std::sqrt(pSquare 451 G4double totalMomentum = std::sqrt(pSquare);
450 452
451 G4int ionizationShell = 0; 453 G4int ionizationShell = 0;
452 454
453 if (!fasterCode) ionizationShell = RandomS 455 if (!fasterCode) ionizationShell = RandomSelect(k,particleName);
454 456
455 // SI: The following protection is necessa 457 // SI: The following protection is necessary to avoid infinite loops :
456 // sigmadiff_ionisation_e_born.dat has no 458 // sigmadiff_ionisation_e_born.dat has non zero partial xs at 18 eV for shell 3 (ionizationShell ==2)
457 // sigmadiff_cumulated_ionisation_e_born. 459 // sigmadiff_cumulated_ionisation_e_born.dat has zero cumulated partial xs at 18 eV for shell 3 (ionizationShell ==2)
458 // this is due to the fact that the max a 460 // this is due to the fact that the max allowed transfered energy is (18+10.79)/2=17.025 eV and only transfered energies
459 // strictly above this value have non zer 461 // strictly above this value have non zero partial xs in sigmadiff_ionisation_e_born.dat (starting at trans = 17.12 eV)
460 462
461 if (fasterCode) 463 if (fasterCode)
462 do 464 do
463 { 465 {
464 ionizationShell = RandomSelect(k,particl 466 ionizationShell = RandomSelect(k,particleName);
465 } while (k<19*eV && ionizationShell==2 && << 467 }while (k<19*eV && ionizationShell==2 && particle->GetDefinition()==G4Electron::ElectronDefinition());
>> 468
>> 469 // AM: sample deexcitation
>> 470 // here we assume that H_{2}O electronic levels are the same as Oxygen.
>> 471 // this can be considered true with a rough 10% error in energy on K-shell,
>> 472
>> 473 G4int secNumberInit = 0;// need to know at a certain point the energy of secondaries
>> 474 G4int secNumberFinal = 0;// So I'll make the diference and then sum the energies
466 475
467 G4double bindingEnergy = 0; 476 G4double bindingEnergy = 0;
468 bindingEnergy = waterStructure.IonisationE 477 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
469 478
470 // SI: additional protection if tcs interp << 479 //SI: additional protection if tcs interpolation method is modified
471 if (k<bindingEnergy) return; 480 if (k<bindingEnergy) return;
472 // 481 //
473 482
>> 483 G4int Z = 8;
>> 484 if(fAtomDeexcitation)
>> 485 {
>> 486 G4AtomicShellEnumerator as = fKShell;
>> 487
>> 488 if (ionizationShell <5 && ionizationShell >1)
>> 489 {
>> 490 as = G4AtomicShellEnumerator(4-ionizationShell);
>> 491 }
>> 492 else if (ionizationShell <2)
>> 493 {
>> 494 as = G4AtomicShellEnumerator(3);
>> 495 }
>> 496
>> 497 // FOR DEBUG ONLY
>> 498 // if (ionizationShell == 4) {
>> 499 //
>> 500 // G4cout << "Z: " << Z << " as: " << as
>> 501 // << " ionizationShell: " << ionizationShell << " bindingEnergy: "<< bindingEnergy/eV << G4endl;
>> 502 // G4cout << "Press <Enter> key to continue..." << G4endl;
>> 503 // G4cin.ignore();
>> 504 // }
>> 505
>> 506 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
>> 507 secNumberInit = fvect->size();
>> 508 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
>> 509 secNumberFinal = fvect->size();
>> 510 }
>> 511
474 G4double secondaryKinetic=-1000*eV; 512 G4double secondaryKinetic=-1000*eV;
475 513
476 if (!fasterCode) << 514 if (fasterCode == false)
477 { 515 {
478 secondaryKinetic = RandomizeEjectedElect 516 secondaryKinetic = RandomizeEjectedElectronEnergy(particle->GetDefinition(),k,ionizationShell);
479 } 517 }
>> 518 // SI - 01/04/2014
480 else 519 else
481 { 520 {
482 secondaryKinetic = RandomizeEjectedElect 521 secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(particle->GetDefinition(),k,ionizationShell);
483 } 522 }
484 // 523 //
485 524
486 G4int Z = 8; <<
487 <<
488 G4ThreeVector deltaDirection = 525 G4ThreeVector deltaDirection =
489 GetAngularDistribution()->SampleDirectionF 526 GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
490 Z, ionizationShell, 527 Z, ionizationShell,
491 couple->GetMaterial()); 528 couple->GetMaterial());
492 529
493 if (secondaryKinetic>0) <<
494 { <<
495 auto dp = new G4DynamicParticle (G4Elec <<
496 fvect->push_back(dp); <<
497 } <<
498 <<
499 if (particle->GetDefinition() == G4Electro 530 if (particle->GetDefinition() == G4Electron::ElectronDefinition())
500 { 531 {
501 G4double deltaTotalMomentum = std::sqrt( 532 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
502 533
503 G4double finalPx = totalMomentum*primary 534 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
504 G4double finalPy = totalMomentum*primary 535 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
505 G4double finalPz = totalMomentum*primary 536 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
506 G4double finalMomentum = std::sqrt(final 537 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
507 finalPx /= finalMomentum; 538 finalPx /= finalMomentum;
508 finalPy /= finalMomentum; 539 finalPy /= finalMomentum;
509 finalPz /= finalMomentum; 540 finalPz /= finalMomentum;
510 541
511 G4ThreeVector direction; 542 G4ThreeVector direction;
512 direction.set(finalPx,finalPy,finalPz); 543 direction.set(finalPx,finalPy,finalPz);
513 544
514 fParticleChangeForGamma->ProposeMomentum 545 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit());
515 } 546 }
516 547
517 else fParticleChangeForGamma->ProposeMomen 548 else fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
518 549
519 // AM: sample deexcitation << 550 // note that secondaryKinetic is the energy of the delta ray, not of all secondaries.
520 // here we assume that H_{2}O electronic l <<
521 // this can be considered true with a roug <<
522 <<
523 std::size_t secNumberInit = 0;// need to k <<
524 std::size_t secNumberFinal = 0;// So I'll <<
525 <<
526 G4double scatteredEnergy = k-bindingEnergy 551 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
527 << 552 G4double deexSecEnergy = 0;
528 // SI: only atomic deexcitation from K she << 553 for (G4int j=secNumberInit; j < secNumberFinal; j++)
529 if((fAtomDeexcitation != nullptr) && ioniz <<
530 { 554 {
531 const G4AtomicShell* shell = << 555 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
532 fAtomDeexcitation->GetAtomicShell(Z, G <<
533 secNumberInit = fvect->size(); <<
534 fAtomDeexcitation->GenerateParticles(fve <<
535 secNumberFinal = fvect->size(); <<
536 <<
537 //TEST <<
538 //G4cout << "ionizationShell=" << ioniza <<
539 //G4cout << "bindingEnergy=" << bindingE <<
540 <<
541 if(secNumberFinal > secNumberInit) <<
542 { <<
543 for (std::size_t i=secNumberInit; i<secNumbe <<
544 { <<
545 //Check if there is enough residual <<
546 if (bindingEnergy >= ((*fvect)[i])-> <<
547 { <<
548 //Ok, this is a valid secondary: <<
549 bindingEnergy -= ((*fvect)[i])->GetKine <<
550 //G4cout << "--deex nrj=" << ((*fvect)[ <<
551 //<< G4endl; <<
552 } <<
553 else <<
554 { <<
555 //Invalid secondary: not enough energy <<
556 //Keep its energy in the local deposit <<
557 delete (*fvect)[i]; <<
558 (*fvect)[i]=nullptr; <<
559 } <<
560 } <<
561 } <<
562 <<
563 //TEST <<
564 //G4cout << "k=" << k/eV<< G4endl; <<
565 //G4cout << "secondaryKinetic=" << seconda <<
566 //G4cout << "scatteredEnergy=" << scattere <<
567 //G4cout << "deposited energy=" << binding <<
568 // <<
569 <<
570 } 556 }
571 557
572 //This should never happen <<
573 if(bindingEnergy < 0.0) <<
574 G4Exception("G4DNABornIonisatioModel1::Sa <<
575 "em2050",FatalException,"Nega <<
576 <<
577 //bindingEnergy has been decreased <<
578 //by the amount of energy taken away by de <<
579 if (!statCode) 558 if (!statCode)
580 { 559 {
581 fParticleChangeForGamma->SetProposedKine 560 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
582 fParticleChangeForGamma->ProposeLocalEne << 561 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy-secondaryKinetic-deexSecEnergy);
583 } 562 }
584 else 563 else
585 { 564 {
586 fParticleChangeForGamma->SetProposedKine 565 fParticleChangeForGamma->SetProposedKineticEnergy(k);
587 fParticleChangeForGamma->ProposeLocalEne 566 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy);
588 } 567 }
589 <<
590 // TEST ////////////////////////// <<
591 // if (secondaryKinetic<0) abort(); <<
592 // if (scatteredEnergy<0) abort(); <<
593 // if (k-scatteredEnergy-secondaryKinetic- <<
594 // if (k-scatteredEnergy<0) abort(); <<
595 ///////////////////////////////// <<
596 568
>> 569 // SI - 01/04/2014
>> 570 if (secondaryKinetic>0)
>> 571 {
>> 572 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic);
>> 573 fvect->push_back(dp);
>> 574 }
>> 575 //
>> 576
597 const G4Track * theIncomingTrack = fPartic 577 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
598 G4DNAChemistryManager::Instance()->CreateW 578 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
599 ionizationShell, 579 ionizationShell,
600 theIncomingTrack); 580 theIncomingTrack);
601 } 581 }
602 } 582 }
603 583
604 //....oooOO0OOooo........oooOO0OOooo........oo 584 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
605 585
606 G4double G4DNABornIonisationModel1::RandomizeE 586 G4double G4DNABornIonisationModel1::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
607 587 G4double k,
608 588 G4int shell)
609 { 589 {
610 // G4cout << "*** SLOW computation for " << 590 // G4cout << "*** SLOW computation for " << " " << particleDefinition->GetParticleName() << G4endl;
611 591
612 if (particleDefinition == G4Electron::Electr 592 if (particleDefinition == G4Electron::ElectronDefinition())
613 { 593 {
614 G4double maximumEnergyTransfer = 0.; 594 G4double maximumEnergyTransfer = 0.;
615 if ((k + waterStructure.IonisationEnergy(s 595 if ((k + waterStructure.IonisationEnergy(shell)) / 2. > k)
616 maximumEnergyTransfer = k; 596 maximumEnergyTransfer = k;
617 else 597 else
618 maximumEnergyTransfer = (k + waterStruct 598 maximumEnergyTransfer = (k + waterStructure.IonisationEnergy(shell)) / 2.;
619 599
620 // SI : original method 600 // SI : original method
621 /* 601 /*
622 G4double crossSectionMaximum = 0.; << 602 G4double crossSectionMaximum = 0.;
623 for(G4double value=waterStructure.Ionisati << 603 for(G4double value=waterStructure.IonisationEnergy(shell); value<=maximumEnergyTransfer; value+=0.1*eV)
624 { << 604 {
625 G4double differentialCrossSection = Diffe 605 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
626 if(differentialCrossSection >= crossSecti 606 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
627 } << 607 }
628 */ << 608 */
629 609
630 // SI : alternative method 610 // SI : alternative method
631 G4double crossSectionMaximum = 0.; 611 G4double crossSectionMaximum = 0.;
632 612
633 G4double minEnergy = waterStructure.Ionisa 613 G4double minEnergy = waterStructure.IonisationEnergy(shell);
634 G4double maxEnergy = maximumEnergyTransfer 614 G4double maxEnergy = maximumEnergyTransfer;
635 G4int nEnergySteps = 50; 615 G4int nEnergySteps = 50;
636 616
637 G4double value(minEnergy); 617 G4double value(minEnergy);
638 G4double stpEnergy(std::pow(maxEnergy / va 618 G4double stpEnergy(std::pow(maxEnergy / value,
639 1. / static_ca 619 1. / static_cast<G4double>(nEnergySteps - 1)));
640 G4int step(nEnergySteps); 620 G4int step(nEnergySteps);
641 while (step > 0) 621 while (step > 0)
642 { 622 {
643 step--; 623 step--;
644 G4double differentialCrossSection = 624 G4double differentialCrossSection =
645 DifferentialCrossSection(particleDef 625 DifferentialCrossSection(particleDefinition,
646 k / eV, 626 k / eV,
647 value / eV, 627 value / eV,
648 shell); 628 shell);
649 if (differentialCrossSection >= crossSec 629 if (differentialCrossSection >= crossSectionMaximum)
650 crossSectionMaximum = differentialCros 630 crossSectionMaximum = differentialCrossSection;
651 value *= stpEnergy; 631 value *= stpEnergy;
652 } 632 }
653 // 633 //
654 634
655 G4double secondaryElectronKineticEnergy = 635 G4double secondaryElectronKineticEnergy = 0.;
656 do 636 do
657 { 637 {
658 secondaryElectronKineticEnergy = G4Unifo 638 secondaryElectronKineticEnergy = G4UniformRand()* (maximumEnergyTransfer-waterStructure.IonisationEnergy(shell));
659 } while(G4UniformRand()*crossSectionMaximu << 639 }while(G4UniformRand()*crossSectionMaximum >
660 DifferentialCrossSection(particleDefin 640 DifferentialCrossSection(particleDefinition, k/eV,
661 (secondaryElectronKineticEnergy+wa 641 (secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
662 642
663 return secondaryElectronKineticEnergy; 643 return secondaryElectronKineticEnergy;
664 644
665 } 645 }
666 646
667 if (particleDefinition == G4Proton::ProtonDe << 647 else if (particleDefinition == G4Proton::ProtonDefinition())
668 { 648 {
669 G4double maximumKineticEnergyTransfer = 4. 649 G4double maximumKineticEnergyTransfer = 4.
670 * (electron_mass_c2 / proton_mass_c2) 650 * (electron_mass_c2 / proton_mass_c2) * k;
671 651
672 G4double crossSectionMaximum = 0.; 652 G4double crossSectionMaximum = 0.;
673 for (G4double value = waterStructure.Ionis 653 for (G4double value = waterStructure.IonisationEnergy(shell);
674 value <= 4. * waterStructure.Ionisatio 654 value <= 4. * waterStructure.IonisationEnergy(shell); value += 0.1 * eV)
675 { 655 {
676 G4double differentialCrossSection = 656 G4double differentialCrossSection =
677 DifferentialCrossSection(particleDef 657 DifferentialCrossSection(particleDefinition,
678 k / eV, 658 k / eV,
679 value / eV, 659 value / eV,
680 shell); 660 shell);
681 if (differentialCrossSection >= crossSec 661 if (differentialCrossSection >= crossSectionMaximum)
682 crossSectionMaximum = differentialCros 662 crossSectionMaximum = differentialCrossSection;
683 } 663 }
684 664
685 G4double secondaryElectronKineticEnergy = 665 G4double secondaryElectronKineticEnergy = 0.;
686 do 666 do
687 { 667 {
688 secondaryElectronKineticEnergy = G4Unifo 668 secondaryElectronKineticEnergy = G4UniformRand()* maximumKineticEnergyTransfer;
689 } while(G4UniformRand()*crossSectionMaximu << 669 }while(G4UniformRand()*crossSectionMaximum >=
690 DifferentialCrossSection(particleDefin 670 DifferentialCrossSection(particleDefinition, k/eV,
691 (secondaryElectronKineticEnergy+wa 671 (secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
692 672
693 return secondaryElectronKineticEnergy; 673 return secondaryElectronKineticEnergy;
694 } 674 }
695 675
696 return 0; 676 return 0;
697 } 677 }
698 678
699 //....oooOO0OOooo........oooOO0OOooo........oo 679 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
700 680
701 // The following section is not used anymore b 681 // The following section is not used anymore but is kept for memory
702 // GetAngularDistribution()->SampleDirectionFo 682 // GetAngularDistribution()->SampleDirectionForShell is used instead
703 683
704 /* 684 /*
705 void G4DNABornIonisationModel1::RandomizeEjec 685 void G4DNABornIonisationModel1::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
706 G4double k, 686 G4double k,
707 G4double secKinetic, 687 G4double secKinetic,
708 G4double & cosTheta, 688 G4double & cosTheta,
709 G4double & phi ) 689 G4double & phi )
710 { 690 {
711 if (particleDefinition == G4Electron::Electro 691 if (particleDefinition == G4Electron::ElectronDefinition())
712 { 692 {
713 phi = twopi * G4UniformRand(); 693 phi = twopi * G4UniformRand();
714 if (secKinetic < 50.*eV) cosTheta = (2.*G4Uni 694 if (secKinetic < 50.*eV) cosTheta = (2.*G4UniformRand())-1.;
715 else if (secKinetic <= 200.*eV) 695 else if (secKinetic <= 200.*eV)
716 { 696 {
717 if (G4UniformRand() <= 0.1) cosTheta = (2.*G4 697 if (G4UniformRand() <= 0.1) cosTheta = (2.*G4UniformRand())-1.;
718 else cosTheta = G4UniformRand()*(std::sqrt(2. 698 else cosTheta = G4UniformRand()*(std::sqrt(2.)/2);
719 } 699 }
720 else 700 else
721 { 701 {
722 G4double sin2O = (1.-secKinetic/k) / (1.+secK 702 G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
723 cosTheta = std::sqrt(1.-sin2O); 703 cosTheta = std::sqrt(1.-sin2O);
724 } 704 }
725 } 705 }
726 706
727 else if (particleDefinition == G4Proton::Prot 707 else if (particleDefinition == G4Proton::ProtonDefinition())
728 { 708 {
729 G4double maxSecKinetic = 4.* (electron_mass_c 709 G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
730 phi = twopi * G4UniformRand(); 710 phi = twopi * G4UniformRand();
731 711
732 // cosTheta = std::sqrt(secKinetic / maxSecKi 712 // cosTheta = std::sqrt(secKinetic / maxSecKinetic);
733 713
734 // Restriction below 100 eV from Emfietzoglou 714 // Restriction below 100 eV from Emfietzoglou (2000)
735 715
736 if (secKinetic>100*eV) cosTheta = std::sqrt(s 716 if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
737 else cosTheta = (2.*G4UniformRand())-1.; 717 else cosTheta = (2.*G4UniformRand())-1.;
738 718
739 } 719 }
740 } 720 }
741 */ 721 */
742 722
743 //....oooOO0OOooo........oooOO0OOooo........oo 723 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
744 G4double G4DNABornIonisationModel1::Differenti << 724 double G4DNABornIonisationModel1::DifferentialCrossSection(G4ParticleDefinition * particleDefinition,
745 725 G4double k,
746 726 G4double energyTransfer,
747 727 G4int ionizationLevelIndex)
748 { 728 {
749 G4double sigma = 0.; 729 G4double sigma = 0.;
750 730
751 if (energyTransfer >= waterStructure.Ionisat << 731 if (energyTransfer >= waterStructure.IonisationEnergy(ionizationLevelIndex))
752 { 732 {
753 G4double valueT1 = 0; 733 G4double valueT1 = 0;
754 G4double valueT2 = 0; 734 G4double valueT2 = 0;
755 G4double valueE21 = 0; 735 G4double valueE21 = 0;
756 G4double valueE22 = 0; 736 G4double valueE22 = 0;
757 G4double valueE12 = 0; 737 G4double valueE12 = 0;
758 G4double valueE11 = 0; 738 G4double valueE11 = 0;
759 739
760 G4double xs11 = 0; 740 G4double xs11 = 0;
761 G4double xs12 = 0; 741 G4double xs12 = 0;
762 G4double xs21 = 0; 742 G4double xs21 = 0;
763 G4double xs22 = 0; 743 G4double xs22 = 0;
764 744
765 if (particleDefinition == G4Electron::Elec 745 if (particleDefinition == G4Electron::ElectronDefinition())
766 { 746 {
767 <<
768 // Protection against out of boundary ac <<
769 if (k==eTdummyVec.back()) k=k*(1.-1e-12) <<
770 // <<
771 <<
772 // k should be in eV and energy transfer 747 // k should be in eV and energy transfer eV also
773 748
774 auto t2 = std::upper_bound(eTdummyVec.be << 749 std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),
775 750 eTdummyVec.end(),
776 751 k);
777 752
778 auto t1 = t2 - 1; << 753 std::vector<double>::iterator t1 = t2 - 1;
779 754
780 // SI : the following condition avoids s 755 // SI : the following condition avoids situations where energyTransfer >last vector element
781 if (energyTransfer <= eVecm[(*t1)].back( 756 if (energyTransfer <= eVecm[(*t1)].back()
782 && energyTransfer <= eVecm[(*t2)].ba 757 && energyTransfer <= eVecm[(*t2)].back())
783 { 758 {
784 auto e12 = << 759 std::vector<double>::iterator e12 =
785 std::upper_bound(eVecm[(*t1)].begi 760 std::upper_bound(eVecm[(*t1)].begin(),
786 eVecm[(*t1)].end( 761 eVecm[(*t1)].end(),
787 energyTransfer); 762 energyTransfer);
788 auto e11 = e12 - 1; << 763 std::vector<double>::iterator e11 = e12 - 1;
789 764
790 auto e22 = << 765 std::vector<double>::iterator e22 =
791 std::upper_bound(eVecm[(*t2)].begi 766 std::upper_bound(eVecm[(*t2)].begin(),
792 eVecm[(*t2)].end( 767 eVecm[(*t2)].end(),
793 energyTransfer); 768 energyTransfer);
794 auto e21 = e22 - 1; << 769 std::vector<double>::iterator e21 = e22 - 1;
795 770
796 valueT1 = *t1; 771 valueT1 = *t1;
797 valueT2 = *t2; 772 valueT2 = *t2;
798 valueE21 = *e21; 773 valueE21 = *e21;
799 valueE22 = *e22; 774 valueE22 = *e22;
800 valueE12 = *e12; 775 valueE12 = *e12;
801 valueE11 = *e11; 776 valueE11 = *e11;
802 777
803 xs11 = eDiffCrossSectionData[ionizatio 778 xs11 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
804 xs12 = eDiffCrossSectionData[ionizatio 779 xs12 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
805 xs21 = eDiffCrossSectionData[ionizatio 780 xs21 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
806 xs22 = eDiffCrossSectionData[ionizatio 781 xs22 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
807 782
808 } 783 }
809 784
810 } 785 }
811 786
812 if (particleDefinition == G4Proton::Proton 787 if (particleDefinition == G4Proton::ProtonDefinition())
813 { 788 {
814 // Protection against out of boundary ac <<
815 if (k==pTdummyVec.back()) k=k*(1.-1e-12) <<
816 // <<
817 <<
818 // k should be in eV and energy transfer 789 // k should be in eV and energy transfer eV also
819 auto t2 = std::upper_bound(pTdummyVec.be << 790 std::vector<double>::iterator t2 = std::upper_bound(pTdummyVec.begin(),
820 791 pTdummyVec.end(),
821 792 k);
822 auto t1 = t2 - 1; << 793 std::vector<double>::iterator t1 = t2 - 1;
823 794
824 auto e12 = std::upper_bound(pVecm[(*t1)] << 795 std::vector<double>::iterator e12 = std::upper_bound(pVecm[(*t1)].begin(),
825 796 pVecm[(*t1)].end(),
826 797 energyTransfer);
827 auto e11 = e12 - 1; << 798 std::vector<double>::iterator e11 = e12 - 1;
828 799
829 auto e22 = std::upper_bound(pVecm[(*t2)] << 800 std::vector<double>::iterator e22 = std::upper_bound(pVecm[(*t2)].begin(),
830 801 pVecm[(*t2)].end(),
831 802 energyTransfer);
832 auto e21 = e22 - 1; << 803 std::vector<double>::iterator e21 = e22 - 1;
833 804
834 valueT1 = *t1; 805 valueT1 = *t1;
835 valueT2 = *t2; 806 valueT2 = *t2;
836 valueE21 = *e21; 807 valueE21 = *e21;
837 valueE22 = *e22; 808 valueE22 = *e22;
838 valueE12 = *e12; 809 valueE12 = *e12;
839 valueE11 = *e11; 810 valueE11 = *e11;
840 811
841 xs11 = pDiffCrossSectionData[ionizationL 812 xs11 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
842 xs12 = pDiffCrossSectionData[ionizationL 813 xs12 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
843 xs21 = pDiffCrossSectionData[ionizationL 814 xs21 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
844 xs22 = pDiffCrossSectionData[ionizationL 815 xs22 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
845 816
846 } 817 }
847 818
848 G4double xsProduct = xs11 * xs12 * xs21 * 819 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
849 if (xsProduct != 0.) 820 if (xsProduct != 0.)
850 { 821 {
851 sigma = QuadInterpolator(valueE11, 822 sigma = QuadInterpolator(valueE11,
852 valueE12, 823 valueE12,
853 valueE21, 824 valueE21,
854 valueE22, 825 valueE22,
855 xs11, 826 xs11,
856 xs12, 827 xs12,
857 xs21, 828 xs21,
858 xs22, 829 xs22,
859 valueT1, 830 valueT1,
860 valueT2, 831 valueT2,
861 k, 832 k,
862 energyTransfer) 833 energyTransfer);
863 } 834 }
864 } 835 }
865 836
866 return sigma; 837 return sigma;
867 } 838 }
868 839
869 //....oooOO0OOooo........oooOO0OOooo........oo 840 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
870 841
871 G4double G4DNABornIonisationModel1::Interpolat 842 G4double G4DNABornIonisationModel1::Interpolate(G4double e1,
872 843 G4double e2,
873 844 G4double e,
874 845 G4double xs1,
875 846 G4double xs2)
876 { 847 {
877 G4double value = 0.; 848 G4double value = 0.;
878 849
879 // Log-log interpolation by default 850 // Log-log interpolation by default
880 851
881 if (e1 != 0 && e2 != 0 && (std::log10(e2) - 852 if (e1 != 0 && e2 != 0 && (std::log10(e2) - std::log10(e1)) != 0
882 && !fasterCode) 853 && !fasterCode)
883 { 854 {
884 G4double a = (std::log10(xs2) - std::log10 855 G4double a = (std::log10(xs2) - std::log10(xs1))
885 / (std::log10(e2) - std::log10(e1)); 856 / (std::log10(e2) - std::log10(e1));
886 G4double b = std::log10(xs2) - a * std::lo 857 G4double b = std::log10(xs2) - a * std::log10(e2);
887 G4double sigma = a * std::log10(e) + b; 858 G4double sigma = a * std::log10(e) + b;
888 value = (std::pow(10., sigma)); 859 value = (std::pow(10., sigma));
889 } 860 }
890 861
891 // Switch to lin-lin interpolation 862 // Switch to lin-lin interpolation
892 /* 863 /*
893 if ((e2-e1)!=0) 864 if ((e2-e1)!=0)
894 { 865 {
895 G4double d1 = xs1; 866 G4double d1 = xs1;
896 G4double d2 = xs2; 867 G4double d2 = xs2;
897 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1) 868 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
898 } 869 }
899 */ << 870 */
900 871
901 // Switch to log-lin interpolation for faste 872 // Switch to log-lin interpolation for faster code
902 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 & 873 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 && fasterCode)
903 { 874 {
904 G4double d1 = std::log10(xs1); 875 G4double d1 = std::log10(xs1);
905 G4double d2 = std::log10(xs2); 876 G4double d2 = std::log10(xs2);
906 value = std::pow(10., (d1 + (d2 - d1) * (e 877 value = std::pow(10., (d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
907 } 878 }
908 879
909 // Switch to lin-lin interpolation for faste 880 // Switch to lin-lin interpolation for faster code
910 // in case one of xs1 or xs2 (=cum proba) va 881 // in case one of xs1 or xs2 (=cum proba) value is zero
911 882
912 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) 883 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) && fasterCode)
913 { 884 {
914 G4double d1 = xs1; 885 G4double d1 = xs1;
915 G4double d2 = xs2; 886 G4double d2 = xs2;
916 value = (d1 + (d2 - d1) * (e - e1) / (e2 - 887 value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
917 } 888 }
918 889
919 /* 890 /*
920 G4cout 891 G4cout
921 << e1 << " " 892 << e1 << " "
922 << e2 << " " 893 << e2 << " "
923 << e << " " 894 << e << " "
924 << xs1 << " " 895 << xs1 << " "
925 << xs2 << " " 896 << xs2 << " "
926 << value 897 << value
927 << G4endl; 898 << G4endl;
928 */ 899 */
929 900
930 return value; 901 return value;
931 } 902 }
932 903
933 //....oooOO0OOooo........oooOO0OOooo........oo 904 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
934 905
935 G4double G4DNABornIonisationModel1::QuadInterp 906 G4double G4DNABornIonisationModel1::QuadInterpolator(G4double e11,
936 907 G4double e12,
937 908 G4double e21,
938 909 G4double e22,
939 910 G4double xs11,
940 911 G4double xs12,
941 912 G4double xs21,
942 913 G4double xs22,
943 914 G4double t1,
944 915 G4double t2,
945 916 G4double t,
946 917 G4double e)
947 { 918 {
948 G4double interpolatedvalue1 = Interpolate(e1 919 G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
949 G4double interpolatedvalue2 = Interpolate(e2 920 G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
950 G4double value = Interpolate(t1, 921 G4double value = Interpolate(t1,
951 t2, 922 t2,
952 t, 923 t,
953 interpolatedval 924 interpolatedvalue1,
954 interpolatedval 925 interpolatedvalue2);
955 926
956 return value; 927 return value;
957 } 928 }
958 929
959 G4double G4DNABornIonisationModel1::GetPartial 930 G4double G4DNABornIonisationModel1::GetPartialCrossSection(const G4Material* /*material*/,
960 931 G4int level,
961 932 const G4ParticleDefinition* particle,
962 933 G4double kineticEnergy)
963 { 934 {
964 std::map<G4String, G4DNACrossSectionDataSet* 935 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
965 pos = tableData.find(particle->GetParticleNa 936 pos = tableData.find(particle->GetParticleName());
966 937
967 if (pos != tableData.end()) 938 if (pos != tableData.end())
968 { 939 {
969 G4DNACrossSectionDataSet* table = pos->sec 940 G4DNACrossSectionDataSet* table = pos->second;
970 return table->GetComponent(level)->FindVal 941 return table->GetComponent(level)->FindValue(kineticEnergy);
971 } 942 }
972 943
973 return 0; 944 return 0;
974 } 945 }
975 946
976 //....oooOO0OOooo........oooOO0OOooo........oo 947 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
977 948
978 G4int G4DNABornIonisationModel1::RandomSelect( 949 G4int G4DNABornIonisationModel1::RandomSelect(G4double k,
979 950 const G4String& particle)
980 { 951 {
981 G4int level = 0; 952 G4int level = 0;
982 953
983 std::map<G4String, G4DNACrossSectionDataSet* 954 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
984 pos = tableData.find(particle); 955 pos = tableData.find(particle);
985 956
986 if (pos != tableData.end()) 957 if (pos != tableData.end())
987 { 958 {
988 G4DNACrossSectionDataSet* table = pos->sec 959 G4DNACrossSectionDataSet* table = pos->second;
989 960
990 if (table != nullptr) << 961 if (table != 0)
991 { 962 {
992 auto valuesBuffer = new G4double[table- << 963 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
993 const auto n = (G4int)table->NumberOfCo << 964 const size_t n(table->NumberOfComponents());
994 G4int i(n); << 965 size_t i(n);
995 G4double value = 0.; 966 G4double value = 0.;
996 967
997 while (i > 0) 968 while (i > 0)
998 { 969 {
999 i--; 970 i--;
1000 valuesBuffer[i] = table->GetComponent 971 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
1001 value += valuesBuffer[i]; 972 value += valuesBuffer[i];
1002 } 973 }
1003 974
1004 value *= G4UniformRand(); 975 value *= G4UniformRand();
1005 976
1006 i = n; 977 i = n;
1007 978
1008 while (i > 0) 979 while (i > 0)
1009 { 980 {
1010 i--; 981 i--;
1011 982
1012 if (valuesBuffer[i] > value) 983 if (valuesBuffer[i] > value)
1013 { 984 {
1014 delete[] valuesBuffer; 985 delete[] valuesBuffer;
1015 return i; 986 return i;
1016 } 987 }
1017 value -= valuesBuffer[i]; 988 value -= valuesBuffer[i];
1018 } 989 }
1019 990
1020 << 991 if (valuesBuffer)
1021 delete[] valuesBuffer; 992 delete[] valuesBuffer;
1022 993
1023 } 994 }
1024 } else 995 } else
1025 { 996 {
1026 G4Exception("G4DNABornIonisationModel1::R 997 G4Exception("G4DNABornIonisationModel1::RandomSelect",
1027 "em0002", 998 "em0002",
1028 FatalException, 999 FatalException,
1029 "Model not applicable to part 1000 "Model not applicable to particle type.");
1030 } 1001 }
1031 1002
1032 return level; 1003 return level;
1033 } 1004 }
1034 1005
1035 //....oooOO0OOooo........oooOO0OOooo........o 1006 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1036 1007
1037 G4double G4DNABornIonisationModel1::Randomize 1008 G4double G4DNABornIonisationModel1::RandomizeEjectedElectronEnergyFromCumulatedDcs(G4ParticleDefinition* particleDefinition,
1038 1009 G4double k,
1039 1010 G4int shell)
1040 { 1011 {
1041 //G4cout << "*** FAST computation for " << 1012 //G4cout << "*** FAST computation for " << " " << particleDefinition->GetParticleName() << G4endl;
1042 1013
1043 G4double secondaryElectronKineticEnergy = 0 1014 G4double secondaryElectronKineticEnergy = 0.;
1044 1015
1045 G4double random = G4UniformRand(); 1016 G4double random = G4UniformRand();
1046 1017
1047 secondaryElectronKineticEnergy = Transfered 1018 secondaryElectronKineticEnergy = TransferedEnergy(particleDefinition,
1048 1019 k / eV,
1049 1020 shell,
1050 1021 random) * eV
1051 - waterStructure.IonisationEnergy(shell 1022 - waterStructure.IonisationEnergy(shell);
1052 1023
1053 //G4cout << RandomTransferedEnergy(particle 1024 //G4cout << RandomTransferedEnergy(particleDefinition, k/eV, shell) << G4endl;
1054 << 1025 // SI - 01/04/2014
1055 if (secondaryElectronKineticEnergy < 0.) 1026 if (secondaryElectronKineticEnergy < 0.)
1056 return 0.; 1027 return 0.;
1057 // 1028 //
1058 1029
1059 return secondaryElectronKineticEnergy; 1030 return secondaryElectronKineticEnergy;
1060 } 1031 }
1061 1032
1062 //....oooOO0OOooo........oooOO0OOooo........o 1033 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1063 1034
1064 G4double G4DNABornIonisationModel1::Transfere 1035 G4double G4DNABornIonisationModel1::TransferedEnergy(G4ParticleDefinition* particleDefinition,
1065 1036 G4double k,
1066 1037 G4int ionizationLevelIndex,
1067 1038 G4double random)
1068 { 1039 {
1069 G4double nrj = 0.; 1040 G4double nrj = 0.;
1070 1041
1071 G4double valueK1 = 0; 1042 G4double valueK1 = 0;
1072 G4double valueK2 = 0; 1043 G4double valueK2 = 0;
1073 G4double valuePROB21 = 0; 1044 G4double valuePROB21 = 0;
1074 G4double valuePROB22 = 0; 1045 G4double valuePROB22 = 0;
1075 G4double valuePROB12 = 0; 1046 G4double valuePROB12 = 0;
1076 G4double valuePROB11 = 0; 1047 G4double valuePROB11 = 0;
1077 1048
1078 G4double nrjTransf11 = 0; 1049 G4double nrjTransf11 = 0;
1079 G4double nrjTransf12 = 0; 1050 G4double nrjTransf12 = 0;
1080 G4double nrjTransf21 = 0; 1051 G4double nrjTransf21 = 0;
1081 G4double nrjTransf22 = 0; 1052 G4double nrjTransf22 = 0;
1082 1053
1083 if (particleDefinition == G4Electron::Elect 1054 if (particleDefinition == G4Electron::ElectronDefinition())
1084 { 1055 {
1085 // Protection against out of boundary acc <<
1086 if (k==eTdummyVec.back()) k=k*(1.-1e-12); <<
1087 // <<
1088 <<
1089 // k should be in eV 1056 // k should be in eV
1090 auto k2 = std::upper_bound(eTdummyVec.beg << 1057 std::vector<double>::iterator k2 = std::upper_bound(eTdummyVec.begin(),
1091 1058 eTdummyVec.end(),
1092 1059 k);
1093 auto k1 = k2 - 1; << 1060 std::vector<double>::iterator k1 = k2 - 1;
1094 1061
1095 /* 1062 /*
1096 G4cout << "----> k=" << k 1063 G4cout << "----> k=" << k
1097 << " " << *k1 1064 << " " << *k1
1098 << " " << *k2 1065 << " " << *k2
1099 << " " << random 1066 << " " << random
1100 << " " << ionizationLevelIndex 1067 << " " << ionizationLevelIndex
1101 << " " << eProbaShellMap[ionizationLevel 1068 << " " << eProbaShellMap[ionizationLevelIndex][(*k1)].back()
1102 << " " << eProbaShellMap[ionizationLevel 1069 << " " << eProbaShellMap[ionizationLevelIndex][(*k2)].back()
1103 << G4endl; 1070 << G4endl;
1104 */ 1071 */
1105 1072
1106 // SI : the following condition avoids si 1073 // SI : the following condition avoids situations where random >last vector element
1107 if (random <= eProbaShellMap[ionizationLe 1074 if (random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
1108 && random <= eProbaShellMap[ionizatio 1075 && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back())
1109 { 1076 {
1110 auto prob12 = << 1077 std::vector<double>::iterator prob12 =
1111 std::upper_bound(eProbaShellMap[ion 1078 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
1112 eProbaShellMap[ion 1079 eProbaShellMap[ionizationLevelIndex][(*k1)].end(),
1113 random); 1080 random);
1114 1081
1115 auto prob11 = prob12 - 1; << 1082 std::vector<double>::iterator prob11 = prob12 - 1;
1116 1083
1117 auto prob22 = << 1084 std::vector<double>::iterator prob22 =
1118 std::upper_bound(eProbaShellMap[ion 1085 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1119 eProbaShellMap[ion 1086 eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1120 random); 1087 random);
1121 1088
1122 auto prob21 = prob22 - 1; << 1089 std::vector<double>::iterator prob21 = prob22 - 1;
1123 1090
1124 valueK1 = *k1; 1091 valueK1 = *k1;
1125 valueK2 = *k2; 1092 valueK2 = *k2;
1126 valuePROB21 = *prob21; 1093 valuePROB21 = *prob21;
1127 valuePROB22 = *prob22; 1094 valuePROB22 = *prob22;
1128 valuePROB12 = *prob12; 1095 valuePROB12 = *prob12;
1129 valuePROB11 = *prob11; 1096 valuePROB11 = *prob11;
1130 1097
1131 /* 1098 /*
1132 G4cout << " " << random << " " 1099 G4cout << " " << random << " " << valuePROB11 << " "
1133 << valuePROB12 << " " << valuePROB21 < 1100 << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1134 */ 1101 */
1135 1102
1136 nrjTransf11 = eNrjTransfData[ionization 1103 nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1137 nrjTransf12 = eNrjTransfData[ionization 1104 nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1138 nrjTransf21 = eNrjTransfData[ionization 1105 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1139 nrjTransf22 = eNrjTransfData[ionization 1106 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1140 1107
1141 /* 1108 /*
1142 G4cout << " " << ionizationLeve 1109 G4cout << " " << ionizationLevelIndex << " "
1143 << random << " " <<valueK1 << " " << v 1110 << random << " " <<valueK1 << " " << valueK2 << G4endl;
1144 1111
1145 G4cout << " " << random << " " 1112 G4cout << " " << random << " " << nrjTransf11 << " "
1146 << nrjTransf12 << " " << nrjTransf21 < 1113 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1147 */ 1114 */
1148 1115
1149 } 1116 }
1150 <<
1151 // Avoids cases where cum xs is zero for 1117 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
1152 if (random > eProbaShellMap[ionizationLev 1118 if (random > eProbaShellMap[ionizationLevelIndex][(*k1)].back())
1153 { 1119 {
1154 auto prob22 = << 1120 std::vector<double>::iterator prob22 =
1155 std::upper_bound(eProbaShellMap[ion 1121 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1156 eProbaShellMap[ion 1122 eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1157 random); 1123 random);
1158 1124
1159 auto prob21 = prob22 - 1; << 1125 std::vector<double>::iterator prob21 = prob22 - 1;
1160 1126
1161 valueK1 = *k1; 1127 valueK1 = *k1;
1162 valueK2 = *k2; 1128 valueK2 = *k2;
1163 valuePROB21 = *prob21; 1129 valuePROB21 = *prob21;
1164 valuePROB22 = *prob22; 1130 valuePROB22 = *prob22;
1165 1131
1166 //G4cout << " " << random << " " 1132 //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1167 1133
1168 nrjTransf21 = eNrjTransfData[ionization 1134 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1169 nrjTransf22 = eNrjTransfData[ionization 1135 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1170 1136
1171 G4double interpolatedvalue2 = Interpola 1137 G4double interpolatedvalue2 = Interpolate(valuePROB21,
1172 1138 valuePROB22,
1173 1139 random,
1174 1140 nrjTransf21,
1175 1141 nrjTransf22);
1176 1142
1177 // zeros are explicitly set << 1143 // zeros are explicitely set
1178 1144
1179 G4double value = Interpolate(valueK1, v 1145 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1180 1146
1181 /* 1147 /*
1182 G4cout << " " << ionizationLeve 1148 G4cout << " " << ionizationLevelIndex << " "
1183 << random << " " <<valueK1 << " " << v 1149 << random << " " <<valueK1 << " " << valueK2 << G4endl;
1184 1150
1185 G4cout << " " << random << " " 1151 G4cout << " " << random << " " << nrjTransf11 << " "
1186 << nrjTransf12 << " " << nrjTransf21 < 1152 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1187 1153
1188 G4cout << "ici" << " " << value << G4e 1154 G4cout << "ici" << " " << value << G4endl;
1189 */ 1155 */
1190 1156
1191 return value; 1157 return value;
1192 } 1158 }
1193 } 1159 }
1194 // 1160 //
1195 else if (particleDefinition == G4Proton::Pr 1161 else if (particleDefinition == G4Proton::ProtonDefinition())
1196 { 1162 {
1197 // Protection against out of boundary acc <<
1198 if (k==pTdummyVec.back()) k=k*(1.-1e-12); <<
1199 // <<
1200 <<
1201 // k should be in eV 1163 // k should be in eV
1202 1164
1203 auto k2 = std::upper_bound(pTdummyVec.beg << 1165 std::vector<double>::iterator k2 = std::upper_bound(pTdummyVec.begin(),
1204 1166 pTdummyVec.end(),
1205 1167 k);
1206 1168
1207 auto k1 = k2 - 1; << 1169 std::vector<double>::iterator k1 = k2 - 1;
1208 1170
1209 /* 1171 /*
1210 G4cout << "----> k=" << k 1172 G4cout << "----> k=" << k
1211 << " " << *k1 1173 << " " << *k1
1212 << " " << *k2 1174 << " " << *k2
1213 << " " << random 1175 << " " << random
1214 << " " << ionizationLevelIndex 1176 << " " << ionizationLevelIndex
1215 << " " << pProbaShellMap[ionizationLevel 1177 << " " << pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1216 << " " << pProbaShellMap[ionizationLevel 1178 << " " << pProbaShellMap[ionizationLevelIndex][(*k2)].back()
1217 << G4endl; 1179 << G4endl;
1218 */ 1180 */
1219 1181
1220 // SI : the following condition avoids si 1182 // SI : the following condition avoids situations where random > last vector element,
1221 // for eg. when the last element is 1183 // for eg. when the last element is zero
1222 if (random <= pProbaShellMap[ionizationLe 1184 if (random <= pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1223 && random <= pProbaShellMap[ionizatio 1185 && random <= pProbaShellMap[ionizationLevelIndex][(*k2)].back())
1224 { 1186 {
1225 auto prob12 = << 1187 std::vector<double>::iterator prob12 =
1226 std::upper_bound(pProbaShellMap[ion 1188 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
1227 pProbaShellMap[ion 1189 pProbaShellMap[ionizationLevelIndex][(*k1)].end(),
1228 random); 1190 random);
1229 1191
1230 auto prob11 = prob12 - 1; << 1192 std::vector<double>::iterator prob11 = prob12 - 1;
1231 1193
1232 auto prob22 = << 1194 std::vector<double>::iterator prob22 =
1233 std::upper_bound(pProbaShellMap[ion 1195 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1234 pProbaShellMap[ion 1196 pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1235 random); 1197 random);
1236 1198
1237 auto prob21 = prob22 - 1; << 1199 std::vector<double>::iterator prob21 = prob22 - 1;
1238 1200
1239 valueK1 = *k1; 1201 valueK1 = *k1;
1240 valueK2 = *k2; 1202 valueK2 = *k2;
1241 valuePROB21 = *prob21; 1203 valuePROB21 = *prob21;
1242 valuePROB22 = *prob22; 1204 valuePROB22 = *prob22;
1243 valuePROB12 = *prob12; 1205 valuePROB12 = *prob12;
1244 valuePROB11 = *prob11; 1206 valuePROB11 = *prob11;
1245 1207
1246 /* 1208 /*
1247 G4cout << " " << random << " " 1209 G4cout << " " << random << " " << valuePROB11 << " "
1248 << valuePROB12 << " " << valuePROB21 < 1210 << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1249 */ 1211 */
1250 1212
1251 nrjTransf11 = pNrjTransfData[ionization 1213 nrjTransf11 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1252 nrjTransf12 = pNrjTransfData[ionization 1214 nrjTransf12 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1253 nrjTransf21 = pNrjTransfData[ionization 1215 nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1254 nrjTransf22 = pNrjTransfData[ionization 1216 nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1255 1217
1256 /* 1218 /*
1257 G4cout << " " << ionizationLeve 1219 G4cout << " " << ionizationLevelIndex << " "
1258 << random << " " <<valueK1 << " " << v 1220 << random << " " <<valueK1 << " " << valueK2 << G4endl;
1259 1221
1260 G4cout << " " << random << " " 1222 G4cout << " " << random << " " << nrjTransf11 << " "
1261 << nrjTransf12 << " " << nrjTransf21 < 1223 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1262 */ 1224 */
1263 } 1225 }
1264 1226
1265 // Avoids cases where cum xs is zero for 1227 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
1266 1228
1267 if (random > pProbaShellMap[ionizationLev 1229 if (random > pProbaShellMap[ionizationLevelIndex][(*k1)].back())
1268 { 1230 {
1269 auto prob22 = << 1231 std::vector<double>::iterator prob22 =
1270 std::upper_bound(pProbaShellMap[ion 1232 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1271 pProbaShellMap[ion 1233 pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1272 random); 1234 random);
1273 1235
1274 auto prob21 = prob22 - 1; << 1236 std::vector<double>::iterator prob21 = prob22 - 1;
1275 1237
1276 valueK1 = *k1; 1238 valueK1 = *k1;
1277 valueK2 = *k2; 1239 valueK2 = *k2;
1278 valuePROB21 = *prob21; 1240 valuePROB21 = *prob21;
1279 valuePROB22 = *prob22; 1241 valuePROB22 = *prob22;
1280 1242
1281 //G4cout << " " << random << " " 1243 //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1282 1244
1283 nrjTransf21 = pNrjTransfData[ionization 1245 nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1284 nrjTransf22 = pNrjTransfData[ionization 1246 nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1285 1247
1286 G4double interpolatedvalue2 = Interpola 1248 G4double interpolatedvalue2 = Interpolate(valuePROB21,
1287 1249 valuePROB22,
1288 1250 random,
1289 1251 nrjTransf21,
1290 1252 nrjTransf22);
1291 1253
1292 // zeros are explicitly set << 1254 // zeros are explicitely set
1293 1255
1294 G4double value = Interpolate(valueK1, v 1256 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1295 1257
1296 /* 1258 /*
1297 G4cout << " " << ionizationLeve 1259 G4cout << " " << ionizationLevelIndex << " "
1298 << random << " " <<valueK1 << " " << v 1260 << random << " " <<valueK1 << " " << valueK2 << G4endl;
1299 1261
1300 G4cout << " " << random << " " 1262 G4cout << " " << random << " " << nrjTransf11 << " "
1301 << nrjTransf12 << " " << nrjTransf21 < 1263 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1302 1264
1303 G4cout << "ici" << " " << value << G4e 1265 G4cout << "ici" << " " << value << G4endl;
1304 */ 1266 */
1305 1267
1306 return value; 1268 return value;
1307 } 1269 }
1308 } 1270 }
1309 // End electron and proton cases 1271 // End electron and proton cases
1310 1272
1311 G4double nrjTransfProduct = nrjTransf11 * n 1273 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21
1312 * nrjTransf22; 1274 * nrjTransf22;
1313 //G4cout << "nrjTransfProduct=" << nrjTrans 1275 //G4cout << "nrjTransfProduct=" << nrjTransfProduct << G4endl;
1314 1276
1315 if (nrjTransfProduct != 0.) 1277 if (nrjTransfProduct != 0.)
1316 { 1278 {
1317 nrj = QuadInterpolator(valuePROB11, 1279 nrj = QuadInterpolator(valuePROB11,
1318 valuePROB12, 1280 valuePROB12,
1319 valuePROB21, 1281 valuePROB21,
1320 valuePROB22, 1282 valuePROB22,
1321 nrjTransf11, 1283 nrjTransf11,
1322 nrjTransf12, 1284 nrjTransf12,
1323 nrjTransf21, 1285 nrjTransf21,
1324 nrjTransf22, 1286 nrjTransf22,
1325 valueK1, 1287 valueK1,
1326 valueK2, 1288 valueK2,
1327 k, 1289 k,
1328 random); 1290 random);
1329 } 1291 }
1330 //G4cout << nrj << endl; 1292 //G4cout << nrj << endl;
1331 1293
1332 return nrj; 1294 return nrj;
1333 } 1295 }
1334 1296