Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/dna/models/src/G4DNABornIonisationModel2.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 
 28 #include "G4DNABornIonisationModel2.hh"
 29 #include "G4PhysicalConstants.hh"
 30 #include "G4SystemOfUnits.hh"
 31 #include "G4UAtomicDeexcitation.hh"
 32 #include "G4LossTableManager.hh"
 33 #include "G4DNAChemistryManager.hh"
 34 #include "G4DNAMolecularMaterial.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 
 45 G4DNABornIonisationModel2::G4DNABornIonisationModel2(const G4ParticleDefinition*,
 46                                                      const G4String& nam) :
 47 G4VEmModel(nam) 
 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);
 65   fAtomDeexcitation = nullptr;
 66   fParticleChangeForGamma = nullptr;
 67   fpMolWaterDensity = nullptr;
 68   fTableData = nullptr;
 69   fLowEnergyLimit = 0;
 70   fHighEnergyLimit = 0;
 71   fParticleDef = nullptr;
 72 
 73   // Define default angular generator
 74   
 75   SetAngularDistribution(new G4DNABornAngle());
 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 
 92 G4DNABornIonisationModel2::~G4DNABornIonisationModel2()
 93 {
 94   // Cross section
 95 
 96   
 97     delete fTableData;
 98 
 99   // Final state
100 
101   fVecm.clear();
102 }
103 
104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
105 
106 void G4DNABornIonisationModel2::Initialise(const G4ParticleDefinition* particle,
107                                            const G4DataVector& /*cuts*/)
108 {
109 
110   if (verboseLevel > 3)
111   {
112     G4cout << "Calling G4DNABornIonisationModel2::Initialise()" << G4endl;
113   }
114 
115   if(fParticleDef != nullptr && 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 initialised "
123     "for particle:" << fParticleDef->GetParticleName() << G4endl;
124     G4Exception("G4DNABornIonisationModel2::Initialise","bornIonInit",
125         FatalException,description);
126   }
127 
128   fParticleDef = particle;
129 
130   // Energy limits
131   const char* path = G4FindDataDir("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;
173   fTableData = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
174   fTableData->LoadData(fTableFile);
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   //
233   SetLowEnergyLimit(fLowEnergyLimit);
234   SetHighEnergyLimit(fHighEnergyLimit);
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   
248   fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->
249   GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
250 
251   // AD
252   
253   fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
254 
255   if (isInitialised)
256   { return;}
257   fParticleChangeForGamma = GetParticleChangeForGamma();
258   isInitialised = true;
259 }
260 
261 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
262 
263 G4double G4DNABornIonisationModel2::CrossSectionPerVolume(const G4Material* material,
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,
315                                                   const G4DynamicParticle* particle,
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)
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       auto  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 
391       fParticleChangeForGamma->ProposeMomentumDirection(direction.unit());
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     std::size_t secNumberInit = 0;
401     std::size_t secNumberFinal = 0;
402 
403     G4double bindingEnergy = 0;
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 != nullptr) && ionizationShell == 4)
414     {
415       const G4AtomicShell* shell = 
416         fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0));
417       secNumberInit = fvect->size();
418       fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
419       secNumberFinal = fvect->size();
420 
421       if(secNumberFinal > secNumberInit) 
422       {
423   for (std::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]=nullptr;
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     {
452       fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
453       fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
454     }
455     else
456     {
457       fParticleChangeForGamma->SetProposedKineticEnergy(k);
458       fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy);
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();
466     /////////////////////////////////
467 
468     const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
469     G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
470         ionizationShell,
471         theIncomingTrack);
472   }
473 }
474 
475 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
476 
477 G4double G4DNABornIonisationModel2::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
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   if (particleDefinition == G4Proton::ProtonDefinition())
539   {
540     G4double maximumKineticEnergyTransfer = 4.
541         * (electron_mass_c2 / proton_mass_c2) * k;
542 
543     G4double crossSectionMaximum = 0.;
544     for (G4double value = waterStructure.IonisationEnergy(shell);
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......
615 G4double G4DNABornIonisationModel2::DifferentialCrossSection(G4ParticleDefinition * /*particleDefinition*/,
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     auto t2 = std::upper_bound(fTdummyVec.begin(),
643                                                         fTdummyVec.end(),
644                                                         k);
645 
646     auto 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       auto e12 = std::upper_bound(fVecm[(*t1)].begin(),
654                                                            fVecm[(*t1)].end(),
655                                                            energyTransfer);
656       auto e11 = e12 - 1;
657 
658       auto e22 = std::upper_bound(fVecm[(*t2)].begin(),
659                                                            fVecm[(*t2)].end(),
660                                                            energyTransfer);
661       auto 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 
700 G4double G4DNABornIonisationModel2::Interpolate(G4double e1,
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 
764 G4double G4DNABornIonisationModel2::QuadInterpolator(G4double e11,
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);
779   G4double value = Interpolate(t1,
780                                t2,
781                                t,
782                                interpolatedvalue1,
783                                interpolatedvalue2);
784 
785   return value;
786 }
787 
788 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
789 
790 G4double G4DNABornIonisationModel2::GetPartialCrossSection(const G4Material*,
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 
800 G4int G4DNABornIonisationModel2::RandomSelect(G4double k)
801 {
802   G4int level = 0;
803 
804   auto  valuesBuffer = new G4double[fTableData->NumberOfComponents()];
805   const auto  n = (G4int)fTableData->NumberOfComponents();
806   G4int 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   
833     delete[] valuesBuffer;
834 
835   return level;
836 }
837 
838 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
839 
840 G4double G4DNABornIonisationModel2::RandomizeEjectedElectronEnergyFromCumulatedDcs(G4ParticleDefinition* particleDefinition,
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
854       - waterStructure.IonisationEnergy(shell);
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 
866 G4double G4DNABornIonisationModel2::TransferedEnergy(G4ParticleDefinition* /*particleDefinition*/,
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   auto k2 = std::upper_bound(fTdummyVec.begin(),
892                                                       fTdummyVec.end(),
893                                                       k);
894   auto 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     auto prob12 =
912         std::upper_bound(fProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
913                          fProbaShellMap[ionizationLevelIndex][(*k1)].end(),
914                          random);
915 
916     auto prob11 = prob12 - 1;
917 
918     auto prob22 =
919         std::upper_bound(fProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
920                          fProbaShellMap[ionizationLevelIndex][(*k2)].end(),
921                          random);
922 
923     auto 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     auto prob22 =
955         std::upper_bound(fProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
956                          fProbaShellMap[ionizationLevelIndex][(*k2)].end(),
957                          random);
958 
959     auto 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 }
1019