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 // Based on the work described in 26 // Based on the work described in 27 // Rad Res 163, 98-111 (2005) 27 // Rad Res 163, 98-111 (2005) 28 // D. Emfietzoglou, H. Nikjoo 28 // D. Emfietzoglou, H. Nikjoo 29 // 29 // 30 // Authors of the class (2014): 30 // Authors of the class (2014): 31 // I. Kyriakou (kyriak@cc.uoi.gr) 31 // I. Kyriakou (kyriak@cc.uoi.gr) 32 // D. Emfietzoglou (demfietz@cc.uoi.gr) 32 // D. Emfietzoglou (demfietz@cc.uoi.gr) 33 // S. Incerti (incerti@cenbg.in2p3.fr) 33 // S. Incerti (incerti@cenbg.in2p3.fr) 34 // 34 // 35 35 36 #include "G4DNAEmfietzoglouIonisationModel.hh" 36 #include "G4DNAEmfietzoglouIonisationModel.hh" 37 #include "G4PhysicalConstants.hh" 37 #include "G4PhysicalConstants.hh" 38 #include "G4SystemOfUnits.hh" 38 #include "G4SystemOfUnits.hh" 39 #include "G4UAtomicDeexcitation.hh" 39 #include "G4UAtomicDeexcitation.hh" 40 #include "G4LossTableManager.hh" 40 #include "G4LossTableManager.hh" 41 #include "G4DNAChemistryManager.hh" 41 #include "G4DNAChemistryManager.hh" 42 #include "G4DNAMolecularMaterial.hh" 42 #include "G4DNAMolecularMaterial.hh" 43 #include "G4DNABornAngle.hh" 43 #include "G4DNABornAngle.hh" 44 #include "G4DeltaAngle.hh" 44 #include "G4DeltaAngle.hh" 45 45 46 //....oooOO0OOooo........oooOO0OOooo........oo 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 47 47 48 using namespace std; 48 using namespace std; 49 49 50 //....oooOO0OOooo........oooOO0OOooo........oo 50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 51 51 52 G4DNAEmfietzoglouIonisationModel::G4DNAEmfietz 52 G4DNAEmfietzoglouIonisationModel::G4DNAEmfietzoglouIonisationModel(const G4ParticleDefinition*, 53 53 const G4String& nam) : 54 G4VEmModel(nam) << 54 G4VEmModel(nam), isInitialised(false) 55 { 55 { 56 verboseLevel = 0; 56 verboseLevel = 0; 57 // Verbosity scale: 57 // Verbosity scale: 58 // 0 = nothing 58 // 0 = nothing 59 // 1 = warning for energy non-conservation 59 // 1 = warning for energy non-conservation 60 // 2 = details of energy budget 60 // 2 = details of energy budget 61 // 3 = calculation of cross sections, file o 61 // 3 = calculation of cross sections, file openings, sampling of atoms 62 // 4 = entering in methods 62 // 4 = entering in methods 63 63 64 if(verboseLevel > 0) 64 if(verboseLevel > 0) 65 { 65 { 66 G4cout << "Emfietzoglou ionisation model i 66 G4cout << "Emfietzoglou ionisation model is constructed " << G4endl; 67 } 67 } 68 68 69 // Mark this model as "applicable" for atomi 69 // Mark this model as "applicable" for atomic deexcitation 70 SetDeexcitationFlag(true); 70 SetDeexcitationFlag(true); 71 fAtomDeexcitation = nullptr; << 71 fAtomDeexcitation = 0; 72 fParticleChangeForGamma = nullptr; << 72 fParticleChangeForGamma = 0; 73 fpMolWaterDensity = nullptr; << 73 fpMolWaterDensity = 0; 74 74 75 // Define default angular generator 75 // Define default angular generator 76 SetAngularDistribution(new G4DNABornAngle()) 76 SetAngularDistribution(new G4DNABornAngle()); 77 77 78 SetLowEnergyLimit(10. * eV); 78 SetLowEnergyLimit(10. * eV); 79 SetHighEnergyLimit(10. * keV); 79 SetHighEnergyLimit(10. * keV); 80 80 81 // Selection of computation method 81 // Selection of computation method 82 82 83 fasterCode = false; 83 fasterCode = false; 84 84 85 // Selection of stationary mode 85 // Selection of stationary mode 86 86 87 statCode = false; 87 statCode = false; 88 } 88 } 89 89 90 //....oooOO0OOooo........oooOO0OOooo........oo 90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 91 91 92 G4DNAEmfietzoglouIonisationModel::~G4DNAEmfiet 92 G4DNAEmfietzoglouIonisationModel::~G4DNAEmfietzoglouIonisationModel() 93 { 93 { 94 // Cross section 94 // Cross section 95 95 96 std::map<G4String, G4DNACrossSectionDataSet* 96 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos; 97 for(pos = tableData.begin(); pos != tableDat 97 for(pos = tableData.begin(); pos != tableData.end(); ++pos) 98 { 98 { 99 G4DNACrossSectionDataSet* table = pos->sec 99 G4DNACrossSectionDataSet* table = pos->second; 100 delete table; 100 delete table; 101 } 101 } 102 102 103 // Final state 103 // Final state 104 104 105 eVecm.clear(); 105 eVecm.clear(); 106 106 107 } 107 } 108 108 109 //....oooOO0OOooo........oooOO0OOooo........oo 109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 110 110 111 void G4DNAEmfietzoglouIonisationModel::Initial 111 void G4DNAEmfietzoglouIonisationModel::Initialise(const G4ParticleDefinition* particle, 112 112 const G4DataVector& /*cuts*/) 113 { 113 { 114 114 115 if(verboseLevel > 3) 115 if(verboseLevel > 3) 116 { 116 { 117 G4cout << "Calling G4DNAEmfietzoglouIonisa 117 G4cout << "Calling G4DNAEmfietzoglouIonisationModel::Initialise()" << G4endl; 118 } 118 } 119 119 120 // Energy limits 120 // Energy limits 121 121 122 G4String fileElectron("dna/sigma_ionisation_ 122 G4String fileElectron("dna/sigma_ionisation_e_emfietzoglou"); 123 123 124 G4ParticleDefinition* electronDef = G4Electr 124 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition(); 125 125 126 G4String electron; 126 G4String electron; 127 127 128 G4double scaleFactor = (1.e-22 / 3.343) * m* 128 G4double scaleFactor = (1.e-22 / 3.343) * m*m; 129 129 130 const char *path = G4FindDataDir("G4LEDATA") << 130 char *path = getenv("G4LEDATA"); 131 131 132 // *** ELECTRON 132 // *** ELECTRON 133 133 134 electron = electronDef->GetParticleName(); 134 electron = electronDef->GetParticleName(); 135 135 136 tableFile[electron] = fileElectron; 136 tableFile[electron] = fileElectron; 137 137 138 // Cross section 138 // Cross section 139 139 140 auto tableE = new G4DNACrossSectionDataSet( << 140 G4DNACrossSectionDataSet* tableE = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor ); 141 tableE->LoadData(fileElectron); 141 tableE->LoadData(fileElectron); 142 142 143 tableData[electron] = tableE; 143 tableData[electron] = tableE; 144 144 145 // Final state 145 // Final state 146 146 147 std::ostringstream eFullFileName; 147 std::ostringstream eFullFileName; 148 148 149 if (fasterCode) eFullFileName << path << "/d 149 if (fasterCode) eFullFileName << path << "/dna/sigmadiff_cumulated_ionisation_e_emfietzoglou.dat"; 150 if (!fasterCode) eFullFileName << path << "/ 150 if (!fasterCode) eFullFileName << path << "/dna/sigmadiff_ionisation_e_emfietzoglou.dat"; 151 151 152 std::ifstream eDiffCrossSection(eFullFileNam 152 std::ifstream eDiffCrossSection(eFullFileName.str().c_str()); 153 153 154 if (!eDiffCrossSection) 154 if (!eDiffCrossSection) 155 { 155 { 156 if (fasterCode) G4Exception("G4DNAEmfietzo 156 if (fasterCode) G4Exception("G4DNAEmfietzoglouIonisationModel::Initialise","em0003", 157 FatalException,"Missing data file:/dna 157 FatalException,"Missing data file:/dna/sigmadiff_cumulated_ionisation_e_emfietzoglou.dat"); 158 158 159 if (!fasterCode) G4Exception("G4DNAEmfietz 159 if (!fasterCode) G4Exception("G4DNAEmfietzoglouIonisationModel::Initialise","em0003", 160 FatalException,"Missing data file:/dna 160 FatalException,"Missing data file:/dna/sigmadiff_ionisation_e_emfietzoglou.dat"); 161 } 161 } 162 162 163 // 163 // 164 164 165 // Clear the arrays for re-initialization ca 165 // Clear the arrays for re-initialization case (MT mode) 166 // March 25th, 2014 - Vaclav Stepan, Sebasti 166 // March 25th, 2014 - Vaclav Stepan, Sebastien Incerti 167 167 168 eTdummyVec.clear(); 168 eTdummyVec.clear(); 169 169 170 eVecm.clear(); 170 eVecm.clear(); 171 171 172 eProbaShellMap->clear(); 172 eProbaShellMap->clear(); 173 173 174 eDiffCrossSectionData->clear(); 174 eDiffCrossSectionData->clear(); 175 175 176 eNrjTransfData->clear(); 176 eNrjTransfData->clear(); 177 177 178 // 178 // 179 179 180 eTdummyVec.push_back(0.); 180 eTdummyVec.push_back(0.); 181 while(!eDiffCrossSection.eof()) 181 while(!eDiffCrossSection.eof()) 182 { 182 { 183 G4double tDummy; 183 G4double tDummy; 184 G4double eDummy; 184 G4double eDummy; 185 eDiffCrossSection>>tDummy>>eDummy; 185 eDiffCrossSection>>tDummy>>eDummy; 186 if (tDummy != eTdummyVec.back()) eTdummyVe 186 if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy); 187 for (G4int j=0; j<5; j++) 187 for (G4int j=0; j<5; j++) 188 { 188 { 189 eDiffCrossSection>>eDiffCrossSectionData 189 eDiffCrossSection>>eDiffCrossSectionData[j][tDummy][eDummy]; 190 190 191 if (fasterCode) 191 if (fasterCode) 192 { 192 { 193 eNrjTransfData[j][tDummy][eDiffCrossSe 193 eNrjTransfData[j][tDummy][eDiffCrossSectionData[j][tDummy][eDummy]]=eDummy; 194 eProbaShellMap[j][tDummy].push_back(eD 194 eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]); 195 } 195 } 196 196 197 // SI - only if eof is not reached 197 // SI - only if eof is not reached 198 if (!eDiffCrossSection.eof() && !fasterC 198 if (!eDiffCrossSection.eof() && !fasterCode) 199 { 199 { 200 eDiffCrossSectionData[j][tDummy][eDumm 200 eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor; 201 } 201 } 202 202 203 if (!fasterCode) eVecm[tDummy].push_back 203 if (!fasterCode) eVecm[tDummy].push_back(eDummy); 204 204 205 } 205 } 206 } 206 } 207 207 208 // 208 // 209 209 210 if( verboseLevel>0 ) 210 if( verboseLevel>0 ) 211 { 211 { 212 G4cout << "Emfietzoglou ionisation model i 212 G4cout << "Emfietzoglou ionisation model is initialized " << G4endl 213 << "Energy range: " 213 << "Energy range: " 214 << LowEnergyLimit() / eV << " eV - " 214 << LowEnergyLimit() / eV << " eV - " 215 << HighEnergyLimit() / keV << " keV for " 215 << HighEnergyLimit() / keV << " keV for " 216 << particle->GetParticleName() 216 << particle->GetParticleName() 217 << G4endl; 217 << G4endl; 218 } 218 } 219 219 220 // Initialize water density pointer 220 // Initialize water density pointer 221 221 222 fpMolWaterDensity = 222 fpMolWaterDensity = 223 G4DNAMolecularMaterial::Instance()-> 223 G4DNAMolecularMaterial::Instance()-> 224 GetNumMolPerVolTableFor(G4Material::Ge 224 GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER")); 225 225 226 // AD 226 // AD 227 227 228 fAtomDeexcitation = G4LossTableManager::Inst 228 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation(); 229 229 230 if (isInitialised) 230 if (isInitialised) 231 { return;} 231 { return;} 232 fParticleChangeForGamma = GetParticleChangeF 232 fParticleChangeForGamma = GetParticleChangeForGamma(); 233 isInitialised = true; 233 isInitialised = true; 234 } 234 } 235 235 236 //....oooOO0OOooo........oooOO0OOooo........oo 236 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 237 237 238 G4double G4DNAEmfietzoglouIonisationModel:: 238 G4double G4DNAEmfietzoglouIonisationModel:: 239 CrossSectionPerVolume(const G4Material* materi 239 CrossSectionPerVolume(const G4Material* material, 240 const G4ParticleDefiniti 240 const G4ParticleDefinition* particleDefinition, 241 G4double ekin, 241 G4double ekin, 242 G4double, 242 G4double, 243 G4double) 243 G4double) 244 { 244 { 245 if(verboseLevel > 3) 245 if(verboseLevel > 3) 246 { 246 { 247 G4cout 247 G4cout 248 << "Calling CrossSectionPerVolume() of 248 << "Calling CrossSectionPerVolume() of G4DNAEmfietzoglouIonisationModel" 249 << G4endl; 249 << G4endl; 250 } 250 } 251 251 252 if (particleDefinition != G4Electron::Electr 252 if (particleDefinition != G4Electron::ElectronDefinition()) return 0; // necessary ?? 253 253 254 // Calculate total cross section for model 254 // Calculate total cross section for model 255 255 256 G4double sigma=0; 256 G4double sigma=0; 257 257 258 G4double waterDensity = (*fpMolWaterDensity) 258 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()]; 259 259 260 const G4String& particleName = particleDefin 260 const G4String& particleName = particleDefinition->GetParticleName(); 261 261 262 if (ekin >= LowEnergyLimit() && ekin <= High 262 if (ekin >= LowEnergyLimit() && ekin <= HighEnergyLimit()) 263 { 263 { 264 std::map< G4String,G4DNACrossSectionDataSe 264 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos; 265 pos = tableData.find(particleName); 265 pos = tableData.find(particleName); 266 266 267 if (pos != tableData.end()) 267 if (pos != tableData.end()) 268 { 268 { 269 G4DNACrossSectionDataSet* table = pos->s 269 G4DNACrossSectionDataSet* table = pos->second; 270 if (table != nullptr) << 270 if (table != 0) 271 { 271 { 272 sigma = table->FindValue(ekin); 272 sigma = table->FindValue(ekin); 273 } 273 } 274 } 274 } 275 else 275 else 276 { 276 { 277 G4Exception("G4DNAEmfietzoglouIonisation 277 G4Exception("G4DNAEmfietzoglouIonisationModel::CrossSectionPerVolume","em0002", 278 FatalException,"Model not applicable 278 FatalException,"Model not applicable to particle type."); 279 } 279 } 280 } 280 } 281 281 282 if (verboseLevel > 2) 282 if (verboseLevel > 2) 283 { 283 { 284 G4cout << "_______________________________ 284 G4cout << "__________________________________" << G4endl; 285 G4cout << "G4DNAEmfietzoglouIonisationMode 285 G4cout << "G4DNAEmfietzoglouIonisationModel - XS INFO START" << G4endl; 286 G4cout << "Kinetic energy(eV)=" << ekin/eV 286 G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl; 287 G4cout << "Cross section per water molecul 287 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl; 288 G4cout << "Cross section per water molecul 288 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl; 289 G4cout << "G4DNAEmfietzoglouIonisationMode 289 G4cout << "G4DNAEmfietzoglouIonisationModel - XS INFO END" << G4endl; 290 } 290 } 291 291 292 return sigma*waterDensity; 292 return sigma*waterDensity; 293 } 293 } 294 294 295 //....oooOO0OOooo........oooOO0OOooo........oo 295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 296 296 297 void G4DNAEmfietzoglouIonisationModel:: 297 void G4DNAEmfietzoglouIonisationModel:: 298 SampleSecondaries(std::vector<G4DynamicParticl 298 SampleSecondaries(std::vector<G4DynamicParticle*>* fvect, 299 const G4MaterialCutsCouple* 299 const G4MaterialCutsCouple* couple, 300 const G4DynamicParticle* par 300 const G4DynamicParticle* particle, 301 G4double, 301 G4double, 302 G4double) 302 G4double) 303 { 303 { 304 304 305 if(verboseLevel > 3) 305 if(verboseLevel > 3) 306 { 306 { 307 G4cout << "Calling SampleSecondaries() of 307 G4cout << "Calling SampleSecondaries() of G4DNAEmfietzoglouIonisationModel" 308 << G4endl; 308 << G4endl; 309 } 309 } 310 310 311 G4double k = particle->GetKineticEnergy(); 311 G4double k = particle->GetKineticEnergy(); 312 312 313 const G4String& particleName = particle->Get 313 const G4String& particleName = particle->GetDefinition()->GetParticleName(); 314 314 315 if (k >= LowEnergyLimit() && k <= HighEnergy 315 if (k >= LowEnergyLimit() && k <= HighEnergyLimit()) 316 { 316 { 317 G4ParticleMomentum primaryDirection = part 317 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection(); 318 G4double particleMass = particle->GetDefin 318 G4double particleMass = particle->GetDefinition()->GetPDGMass(); 319 G4double totalEnergy = k + particleMass; 319 G4double totalEnergy = k + particleMass; 320 G4double pSquare = k * (totalEnergy + part 320 G4double pSquare = k * (totalEnergy + particleMass); 321 G4double totalMomentum = std::sqrt(pSquare 321 G4double totalMomentum = std::sqrt(pSquare); 322 322 323 G4int ionizationShell = 0; 323 G4int ionizationShell = 0; 324 324 325 ionizationShell = RandomSelect(k,particleN 325 ionizationShell = RandomSelect(k,particleName); 326 326 327 G4double bindingEnergy = 0; 327 G4double bindingEnergy = 0; 328 bindingEnergy = waterStructure.IonisationE 328 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell); 329 329 330 // SI : additional protection if tcs inter 330 // SI : additional protection if tcs interpolation method is modified 331 if (k<bindingEnergy) return; 331 if (k<bindingEnergy) return; 332 // 332 // 333 333 334 G4double secondaryKinetic=-1000*eV; 334 G4double secondaryKinetic=-1000*eV; 335 335 336 if (!fasterCode) secondaryKinetic = Random 336 if (!fasterCode) secondaryKinetic = RandomizeEjectedElectronEnergy(particle->GetDefinition(),k,ionizationShell); 337 337 338 if (fasterCode) 338 if (fasterCode) 339 secondaryKinetic = RandomizeEjectedElectro 339 secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(particle->GetDefinition(),k,ionizationShell); 340 340 341 // SI - For atom. deexc. tagging - 23/05/2 341 // SI - For atom. deexc. tagging - 23/05/2017 342 342 343 G4int Z = 8; 343 G4int Z = 8; 344 344 345 G4ThreeVector deltaDirection = 345 G4ThreeVector deltaDirection = 346 GetAngularDistribution()->SampleDirectionF 346 GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic, 347 Z, ionizationShell, 347 Z, ionizationShell, 348 couple->GetMaterial()); 348 couple->GetMaterial()); 349 349 350 if (secondaryKinetic>0) 350 if (secondaryKinetic>0) 351 { 351 { 352 auto dp = new G4DynamicParticle (G4Elec << 352 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic); 353 fvect->push_back(dp); 353 fvect->push_back(dp); 354 } 354 } 355 355 356 G4double deltaTotalMomentum = std::sqrt(se 356 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 )); 357 357 358 G4double finalPx = totalMomentum*primaryDi 358 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x(); 359 G4double finalPy = totalMomentum*primaryDi 359 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y(); 360 G4double finalPz = totalMomentum*primaryDi 360 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z(); 361 G4double finalMomentum = std::sqrt(finalPx 361 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz); 362 finalPx /= finalMomentum; 362 finalPx /= finalMomentum; 363 finalPy /= finalMomentum; 363 finalPy /= finalMomentum; 364 finalPz /= finalMomentum; 364 finalPz /= finalMomentum; 365 365 366 G4ThreeVector direction; 366 G4ThreeVector direction; 367 direction.set(finalPx,finalPy,finalPz); 367 direction.set(finalPx,finalPy,finalPz); 368 368 369 fParticleChangeForGamma->ProposeMomentumDi 369 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()); 370 370 371 // AM: sample deexcitation 371 // AM: sample deexcitation 372 // here we assume that H_{2}O electronic l 372 // here we assume that H_{2}O electronic levels are the same as Oxygen. 373 // this can be considered true with a roug 373 // this can be considered true with a rough 10% error in energy on K-shell, 374 374 375 size_t secNumberInit = 0;// need to know a 375 size_t secNumberInit = 0;// need to know at a certain point the energy of secondaries 376 size_t secNumberFinal = 0;// So I'll make 376 size_t secNumberFinal = 0;// So I'll make the diference and then sum the energies 377 377 378 G4double scatteredEnergy = k-bindingEnergy 378 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic; 379 379 380 // SI: only atomic deexcitation from K she 380 // SI: only atomic deexcitation from K shell is considered 381 if((fAtomDeexcitation != nullptr) && ioniz << 381 if(fAtomDeexcitation && ionizationShell == 4) 382 { 382 { 383 const G4AtomicShell* shell 383 const G4AtomicShell* shell 384 = fAtomDeexcitation->GetAtomicShell(Z, 384 = fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0)); 385 secNumberInit = fvect->size(); 385 secNumberInit = fvect->size(); 386 fAtomDeexcitation->GenerateParticles(fve 386 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0); 387 secNumberFinal = fvect->size(); 387 secNumberFinal = fvect->size(); 388 388 389 if(secNumberFinal > secNumberInit) { 389 if(secNumberFinal > secNumberInit) { 390 for (size_t i=secNumberInit; i<secNumberFina 390 for (size_t i=secNumberInit; i<secNumberFinal; ++i) { 391 //Check if there is enough residual 391 //Check if there is enough residual energy 392 if (bindingEnergy >= ((*fvect)[i])-> 392 if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy()) 393 { 393 { 394 //Ok, this is a valid secondary: 394 //Ok, this is a valid secondary: keep it 395 bindingEnergy -= ((*fvect)[i])->GetKine 395 bindingEnergy -= ((*fvect)[i])->GetKineticEnergy(); 396 } 396 } 397 else 397 else 398 { 398 { 399 //Invalid secondary: not enough energy 399 //Invalid secondary: not enough energy to create it! 400 //Keep its energy in the local deposit 400 //Keep its energy in the local deposit 401 delete (*fvect)[i]; 401 delete (*fvect)[i]; 402 (*fvect)[i]=nullptr; << 402 (*fvect)[i]=0; 403 } 403 } 404 } 404 } 405 } 405 } 406 406 407 } 407 } 408 408 409 //This should never happen 409 //This should never happen 410 if(bindingEnergy < 0.0) 410 if(bindingEnergy < 0.0) 411 G4Exception("G4DNAEmfietzoglouIonisatioMo 411 G4Exception("G4DNAEmfietzoglouIonisatioModel1::SampleSecondaries()", 412 "em2050",FatalException,"Nega 412 "em2050",FatalException,"Negative local energy deposit"); 413 413 414 //bindingEnergy has been decreased 414 //bindingEnergy has been decreased 415 //by the amount of energy taken away by de 415 //by the amount of energy taken away by deexc. products 416 if (!statCode) 416 if (!statCode) 417 { 417 { 418 fParticleChangeForGamma->SetProposedKine 418 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy); 419 fParticleChangeForGamma->ProposeLocalEne 419 fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy); 420 } 420 } 421 else 421 else 422 { 422 { 423 fParticleChangeForGamma->SetProposedKine 423 fParticleChangeForGamma->SetProposedKineticEnergy(k); 424 fParticleChangeForGamma->ProposeLocalEne 424 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy); 425 } 425 } 426 // TEST ////////////////////////// 426 // TEST ////////////////////////// 427 // if (secondaryKinetic<0) abort(); 427 // if (secondaryKinetic<0) abort(); 428 // if (scatteredEnergy<0) abort(); 428 // if (scatteredEnergy<0) abort(); 429 // if (k-scatteredEnergy-secondaryKinetic- 429 // if (k-scatteredEnergy-secondaryKinetic-deexSecEnergy<0) abort(); 430 // if (k-scatteredEnergy<0) abort(); 430 // if (k-scatteredEnergy<0) abort(); 431 ///////////////////////////////// 431 ///////////////////////////////// 432 432 433 const G4Track * theIncomingTrack = fPartic 433 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack(); 434 G4DNAChemistryManager::Instance()->CreateW 434 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule, 435 ionizationShell, 435 ionizationShell, 436 theIncomingTrack); 436 theIncomingTrack); 437 } 437 } 438 438 439 } 439 } 440 440 441 //....oooOO0OOooo........oooOO0OOooo........oo 441 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 442 442 443 G4double 443 G4double 444 G4DNAEmfietzoglouIonisationModel:: 444 G4DNAEmfietzoglouIonisationModel:: 445 RandomizeEjectedElectronEnergy(G4ParticleDefin 445 RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition, 446 G4double k, 446 G4double k, 447 G4int shell) 447 G4int shell) 448 { 448 { 449 // G4cout << "*** SLOW computation for " 449 // G4cout << "*** SLOW computation for " 450 // << " " << particleDefinition->GetP 450 // << " " << particleDefinition->GetParticleName() << G4endl; 451 451 452 if(particleDefinition == G4Electron::Electro 452 if(particleDefinition == G4Electron::ElectronDefinition()) 453 { 453 { 454 G4double maximumEnergyTransfer = 0.; 454 G4double maximumEnergyTransfer = 0.; 455 if((k + waterStructure.IonisationEnergy(sh 455 if((k + waterStructure.IonisationEnergy(shell)) / 2. > k) 456 maximumEnergyTransfer = k; 456 maximumEnergyTransfer = k; 457 else 457 else 458 maximumEnergyTransfer = (k + waterStruct 458 maximumEnergyTransfer = (k + waterStructure.IonisationEnergy(shell))/ 2.; 459 459 460 // SI : original method 460 // SI : original method 461 /* 461 /* 462 G4double crossSectionMaximum = 0.; 462 G4double crossSectionMaximum = 0.; 463 for(G4double value=waterStructure.Ionisat 463 for(G4double value=waterStructure.IonisationEnergy(shell); value<=maximumEnergyTransfer; value+=0.1*eV) 464 { 464 { 465 G4double differentialCrossSection = Diffe 465 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell); 466 if(differentialCrossSection >= crossSecti 466 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection; 467 } 467 } 468 */ 468 */ 469 469 470 // SI : alternative method 470 // SI : alternative method 471 G4double crossSectionMaximum = 0.; 471 G4double crossSectionMaximum = 0.; 472 472 473 G4double minEnergy = waterStructure.Ionisa 473 G4double minEnergy = waterStructure.IonisationEnergy(shell); 474 G4double maxEnergy = maximumEnergyTransfer 474 G4double maxEnergy = maximumEnergyTransfer; 475 G4int nEnergySteps = 50; 475 G4int nEnergySteps = 50; 476 476 477 G4double value(minEnergy); 477 G4double value(minEnergy); 478 G4double stpEnergy(std::pow(maxEnergy / va 478 G4double stpEnergy(std::pow(maxEnergy / value, 479 1. / static_ca 479 1. / static_cast<G4double>(nEnergySteps - 1))); 480 G4int step(nEnergySteps); 480 G4int step(nEnergySteps); 481 while(step > 0) 481 while(step > 0) 482 { 482 { 483 step--; 483 step--; 484 G4double differentialCrossSection = 484 G4double differentialCrossSection = 485 DifferentialCrossSection(particleDef 485 DifferentialCrossSection(particleDefinition, 486 k / eV, 486 k / eV, 487 value / eV, 487 value / eV, 488 shell); 488 shell); 489 if(differentialCrossSection >= crossSect 489 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = 490 differentialCrossSection; 490 differentialCrossSection; 491 value *= stpEnergy; 491 value *= stpEnergy; 492 } 492 } 493 // 493 // 494 494 495 G4double secondaryElectronKineticEnergy = 495 G4double secondaryElectronKineticEnergy = 0.; 496 do 496 do 497 { 497 { 498 secondaryElectronKineticEnergy = G4Unifo 498 secondaryElectronKineticEnergy = G4UniformRand()* (maximumEnergyTransfer-waterStructure.IonisationEnergy(shell)); 499 }while(G4UniformRand()*crossSectionMaximum 499 }while(G4UniformRand()*crossSectionMaximum > 500 DifferentialCrossSection(particleDefin 500 DifferentialCrossSection(particleDefinition, k/eV, 501 (secondaryElectronKineticEnergy+wa 501 (secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell)); 502 502 503 return secondaryElectronKineticEnergy; 503 return secondaryElectronKineticEnergy; 504 504 505 } 505 } 506 506 507 return 0; 507 return 0; 508 } 508 } 509 509 510 //....oooOO0OOooo........oooOO0OOooo........oo 510 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 511 511 512 // The following section is not used anymore b 512 // The following section is not used anymore but is kept for memory 513 // GetAngularDistribution()->SampleDirectionFo 513 // GetAngularDistribution()->SampleDirectionForShell is used instead 514 514 515 /* 515 /* 516 void G4DNAEmfietzoglouIonisationModel::Random 516 void G4DNAEmfietzoglouIonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition, 517 G4double k, 517 G4double k, 518 G4double secKinetic, 518 G4double secKinetic, 519 G4double & cosTheta, 519 G4double & cosTheta, 520 G4double & phi ) 520 G4double & phi ) 521 { 521 { 522 if (particleDefinition == G4Electron::Electro 522 if (particleDefinition == G4Electron::ElectronDefinition()) 523 { 523 { 524 phi = twopi * G4UniformRand(); 524 phi = twopi * G4UniformRand(); 525 if (secKinetic < 50.*eV) cosTheta = (2.*G4Uni 525 if (secKinetic < 50.*eV) cosTheta = (2.*G4UniformRand())-1.; 526 else if (secKinetic <= 200.*eV) 526 else if (secKinetic <= 200.*eV) 527 { 527 { 528 if (G4UniformRand() <= 0.1) cosTheta = (2.*G4 528 if (G4UniformRand() <= 0.1) cosTheta = (2.*G4UniformRand())-1.; 529 else cosTheta = G4UniformRand()*(std::sqrt(2. 529 else cosTheta = G4UniformRand()*(std::sqrt(2.)/2); 530 } 530 } 531 else 531 else 532 { 532 { 533 G4double sin2O = (1.-secKinetic/k) / (1.+secK 533 G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2)); 534 cosTheta = std::sqrt(1.-sin2O); 534 cosTheta = std::sqrt(1.-sin2O); 535 } 535 } 536 } 536 } 537 537 538 else if (particleDefinition == G4Proton::Prot 538 else if (particleDefinition == G4Proton::ProtonDefinition()) 539 { 539 { 540 G4double maxSecKinetic = 4.* (electron_mass_c 540 G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k; 541 phi = twopi * G4UniformRand(); 541 phi = twopi * G4UniformRand(); 542 542 543 // cosTheta = std::sqrt(secKinetic / maxSecKi 543 // cosTheta = std::sqrt(secKinetic / maxSecKinetic); 544 544 545 // Restriction below 100 eV from Emfietzoglou 545 // Restriction below 100 eV from Emfietzoglou (2000) 546 546 547 if (secKinetic>100*eV) cosTheta = std::sqrt(s 547 if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic); 548 else cosTheta = (2.*G4UniformRand())-1.; 548 else cosTheta = (2.*G4UniformRand())-1.; 549 549 550 } 550 } 551 } 551 } 552 */ 552 */ 553 553 554 //....oooOO0OOooo........oooOO0OOooo........oo 554 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 555 G4double G4DNAEmfietzoglouIonisationModel::Dif 555 G4double G4DNAEmfietzoglouIonisationModel::DifferentialCrossSection(G4ParticleDefinition * particleDefinition, 556 556 G4double k, 557 557 G4double energyTransfer, 558 558 G4int ionizationLevelIndex) 559 { 559 { 560 G4double sigma = 0.; 560 G4double sigma = 0.; 561 561 562 if(energyTransfer >= waterStructure.Ionisati 562 if(energyTransfer >= waterStructure.IonisationEnergy(ionizationLevelIndex)/eV) 563 { 563 { 564 G4double valueT1 = 0; 564 G4double valueT1 = 0; 565 G4double valueT2 = 0; 565 G4double valueT2 = 0; 566 G4double valueE21 = 0; 566 G4double valueE21 = 0; 567 G4double valueE22 = 0; 567 G4double valueE22 = 0; 568 G4double valueE12 = 0; 568 G4double valueE12 = 0; 569 G4double valueE11 = 0; 569 G4double valueE11 = 0; 570 570 571 G4double xs11 = 0; 571 G4double xs11 = 0; 572 G4double xs12 = 0; 572 G4double xs12 = 0; 573 G4double xs21 = 0; 573 G4double xs21 = 0; 574 G4double xs22 = 0; 574 G4double xs22 = 0; 575 575 576 if(particleDefinition == G4Electron::Elect 576 if(particleDefinition == G4Electron::ElectronDefinition()) 577 { 577 { 578 // Protection against out of boundary ac 578 // Protection against out of boundary access 579 if (k==eTdummyVec.back()) k=k*(1.-1e-12) 579 if (k==eTdummyVec.back()) k=k*(1.-1e-12); 580 // 580 // 581 581 582 // k should be in eV and energy transfer 582 // k should be in eV and energy transfer eV also 583 583 584 auto t2 = std::upper_bound(eTdummyVec.be << 584 std::vector<G4double>::iterator t2 = std::upper_bound(eTdummyVec.begin(), 585 585 eTdummyVec.end(), 586 586 k); 587 587 588 auto t1 = t2 - 1; << 588 std::vector<G4double>::iterator t1 = t2 - 1; 589 589 590 // SI : the following condition avoids s 590 // SI : the following condition avoids situations where energyTransfer >last vector element 591 // added strict limitations (09/08/2017) 591 // added strict limitations (09/08/2017) 592 if(energyTransfer < eVecm[(*t1)].back() 592 if(energyTransfer < eVecm[(*t1)].back() && 593 energyTransfer < eVecm[(*t2)].back()) 593 energyTransfer < eVecm[(*t2)].back()) 594 { 594 { 595 auto e12 = << 595 std::vector<G4double>::iterator e12 = 596 std::upper_bound(eVecm[(*t1)].begi 596 std::upper_bound(eVecm[(*t1)].begin(), 597 eVecm[(*t1)].end( 597 eVecm[(*t1)].end(), 598 energyTransfer); 598 energyTransfer); 599 auto e11 = e12 - 1; << 599 std::vector<G4double>::iterator e11 = e12 - 1; 600 600 601 auto e22 = << 601 std::vector<G4double>::iterator e22 = 602 std::upper_bound(eVecm[(*t2)].begi 602 std::upper_bound(eVecm[(*t2)].begin(), 603 eVecm[(*t2)].end( 603 eVecm[(*t2)].end(), 604 energyTransfer); 604 energyTransfer); 605 auto e21 = e22 - 1; << 605 std::vector<G4double>::iterator e21 = e22 - 1; 606 606 607 valueT1 = *t1; 607 valueT1 = *t1; 608 valueT2 = *t2; 608 valueT2 = *t2; 609 valueE21 = *e21; 609 valueE21 = *e21; 610 valueE22 = *e22; 610 valueE22 = *e22; 611 valueE12 = *e12; 611 valueE12 = *e12; 612 valueE11 = *e11; 612 valueE11 = *e11; 613 613 614 xs11 = eDiffCrossSectionData[ionizatio 614 xs11 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11]; 615 xs12 = eDiffCrossSectionData[ionizatio 615 xs12 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12]; 616 xs21 = eDiffCrossSectionData[ionizatio 616 xs21 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21]; 617 xs22 = eDiffCrossSectionData[ionizatio 617 xs22 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22]; 618 618 619 //G4cout << "-------------------" << G 619 //G4cout << "-------------------" << G4endl; 620 //G4cout << "ionizationLevelIndex=" << 620 //G4cout << "ionizationLevelIndex=" << ionizationLevelIndex << G4endl; 621 //G4cout << "valueT1/eV=" << valueT1 < 621 //G4cout << "valueT1/eV=" << valueT1 << " valueT2/eV=" << valueT2 << G4endl; 622 //G4cout << "valueE11/eV=" << valueE11 622 //G4cout << "valueE11/eV=" << valueE11 << " valueE12/eV=" << valueE12 623 // << " valueE21/eV=" << valueE2 623 // << " valueE21/eV=" << valueE21 << " valueE22/eV=" << valueE22 << G4endl; 624 //G4cout << "xs11=" << xs11 / ((1.e-22 624 //G4cout << "xs11=" << xs11 / ((1.e-22 / 3.343) * m*m) << G4endl; 625 //G4cout << "xs12=" << xs12 / ((1.e-22 625 //G4cout << "xs12=" << xs12 / ((1.e-22 / 3.343) * m*m) << G4endl; 626 //G4cout << "xs21=" << xs21 / ((1.e-22 626 //G4cout << "xs21=" << xs21 / ((1.e-22 / 3.343) * m*m) << G4endl; 627 //G4cout << "xs22=" << xs22 / ((1.e-22 627 //G4cout << "xs22=" << xs22 / ((1.e-22 / 3.343) * m*m) << G4endl; 628 //G4cout << "###################" << G 628 //G4cout << "###################" << G4endl; 629 629 630 } 630 } 631 631 632 } 632 } 633 633 634 G4double xsProduct = xs11 * xs12 * xs21 * 634 G4double xsProduct = xs11 * xs12 * xs21 * xs22; 635 if(xsProduct != 0.) 635 if(xsProduct != 0.) 636 { 636 { 637 sigma = QuadInterpolator(valueE11, 637 sigma = QuadInterpolator(valueE11, 638 valueE12, 638 valueE12, 639 valueE21, 639 valueE21, 640 valueE22, 640 valueE22, 641 xs11, 641 xs11, 642 xs12, 642 xs12, 643 xs21, 643 xs21, 644 xs22, 644 xs22, 645 valueT1, 645 valueT1, 646 valueT2, 646 valueT2, 647 k, 647 k, 648 energyTransfer) 648 energyTransfer); 649 } 649 } 650 650 651 } 651 } 652 652 653 return sigma; 653 return sigma; 654 } 654 } 655 655 656 //....oooOO0OOooo........oooOO0OOooo........oo 656 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 657 657 658 G4double G4DNAEmfietzoglouIonisationModel::Int 658 G4double G4DNAEmfietzoglouIonisationModel::Interpolate(G4double e1, 659 659 G4double e2, 660 660 G4double e, 661 661 G4double xs1, 662 662 G4double xs2) 663 { 663 { 664 664 665 G4double value = 0.; 665 G4double value = 0.; 666 666 667 // Log-log interpolation by default 667 // Log-log interpolation by default 668 668 669 if(e1 != 0 && e2 != 0 && (std::log10(e2) - s 669 if(e1 != 0 && e2 != 0 && (std::log10(e2) - std::log10(e1)) != 0 670 && !fasterCode) 670 && !fasterCode) 671 { 671 { 672 G4double a = (std::log10(xs2) - std::log10 672 G4double a = (std::log10(xs2) - std::log10(xs1)) 673 / (std::log10(e2) - std::log10(e1)); 673 / (std::log10(e2) - std::log10(e1)); 674 G4double b = std::log10(xs2) - a * std::lo 674 G4double b = std::log10(xs2) - a * std::log10(e2); 675 G4double sigma = a * std::log10(e) + b; 675 G4double sigma = a * std::log10(e) + b; 676 value = (std::pow(10., sigma)); 676 value = (std::pow(10., sigma)); 677 } 677 } 678 678 679 // Switch to lin-lin interpolation 679 // Switch to lin-lin interpolation 680 /* 680 /* 681 if ((e2-e1)!=0) 681 if ((e2-e1)!=0) 682 { 682 { 683 G4double d1 = xs1; 683 G4double d1 = xs1; 684 G4double d2 = xs2; 684 G4double d2 = xs2; 685 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1) 685 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1)); 686 } 686 } 687 */ 687 */ 688 688 689 // Switch to log-lin interpolation for faste 689 // Switch to log-lin interpolation for faster code 690 if((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 && 690 if((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 && fasterCode) 691 { 691 { 692 G4double d1 = std::log10(xs1); 692 G4double d1 = std::log10(xs1); 693 G4double d2 = std::log10(xs2); 693 G4double d2 = std::log10(xs2); 694 value = std::pow(10., (d1 + (d2 - d1) * (e 694 value = std::pow(10., (d1 + (d2 - d1) * (e - e1) / (e2 - e1))); 695 } 695 } 696 696 697 // Switch to lin-lin interpolation for faste 697 // Switch to lin-lin interpolation for faster code 698 // in case one of xs1 or xs2 (=cum proba) va 698 // in case one of xs1 or xs2 (=cum proba) value is zero 699 699 700 if((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) 700 if((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) && fasterCode) 701 { 701 { 702 G4double d1 = xs1; 702 G4double d1 = xs1; 703 G4double d2 = xs2; 703 G4double d2 = xs2; 704 value = (d1 + (d2 - d1) * (e - e1) / (e2 - 704 value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1)); 705 } 705 } 706 706 707 /* 707 /* 708 G4cout 708 G4cout 709 << e1 << " " 709 << e1 << " " 710 << e2 << " " 710 << e2 << " " 711 << e << " " 711 << e << " " 712 << xs1 << " " 712 << xs1 << " " 713 << xs2 << " " 713 << xs2 << " " 714 << value 714 << value 715 << G4endl; 715 << G4endl; 716 */ 716 */ 717 717 718 return value; 718 return value; 719 } 719 } 720 720 721 //....oooOO0OOooo........oooOO0OOooo........oo 721 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 722 722 723 G4double G4DNAEmfietzoglouIonisationModel::Qua 723 G4double G4DNAEmfietzoglouIonisationModel::QuadInterpolator(G4double e11, 724 724 G4double e12, 725 725 G4double e21, 726 726 G4double e22, 727 727 G4double xs11, 728 728 G4double xs12, 729 729 G4double xs21, 730 730 G4double xs22, 731 731 G4double t1, 732 732 G4double t2, 733 733 G4double t, 734 734 G4double e) 735 { 735 { 736 G4double interpolatedvalue1 = Interpolate(e1 736 G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12); 737 G4double interpolatedvalue2 = Interpolate(e2 737 G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22); 738 G4double value = Interpolate(t1, 738 G4double value = Interpolate(t1, 739 t2, 739 t2, 740 t, 740 t, 741 interpolatedval 741 interpolatedvalue1, 742 interpolatedval 742 interpolatedvalue2); 743 743 744 return value; 744 return value; 745 } 745 } 746 746 747 //....oooOO0OOooo........oooOO0OOooo........oo 747 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 748 748 749 G4int G4DNAEmfietzoglouIonisationModel::Random 749 G4int G4DNAEmfietzoglouIonisationModel::RandomSelect(G4double k, 750 750 const G4String& particle) 751 { 751 { 752 G4int level = 0; 752 G4int level = 0; 753 753 754 auto pos = tableData.find(particle); << 754 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos; >> 755 pos = tableData.find(particle); 755 756 756 if(pos != tableData.cend()) << 757 if(pos != tableData.end()) 757 { 758 { 758 G4DNACrossSectionDataSet* table = pos->sec 759 G4DNACrossSectionDataSet* table = pos->second; 759 760 760 if(table != nullptr) << 761 if(table != 0) 761 { 762 { 762 auto valuesBuffer = new G4double[table- << 763 G4double* valuesBuffer = new G4double[table->NumberOfComponents()]; 763 const auto n = (G4int)table->NumberOfCo << 764 const size_t n(table->NumberOfComponents()); 764 G4int i(n); << 765 size_t i(n); 765 G4double value = 0.; 766 G4double value = 0.; 766 767 767 while(i > 0) 768 while(i > 0) 768 { 769 { 769 i--; 770 i--; 770 valuesBuffer[i] = table->GetComponent( 771 valuesBuffer[i] = table->GetComponent(i)->FindValue(k); 771 value += valuesBuffer[i]; 772 value += valuesBuffer[i]; 772 } 773 } 773 774 774 value *= G4UniformRand(); 775 value *= G4UniformRand(); 775 776 776 i = n; 777 i = n; 777 778 778 while(i > 0) 779 while(i > 0) 779 { 780 { 780 i--; 781 i--; 781 782 782 if(valuesBuffer[i] > value) 783 if(valuesBuffer[i] > value) 783 { 784 { 784 delete[] valuesBuffer; 785 delete[] valuesBuffer; 785 return i; 786 return i; 786 } 787 } 787 value -= valuesBuffer[i]; 788 value -= valuesBuffer[i]; 788 } 789 } 789 790 790 delete[] valuesBuffer; << 791 if(valuesBuffer) delete[] valuesBuffer; 791 792 792 } 793 } 793 } 794 } 794 else 795 else 795 { 796 { 796 G4Exception("G4DNAEmfietzoglouIonisationMo 797 G4Exception("G4DNAEmfietzoglouIonisationModel::RandomSelect", 797 "em0002", 798 "em0002", 798 FatalException, 799 FatalException, 799 "Model not applicable to parti 800 "Model not applicable to particle type."); 800 } 801 } 801 802 802 return level; 803 return level; 803 } 804 } 804 805 805 //....oooOO0OOooo........oooOO0OOooo........oo 806 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 806 807 807 G4double G4DNAEmfietzoglouIonisationModel::Ran 808 G4double G4DNAEmfietzoglouIonisationModel::RandomizeEjectedElectronEnergyFromCumulatedDcs(G4ParticleDefinition* particleDefinition, 808 809 G4double k, 809 810 G4int shell) 810 { 811 { 811 //G4cout << "*** FAST computation for " << " 812 //G4cout << "*** FAST computation for " << " " << particleDefinition->GetParticleName() << G4endl; 812 813 813 G4double secondaryElectronKineticEnergy = 0. 814 G4double secondaryElectronKineticEnergy = 0.; 814 815 815 secondaryElectronKineticEnergy = RandomTrans 816 secondaryElectronKineticEnergy = RandomTransferedEnergy(particleDefinition, 816 817 k / eV, 817 818 shell) 818 * eV 819 * eV 819 - waterStru 820 - waterStructure.IonisationEnergy(shell); 820 821 821 //G4cout << RandomTransferedEnergy(particleD 822 //G4cout << RandomTransferedEnergy(particleDefinition, k/eV, shell) << G4endl; 822 if(secondaryElectronKineticEnergy < 0.) retu 823 if(secondaryElectronKineticEnergy < 0.) return 0.; 823 824 824 return secondaryElectronKineticEnergy; 825 return secondaryElectronKineticEnergy; 825 } 826 } 826 827 827 //....oooOO0OOooo........oooOO0OOooo........oo 828 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 828 829 829 G4double G4DNAEmfietzoglouIonisationModel::Ran 830 G4double G4DNAEmfietzoglouIonisationModel::RandomTransferedEnergy(G4ParticleDefinition* particleDefinition, 830 831 G4double k, 831 832 G4int ionizationLevelIndex) 832 { 833 { 833 834 834 G4double random = G4UniformRand(); 835 G4double random = G4UniformRand(); 835 836 836 G4double nrj = 0.; 837 G4double nrj = 0.; 837 838 838 G4double valueK1 = 0; 839 G4double valueK1 = 0; 839 G4double valueK2 = 0; 840 G4double valueK2 = 0; 840 G4double valuePROB21 = 0; 841 G4double valuePROB21 = 0; 841 G4double valuePROB22 = 0; 842 G4double valuePROB22 = 0; 842 G4double valuePROB12 = 0; 843 G4double valuePROB12 = 0; 843 G4double valuePROB11 = 0; 844 G4double valuePROB11 = 0; 844 845 845 G4double nrjTransf11 = 0; 846 G4double nrjTransf11 = 0; 846 G4double nrjTransf12 = 0; 847 G4double nrjTransf12 = 0; 847 G4double nrjTransf21 = 0; 848 G4double nrjTransf21 = 0; 848 G4double nrjTransf22 = 0; 849 G4double nrjTransf22 = 0; 849 850 850 if (particleDefinition == G4Electron::Electr 851 if (particleDefinition == G4Electron::ElectronDefinition()) 851 { 852 { 852 // Protection against out of boundary acce 853 // Protection against out of boundary access 853 if (k==eTdummyVec.back()) k=k*(1.-1e-12); 854 if (k==eTdummyVec.back()) k=k*(1.-1e-12); 854 // 855 // 855 856 856 // k should be in eV 857 // k should be in eV 857 auto k2 = std::upper_bound(eTdummyVec.begi << 858 std::vector<G4double>::iterator k2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k); 858 859 859 auto k1 = k2-1; << 860 std::vector<G4double>::iterator k1 = k2-1; 860 861 861 /* 862 /* 862 G4cout << "----> k=" << k 863 G4cout << "----> k=" << k 863 << " " << *k1 864 << " " << *k1 864 << " " << *k2 865 << " " << *k2 865 << " " << random 866 << " " << random 866 << " " << ionizationLevelIndex 867 << " " << ionizationLevelIndex 867 << " " << eProbaShellMap[ionizationLevelI 868 << " " << eProbaShellMap[ionizationLevelIndex][(*k1)].back() 868 << " " << eProbaShellMap[ionizationLevelI 869 << " " << eProbaShellMap[ionizationLevelIndex][(*k2)].back() 869 << G4endl; 870 << G4endl; 870 */ 871 */ 871 872 872 // SI : the following condition avoids sit 873 // SI : the following condition avoids situations where random >last vector element 873 if ( random <= eProbaShellMap[ionizationLe 874 if ( random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back() 874 && random <= eProbaShellMap[ionization 875 && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back() ) 875 876 876 { 877 { 877 auto prob12 = std::upper_bound(eProbaShe << 878 std::vector<G4double>::iterator prob12 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(), 878 eProbaShellMap[ionizationLevelIndex] 879 eProbaShellMap[ionizationLevelIndex][(*k1)].end(), random); 879 880 880 auto prob11 = prob12-1; << 881 std::vector<G4double>::iterator prob11 = prob12-1; 881 882 882 auto prob22 = std::upper_bound(eProbaShe << 883 std::vector<G4double>::iterator prob22 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(), 883 eProbaShellMap[ionizationLevelIndex] 884 eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random); 884 885 885 auto prob21 = prob22-1; << 886 std::vector<G4double>::iterator prob21 = prob22-1; 886 887 887 valueK1 =*k1; 888 valueK1 =*k1; 888 valueK2 =*k2; 889 valueK2 =*k2; 889 valuePROB21 =*prob21; 890 valuePROB21 =*prob21; 890 valuePROB22 =*prob22; 891 valuePROB22 =*prob22; 891 valuePROB12 =*prob12; 892 valuePROB12 =*prob12; 892 valuePROB11 =*prob11; 893 valuePROB11 =*prob11; 893 894 894 /* 895 /* 895 G4cout << " " << random << " " < 896 G4cout << " " << random << " " << valuePROB11 << " " 896 << valuePROB12 << " " << valuePROB21 << 897 << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl; 897 */ 898 */ 898 899 899 nrjTransf11 = eNrjTransfData[ionizationL 900 nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11]; 900 nrjTransf12 = eNrjTransfData[ionizationL 901 nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12]; 901 nrjTransf21 = eNrjTransfData[ionizationL 902 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21]; 902 nrjTransf22 = eNrjTransfData[ionizationL 903 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22]; 903 904 904 /* 905 /* 905 G4cout << " " << ionizationLevel 906 G4cout << " " << ionizationLevelIndex << " " 906 << random << " " <<valueK1 << " " << va 907 << random << " " <<valueK1 << " " << valueK2 << G4endl; 907 908 908 G4cout << " " << random << " " < 909 G4cout << " " << random << " " << nrjTransf11 << " " 909 << nrjTransf12 << " " << nrjTransf21 << 910 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl; 910 */ 911 */ 911 912 912 } 913 } 913 914 914 // Avoids cases where cum xs is zero for k 915 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2) 915 916 916 if ( random > eProbaShellMap[ionizationLev 917 if ( random > eProbaShellMap[ionizationLevelIndex][(*k1)].back() ) 917 918 918 { 919 { 919 auto prob22 = std::upper_bound(eProbaShe << 920 std::vector<G4double>::iterator prob22 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(), 920 eProbaShellMap[ionizationLevelIndex] 921 eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random); 921 922 922 auto prob21 = prob22-1; << 923 std::vector<G4double>::iterator prob21 = prob22-1; 923 924 924 valueK1 =*k1; 925 valueK1 =*k1; 925 valueK2 =*k2; 926 valueK2 =*k2; 926 valuePROB21 =*prob21; 927 valuePROB21 =*prob21; 927 valuePROB22 =*prob22; 928 valuePROB22 =*prob22; 928 929 929 //G4cout << " " << random << " " 930 //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl; 930 931 931 nrjTransf21 = eNrjTransfData[ionizationL 932 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21]; 932 nrjTransf22 = eNrjTransfData[ionizationL 933 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22]; 933 934 934 G4double interpolatedvalue2 = Interpolat 935 G4double interpolatedvalue2 = Interpolate(valuePROB21, valuePROB22, random, nrjTransf21, nrjTransf22); 935 936 936 // zeros are explicitly set << 937 // zeros are explicitely set 937 938 938 G4double value = Interpolate(valueK1, va 939 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2); 939 940 940 /* 941 /* 941 G4cout << " " << ionizationLevel 942 G4cout << " " << ionizationLevelIndex << " " 942 << random << " " <<valueK1 << " " << va 943 << random << " " <<valueK1 << " " << valueK2 << G4endl; 943 944 944 G4cout << " " << random << " " < 945 G4cout << " " << random << " " << nrjTransf11 << " " 945 << nrjTransf12 << " " << nrjTransf21 << 946 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl; 946 947 947 G4cout << "ici" << " " << value << G4en 948 G4cout << "ici" << " " << value << G4endl; 948 */ 949 */ 949 950 950 return value; 951 return value; 951 } 952 } 952 953 953 } 954 } 954 955 955 // End electron 956 // End electron 956 957 957 G4double nrjTransfProduct = nrjTransf11 * nr 958 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21 * nrjTransf22; 958 959 959 //G4cout << "nrjTransfProduct=" << nrjTransf 960 //G4cout << "nrjTransfProduct=" << nrjTransfProduct << G4endl; 960 961 961 if (nrjTransfProduct != 0.) 962 if (nrjTransfProduct != 0.) 962 { 963 { 963 nrj = QuadInterpolator( valuePROB11, value 964 nrj = QuadInterpolator( valuePROB11, valuePROB12, 964 valuePROB21, valuePROB22, 965 valuePROB21, valuePROB22, 965 nrjTransf11, nrjTransf12, 966 nrjTransf11, nrjTransf12, 966 nrjTransf21, nrjTransf22, 967 nrjTransf21, nrjTransf22, 967 valueK1, valueK2, 968 valueK1, valueK2, 968 k, random); 969 k, random); 969 } 970 } 970 971 971 //G4cout << nrj << endl; 972 //G4cout << nrj << endl; 972 973 973 return nrj; 974 return nrj; 974 } 975 } 975 976