ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNAEmfietzoglouIonisationModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNAEmfietzoglouIonisationModel.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Based on the work described in
27 // Rad Res 163, 98-111 (2005)
28 // D. Emfietzoglou, H. Nikjoo
29 //
30 // Authors of the class (2014):
31 // I. Kyriakou (kyriak@cc.uoi.gr)
32 // D. Emfietzoglou (demfietz@cc.uoi.gr)
33 // S. Incerti (incerti@cenbg.in2p3.fr)
34 //
35 
36 #include "G4DNAEmfietzoglouIonisationModel.hh"
37 #include "G4PhysicalConstants.hh"
38 #include "G4SystemOfUnits.hh"
39 #include "G4UAtomicDeexcitation.hh"
40 #include "G4LossTableManager.hh"
41 #include "G4DNAChemistryManager.hh"
42 #include "G4DNAMolecularMaterial.hh"
43 #include "G4DNABornAngle.hh"
44 #include "G4DeltaAngle.hh"
45 
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47 
48 using namespace std;
49 
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51 
52 G4DNAEmfietzoglouIonisationModel::G4DNAEmfietzoglouIonisationModel(const G4ParticleDefinition*,
53  const G4String& nam) :
54 G4VEmModel(nam), isInitialised(false)
55 {
56  verboseLevel = 0;
57  // Verbosity scale:
58  // 0 = nothing
59  // 1 = warning for energy non-conservation
60  // 2 = details of energy budget
61  // 3 = calculation of cross sections, file openings, sampling of atoms
62  // 4 = entering in methods
63 
64  if(verboseLevel > 0)
65  {
66  G4cout << "Emfietzoglou ionisation model is constructed " << G4endl;
67  }
68 
69  // Mark this model as "applicable" for atomic deexcitation
70  SetDeexcitationFlag(true);
71  fAtomDeexcitation = 0;
72  fParticleChangeForGamma = 0;
73  fpMolWaterDensity = 0;
74 
75  // Define default angular generator
76  SetAngularDistribution(new G4DNABornAngle());
77 
78  SetLowEnergyLimit(10. * eV);
79  SetHighEnergyLimit(10. * keV);
80 
81  // Selection of computation method
82 
83  fasterCode = false;
84 
85  // Selection of stationary mode
86 
87  statCode = false;
88 }
89 
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91 
92 G4DNAEmfietzoglouIonisationModel::~G4DNAEmfietzoglouIonisationModel()
93 {
94  // Cross section
95 
96  std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
97  for(pos = tableData.begin(); pos != tableData.end(); ++pos)
98  {
99  G4DNACrossSectionDataSet* table = pos->second;
100  delete table;
101  }
102 
103  // Final state
104 
105  eVecm.clear();
106 
107 }
108 
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
110 
111 void G4DNAEmfietzoglouIonisationModel::Initialise(const G4ParticleDefinition* particle,
112  const G4DataVector& /*cuts*/)
113 {
114 
115  if(verboseLevel > 3)
116  {
117  G4cout << "Calling G4DNAEmfietzoglouIonisationModel::Initialise()" << G4endl;
118  }
119 
120  // Energy limits
121 
122  G4String fileElectron("dna/sigma_ionisation_e_emfietzoglou");
123 
124  G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
125 
126  G4String electron;
127 
128  G4double scaleFactor = (1.e-22 / 3.343) * m*m;
129 
130  char *path = getenv("G4LEDATA");
131 
132  // *** ELECTRON
133 
134  electron = electronDef->GetParticleName();
135 
136  tableFile[electron] = fileElectron;
137 
138  // Cross section
139 
140  G4DNACrossSectionDataSet* tableE = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
141  tableE->LoadData(fileElectron);
142 
143  tableData[electron] = tableE;
144 
145  // Final state
146 
147  std::ostringstream eFullFileName;
148 
149  if (fasterCode) eFullFileName << path << "/dna/sigmadiff_cumulated_ionisation_e_emfietzoglou.dat";
150  if (!fasterCode) eFullFileName << path << "/dna/sigmadiff_ionisation_e_emfietzoglou.dat";
151 
152  std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
153 
154  if (!eDiffCrossSection)
155  {
156  if (fasterCode) G4Exception("G4DNAEmfietzoglouIonisationModel::Initialise","em0003",
157  FatalException,"Missing data file:/dna/sigmadiff_cumulated_ionisation_e_emfietzoglou.dat");
158 
159  if (!fasterCode) G4Exception("G4DNAEmfietzoglouIonisationModel::Initialise","em0003",
160  FatalException,"Missing data file:/dna/sigmadiff_ionisation_e_emfietzoglou.dat");
161  }
162 
163  //
164 
165  // Clear the arrays for re-initialization case (MT mode)
166  // March 25th, 2014 - Vaclav Stepan, Sebastien Incerti
167 
168  eTdummyVec.clear();
169 
170  eVecm.clear();
171 
172  eProbaShellMap->clear();
173 
174  eDiffCrossSectionData->clear();
175 
176  eNrjTransfData->clear();
177 
178  //
179 
180  eTdummyVec.push_back(0.);
181  while(!eDiffCrossSection.eof())
182  {
183  G4double tDummy;
184  G4double eDummy;
185  eDiffCrossSection>>tDummy>>eDummy;
186  if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
187  for (G4int j=0; j<5; j++)
188  {
189  eDiffCrossSection>>eDiffCrossSectionData[j][tDummy][eDummy];
190 
191  if (fasterCode)
192  {
193  eNrjTransfData[j][tDummy][eDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
194  eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]);
195  }
196 
197  // SI - only if eof is not reached
198  if (!eDiffCrossSection.eof() && !fasterCode)
199  {
200  eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
201  }
202 
203  if (!fasterCode) eVecm[tDummy].push_back(eDummy);
204 
205  }
206  }
207 
208  //
209 
210  if( verboseLevel>0 )
211  {
212  G4cout << "Emfietzoglou ionisation model is initialized " << G4endl
213  << "Energy range: "
214  << LowEnergyLimit() / eV << " eV - "
215  << HighEnergyLimit() / keV << " keV for "
216  << particle->GetParticleName()
217  << G4endl;
218  }
219 
220  // Initialize water density pointer
221 
222  fpMolWaterDensity =
223  G4DNAMolecularMaterial::Instance()->
224  GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
225 
226  // AD
227 
228  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
229 
230  if (isInitialised)
231  { return;}
232  fParticleChangeForGamma = GetParticleChangeForGamma();
233  isInitialised = true;
234 }
235 
236 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
237 
238 G4double G4DNAEmfietzoglouIonisationModel::
239 CrossSectionPerVolume(const G4Material* material,
240  const G4ParticleDefinition* particleDefinition,
241  G4double ekin,
242  G4double,
243  G4double)
244 {
245  if(verboseLevel > 3)
246  {
247  G4cout
248  << "Calling CrossSectionPerVolume() of G4DNAEmfietzoglouIonisationModel"
249  << G4endl;
250  }
251 
252  if (particleDefinition != G4Electron::ElectronDefinition()) return 0; // necessary ??
253 
254  // Calculate total cross section for model
255 
256  G4double sigma=0;
257 
258  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
259 
260  const G4String& particleName = particleDefinition->GetParticleName();
261 
262  if (ekin >= LowEnergyLimit() && ekin <= HighEnergyLimit())
263  {
264  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
265  pos = tableData.find(particleName);
266 
267  if (pos != tableData.end())
268  {
269  G4DNACrossSectionDataSet* table = pos->second;
270  if (table != 0)
271  {
272  sigma = table->FindValue(ekin);
273  }
274  }
275  else
276  {
277  G4Exception("G4DNAEmfietzoglouIonisationModel::CrossSectionPerVolume","em0002",
278  FatalException,"Model not applicable to particle type.");
279  }
280  }
281 
282  if (verboseLevel > 2)
283  {
284  G4cout << "__________________________________" << G4endl;
285  G4cout << "G4DNAEmfietzoglouIonisationModel - XS INFO START" << G4endl;
286  G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
287  G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
288  G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
289  G4cout << "G4DNAEmfietzoglouIonisationModel - XS INFO END" << G4endl;
290  }
291 
292  return sigma*waterDensity;
293 }
294 
295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
296 
297 void G4DNAEmfietzoglouIonisationModel::
298 SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
299  const G4MaterialCutsCouple* couple,
300  const G4DynamicParticle* particle,
301  G4double,
302  G4double)
303 {
304 
305  if(verboseLevel > 3)
306  {
307  G4cout << "Calling SampleSecondaries() of G4DNAEmfietzoglouIonisationModel"
308  << G4endl;
309  }
310 
311  G4double k = particle->GetKineticEnergy();
312 
313  const G4String& particleName = particle->GetDefinition()->GetParticleName();
314 
315  if (k >= LowEnergyLimit() && k <= HighEnergyLimit())
316  {
317  G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
318  G4double particleMass = particle->GetDefinition()->GetPDGMass();
319  G4double totalEnergy = k + particleMass;
320  G4double pSquare = k * (totalEnergy + particleMass);
321  G4double totalMomentum = std::sqrt(pSquare);
322 
323  G4int ionizationShell = 0;
324 
325  ionizationShell = RandomSelect(k,particleName);
326 
327  G4double bindingEnergy = 0;
328  bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
329 
330  // SI : additional protection if tcs interpolation method is modified
331  if (k<bindingEnergy) return;
332  //
333 
334  G4double secondaryKinetic=-1000*eV;
335 
336  if (!fasterCode) secondaryKinetic = RandomizeEjectedElectronEnergy(particle->GetDefinition(),k,ionizationShell);
337 
338  if (fasterCode)
339  secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(particle->GetDefinition(),k,ionizationShell);
340 
341  // SI - For atom. deexc. tagging - 23/05/2017
342 
343  G4int Z = 8;
344 
345  G4ThreeVector deltaDirection =
346  GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
347  Z, ionizationShell,
348  couple->GetMaterial());
349 
350  if (secondaryKinetic>0)
351  {
352  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic);
353  fvect->push_back(dp);
354  }
355 
356  G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
357 
358  G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
359  G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
360  G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
361  G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
362  finalPx /= finalMomentum;
363  finalPy /= finalMomentum;
364  finalPz /= finalMomentum;
365 
366  G4ThreeVector direction;
367  direction.set(finalPx,finalPy,finalPz);
368 
369  fParticleChangeForGamma->ProposeMomentumDirection(direction.unit());
370 
371  // AM: sample deexcitation
372  // here we assume that H_{2}O electronic levels are the same as Oxygen.
373  // this can be considered true with a rough 10% error in energy on K-shell,
374 
375  size_t secNumberInit = 0;// need to know at a certain point the energy of secondaries
376  size_t secNumberFinal = 0;// So I'll make the diference and then sum the energies
377 
378  G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
379 
380  // SI: only atomic deexcitation from K shell is considered
381  if(fAtomDeexcitation && ionizationShell == 4)
382  {
383  const G4AtomicShell* shell
384  = fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0));
385  secNumberInit = fvect->size();
386  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
387  secNumberFinal = fvect->size();
388 
389  if(secNumberFinal > secNumberInit) {
390  for (size_t i=secNumberInit; i<secNumberFinal; ++i) {
391  //Check if there is enough residual energy
392  if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy())
393  {
394  //Ok, this is a valid secondary: keep it
395  bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
396  }
397  else
398  {
399  //Invalid secondary: not enough energy to create it!
400  //Keep its energy in the local deposit
401  delete (*fvect)[i];
402  (*fvect)[i]=0;
403  }
404  }
405  }
406 
407  }
408 
409  //This should never happen
410  if(bindingEnergy < 0.0)
411  G4Exception("G4DNAEmfietzoglouIonisatioModel1::SampleSecondaries()",
412  "em2050",FatalException,"Negative local energy deposit");
413 
414  //bindingEnergy has been decreased
415  //by the amount of energy taken away by deexc. products
416  if (!statCode)
417  {
418  fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
419  fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
420  }
421  else
422  {
423  fParticleChangeForGamma->SetProposedKineticEnergy(k);
424  fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy);
425  }
426  // TEST //////////////////////////
427  // if (secondaryKinetic<0) abort();
428  // if (scatteredEnergy<0) abort();
429  // if (k-scatteredEnergy-secondaryKinetic-deexSecEnergy<0) abort();
430  // if (k-scatteredEnergy<0) abort();
432 
433  const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
434  G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
435  ionizationShell,
436  theIncomingTrack);
437  }
438 
439 }
440 
441 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
442 
443 G4double
444 G4DNAEmfietzoglouIonisationModel::
445 RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
446  G4double k,
447  G4int shell)
448 {
449  // G4cout << "*** SLOW computation for "
450  // << " " << particleDefinition->GetParticleName() << G4endl;
451 
452  if(particleDefinition == G4Electron::ElectronDefinition())
453  {
454  G4double maximumEnergyTransfer = 0.;
455  if((k + waterStructure.IonisationEnergy(shell)) / 2. > k)
456  maximumEnergyTransfer = k;
457  else
458  maximumEnergyTransfer = (k + waterStructure.IonisationEnergy(shell))/ 2.;
459 
460  // SI : original method
461  /*
462  G4double crossSectionMaximum = 0.;
463  for(G4double value=waterStructure.IonisationEnergy(shell); value<=maximumEnergyTransfer; value+=0.1*eV)
464  {
465  G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
466  if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
467  }
468  */
469 
470  // SI : alternative method
471  G4double crossSectionMaximum = 0.;
472 
473  G4double minEnergy = waterStructure.IonisationEnergy(shell);
474  G4double maxEnergy = maximumEnergyTransfer;
475  G4int nEnergySteps = 50;
476 
477  G4double value(minEnergy);
478  G4double stpEnergy(std::pow(maxEnergy / value,
479  1. / static_cast<G4double>(nEnergySteps - 1)));
480  G4int step(nEnergySteps);
481  while(step > 0)
482  {
483  step--;
484  G4double differentialCrossSection =
485  DifferentialCrossSection(particleDefinition,
486  k / eV,
487  value / eV,
488  shell);
489  if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum =
490  differentialCrossSection;
491  value *= stpEnergy;
492  }
493  //
494 
495  G4double secondaryElectronKineticEnergy = 0.;
496  do
497  {
498  secondaryElectronKineticEnergy = G4UniformRand()* (maximumEnergyTransfer-waterStructure.IonisationEnergy(shell));
499  }while(G4UniformRand()*crossSectionMaximum >
500  DifferentialCrossSection(particleDefinition, k/eV,
501  (secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
502 
503  return secondaryElectronKineticEnergy;
504 
505  }
506 
507  return 0;
508 }
509 
510 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
511 
512 // The following section is not used anymore but is kept for memory
513 // GetAngularDistribution()->SampleDirectionForShell is used instead
514 
515 /*
516  void G4DNAEmfietzoglouIonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
517  G4double k,
518  G4double secKinetic,
519  G4double & cosTheta,
520  G4double & phi )
521  {
522  if (particleDefinition == G4Electron::ElectronDefinition())
523  {
524  phi = twopi * G4UniformRand();
525  if (secKinetic < 50.*eV) cosTheta = (2.*G4UniformRand())-1.;
526  else if (secKinetic <= 200.*eV)
527  {
528  if (G4UniformRand() <= 0.1) cosTheta = (2.*G4UniformRand())-1.;
529  else cosTheta = G4UniformRand()*(std::sqrt(2.)/2);
530  }
531  else
532  {
533  G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
534  cosTheta = std::sqrt(1.-sin2O);
535  }
536  }
537 
538  else if (particleDefinition == G4Proton::ProtonDefinition())
539  {
540  G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
541  phi = twopi * G4UniformRand();
542 
543  // cosTheta = std::sqrt(secKinetic / maxSecKinetic);
544 
545  // Restriction below 100 eV from Emfietzoglou (2000)
546 
547  if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
548  else cosTheta = (2.*G4UniformRand())-1.;
549 
550  }
551  }
552  */
553 
554 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
555 G4double G4DNAEmfietzoglouIonisationModel::DifferentialCrossSection(G4ParticleDefinition * particleDefinition,
556  G4double k,
557  G4double energyTransfer,
558  G4int ionizationLevelIndex)
559 {
560  G4double sigma = 0.;
561 
562  if(energyTransfer >= waterStructure.IonisationEnergy(ionizationLevelIndex)/eV)
563  {
564  G4double valueT1 = 0;
565  G4double valueT2 = 0;
566  G4double valueE21 = 0;
567  G4double valueE22 = 0;
568  G4double valueE12 = 0;
569  G4double valueE11 = 0;
570 
571  G4double xs11 = 0;
572  G4double xs12 = 0;
573  G4double xs21 = 0;
574  G4double xs22 = 0;
575 
576  if(particleDefinition == G4Electron::ElectronDefinition())
577  {
578  // Protection against out of boundary access
579  if (k==eTdummyVec.back()) k=k*(1.-1e-12);
580  //
581 
582  // k should be in eV and energy transfer eV also
583 
584  std::vector<G4double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),
585  eTdummyVec.end(),
586  k);
587 
588  std::vector<G4double>::iterator t1 = t2 - 1;
589 
590  // SI : the following condition avoids situations where energyTransfer >last vector element
591  // added strict limitations (09/08/2017)
592  if(energyTransfer < eVecm[(*t1)].back() &&
593  energyTransfer < eVecm[(*t2)].back())
594  {
595  std::vector<G4double>::iterator e12 =
596  std::upper_bound(eVecm[(*t1)].begin(),
597  eVecm[(*t1)].end(),
598  energyTransfer);
599  std::vector<G4double>::iterator e11 = e12 - 1;
600 
601  std::vector<G4double>::iterator e22 =
602  std::upper_bound(eVecm[(*t2)].begin(),
603  eVecm[(*t2)].end(),
604  energyTransfer);
605  std::vector<G4double>::iterator e21 = e22 - 1;
606 
607  valueT1 = *t1;
608  valueT2 = *t2;
609  valueE21 = *e21;
610  valueE22 = *e22;
611  valueE12 = *e12;
612  valueE11 = *e11;
613 
614  xs11 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
615  xs12 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
616  xs21 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
617  xs22 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
618 
619  //G4cout << "-------------------" << G4endl;
620  //G4cout << "ionizationLevelIndex=" << ionizationLevelIndex << G4endl;
621  //G4cout << "valueT1/eV=" << valueT1 << " valueT2/eV=" << valueT2 << G4endl;
622  //G4cout << "valueE11/eV=" << valueE11 << " valueE12/eV=" << valueE12
623  // << " valueE21/eV=" << valueE21 << " valueE22/eV=" << valueE22 << G4endl;
624  //G4cout << "xs11=" << xs11 / ((1.e-22 / 3.343) * m*m) << G4endl;
625  //G4cout << "xs12=" << xs12 / ((1.e-22 / 3.343) * m*m) << G4endl;
626  //G4cout << "xs21=" << xs21 / ((1.e-22 / 3.343) * m*m) << G4endl;
627  //G4cout << "xs22=" << xs22 / ((1.e-22 / 3.343) * m*m) << G4endl;
628  //G4cout << "###################" << G4endl;
629 
630  }
631 
632  }
633 
634  G4double xsProduct = xs11 * xs12 * xs21 * xs22;
635  if(xsProduct != 0.)
636  {
637  sigma = QuadInterpolator(valueE11,
638  valueE12,
639  valueE21,
640  valueE22,
641  xs11,
642  xs12,
643  xs21,
644  xs22,
645  valueT1,
646  valueT2,
647  k,
648  energyTransfer);
649  }
650 
651  }
652 
653  return sigma;
654 }
655 
656 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
657 
658 G4double G4DNAEmfietzoglouIonisationModel::Interpolate(G4double e1,
659  G4double e2,
660  G4double e,
661  G4double xs1,
662  G4double xs2)
663 {
664 
665  G4double value = 0.;
666 
667  // Log-log interpolation by default
668 
669  if(e1 != 0 && e2 != 0 && (std::log10(e2) - std::log10(e1)) != 0
670  && !fasterCode)
671  {
672  G4double a = (std::log10(xs2) - std::log10(xs1))
673  / (std::log10(e2) - std::log10(e1));
674  G4double b = std::log10(xs2) - a * std::log10(e2);
675  G4double sigma = a * std::log10(e) + b;
676  value = (std::pow(10., sigma));
677  }
678 
679  // Switch to lin-lin interpolation
680  /*
681  if ((e2-e1)!=0)
682  {
683  G4double d1 = xs1;
684  G4double d2 = xs2;
685  value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
686  }
687  */
688 
689  // Switch to log-lin interpolation for faster code
690  if((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 && fasterCode)
691  {
692  G4double d1 = std::log10(xs1);
693  G4double d2 = std::log10(xs2);
694  value = std::pow(10., (d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
695  }
696 
697  // Switch to lin-lin interpolation for faster code
698  // in case one of xs1 or xs2 (=cum proba) value is zero
699 
700  if((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) && fasterCode)
701  {
702  G4double d1 = xs1;
703  G4double d2 = xs2;
704  value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
705  }
706 
707  /*
708  G4cout
709  << e1 << " "
710  << e2 << " "
711  << e << " "
712  << xs1 << " "
713  << xs2 << " "
714  << value
715  << G4endl;
716  */
717 
718  return value;
719 }
720 
721 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
722 
723 G4double G4DNAEmfietzoglouIonisationModel::QuadInterpolator(G4double e11,
724  G4double e12,
725  G4double e21,
726  G4double e22,
727  G4double xs11,
728  G4double xs12,
729  G4double xs21,
730  G4double xs22,
731  G4double t1,
732  G4double t2,
733  G4double t,
734  G4double e)
735 {
736  G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
737  G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
738  G4double value = Interpolate(t1,
739  t2,
740  t,
741  interpolatedvalue1,
742  interpolatedvalue2);
743 
744  return value;
745 }
746 
747 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
748 
749 G4int G4DNAEmfietzoglouIonisationModel::RandomSelect(G4double k,
750  const G4String& particle)
751 {
752  G4int level = 0;
753 
754  std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
755  pos = tableData.find(particle);
756 
757  if(pos != tableData.end())
758  {
759  G4DNACrossSectionDataSet* table = pos->second;
760 
761  if(table != 0)
762  {
763  G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
764  const size_t n(table->NumberOfComponents());
765  size_t i(n);
766  G4double value = 0.;
767 
768  while(i > 0)
769  {
770  i--;
771  valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
772  value += valuesBuffer[i];
773  }
774 
775  value *= G4UniformRand();
776 
777  i = n;
778 
779  while(i > 0)
780  {
781  i--;
782 
783  if(valuesBuffer[i] > value)
784  {
785  delete[] valuesBuffer;
786  return i;
787  }
788  value -= valuesBuffer[i];
789  }
790 
791  if(valuesBuffer) delete[] valuesBuffer;
792 
793  }
794  }
795  else
796  {
797  G4Exception("G4DNAEmfietzoglouIonisationModel::RandomSelect",
798  "em0002",
799  FatalException,
800  "Model not applicable to particle type.");
801  }
802 
803  return level;
804 }
805 
806 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
807 
808 G4double G4DNAEmfietzoglouIonisationModel::RandomizeEjectedElectronEnergyFromCumulatedDcs(G4ParticleDefinition* particleDefinition,
809  G4double k,
810  G4int shell)
811 {
812  //G4cout << "*** FAST computation for " << " " << particleDefinition->GetParticleName() << G4endl;
813 
814  G4double secondaryElectronKineticEnergy = 0.;
815 
816  secondaryElectronKineticEnergy = RandomTransferedEnergy(particleDefinition,
817  k / eV,
818  shell)
819  * eV
820  - waterStructure.IonisationEnergy(shell);
821 
822  //G4cout << RandomTransferedEnergy(particleDefinition, k/eV, shell) << G4endl;
823  if(secondaryElectronKineticEnergy < 0.) return 0.;
824 
825  return secondaryElectronKineticEnergy;
826 }
827 
828 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
829 
830 G4double G4DNAEmfietzoglouIonisationModel::RandomTransferedEnergy(G4ParticleDefinition* particleDefinition,
831  G4double k,
832  G4int ionizationLevelIndex)
833 {
834 
835  G4double random = G4UniformRand();
836 
837  G4double nrj = 0.;
838 
839  G4double valueK1 = 0;
840  G4double valueK2 = 0;
841  G4double valuePROB21 = 0;
842  G4double valuePROB22 = 0;
843  G4double valuePROB12 = 0;
844  G4double valuePROB11 = 0;
845 
846  G4double nrjTransf11 = 0;
847  G4double nrjTransf12 = 0;
848  G4double nrjTransf21 = 0;
849  G4double nrjTransf22 = 0;
850 
851  if (particleDefinition == G4Electron::ElectronDefinition())
852  {
853  // Protection against out of boundary access
854  if (k==eTdummyVec.back()) k=k*(1.-1e-12);
855  //
856 
857  // k should be in eV
858  std::vector<G4double>::iterator k2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
859 
860  std::vector<G4double>::iterator k1 = k2-1;
861 
862  /*
863  G4cout << "----> k=" << k
864  << " " << *k1
865  << " " << *k2
866  << " " << random
867  << " " << ionizationLevelIndex
868  << " " << eProbaShellMap[ionizationLevelIndex][(*k1)].back()
869  << " " << eProbaShellMap[ionizationLevelIndex][(*k2)].back()
870  << G4endl;
871  */
872 
873  // SI : the following condition avoids situations where random >last vector element
874  if ( random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
875  && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back() )
876 
877  {
878  std::vector<G4double>::iterator prob12 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
879  eProbaShellMap[ionizationLevelIndex][(*k1)].end(), random);
880 
881  std::vector<G4double>::iterator prob11 = prob12-1;
882 
883  std::vector<G4double>::iterator prob22 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
884  eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
885 
886  std::vector<G4double>::iterator prob21 = prob22-1;
887 
888  valueK1 =*k1;
889  valueK2 =*k2;
890  valuePROB21 =*prob21;
891  valuePROB22 =*prob22;
892  valuePROB12 =*prob12;
893  valuePROB11 =*prob11;
894 
895  /*
896  G4cout << " " << random << " " << valuePROB11 << " "
897  << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
898  */
899 
900  nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
901  nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
902  nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
903  nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
904 
905  /*
906  G4cout << " " << ionizationLevelIndex << " "
907  << random << " " <<valueK1 << " " << valueK2 << G4endl;
908 
909  G4cout << " " << random << " " << nrjTransf11 << " "
910  << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
911  */
912 
913  }
914 
915  // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
916 
917  if ( random > eProbaShellMap[ionizationLevelIndex][(*k1)].back() )
918 
919  {
920  std::vector<G4double>::iterator prob22 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
921  eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
922 
923  std::vector<G4double>::iterator prob21 = prob22-1;
924 
925  valueK1 =*k1;
926  valueK2 =*k2;
927  valuePROB21 =*prob21;
928  valuePROB22 =*prob22;
929 
930  //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
931 
932  nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
933  nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
934 
935  G4double interpolatedvalue2 = Interpolate(valuePROB21, valuePROB22, random, nrjTransf21, nrjTransf22);
936 
937  // zeros are explicitely set
938 
939  G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
940 
941  /*
942  G4cout << " " << ionizationLevelIndex << " "
943  << random << " " <<valueK1 << " " << valueK2 << G4endl;
944 
945  G4cout << " " << random << " " << nrjTransf11 << " "
946  << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
947 
948  G4cout << "ici" << " " << value << G4endl;
949  */
950 
951  return value;
952  }
953 
954  }
955 
956  // End electron
957 
958  G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21 * nrjTransf22;
959 
960  //G4cout << "nrjTransfProduct=" << nrjTransfProduct << G4endl;
961 
962  if (nrjTransfProduct != 0.)
963  {
964  nrj = QuadInterpolator( valuePROB11, valuePROB12,
965  valuePROB21, valuePROB22,
966  nrjTransf11, nrjTransf12,
967  nrjTransf21, nrjTransf22,
968  valueK1, valueK2,
969  k, random);
970  }
971 
972  //G4cout << nrj << endl;
973 
974  return nrj;
975 }