Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // Authors: S. Meylan and C. Villagrasa (IRSN, 27 // Models come from 28 // M. Bug et al, Rad. Phys and Chem. 130, 459- 29 // 30 31 #include "G4DNAPTBExcitationModel.hh" 32 33 #include "G4DNAChemistryManager.hh" 34 #include "G4DNAMaterialManager.hh" 35 #include "G4DNAMolecularMaterial.hh" 36 #include "G4SystemOfUnits.hh" 37 38 G4DNAPTBExcitationModel::G4DNAPTBExcitationMod 39 40 : G4VDNAModel(nam, applyToMaterial) 41 { 42 fpTHF = G4Material::GetMaterial("THF", false 43 fpPY = G4Material::GetMaterial("PY", false); 44 fpPU = G4Material::GetMaterial("PU", false); 45 fpTMP = G4Material::GetMaterial("TMP", false 46 fpG4_WATER = G4Material::GetMaterial("G4_WAT 47 fpBackbone_THF = G4Material::GetMaterial("ba 48 fpCytosine_PY = G4Material::GetMaterial("cyt 49 fpThymine_PY = G4Material::GetMaterial("thym 50 fpAdenine_PU = G4Material::GetMaterial("aden 51 fpBackbone_TMP = G4Material::GetMaterial("ba 52 fpGuanine_PU = G4Material::GetMaterial("guan 53 fpN2 = G4Material::GetMaterial("N2", false); 54 // initialisation of mean energy loss for ea 55 56 if (fpTHF != nullptr) { 57 fTableMeanEnergyPTB[fpTHF->GetIndex()] = 8 58 } 59 60 if (fpPY != nullptr) { 61 fTableMeanEnergyPTB[fpPY->GetIndex()] = 7. 62 } 63 64 if (fpPU != nullptr) { 65 fTableMeanEnergyPTB[fpPU->GetIndex()] = 7. 66 } 67 if (fpTMP != nullptr) { 68 fTableMeanEnergyPTB[fpTMP->GetIndex()] = 8 69 } 70 71 if (verboseLevel > 0) { 72 G4cout << "PTB excitation model is constru 73 } 74 } 75 76 //....oooOO0OOooo........oooOO0OOooo........oo 77 78 void G4DNAPTBExcitationModel::Initialise(const 79 const 80 { 81 if (isInitialised) { 82 return; 83 } 84 if (verboseLevel > 3) 85 { 86 G4cout << "Calling G4DNAPTBExcitationModel 87 } 88 89 if (particle != G4Electron::ElectronDefiniti 90 std::ostringstream oss; 91 oss << " Model is not applied for this par 92 G4Exception("G4DNAPTBExcitationModel::Init 93 } 94 95 G4double scaleFactor = 1e-16 * cm * cm; 96 G4double scaleFactorBorn = (1.e-22 / 3.343) 97 98 //****************************************** 99 // Cross section data 100 //****************************************** 101 std::size_t index; 102 if (fpTHF != nullptr) { 103 index = fpTHF->GetIndex(); 104 AddCrossSectionData(index, particle, "dna/ 105 SetLowELimit(index, particle, 9. * eV); 106 SetHighELimit(index, particle, 1. * keV); 107 } 108 if (fpPY != nullptr) { 109 index = fpPY->GetIndex(); 110 AddCrossSectionData(index, particle, "dna/ 111 SetLowELimit(index, particle, 9. * eV); 112 SetHighELimit(index, particle, 1. * keV); 113 } 114 115 if (fpPU != nullptr) { 116 index = fpPU->GetIndex(); 117 AddCrossSectionData(index, particle, "dna/ 118 SetLowELimit(index, particle, 9. * eV); 119 SetHighELimit(index, particle, 1. * keV); 120 } 121 122 if (fpTMP != nullptr) { 123 index = fpTMP->GetIndex(); 124 AddCrossSectionData(index, particle, "dna/ 125 SetLowELimit(index, particle, 9. * eV); 126 SetHighELimit(index, particle, 1. * keV); 127 } 128 if (fpG4_WATER != nullptr) { 129 index = fpG4_WATER->GetIndex(); 130 AddCrossSectionData(index, particle, "dna/ 131 SetLowELimit(index, particle, 9. * eV); 132 SetHighELimit(index, particle, 1. * keV); 133 } 134 // DNA materials 135 // 136 if (fpBackbone_THF != nullptr) { 137 index = fpBackbone_THF->GetIndex(); 138 AddCrossSectionData(index, particle, "dna/ 139 SetLowELimit(index, particle, 9. * eV); 140 SetHighELimit(index, particle, 1. * keV); 141 } 142 if (fpCytosine_PY != nullptr) { 143 index = fpCytosine_PY->GetIndex(); 144 AddCrossSectionData(index, particle, "dna/ 145 SetLowELimit(index, particle, 9. * eV); 146 SetHighELimit(index, particle, 1. * keV); 147 } 148 if (fpThymine_PY != nullptr) { 149 index = fpThymine_PY->GetIndex(); 150 AddCrossSectionData(index, particle, "dna/ 151 SetLowELimit(index, particle, 9. * eV); 152 SetHighELimit(index, particle, 1. * keV); 153 } 154 if (fpAdenine_PU != nullptr) { 155 index = fpAdenine_PU->GetIndex(); 156 AddCrossSectionData(index, particle, "dna/ 157 SetLowELimit(index, particle, 9. * eV); 158 SetHighELimit(index, particle, 1. * keV); 159 } 160 if (fpGuanine_PU != nullptr) { 161 index = fpGuanine_PU->GetIndex(); 162 AddCrossSectionData(index, particle, "dna/ 163 SetLowELimit(index, particle, 9. * eV); 164 SetHighELimit(index, particle, 1. * keV); 165 } 166 if (fpBackbone_TMP != nullptr) { 167 index = fpBackbone_TMP->GetIndex(); 168 AddCrossSectionData(index, particle, "dna/ 169 SetLowELimit(index, particle, 9. * eV); 170 SetHighELimit(index, particle, 1. * keV); 171 } 172 // MPietrzak, adding paths for N2 173 if (fpN2 != nullptr) { 174 index = fpN2->GetIndex(); 175 AddCrossSectionData(index, particle, "dna/ 176 SetLowELimit(index, particle, 13. * eV); 177 SetHighELimit(index, particle, 1.02 * MeV) 178 } 179 if (!G4DNAMaterialManager::Instance()->IsLoc 180 // Load data 181 LoadCrossSectionData(particle); 182 G4DNAMaterialManager::Instance()->SetMaste 183 fpModelData = this; 184 } 185 else { 186 auto dataModel = dynamic_cast<G4DNAPTBExci 187 G4DNAMaterialManager::Instance()->GetMod 188 if (dataModel == nullptr) { 189 G4cout << "G4DNAPTBExcitationModel::Init 190 G4Exception("G4DNAPTBExcitationModel::In 191 "not good modelData"); 192 } 193 else { 194 fpModelData = dataModel; 195 } 196 } 197 198 fParticleChangeForGamma = GetParticleChangeF 199 isInitialised = true; 200 } 201 202 //....oooOO0OOooo........oooOO0OOooo........oo 203 204 G4double G4DNAPTBExcitationModel::CrossSection 205 206 207 208 { 209 // Get the name of the current particle 210 G4String particleName = p->GetParticleName() 211 const std::size_t& MatID = material->GetInde 212 // initialise variables 213 G4double lowLim; 214 G4double highLim; 215 G4double sigma = 0; 216 217 // Get the low energy limit for the current 218 lowLim = fpModelData->GetLowELimit(MatID, p) 219 220 // Get the high energy limit for the current 221 highLim = fpModelData->GetHighELimit(MatID, 222 223 // Check that we are in the correct energy r 224 if (ekin >= lowLim && ekin < highLim) { 225 // Get the map with all the data tables 226 auto Data = fpModelData->GetData(); 227 228 if ((*Data)[MatID][p] == nullptr) { 229 G4Exception("G4DNAPTBExcitationModel::Cr 230 "No model is registered"); 231 } 232 // Retrieve the cross section value 233 sigma = (*Data)[MatID][p]->FindValue(ekin) 234 235 if (verboseLevel > 2) { 236 G4cout << "_____________________________ 237 G4cout << "°°° G4DNAPTBExcitationMode 238 G4cout << "°°° Kinetic energy(eV)=" < 239 G4cout << "°°° Cross section per " << 240 << G4endl; 241 G4cout << "°°° G4DNAPTBExcitationMode 242 } 243 } 244 245 // Return the cross section value 246 auto MolDensity = 247 (*G4DNAMolecularMaterial::Instance()->GetN 248 return sigma * MolDensity; 249 } 250 251 //....oooOO0OOooo........oooOO0OOooo........oo 252 253 void G4DNAPTBExcitationModel::SampleSecondarie 254 255 256 257 { 258 const std::size_t& materialID = (std::size_t 259 260 // Get the incident particle kinetic energy 261 G4double k = aDynamicParticle->GetKineticEne 262 // Get the particle name 263 const auto& particle = aDynamicParticle->Get 264 // Get the energy limits 265 G4double lowLim = fpModelData->GetLowELimit( 266 G4double highLim = fpModelData->GetHighELimi 267 268 // Check if we are in the correct energy ran 269 if (k >= lowLim && k < highLim) { 270 if (fpN2 != nullptr && materialID == fpN2- 271 // Retrieve the excitation energy for th 272 G4int level = fpModelData->RandomSelectS 273 G4double excitationEnergy = ptbExcitatio 274 275 // Calculate the new energy of the parti 276 G4double newEnergy = k - excitationEnerg 277 278 // Check that the new energy is above ze 279 // Otherwise, do nothing. 280 if (newEnergy > 0) { 281 fParticleChangeForGamma->ProposeMoment 282 fParticleChangeForGamma->SetProposedKi 283 fParticleChangeForGamma->ProposeLocalE 284 G4double ioniThres = ptbIonisationStru 285 // if excitation energy greater than i 286 if ((excitationEnergy > ioniThres) && 287 fParticleChangeForGamma->ProposeLoca 288 // energy of ejected electron 289 G4double secondaryKinetic = excitati 290 // random direction 291 G4double cosTheta = 2 * G4UniformRan 292 G4double sinTheta = std::sqrt(1. - c 293 G4double ux = sinTheta * std::cos(ph 294 G4ThreeVector deltaDirection(ux, uy, 295 // Create the new particle with its 296 auto dp = new G4DynamicParticle(G4El 297 fvect->push_back(dp); 298 } 299 } 300 else { 301 G4ExceptionDescription description; 302 description << "Kinetic energy <= 0 at 303 G4Exception("G4DNAPTBExcitationModel:: 304 } 305 } 306 else if (fpG4_WATER == nullptr || material 307 // Retrieve the excitation energy for th 308 G4double excitationEnergy = fTableMeanEn 309 // Calculate the new energy of the parti 310 G4double newEnergy = k - excitationEnerg 311 // Check that the new energy is above ze 312 // Otherwise, do nothing. 313 if (newEnergy > 0) { 314 fParticleChangeForGamma->ProposeMoment 315 fParticleChangeForGamma->SetProposedKi 316 fParticleChangeForGamma->ProposeLocalE 317 } 318 else { 319 G4ExceptionDescription description; 320 description << "Kinetic energy <= 0 at 321 G4Exception("G4DNAPTBExcitationModel:: 322 } 323 } 324 else { 325 G4int level = RandomSelectShell(k, parti 326 G4double excitationEnergy = waterStructu 327 G4double newEnergy = k - excitationEnerg 328 329 if (newEnergy > 0) { 330 fParticleChangeForGamma->ProposeMoment 331 fParticleChangeForGamma->SetProposedKi 332 fParticleChangeForGamma->ProposeLocalE 333 const G4Track* theIncomingTrack = fPar 334 G4DNAChemistryManager::Instance()->Cre 335 336 } 337 else { 338 G4ExceptionDescription description; 339 description << "Kinetic energy <= 0 at 340 G4Exception("G4DNAPTBExcitationModel:: 341 } 342 } 343 } 344 } 345