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