Geant4 Cross Reference |
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