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