Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/dna/models/src/G4DNABornIonisationModel1.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 ]

Diff markup

Differences between /processes/electromagnetic/dna/models/src/G4DNABornIonisationModel1.cc (Version 11.3.0) and /processes/electromagnetic/dna/models/src/G4DNABornIonisationModel1.cc (Version 10.5.p1)


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