Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/dna/models/src/G4DNACPA100IonisationModel.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/G4DNACPA100IonisationModel.cc (Version 11.3.0) and /processes/electromagnetic/dna/models/src/G4DNACPA100IonisationModel.cc (Version 11.2.1)


  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 // CPA100 ionisation model class for electrons     26 // CPA100 ionisation model class for electrons
 27 //                                                 27 //
 28 // Based on the work of M. Terrissol and M. C.     28 // Based on the work of M. Terrissol and M. C. Bordage
 29 //                                                 29 //
 30 // Users are requested to cite the following p     30 // Users are requested to cite the following papers:
 31 // - M. Terrissol, A. Baudre, Radiat. Prot. Do     31 // - M. Terrissol, A. Baudre, Radiat. Prot. Dosim. 31 (1990) 175-177
 32 // - M.C. Bordage, J. Bordes, S. Edel, M. Terr     32 // - M.C. Bordage, J. Bordes, S. Edel, M. Terrissol, X. Franceries,
 33 //   M. Bardies, N. Lampe, S. Incerti, Phys. M     33 //   M. Bardies, N. Lampe, S. Incerti, Phys. Med. 32 (2016) 1833-1840
 34 //                                                 34 //
 35 // Authors of this class:                          35 // Authors of this class:
 36 // M.C. Bordage, M. Terrissol, S. Edel, J. Bor     36 // M.C. Bordage, M. Terrissol, S. Edel, J. Bordes, S. Incerti
 37 //                                                 37 //
 38 // 15.01.2014: creation                            38 // 15.01.2014: creation
 39 //                                                 39 //
 40 // Based on the study by S. Zein et. al. Nucl.     40 // Based on the study by S. Zein et. al. Nucl. Inst. Meth. B 488 (2021) 70-82
 41 // 1/2/2023 : Hoang added modification             41 // 1/2/2023 : Hoang added modification
 42                                                    42 
 43 #include "G4DNACPA100IonisationModel.hh"           43 #include "G4DNACPA100IonisationModel.hh"
 44                                                    44 
 45 #include "G4DNAChemistryManager.hh"                45 #include "G4DNAChemistryManager.hh"
 46 #include "G4DNAMaterialManager.hh"                 46 #include "G4DNAMaterialManager.hh"
 47 #include "G4DNAMolecularMaterial.hh"               47 #include "G4DNAMolecularMaterial.hh"
 48 #include "G4EnvironmentUtils.hh"                   48 #include "G4EnvironmentUtils.hh"
 49 #include "G4LossTableManager.hh"                   49 #include "G4LossTableManager.hh"
 50 #include "G4PhysicalConstants.hh"                  50 #include "G4PhysicalConstants.hh"
 51 #include "G4SystemOfUnits.hh"                      51 #include "G4SystemOfUnits.hh"
 52 #include "G4UAtomicDeexcitation.hh"                52 #include "G4UAtomicDeexcitation.hh"
 53                                                    53 
 54 #include <fstream>                                 54 #include <fstream>
 55                                                    55 
 56 //....oooOO0OOooo........oooOO0OOooo........oo     56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 57                                                    57 
 58 using namespace std;                               58 using namespace std;
 59                                                    59 
 60 //....oooOO0OOooo........oooOO0OOooo........oo     60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 61                                                    61 
 62 G4DNACPA100IonisationModel::G4DNACPA100Ionisat     62 G4DNACPA100IonisationModel::G4DNACPA100IonisationModel(const G4ParticleDefinition*,
 63                                                    63                                                        const G4String& nam)
 64   : G4VDNAModel(nam, "all")                        64   : G4VDNAModel(nam, "all")
 65 {                                                  65 {
 66   fpGuanine = G4Material::GetMaterial("G4_GUAN     66   fpGuanine = G4Material::GetMaterial("G4_GUANINE", false);
 67   fpG4_WATER = G4Material::GetMaterial("G4_WAT     67   fpG4_WATER = G4Material::GetMaterial("G4_WATER", false);
 68   fpDeoxyribose = G4Material::GetMaterial("G4_     68   fpDeoxyribose = G4Material::GetMaterial("G4_DEOXYRIBOSE", false);
 69   fpCytosine = G4Material::GetMaterial("G4_CYT     69   fpCytosine = G4Material::GetMaterial("G4_CYTOSINE", false);
 70   fpThymine = G4Material::GetMaterial("G4_THYM     70   fpThymine = G4Material::GetMaterial("G4_THYMINE", false);
 71   fpAdenine = G4Material::GetMaterial("G4_ADEN     71   fpAdenine = G4Material::GetMaterial("G4_ADENINE", false);
 72   fpPhosphate = G4Material::GetMaterial("G4_PH     72   fpPhosphate = G4Material::GetMaterial("G4_PHOSPHORIC_ACID", false);
 73   fpParticle = G4Electron::ElectronDefinition(     73   fpParticle = G4Electron::ElectronDefinition();
 74 }                                                  74 }
 75                                                    75 
 76 //....oooOO0OOooo........oooOO0OOooo........oo     76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 77                                                    77 
 78 void G4DNACPA100IonisationModel::Initialise(co     78 void G4DNACPA100IonisationModel::Initialise(const G4ParticleDefinition* p,
 79                                             co     79                                             const G4DataVector& /*cuts*/)
 80 {                                                  80 {
 81   if (isInitialised) {                             81   if (isInitialised) {
 82     return;                                        82     return;
 83   }                                                83   }
 84   if (verboseLevel > 3) {                          84   if (verboseLevel > 3) {
 85     G4cout << "Calling G4DNACPA100IonisationMo     85     G4cout << "Calling G4DNACPA100IonisationModel::Initialise()" << G4endl;
 86   }                                                86   }
 87                                                    87 
 88   if (!G4DNAMaterialManager::Instance()->IsLoc     88   if (!G4DNAMaterialManager::Instance()->IsLocked()) {
 89     if (p != fpParticle) {                         89     if (p != fpParticle) {
 90       std::ostringstream oss;                      90       std::ostringstream oss;
 91       oss << " Model is not applied for this p     91       oss << " Model is not applied for this particle " << p->GetParticleName();
 92       G4Exception("G4DNACPA100IonisationModel:     92       G4Exception("G4DNACPA100IonisationModel::G4DNACPA100IonisationModel", "CPA001",
 93                   FatalException, oss.str().c_     93                   FatalException, oss.str().c_str());
 94     }                                              94     }
 95                                                    95 
 96     const char* path = G4FindDataDir("G4LEDATA     96     const char* path = G4FindDataDir("G4LEDATA");
 97                                                    97 
 98     if (path == nullptr) {                         98     if (path == nullptr) {
 99       G4Exception("G4DNACPA100IonisationModel:     99       G4Exception("G4DNACPA100IonisationModel::Initialise", "em0006", FatalException,
100                   "G4LEDATA environment variab    100                   "G4LEDATA environment variable not set.");
101       return;                                     101       return;
102     }                                             102     }
103                                                   103 
104     std::size_t index;                            104     std::size_t index;
105     if (fpG4_WATER != nullptr) {                  105     if (fpG4_WATER != nullptr) {
106       index = fpG4_WATER->GetIndex();             106       index = fpG4_WATER->GetIndex();
107       G4String eFullFileName = "";                107       G4String eFullFileName = "";
108       fasterCode ? eFullFileName = "/dna/sigma    108       fasterCode ? eFullFileName = "/dna/sigmadiff_cumulated_ionisation_e_cpa100_rel"
109                  : eFullFileName = "/dna/sigma    109                  : eFullFileName = "/dna/sigmadiff_ionisation_e_cpa100_rel";
110       AddCrossSectionData(index, p, "dna/sigma    110       AddCrossSectionData(index, p, "dna/sigma_ionisation_e_cpa100_form_rel", eFullFileName,
111                           1.e-20 * m * m);        111                           1.e-20 * m * m);
112       SetLowELimit(index, p, 11 * eV);            112       SetLowELimit(index, p, 11 * eV);
113       SetHighELimit(index, p, 255955 * eV);       113       SetHighELimit(index, p, 255955 * eV);
114     }                                             114     }
115     if (fpGuanine != nullptr) {                   115     if (fpGuanine != nullptr) {
116       index = fpGuanine->GetIndex();              116       index = fpGuanine->GetIndex();
117       G4String eFullFileName = "";                117       G4String eFullFileName = "";
118       if(useDcs) {                                118       if(useDcs) {
119         fasterCode ? eFullFileName = "/dna/sig    119         fasterCode ? eFullFileName = "/dna/sigmadiff_cumulated_elastic_e_cpa100_guanine"
120                    : eFullFileName = "/dna/sig    120                    : eFullFileName = "/dna/sigmadiff_ionisation_e_cpa100_guanine";
121       }                                           121       }
122       AddCrossSectionData(index, p, "dna/sigma    122       AddCrossSectionData(index, p, "dna/sigma_ionisation_e_cpa100_guanine", eFullFileName,
123                           1. * cm * cm);          123                           1. * cm * cm);
124       SetLowELimit(index, p, 11 * eV);            124       SetLowELimit(index, p, 11 * eV);
125       SetHighELimit(index, p, 1 * MeV);           125       SetHighELimit(index, p, 1 * MeV);
126     }                                             126     }
127     if (fpDeoxyribose != nullptr) {               127     if (fpDeoxyribose != nullptr) {
128       index = fpDeoxyribose->GetIndex();          128       index = fpDeoxyribose->GetIndex();
129       G4String eFullFileName = "";                129       G4String eFullFileName = "";
130       if(useDcs) {                                130       if(useDcs) {
131         eFullFileName = "/dna/sigmadiff_cumula    131         eFullFileName = "/dna/sigmadiff_cumulated_ionisation_e_cpa100_deoxyribose";
132       }                                           132       }
133       AddCrossSectionData(index, p, "dna/sigma    133       AddCrossSectionData(index, p, "dna/sigma_ionisation_e_cpa100_deoxyribose", eFullFileName,
134                           1. * cm * cm);          134                           1. * cm * cm);
135       SetLowELimit(index, p, 11 * eV);            135       SetLowELimit(index, p, 11 * eV);
136       SetHighELimit(index, p, 1 * MeV);           136       SetHighELimit(index, p, 1 * MeV);
137     }                                             137     }
138     if (fpCytosine != nullptr) {                  138     if (fpCytosine != nullptr) {
139       index = fpCytosine->GetIndex();             139       index = fpCytosine->GetIndex();
140       G4String eFullFileName = "";                140       G4String eFullFileName = "";
141       if(useDcs) {                                141       if(useDcs) {
142         fasterCode ? eFullFileName = "/dna/sig    142         fasterCode ? eFullFileName = "/dna/sigmadiff_cumulated_ionisation_e_cpa100_cytosine"
143                    : eFullFileName = "/dna/sig    143                    : eFullFileName = "/dna/sigmadiff_ionisation_e_cpa100_cytosine";
144       }                                           144       }
145       AddCrossSectionData(index, p, "dna/sigma    145       AddCrossSectionData(index, p, "dna/sigma_ionisation_e_cpa100_cytosine", eFullFileName,
146                           1. * cm * cm);          146                           1. * cm * cm);
147       SetLowELimit(index, p, 11 * eV);            147       SetLowELimit(index, p, 11 * eV);
148       SetHighELimit(index, p, 1 * MeV);           148       SetHighELimit(index, p, 1 * MeV);
149     }                                             149     }
150     if (fpThymine != nullptr) {                   150     if (fpThymine != nullptr) {
151       index = fpThymine->GetIndex();              151       index = fpThymine->GetIndex();
152       G4String eFullFileName = "";                152       G4String eFullFileName = "";
153       if(useDcs) {                                153       if(useDcs) {
154         fasterCode ? eFullFileName = "/dna/sig    154         fasterCode ? eFullFileName = "/dna/sigmadiff_cumulated_ionisation_e_cpa100_thymine"
155                    : eFullFileName = "/dna/sig    155                    : eFullFileName = "/dna/sigmadiff_ionisation_e_cpa100_thymine";
156       }                                           156       }
157       AddCrossSectionData(index, p, "dna/sigma    157       AddCrossSectionData(index, p, "dna/sigma_ionisation_e_cpa100_thymine", eFullFileName,
158                           1. * cm * cm);          158                           1. * cm * cm);
159       SetLowELimit(index, p, 11 * eV);            159       SetLowELimit(index, p, 11 * eV);
160       SetHighELimit(index, p, 1 * MeV);           160       SetHighELimit(index, p, 1 * MeV);
161     }                                             161     }
162     if (fpAdenine != nullptr) {                   162     if (fpAdenine != nullptr) {
163       index = fpAdenine->GetIndex();              163       index = fpAdenine->GetIndex();
164       G4String eFullFileName = "";                164       G4String eFullFileName = "";
165       if(useDcs) {                                165       if(useDcs) {
166         fasterCode ? eFullFileName = "/dna/sig    166         fasterCode ? eFullFileName = "/dna/sigmadiff_cumulated_ionisation_e_cpa100_adenine"
167                    : eFullFileName = "/dna/sig    167                    : eFullFileName = "/dna/sigmadiff_ionisation_e_cpa100_adenine";
168       }                                           168       }
169       AddCrossSectionData(index, p, "dna/sigma    169       AddCrossSectionData(index, p, "dna/sigma_ionisation_e_cpa100_adenine", eFullFileName,
170                           1. * cm * cm);          170                           1. * cm * cm);
171       SetLowELimit(index, p, 11 * eV);            171       SetLowELimit(index, p, 11 * eV);
172       SetHighELimit(index, p, 1 * MeV);           172       SetHighELimit(index, p, 1 * MeV);
173     }                                             173     }
174     if (fpPhosphate != nullptr) {                 174     if (fpPhosphate != nullptr) {
175       index = fpPhosphate->GetIndex();            175       index = fpPhosphate->GetIndex();
176       G4String eFullFileName = "";                176       G4String eFullFileName = "";
177       if(useDcs) {                                177       if(useDcs) {
178         eFullFileName = "dna/sigmadiff_cumulat    178         eFullFileName = "dna/sigmadiff_cumulated_ionisation_e_cpa100_phosphoric_acid";
179       }                                           179       }
180       AddCrossSectionData(index, p, "dna/sigma    180       AddCrossSectionData(index, p, "dna/sigma_ionisation_e_cpa100_phosphoric_acid",eFullFileName,
181                           1. * cm * cm);          181                           1. * cm * cm);
182       SetLowELimit(index, p, 11 * eV);            182       SetLowELimit(index, p, 11 * eV);
183       SetHighELimit(index, p, 1 * MeV);           183       SetHighELimit(index, p, 1 * MeV);
184     }                                             184     }
185     LoadCrossSectionData(p);                      185     LoadCrossSectionData(p);
186     G4DNAMaterialManager::Instance()->SetMaste    186     G4DNAMaterialManager::Instance()->SetMasterDataModel(DNAModelType::fDNAIonisation, this);
187     fpModelData = this;                           187     fpModelData = this;
188   }                                               188   }
189   else {                                          189   else {
190     auto dataModel = dynamic_cast<G4DNACPA100I    190     auto dataModel = dynamic_cast<G4DNACPA100IonisationModel*>(
191       G4DNAMaterialManager::Instance()->GetMod    191       G4DNAMaterialManager::Instance()->GetModel(DNAModelType::fDNAIonisation));
192     if (dataModel == nullptr) {                   192     if (dataModel == nullptr) {
193       G4cout << "G4DNACPA100IonisationModel::C    193       G4cout << "G4DNACPA100IonisationModel::CrossSectionPerVolume:: not good modelData" << G4endl;
194       throw;                                      194       throw;
195     }                                             195     }
196     fpModelData = dataModel;                      196     fpModelData = dataModel;
197   }                                               197   }
198                                                   198 
199   fAtomDeexcitation = G4LossTableManager::Inst    199   fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
200                                                   200 
201   fParticleChangeForGamma = GetParticleChangeF    201   fParticleChangeForGamma = GetParticleChangeForGamma();
202   isInitialised = true;                           202   isInitialised = true;
203 }                                                 203 }
204                                                   204 
205 G4double G4DNACPA100IonisationModel::CrossSect    205 G4double G4DNACPA100IonisationModel::CrossSectionPerVolume(const G4Material* material,
206                                                   206                                                            const G4ParticleDefinition* p,
207                                                   207                                                            G4double ekin, G4double, G4double)
208 {                                                 208 {
209   // initialise the cross section value (outpu    209   // initialise the cross section value (output value)
210   G4double sigma(0);                              210   G4double sigma(0);
211                                                   211 
212   // Get the current particle name                212   // Get the current particle name
213   const G4String& particleName = p->GetParticl    213   const G4String& particleName = p->GetParticleName();
214                                                   214 
215   if (p != fpParticle) {                          215   if (p != fpParticle) {
216     G4Exception("G4DNACPA100IonisationModel::C    216     G4Exception("G4DNACPA100IonisationModel::CrossSectionPerVolume", "em00223", FatalException,
217                 "No model is registered for th    217                 "No model is registered for this particle");
218   }                                               218   }
219                                                   219 
220   auto matID = material->GetIndex();              220   auto matID = material->GetIndex();
221                                                   221 
222   // Set the low and high energy limits           222   // Set the low and high energy limits
223   G4double lowLim = fpModelData->GetLowELimit(    223   G4double lowLim = fpModelData->GetLowELimit(matID, p);
224   G4double highLim = fpModelData->GetHighELimi    224   G4double highLim = fpModelData->GetHighELimit(matID, p);
225                                                   225 
226   // Check that we are in the correct energy r    226   // Check that we are in the correct energy range
227   if (ekin >= lowLim && ekin < highLim) {         227   if (ekin >= lowLim && ekin < highLim) {
228     // Get the map with all the model data tab    228     // Get the map with all the model data tables
229     auto tableData = fpModelData->GetData();      229     auto tableData = fpModelData->GetData();
230                                                   230 
231     if ((*tableData)[matID][p] == nullptr) {      231     if ((*tableData)[matID][p] == nullptr) {
232       G4Exception("G4DNACPA100IonisationModel:    232       G4Exception("G4DNACPA100IonisationModel::CrossSectionPerVolume", "em00236", FatalException,
233                   "No model is registered");      233                   "No model is registered");
234     }                                             234     }
235     else {                                        235     else {
236       sigma = (*tableData)[matID][p]->FindValu    236       sigma = (*tableData)[matID][p]->FindValue(ekin);
237     }                                             237     }
238                                                   238 
239     if (verboseLevel > 2) {                       239     if (verboseLevel > 2) {
240       auto MolDensity =                           240       auto MolDensity =
241         (*G4DNAMolecularMaterial::Instance()->    241         (*G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(material))[matID];
242       G4cout << "_____________________________    242       G4cout << "__________________________________" << G4endl;
243       G4cout << "°°° G4DNACPA100IonisationM    243       G4cout << "°°° G4DNACPA100IonisationModel - XS INFO START" << G4endl;
244       G4cout << "°°° Kinetic energy(eV)=" <    244       G4cout << "°°° Kinetic energy(eV)=" << ekin / eV << " particle : " << particleName << G4endl;
245       G4cout << "°°° lowLim (eV) = " << low    245       G4cout << "°°° lowLim (eV) = " << lowLim / eV << " highLim (eV) : " << highLim / eV << G4endl;
246       G4cout << "°°° Materials = " << (*G4M    246       G4cout << "°°° Materials = " << (*G4Material::GetMaterialTable())[matID]->GetName() << G4endl;
247       G4cout << "°°° Cross section per " <<    247       G4cout << "°°° Cross section per " << matID << " index molecule (cm^2)=" << sigma / cm / cm
248              << G4endl;                           248              << G4endl;
249       G4cout << "°°° Cross section per Phos    249       G4cout << "°°° Cross section per Phosphate molecule (cm^-1)="
250              << sigma * MolDensity / (1. / cm)    250              << sigma * MolDensity / (1. / cm) << G4endl;
251       G4cout << "°°° G4DNACPA100IonisationM    251       G4cout << "°°° G4DNACPA100IonisationModel - XS INFO END" << G4endl;
252     }                                             252     }
253   }                                               253   }
254                                                   254 
255   auto MolDensity = (*G4DNAMolecularMaterial::    255   auto MolDensity = (*G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(material))[matID];
256   return sigma * MolDensity;                      256   return sigma * MolDensity;
257 }                                                 257 }
258                                                   258 
259 //....oooOO0OOooo........oooOO0OOooo........oo    259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
260                                                   260 
261 void G4DNACPA100IonisationModel::SampleSeconda    261 void G4DNACPA100IonisationModel::SampleSecondaries(
262   std::vector<G4DynamicParticle*>* fvect,         262   std::vector<G4DynamicParticle*>* fvect,
263   const G4MaterialCutsCouple* couple,  // must    263   const G4MaterialCutsCouple* couple,  // must be set!
264   const G4DynamicParticle* particle, G4double,    264   const G4DynamicParticle* particle, G4double, G4double)
265 {                                                 265 {
266   if (verboseLevel > 3) {                         266   if (verboseLevel > 3) {
267     G4cout << "Calling SampleSecondaries() of     267     G4cout << "Calling SampleSecondaries() of G4DNACPA100IonisationModel" << G4endl;
268   }                                               268   }
269   auto k = particle->GetKineticEnergy();          269   auto k = particle->GetKineticEnergy();
270                                                   270 
271   const G4Material* material = couple->GetMate    271   const G4Material* material = couple->GetMaterial();
272                                                   272 
273   auto MatID = material->GetIndex();              273   auto MatID = material->GetIndex();
274                                                   274 
275   auto p = particle->GetDefinition();             275   auto p = particle->GetDefinition();
276                                                   276 
277   auto lowLim = fpModelData->GetLowELimit(MatI    277   auto lowLim = fpModelData->GetLowELimit(MatID, p);
278   auto highLim = fpModelData->GetHighELimit(Ma    278   auto highLim = fpModelData->GetHighELimit(MatID, p);
279                                                   279 
280   // Check if we are in the correct energy ran    280   // Check if we are in the correct energy range
281   if (k >= lowLim && k < highLim) {               281   if (k >= lowLim && k < highLim) {
282     const auto& primaryDirection = particle->G    282     const auto& primaryDirection = particle->GetMomentumDirection();
283     auto particleMass = particle->GetDefinitio    283     auto particleMass = particle->GetDefinition()->GetPDGMass();
284     auto totalEnergy = k + particleMass;          284     auto totalEnergy = k + particleMass;
285     auto pSquare = k * (totalEnergy + particle    285     auto pSquare = k * (totalEnergy + particleMass);
286     auto totalMomentum = std::sqrt(pSquare);      286     auto totalMomentum = std::sqrt(pSquare);
287     G4int shell = -1;                             287     G4int shell = -1;
288     G4double bindingEnergy, secondaryKinetic;     288     G4double bindingEnergy, secondaryKinetic;
289     shell = fpModelData->RandomSelectShell(k,     289     shell = fpModelData->RandomSelectShell(k, p, MatID);
290     bindingEnergy = iStructure.IonisationEnerg    290     bindingEnergy = iStructure.IonisationEnergy(shell, MatID);
291                                                   291 
292     if (k < bindingEnergy) {                      292     if (k < bindingEnergy) {
293       return;                                     293       return;
294     }                                             294     }
295                                                   295 
296     auto info = std::make_tuple(MatID, k, shel    296     auto info = std::make_tuple(MatID, k, shell);
297                                                   297 
298     secondaryKinetic = -1000 * eV;                298     secondaryKinetic = -1000 * eV;
299     if (fpG4_WATER->GetIndex() != MatID) {//fo    299     if (fpG4_WATER->GetIndex() != MatID) {//for DNA material useDcs = false
300       secondaryKinetic = fpModelData->Randomiz    300       secondaryKinetic = fpModelData->RandomizeEjectedElectronEnergyFromanalytical(info);
301     }else if(fasterCode){                         301     }else if(fasterCode){
302         secondaryKinetic = fpModelData->Random    302         secondaryKinetic = fpModelData->RandomizeEjectedElectronEnergyFromCumulatedDcs(info);
303       }else{                                      303       }else{
304         secondaryKinetic = fpModelData->Random    304         secondaryKinetic = fpModelData->RandomizeEjectedElectronEnergy(info);
305       }                                           305       }
306                                                   306 
307     G4double cosTheta = 0.;                       307     G4double cosTheta = 0.;
308     G4double phi = 0.;                            308     G4double phi = 0.;
309     RandomizeEjectedElectronDirection(particle    309     RandomizeEjectedElectronDirection(particle->GetDefinition(), k, secondaryKinetic, cosTheta,
310                                       phi);       310                                       phi);
311                                                   311 
312     G4double sinTheta = std::sqrt(1. - cosThet    312     G4double sinTheta = std::sqrt(1. - cosTheta * cosTheta);
313     G4double dirX = sinTheta * std::cos(phi);     313     G4double dirX = sinTheta * std::cos(phi);
314     G4double dirY = sinTheta * std::sin(phi);     314     G4double dirY = sinTheta * std::sin(phi);
315     G4double dirZ = cosTheta;                     315     G4double dirZ = cosTheta;
316     G4ThreeVector deltaDirection(dirX, dirY, d    316     G4ThreeVector deltaDirection(dirX, dirY, dirZ);
317     deltaDirection.rotateUz(primaryDirection);    317     deltaDirection.rotateUz(primaryDirection);
318                                                   318 
319     // SI - For atom. deexc. tagging - 23/05/2    319     // SI - For atom. deexc. tagging - 23/05/2017
320     if (secondaryKinetic > 0) {                   320     if (secondaryKinetic > 0) {
321       auto dp = new G4DynamicParticle(G4Electr    321       auto dp = new G4DynamicParticle(G4Electron::Electron(), deltaDirection, secondaryKinetic);
322       fvect->push_back(dp);                       322       fvect->push_back(dp);
323     }                                             323     }
324                                                   324 
325     if (particle->GetDefinition() != fpParticl    325     if (particle->GetDefinition() != fpParticle) {
326       fParticleChangeForGamma->ProposeMomentum    326       fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
327     }                                             327     }
328     else {                                        328     else {
329       G4double deltaTotalMomentum =               329       G4double deltaTotalMomentum =
330         std::sqrt(secondaryKinetic * (secondar    330         std::sqrt(secondaryKinetic * (secondaryKinetic + 2. * electron_mass_c2));
331       G4double finalPx =                          331       G4double finalPx =
332         totalMomentum * primaryDirection.x() -    332         totalMomentum * primaryDirection.x() - deltaTotalMomentum * deltaDirection.x();
333       G4double finalPy =                          333       G4double finalPy =
334         totalMomentum * primaryDirection.y() -    334         totalMomentum * primaryDirection.y() - deltaTotalMomentum * deltaDirection.y();
335       G4double finalPz =                          335       G4double finalPz =
336         totalMomentum * primaryDirection.z() -    336         totalMomentum * primaryDirection.z() - deltaTotalMomentum * deltaDirection.z();
337       G4double finalMomentum = std::sqrt(final    337       G4double finalMomentum = std::sqrt(finalPx * finalPx + finalPy * finalPy + finalPz * finalPz);
338       finalPx /= finalMomentum;                   338       finalPx /= finalMomentum;
339       finalPy /= finalMomentum;                   339       finalPy /= finalMomentum;
340       finalPz /= finalMomentum;                   340       finalPz /= finalMomentum;
341                                                   341 
342       G4ThreeVector direction;                    342       G4ThreeVector direction;
343       direction.set(finalPx, finalPy, finalPz)    343       direction.set(finalPx, finalPy, finalPz);
344                                                   344 
345       fParticleChangeForGamma->ProposeMomentum    345       fParticleChangeForGamma->ProposeMomentumDirection(direction.unit());
346     }                                             346     }
347                                                   347 
348     // SI - For atom. deexc. tagging - 23/05/2    348     // SI - For atom. deexc. tagging - 23/05/2017
349                                                   349 
350     // AM: sample deexcitation                    350     // AM: sample deexcitation
351     // here we assume that H_{2}O electronic l    351     // here we assume that H_{2}O electronic levels are the same of Oxigen.
352     // this can be considered true with a roug    352     // this can be considered true with a rough 10% error in energy on K-shell,
353                                                   353 
354     G4double scatteredEnergy = k - bindingEner    354     G4double scatteredEnergy = k - bindingEnergy - secondaryKinetic;
355                                                   355 
356     // SI: only atomic deexcitation from K she    356     // SI: only atomic deexcitation from K shell is considered
357     // Hoang: only for water                      357     // Hoang: only for water
358     if (material == G4Material::GetMaterial("G << 358     if (fpG4_WATER != nullptr && material == G4Material::GetMaterial("G4_WATER")) {
359       std::size_t secNumberInit = 0;  // need     359       std::size_t secNumberInit = 0;  // need to know at a certain point the energy of secondaries
360       std::size_t secNumberFinal = 0;  // So I    360       std::size_t secNumberFinal = 0;  // So I'll make the diference and then sum the energies
361       if ((fAtomDeexcitation != nullptr) && sh    361       if ((fAtomDeexcitation != nullptr) && shell == 4) {
362         G4int Z = 8;                              362         G4int Z = 8;
363         auto Kshell = fAtomDeexcitation->GetAt    363         auto Kshell = fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0));
364         secNumberInit = fvect->size();            364         secNumberInit = fvect->size();
365         fAtomDeexcitation->GenerateParticles(f    365         fAtomDeexcitation->GenerateParticles(fvect, Kshell, Z, 0, 0);
366         secNumberFinal = fvect->size();           366         secNumberFinal = fvect->size();
367         if (secNumberFinal > secNumberInit) {     367         if (secNumberFinal > secNumberInit) {
368           for (std::size_t i = secNumberInit;     368           for (std::size_t i = secNumberInit; i < secNumberFinal; ++i) {
369             // Check if there is enough residu    369             // Check if there is enough residual energy
370             if (bindingEnergy >= ((*fvect)[i])    370             if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy()) {
371               // Ok, this is a valid secondary    371               // Ok, this is a valid secondary: keep it
372               bindingEnergy -= ((*fvect)[i])->    372               bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
373             }                                     373             }
374             else {                                374             else {
375               // Invalid secondary: not enough    375               // Invalid secondary: not enough energy to create it!
376               // Keep its energy in the local     376               // Keep its energy in the local deposit
377               delete (*fvect)[i];                 377               delete (*fvect)[i];
378               (*fvect)[i] = nullptr;              378               (*fvect)[i] = nullptr;
379             }                                     379             }
380           }                                       380           }
381         }                                         381         }
382       }                                           382       }
383     }                                             383     }
384                                                   384 
385     // This should never happen                   385     // This should never happen
386     if (bindingEnergy < 0.0) {                    386     if (bindingEnergy < 0.0) {
387       G4Exception("G4DNACPA100IonisatioModel1:    387       G4Exception("G4DNACPA100IonisatioModel1::SampleSecondaries()", "em2050", FatalException,
388                   "Negative local energy depos    388                   "Negative local energy deposit");
389     }                                             389     }
390     if (!statCode) {                              390     if (!statCode) {
391       fParticleChangeForGamma->SetProposedKine    391       fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
392       fParticleChangeForGamma->ProposeLocalEne    392       fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
393     }                                             393     }
394     else {                                        394     else {
395       fParticleChangeForGamma->SetProposedKine    395       fParticleChangeForGamma->SetProposedKineticEnergy(k);
396       fParticleChangeForGamma->ProposeLocalEne    396       fParticleChangeForGamma->ProposeLocalEnergyDeposit(k - scatteredEnergy);
397     }                                             397     }
398                                                   398 
399     // only water for chemistry                   399     // only water for chemistry
400     if (fpG4_WATER != nullptr && material == G    400     if (fpG4_WATER != nullptr && material == G4Material::GetMaterial("G4_WATER")) {
401       const G4Track* theIncomingTrack = fParti    401       const G4Track* theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
402       G4DNAChemistryManager::Instance()->Creat    402       G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule, shell,
403                                                   403                                                              theIncomingTrack);
404     }                                             404     }
405   }                                               405   }
406 }                                                 406 }
407                                                   407 
408 //....oooOO0OOooo........oooOO0OOooo........oo    408 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
409                                                   409 
410 G4double G4DNACPA100IonisationModel::Randomize    410 G4double G4DNACPA100IonisationModel::RandomizeEjectedElectronEnergy(PartKineticInMat info)
411 {                                                 411 {
412   auto MatID = std::get<0>(info);                 412   auto MatID = std::get<0>(info);
413   auto k = std::get<1>(info);                     413   auto k = std::get<1>(info);
414   auto shell = std::get<2>(info);                 414   auto shell = std::get<2>(info);
415   G4double maximumEnergyTransfer = 0.;            415   G4double maximumEnergyTransfer = 0.;
416   auto IonLevel = iStructure.IonisationEnergy(    416   auto IonLevel = iStructure.IonisationEnergy(shell, MatID);
417   (k + IonLevel) / 2. > k ? maximumEnergyTrans    417   (k + IonLevel) / 2. > k ? maximumEnergyTransfer = k : maximumEnergyTransfer = (k + IonLevel) / 2.;
418                                                   418 
419   G4double crossSectionMaximum = 0.;              419   G4double crossSectionMaximum = 0.;
420                                                   420 
421   G4double minEnergy = IonLevel;                  421   G4double minEnergy = IonLevel;
422   G4double maxEnergy = maximumEnergyTransfer;     422   G4double maxEnergy = maximumEnergyTransfer;
423                                                   423 
424   // nEnergySteps can be optimized - 100 by de    424   // nEnergySteps can be optimized - 100 by default
425   G4int nEnergySteps = 50;                        425   G4int nEnergySteps = 50;
426                                                   426 
427   G4double value(minEnergy);                      427   G4double value(minEnergy);
428   G4double stpEnergy(std::pow(maxEnergy / valu    428   G4double stpEnergy(std::pow(maxEnergy / value, 1. / static_cast<G4double>(nEnergySteps - 1)));
429   G4int step(nEnergySteps);                       429   G4int step(nEnergySteps);
430   G4double differentialCrossSection = 0.;         430   G4double differentialCrossSection = 0.;
431   while (step > 0) {                              431   while (step > 0) {
432     step--;                                       432     step--;
433     differentialCrossSection = DifferentialCro    433     differentialCrossSection = DifferentialCrossSection(info, value / eV);
434                                                   434 
435     if (differentialCrossSection > 0) {           435     if (differentialCrossSection > 0) {
436       crossSectionMaximum = differentialCrossS    436       crossSectionMaximum = differentialCrossSection;
437       break;                                      437       break;
438     }                                             438     }
439     value *= stpEnergy;                           439     value *= stpEnergy;
440   }                                               440   }
441                                                   441 
442   G4double secondaryElectronKineticEnergy = 0.    442   G4double secondaryElectronKineticEnergy = 0.;
443   do {                                            443   do {
444     secondaryElectronKineticEnergy = G4Uniform    444     secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer - IonLevel);
445   } while (G4UniformRand() * crossSectionMaxim    445   } while (G4UniformRand() * crossSectionMaximum
446            > DifferentialCrossSection(info, (s    446            > DifferentialCrossSection(info, (secondaryElectronKineticEnergy + IonLevel) / eV));
447                                                   447 
448   return secondaryElectronKineticEnergy;          448   return secondaryElectronKineticEnergy;
449 }                                                 449 }
450                                                   450 
451 //....oooOO0OOooo........oooOO0OOooo........oo    451 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
452                                                   452 
453 void G4DNACPA100IonisationModel::RandomizeEjec    453 void G4DNACPA100IonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition*,
454                                                   454                                                                    G4double k, G4double secKinetic,
455                                                   455                                                                    G4double& cosTheta,
456                                                   456                                                                    G4double& phi)
457 {                                                 457 {
458   phi = twopi * G4UniformRand();                  458   phi = twopi * G4UniformRand();
459   G4double sin2O = (1. - secKinetic / k) / (1.    459   G4double sin2O = (1. - secKinetic / k) / (1. + secKinetic / (2. * electron_mass_c2));
460   cosTheta = std::sqrt(1. - sin2O);               460   cosTheta = std::sqrt(1. - sin2O);
461 }                                                 461 }
462                                                   462 
463 //....oooOO0OOooo........oooOO0OOooo........oo    463 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
464                                                   464 
465 G4double G4DNACPA100IonisationModel::Different    465 G4double G4DNACPA100IonisationModel::DifferentialCrossSection(PartKineticInMat info,
466                                                   466                                                               const G4double& energyTransfer)
467 {                                                 467 {
468   auto MatID = std::get<0>(info);                 468   auto MatID = std::get<0>(info);
469   auto k = std::get<1>(info) / eV;  // in eV u    469   auto k = std::get<1>(info) / eV;  // in eV unit
470   auto shell = std::get<2>(info);                 470   auto shell = std::get<2>(info);
471   G4double sigma = 0.;                            471   G4double sigma = 0.;
472   G4double shellEnergy = iStructure.Ionisation    472   G4double shellEnergy = iStructure.IonisationEnergy(shell, MatID);
473   G4double kSE(energyTransfer - shellEnergy);     473   G4double kSE(energyTransfer - shellEnergy);
474                                                   474 
475   if (energyTransfer >= shellEnergy) {            475   if (energyTransfer >= shellEnergy) {
476     G4double valueT1 = 0;                         476     G4double valueT1 = 0;
477     G4double valueT2 = 0;                         477     G4double valueT2 = 0;
478     G4double valueE21 = 0;                        478     G4double valueE21 = 0;
479     G4double valueE22 = 0;                        479     G4double valueE22 = 0;
480     G4double valueE12 = 0;                        480     G4double valueE12 = 0;
481     G4double valueE11 = 0;                        481     G4double valueE11 = 0;
482                                                   482 
483     G4double xs11 = 0;                            483     G4double xs11 = 0;
484     G4double xs12 = 0;                            484     G4double xs12 = 0;
485     G4double xs21 = 0;                            485     G4double xs21 = 0;
486     G4double xs22 = 0;                            486     G4double xs22 = 0;
487                                                   487 
488     auto t2 = std::upper_bound(fTMapWithVec[Ma    488     auto t2 = std::upper_bound(fTMapWithVec[MatID][fpParticle].begin(),
489                                fTMapWithVec[Ma    489                                fTMapWithVec[MatID][fpParticle].end(), k);
490     auto t1 = t2 - 1;                             490     auto t1 = t2 - 1;
491                                                   491 
492     if (kSE <= fEMapWithVector[MatID][fpPartic    492     if (kSE <= fEMapWithVector[MatID][fpParticle][(*t1)].back()
493         && kSE <= fEMapWithVector[MatID][fpPar    493         && kSE <= fEMapWithVector[MatID][fpParticle][(*t2)].back())
494     {                                             494     {
495       auto e12 = std::upper_bound(fEMapWithVec    495       auto e12 = std::upper_bound(fEMapWithVector[MatID][fpParticle][(*t1)].begin(),
496                                   fEMapWithVec    496                                   fEMapWithVector[MatID][fpParticle][(*t1)].end(), kSE);
497       auto e11 = e12 - 1;                         497       auto e11 = e12 - 1;
498                                                   498 
499       auto e22 = std::upper_bound(fEMapWithVec    499       auto e22 = std::upper_bound(fEMapWithVector[MatID][fpParticle][(*t2)].begin(),
500                                   fEMapWithVec    500                                   fEMapWithVector[MatID][fpParticle][(*t2)].end(), kSE);
501       auto e21 = e22 - 1;                         501       auto e21 = e22 - 1;
502                                                   502 
503       valueT1 = *t1;                              503       valueT1 = *t1;
504       valueT2 = *t2;                              504       valueT2 = *t2;
505       valueE21 = *e21;                            505       valueE21 = *e21;
506       valueE22 = *e22;                            506       valueE22 = *e22;
507       valueE12 = *e12;                            507       valueE12 = *e12;
508       valueE11 = *e11;                            508       valueE11 = *e11;
509                                                   509 
510       xs11 = diffCrossSectionData[MatID][fpPar    510       xs11 = diffCrossSectionData[MatID][fpParticle][shell][valueT1][valueE11];
511       xs12 = diffCrossSectionData[MatID][fpPar    511       xs12 = diffCrossSectionData[MatID][fpParticle][shell][valueT1][valueE12];
512       xs21 = diffCrossSectionData[MatID][fpPar    512       xs21 = diffCrossSectionData[MatID][fpParticle][shell][valueT2][valueE21];
513       xs22 = diffCrossSectionData[MatID][fpPar    513       xs22 = diffCrossSectionData[MatID][fpParticle][shell][valueT2][valueE22];
514     }                                             514     }
515                                                   515 
516     G4double xsProduct = xs11 * xs12 * xs21 *     516     G4double xsProduct = xs11 * xs12 * xs21 * xs22;
517                                                   517 
518     if (xsProduct != 0.) {                        518     if (xsProduct != 0.) {
519       sigma = QuadInterpolator(valueE11, value    519       sigma = QuadInterpolator(valueE11, valueE12, valueE21, valueE22, xs11, xs12, xs21, xs22,
520                                valueT1, valueT    520                                valueT1, valueT2, k, kSE);
521     }                                             521     }
522   }                                               522   }
523                                                   523 
524   return sigma;                                   524   return sigma;
525 }                                                 525 }
526                                                   526 
527 //....oooOO0OOooo........oooOO0OOooo........oo    527 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
528                                                   528 
529 G4double G4DNACPA100IonisationModel::Interpola    529 G4double G4DNACPA100IonisationModel::Interpolate(G4double e1, G4double e2, G4double e, G4double xs1,
530                                                   530                                                  G4double xs2)
531 {                                                 531 {
532   G4double value = 0.;                            532   G4double value = 0.;
533                                                   533 
534   // Log-log interpolation by default             534   // Log-log interpolation by default
535                                                   535 
536   if (e1 != 0 && e2 != 0 && (std::log10(e2) -     536   if (e1 != 0 && e2 != 0 && (std::log10(e2) - std::log10(e1)) != 0 && !fasterCode) {
537     G4double a = (std::log10(xs2) - std::log10    537     G4double a = (std::log10(xs2) - std::log10(xs1)) / (std::log10(e2) - std::log10(e1));
538     G4double b = std::log10(xs2) - a * std::lo    538     G4double b = std::log10(xs2) - a * std::log10(e2);
539     G4double sigma = a * std::log10(e) + b;       539     G4double sigma = a * std::log10(e) + b;
540     value = (std::pow(10., sigma));               540     value = (std::pow(10., sigma));
541   }                                               541   }
542                                                   542 
543   // Switch to lin-lin interpolation              543   // Switch to lin-lin interpolation
544   /*                                              544   /*
545   if ((e2-e1)!=0)                                 545   if ((e2-e1)!=0)
546   {                                               546   {
547     G4double d1 = xs1;                            547     G4double d1 = xs1;
548     G4double d2 = xs2;                            548     G4double d2 = xs2;
549     value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1    549     value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
550   }                                               550   }
551   */                                              551   */
552                                                   552 
553   // Switch to log-lin interpolation for faste    553   // Switch to log-lin interpolation for faster code
554                                                   554 
555   if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 &    555   if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 && fasterCode) {
556     G4double d1 = std::log10(xs1);                556     G4double d1 = std::log10(xs1);
557     G4double d2 = std::log10(xs2);                557     G4double d2 = std::log10(xs2);
558     value = std::pow(10., (d1 + (d2 - d1) * (e    558     value = std::pow(10., (d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
559   }                                               559   }
560                                                   560 
561   // Switch to lin-lin interpolation for faste    561   // Switch to lin-lin interpolation for faster code
562   // in case one of xs1 or xs2 (=cum proba) va    562   // in case one of xs1 or xs2 (=cum proba) value is zero
563                                                   563 
564   if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0)    564   if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) && fasterCode) {
565     G4double d1 = xs1;                            565     G4double d1 = xs1;
566     G4double d2 = xs2;                            566     G4double d2 = xs2;
567     value = (d1 + (d2 - d1) * (e - e1) / (e2 -    567     value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
568   }                                               568   }
569   return value;                                   569   return value;
570 }                                                 570 }
571                                                   571 
572 //....oooOO0OOooo........oooOO0OOooo........oo    572 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
573                                                   573 
574 G4double G4DNACPA100IonisationModel::QuadInter    574 G4double G4DNACPA100IonisationModel::QuadInterpolator(G4double e11, G4double e12, G4double e21,
575                                                   575                                                       G4double e22, G4double xs11, G4double xs12,
576                                                   576                                                       G4double xs21, G4double xs22, G4double t1,
577                                                   577                                                       G4double t2, G4double t, G4double e)
578 {                                                 578 {
579   G4double interpolatedvalue1 = Interpolate(e1    579   G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
580   G4double interpolatedvalue2 = Interpolate(e2    580   G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
581   G4double value = Interpolate(t1, t2, t, inte    581   G4double value = Interpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
582                                                   582 
583   return value;                                   583   return value;
584 }                                                 584 }
585                                                   585 
586 G4double                                          586 G4double
587 G4DNACPA100IonisationModel::RandomizeEjectedEl    587 G4DNACPA100IonisationModel::RandomizeEjectedElectronEnergyFromCumulatedDcs(PartKineticInMat info)
588 {                                                 588 {
589   auto MatID = std::get<0>(info);                 589   auto MatID = std::get<0>(info);
590   auto shell = std::get<2>(info);                 590   auto shell = std::get<2>(info);
591   G4double secondaryElectronKineticEnergy =       591   G4double secondaryElectronKineticEnergy =
592     RandomTransferedEnergy(info) * eV - iStruc    592     RandomTransferedEnergy(info) * eV - iStructure.IonisationEnergy(shell, MatID);
593   if (secondaryElectronKineticEnergy < 0.) {      593   if (secondaryElectronKineticEnergy < 0.) {
594     return 0.;                                    594     return 0.;
595   }                                               595   }
596   return secondaryElectronKineticEnergy;          596   return secondaryElectronKineticEnergy;
597 }                                                 597 }
598                                                   598 
599 G4double G4DNACPA100IonisationModel::RandomTra    599 G4double G4DNACPA100IonisationModel::RandomTransferedEnergy(PartKineticInMat info)
600 {                                                 600 {
601   auto materialID = std::get<0>(info);            601   auto materialID = std::get<0>(info);
602   auto k = std::get<1>(info) / eV;  // data ta    602   auto k = std::get<1>(info) / eV;  // data table in eV
603   auto shell = std::get<2>(info);                 603   auto shell = std::get<2>(info);
604   G4double ejectedElectronEnergy = 0.;            604   G4double ejectedElectronEnergy = 0.;
605   G4double valueK1 = 0;                           605   G4double valueK1 = 0;
606   G4double valueK2 = 0;                           606   G4double valueK2 = 0;
607   G4double valueCumulCS21 = 0;                    607   G4double valueCumulCS21 = 0;
608   G4double valueCumulCS22 = 0;                    608   G4double valueCumulCS22 = 0;
609   G4double valueCumulCS12 = 0;                    609   G4double valueCumulCS12 = 0;
610   G4double valueCumulCS11 = 0;                    610   G4double valueCumulCS11 = 0;
611   G4double secElecE11 = 0;                        611   G4double secElecE11 = 0;
612   G4double secElecE12 = 0;                        612   G4double secElecE12 = 0;
613   G4double secElecE21 = 0;                        613   G4double secElecE21 = 0;
614   G4double secElecE22 = 0;                        614   G4double secElecE22 = 0;
615                                                   615 
616   if (k == fTMapWithVec[materialID][fpParticle    616   if (k == fTMapWithVec[materialID][fpParticle].back()) {
617     k = k * (1. - 1e-12);                         617     k = k * (1. - 1e-12);
618   }                                               618   }
619                                                   619 
620   G4double random = G4UniformRand();              620   G4double random = G4UniformRand();
621   auto k2 = std::upper_bound(fTMapWithVec[mate    621   auto k2 = std::upper_bound(fTMapWithVec[materialID][fpParticle].begin(),
622                              fTMapWithVec[mate    622                              fTMapWithVec[materialID][fpParticle].end(), k);
623   auto k1 = k2 - 1;                               623   auto k1 = k2 - 1;
624                                                   624 
625   if (random <= fProbaShellMap[materialID][fpP    625   if (random <= fProbaShellMap[materialID][fpParticle][shell][(*k1)].back()
626       && random <= fProbaShellMap[materialID][    626       && random <= fProbaShellMap[materialID][fpParticle][shell][(*k2)].back())
627   {                                               627   {
628     auto cumulCS12 =                              628     auto cumulCS12 =
629       std::upper_bound(fProbaShellMap[material    629       std::upper_bound(fProbaShellMap[materialID][fpParticle][shell][(*k1)].begin(),
630                        fProbaShellMap[material    630                        fProbaShellMap[materialID][fpParticle][shell][(*k1)].end(), random);
631     auto cumulCS11 = cumulCS12 - 1;               631     auto cumulCS11 = cumulCS12 - 1;
632     // Second one.                                632     // Second one.
633     auto cumulCS22 =                              633     auto cumulCS22 =
634       std::upper_bound(fProbaShellMap[material    634       std::upper_bound(fProbaShellMap[materialID][fpParticle][shell][(*k2)].begin(),
635                        fProbaShellMap[material    635                        fProbaShellMap[materialID][fpParticle][shell][(*k2)].end(), random);
636     auto cumulCS21 = cumulCS22 - 1;               636     auto cumulCS21 = cumulCS22 - 1;
637                                                   637 
638     valueK1 = *k1;                                638     valueK1 = *k1;
639     valueK2 = *k2;                                639     valueK2 = *k2;
640     valueCumulCS11 = *cumulCS11;                  640     valueCumulCS11 = *cumulCS11;
641     valueCumulCS12 = *cumulCS12;                  641     valueCumulCS12 = *cumulCS12;
642     valueCumulCS21 = *cumulCS21;                  642     valueCumulCS21 = *cumulCS21;
643     valueCumulCS22 = *cumulCS22;                  643     valueCumulCS22 = *cumulCS22;
644                                                   644 
645     secElecE11 = fEnergySecondaryData[material    645     secElecE11 = fEnergySecondaryData[materialID][fpParticle][shell][valueK1][valueCumulCS11];
646     secElecE12 = fEnergySecondaryData[material    646     secElecE12 = fEnergySecondaryData[materialID][fpParticle][shell][valueK1][valueCumulCS12];
647     secElecE21 = fEnergySecondaryData[material    647     secElecE21 = fEnergySecondaryData[materialID][fpParticle][shell][valueK2][valueCumulCS21];
648     secElecE22 = fEnergySecondaryData[material    648     secElecE22 = fEnergySecondaryData[materialID][fpParticle][shell][valueK2][valueCumulCS22];
649                                                   649 
650     if (valueCumulCS11 == 0. && valueCumulCS12    650     if (valueCumulCS11 == 0. && valueCumulCS12 == 1.) {
651       auto interpolatedvalue2 =                   651       auto interpolatedvalue2 =
652         Interpolate(valueCumulCS21, valueCumul    652         Interpolate(valueCumulCS21, valueCumulCS22, random, secElecE21, secElecE22);
653       G4double valueNrjTransf = Interpolate(va    653       G4double valueNrjTransf = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
654       return valueNrjTransf;                      654       return valueNrjTransf;
655     }                                             655     }
656   }                                               656   }
657                                                   657 
658   if (random > fProbaShellMap[materialID][fpPa    658   if (random > fProbaShellMap[materialID][fpParticle][shell][(*k1)].back()) {
659     auto cumulCS22 =                              659     auto cumulCS22 =
660       std::upper_bound(fProbaShellMap[material    660       std::upper_bound(fProbaShellMap[materialID][fpParticle][shell][(*k2)].begin(),
661                        fProbaShellMap[material    661                        fProbaShellMap[materialID][fpParticle][shell][(*k2)].end(), random);
662     auto cumulCS21 = cumulCS22 - 1;               662     auto cumulCS21 = cumulCS22 - 1;
663     valueK1 = *k1;                                663     valueK1 = *k1;
664     valueK2 = *k2;                                664     valueK2 = *k2;
665     valueCumulCS21 = *cumulCS21;                  665     valueCumulCS21 = *cumulCS21;
666     valueCumulCS22 = *cumulCS22;                  666     valueCumulCS22 = *cumulCS22;
667                                                   667 
668     secElecE21 = fEnergySecondaryData[material    668     secElecE21 = fEnergySecondaryData[materialID][fpParticle][shell][valueK2][valueCumulCS21];
669     secElecE22 = fEnergySecondaryData[material    669     secElecE22 = fEnergySecondaryData[materialID][fpParticle][shell][valueK2][valueCumulCS22];
670                                                   670 
671     G4double interpolatedvalue2 =                 671     G4double interpolatedvalue2 =
672       Interpolate(valueCumulCS21, valueCumulCS    672       Interpolate(valueCumulCS21, valueCumulCS22, random, secElecE21, secElecE22);
673                                                   673 
674     G4double value = Interpolate(valueK1, valu    674     G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
675     return value;                                 675     return value;
676   }                                               676   }
677   G4double nrjTransfProduct = secElecE11 * sec    677   G4double nrjTransfProduct = secElecE11 * secElecE12 * secElecE21 * secElecE22;
678                                                   678 
679   if (nrjTransfProduct != 0.) {                   679   if (nrjTransfProduct != 0.) {
680     ejectedElectronEnergy =                       680     ejectedElectronEnergy =
681       QuadInterpolator(valueCumulCS11, valueCu    681       QuadInterpolator(valueCumulCS11, valueCumulCS12, valueCumulCS21, valueCumulCS22, secElecE11,
682                        secElecE12, secElecE21,    682                        secElecE12, secElecE21, secElecE22, valueK1, valueK2, k, random);
683   }                                               683   }
684   return ejectedElectronEnergy;                   684   return ejectedElectronEnergy;
685 }                                                 685 }
686                                                   686 
687 //....oooOO0OOooo........oooOO0OOooo........oo    687 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
688                                                   688 
689 G4double                                          689 G4double
690 G4DNACPA100IonisationModel::RandomizeEjectedEl    690 G4DNACPA100IonisationModel::RandomizeEjectedElectronEnergyFromanalytical(PartKineticInMat info)
691 {                                                 691 {
692   auto MatID = std::get<0>(info);                 692   auto MatID = std::get<0>(info);
693   auto tt = std::get<1>(info);                    693   auto tt = std::get<1>(info);
694   auto shell = std::get<2>(info);                 694   auto shell = std::get<2>(info);
695   // ***** METHOD by M. C. Bordage ***** (opti    695   // ***** METHOD by M. C. Bordage ***** (optimized)
696   //  Composition sampling method based on eq     696   //  Composition sampling method based on eq 7 in (Guerra et al. 2015) the RBEBV
697                                                   697 
698   //// Defining constants                         698   //// Defining constants
699   G4double alfa = 1. / 137;  // fine structure    699   G4double alfa = 1. / 137;  // fine structure constant
700   G4double e_charge = 1.6e-19;  // electron ch    700   G4double e_charge = 1.6e-19;  // electron charge
701   G4double e_mass = 9.1e-31;  // electron mass    701   G4double e_mass = 9.1e-31;  // electron mass in kg
702   G4double c = 3e8;  // speed of light in vacu    702   G4double c = 3e8;  // speed of light in vacuum constant c (m/s)
703   G4double mc2 = e_mass * c * c / e_charge;  /    703   G4double mc2 = e_mass * c * c / e_charge;  //
704                                                   704 
705   G4double BB = iStructure.IonisationEnergy(sh    705   G4double BB = iStructure.IonisationEnergy(shell, MatID);  // binding energy of the shell (eV)
706                                                   706 
707   if (tt <= BB) return 0.;                        707   if (tt <= BB) return 0.;
708                                                   708 
709   G4double b_prime = BB / mc2;  // binding ene    709   G4double b_prime = BB / mc2;  // binding energy divided by mc2
710   G4double beta_b2 = 1. - 1. / ((1 + b_prime)     710   G4double beta_b2 = 1. - 1. / ((1 + b_prime) * (1 + b_prime));  // binding energy Beta
711                                                   711 
712   //// Indicent energy                            712   //// Indicent energy
713   //// tt is the incident electron energy         713   //// tt is the incident electron energy
714                                                   714 
715   G4double t_prime = tt / mc2;  // incident en    715   G4double t_prime = tt / mc2;  // incident energy divided by mc2
716   G4double t = tt / BB;  // reduced incident e    716   G4double t = tt / BB;  // reduced incident energy by binding energy
717                                                   717 
718   G4double D = (1 + 2 * t_prime) / ((1 + t_pri    718   G4double D = (1 + 2 * t_prime) / ((1 + t_prime / 2) * (1 + t_prime / 2));
719   G4double F = b_prime * b_prime / ((1 + t_pri    719   G4double F = b_prime * b_prime / ((1 + t_prime / 2) * (1 + t_prime / 2));
720                                                   720 
721   G4double beta_t2 = 1 - 1 / ((1 + t_prime) *     721   G4double beta_t2 = 1 - 1 / ((1 + t_prime) * (1 + t_prime));  // incident energy Beta
722                                                   722 
723   G4double PHI_R = std::cos(std::sqrt(alfa * a    723   G4double PHI_R = std::cos(std::sqrt(alfa * alfa / (beta_t2 + beta_b2))
724                             * std::log(beta_t2    724                             * std::log(beta_t2 / beta_b2));  // relativistic Vriens function phi
725   G4double G_R = std::log(beta_t2 / (1 - beta_    725   G4double G_R = std::log(beta_t2 / (1 - beta_t2)) - beta_t2 - std::log(2 * b_prime);
726                                                   726 
727   G4double tplus1 = t + 1;                        727   G4double tplus1 = t + 1;
728   G4double tminus1 = t - 1;                       728   G4double tminus1 = t - 1;
729   G4double tplus12 = tplus1 * tplus1;             729   G4double tplus12 = tplus1 * tplus1;
730   G4double ZH1max = 1 + F - (PHI_R * D * (2 *     730   G4double ZH1max = 1 + F - (PHI_R * D * (2 * t + 1) / (2 * t * tplus1));
731   G4double ZH2max = 1 - PHI_R * D / 4;            731   G4double ZH2max = 1 - PHI_R * D / 4;
732                                                   732 
733   G4double A1_p = ZH1max * tminus1 / tplus1;      733   G4double A1_p = ZH1max * tminus1 / tplus1;  // A1'
734   G4double A2_p = ZH2max * tminus1 / (t * tplu    734   G4double A2_p = ZH2max * tminus1 / (t * tplus1);  // A2'
735   G4double A3_p = ((tplus12 - 4) / tplus12) *     735   G4double A3_p = ((tplus12 - 4) / tplus12) * G_R;  // A3'
736                                                   736 
737   G4double AAA = A1_p + A2_p + A3_p;              737   G4double AAA = A1_p + A2_p + A3_p;
738                                                   738 
739   G4double AA1_R = A1_p / AAA;                    739   G4double AA1_R = A1_p / AAA;
740   G4double AA2_R = (A1_p + A2_p) / AAA;           740   G4double AA2_R = (A1_p + A2_p) / AAA;
741                                                   741 
742   G4int FF = 0;                                   742   G4int FF = 0;
743   G4double fx = 0;                                743   G4double fx = 0;
744   G4double gx = 0;                                744   G4double gx = 0;
745   G4double gg = 0;                                745   G4double gg = 0;
746   G4double wx = 0;                                746   G4double wx = 0;
747                                                   747 
748   G4double r1 = 0;                                748   G4double r1 = 0;
749   G4double r2 = 0;                                749   G4double r2 = 0;
750   G4double r3 = 0;                                750   G4double r3 = 0;
751                                                   751 
752   //                                              752   //
753                                                   753 
754   do {                                            754   do {
755     r1 = G4UniformRand();                         755     r1 = G4UniformRand();
756     r2 = G4UniformRand();                         756     r2 = G4UniformRand();
757     r3 = G4UniformRand();                         757     r3 = G4UniformRand();
758                                                   758 
759     if (r1 > AA2_R)                               759     if (r1 > AA2_R)
760       FF = 3;                                     760       FF = 3;
761     else if ((r1 > AA1_R) && (r1 < AA2_R))        761     else if ((r1 > AA1_R) && (r1 < AA2_R))
762       FF = 2;                                     762       FF = 2;
763     else                                          763     else
764       FF = 1;                                     764       FF = 1;
765                                                   765 
766     switch (FF) {                                 766     switch (FF) {
767       case 1: {                                   767       case 1: {
768         fx = r2 * tminus1 / tplus1;               768         fx = r2 * tminus1 / tplus1;
769         wx = 1 / (1 - fx) - 1;                    769         wx = 1 / (1 - fx) - 1;
770         gg = PHI_R * D * (wx + 1) / tplus1;       770         gg = PHI_R * D * (wx + 1) / tplus1;
771         gx = 1 - gg;                              771         gx = 1 - gg;
772         gx = gx - gg * (wx + 1) / (2 * (t - wx    772         gx = gx - gg * (wx + 1) / (2 * (t - wx));
773         gx = gx + F * (wx + 1) * (wx + 1);        773         gx = gx + F * (wx + 1) * (wx + 1);
774         gx = gx / ZH1max;                         774         gx = gx / ZH1max;
775         break;                                    775         break;
776       }                                           776       }
777                                                   777 
778       case 2: {                                   778       case 2: {
779         fx = tplus1 + r2 * tminus1;               779         fx = tplus1 + r2 * tminus1;
780         wx = t * tminus1 * r2 / fx;               780         wx = t * tminus1 * r2 / fx;
781         gx = 1 - (PHI_R * D * (t - wx) / (2 *     781         gx = 1 - (PHI_R * D * (t - wx) / (2 * tplus1));
782         gx = gx / ZH2max;                         782         gx = gx / ZH2max;
783         break;                                    783         break;
784       }                                           784       }
785                                                   785 
786       case 3: {                                   786       case 3: {
787         fx = 1 - r2 * (tplus12 - 4) / tplus12;    787         fx = 1 - r2 * (tplus12 - 4) / tplus12;
788         wx = std::sqrt(1 / fx) - 1;               788         wx = std::sqrt(1 / fx) - 1;
789         gg = (wx + 1) / (t - wx);                 789         gg = (wx + 1) / (t - wx);
790         gx = (1 + gg * gg * gg) / 2;              790         gx = (1 + gg * gg * gg) / 2;
791         break;                                    791         break;
792       }                                           792       }
793     }  // switch                                  793     }  // switch
794                                                   794 
795   } while (r3 > gx);                              795   } while (r3 > gx);
796                                                   796 
797   return wx * BB;                                 797   return wx * BB;
798 }                                                 798 }
799                                                   799 
800 //....oooOO0OOooo........oooOO0OOooo........oo    800 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
801                                                   801 
802 void G4DNACPA100IonisationModel::ReadDiffCSFil    802 void G4DNACPA100IonisationModel::ReadDiffCSFile(const std::size_t& materialID,
803                                                   803                                                 const G4ParticleDefinition* p, const G4String& file,
804                                                   804                                                 const G4double& scaleFactor)
805 {                                                 805 {
806   const char* path = G4FindDataDir("G4LEDATA")    806   const char* path = G4FindDataDir("G4LEDATA");
807   if (path == nullptr) {                          807   if (path == nullptr) {
808     G4Exception("G4DNACPA100IonisationModel::R    808     G4Exception("G4DNACPA100IonisationModel::ReadAllDiffCSFiles", "em0006", FatalException,
809                 "G4LEDATA environment variable    809                 "G4LEDATA environment variable was not set.");
810     return;                                       810     return;
811   }                                               811   }
812                                                   812 
813   std::ostringstream fullFileName;                813   std::ostringstream fullFileName;
814   fullFileName << path << "/" << file << ".dat    814   fullFileName << path << "/" << file << ".dat";
815                                                   815 
816   std::ifstream diffCrossSection(fullFileName.    816   std::ifstream diffCrossSection(fullFileName.str().c_str());
817   std::stringstream endPath;                      817   std::stringstream endPath;
818   if (!diffCrossSection) {                        818   if (!diffCrossSection) {
819     endPath << "Missing data file: " << file;     819     endPath << "Missing data file: " << file;
820     G4Exception("G4DNACPA100IonisationModel::I    820     G4Exception("G4DNACPA100IonisationModel::Initialise", "em0003", FatalException,
821                 endPath.str().c_str());           821                 endPath.str().c_str());
822   }                                               822   }
823                                                   823 
824   // load data from the file                      824   // load data from the file
825   fTMapWithVec[materialID][p].push_back(0.);      825   fTMapWithVec[materialID][p].push_back(0.);
826                                                   826 
827   G4String line;                                  827   G4String line;
828                                                   828 
829   while (!diffCrossSection.eof()) {               829   while (!diffCrossSection.eof()) {
830     G4double T, E;                                830     G4double T, E;
831     diffCrossSection >> T >> E;                   831     diffCrossSection >> T >> E;
832                                                   832 
833     if (T != fTMapWithVec[materialID][p].back(    833     if (T != fTMapWithVec[materialID][p].back()) {
834       fTMapWithVec[materialID][p].push_back(T)    834       fTMapWithVec[materialID][p].push_back(T);
835     }                                             835     }
836                                                   836 
837     // T is incident energy, E is the energy t    837     // T is incident energy, E is the energy transferred
838     if (T != fTMapWithVec[materialID][p].back(    838     if (T != fTMapWithVec[materialID][p].back()) {
839       fTMapWithVec[materialID][p].push_back(T)    839       fTMapWithVec[materialID][p].push_back(T);
840     }                                             840     }
841                                                   841 
842     auto eshell = (G4int)iStructure.NumberOfLe    842     auto eshell = (G4int)iStructure.NumberOfLevels(materialID);
843     for (G4int shell = 0; shell < eshell; ++sh    843     for (G4int shell = 0; shell < eshell; ++shell) {
844       diffCrossSection >> diffCrossSectionData    844       diffCrossSection >> diffCrossSectionData[materialID][p][shell][T][E];
845       if (fasterCode) {                           845       if (fasterCode) {
846         fEnergySecondaryData[materialID][p][sh    846         fEnergySecondaryData[materialID][p][shell][T]
847                             [diffCrossSectionD    847                             [diffCrossSectionData[materialID][p][shell][T][E]] = E;
848                                                   848 
849         fProbaShellMap[materialID][p][shell][T    849         fProbaShellMap[materialID][p][shell][T].push_back(
850           diffCrossSectionData[materialID][p][    850           diffCrossSectionData[materialID][p][shell][T][E]);
851       }                                           851       }
852       else {                                      852       else {
853         diffCrossSectionData[materialID][p][sh    853         diffCrossSectionData[materialID][p][shell][T][E] *= scaleFactor;
854         fEMapWithVector[materialID][p][T].push    854         fEMapWithVector[materialID][p][T].push_back(E);
855       }                                           855       }
856     }                                             856     }
857   }                                               857   }
858 }                                                 858 }
859                                                   859