Geant4 Cross Reference

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

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

Diff markup

Differences between /processes/electromagnetic/dna/models/src/G4DNABornIonisationModel2.cc (Version 11.3.0) and /processes/electromagnetic/dna/models/src/G4DNABornIonisationModel2.cc (Version 10.6.p3)


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