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 // G4DNATripleIonisationModel.cc 27 // 28 // Created at 2024/04/03 (Thu.) 29 // Author: Shogo OKADA @KEK-CRC (shogo.okada@ 30 // 31 // Reference: J.Meesungnoen et. al, DOI: 10.1 32 // 33 34 #include "G4DNATripleIonisationModel.hh" 35 36 #include "G4PhysicalConstants.hh" 37 #include "G4SystemOfUnits.hh" 38 #include "G4UAtomicDeexcitation.hh" 39 #include "G4LossTableManager.hh" 40 41 #include "G4SystemOfUnits.hh" 42 #include "G4DNAMolecularMaterial.hh" 43 #include "G4IonTable.hh" 44 #include "G4GenericIon.hh" 45 #include "G4DNARuddAngle.hh" 46 47 #include <sstream> 48 49 namespace { 50 51 G4DNAWaterIonisationStructure water_structure; 52 53 } // end of anonymous namespace 54 55 //============================================ 56 57 // constructor 58 G4DNATripleIonisationModel::G4DNATripleIonisat 59 const G4ParticleDefinition* p, const G4Strin 60 : G4DNADoubleIonisationModel(p, model_name 61 { 62 // Triple-ionisation energy 63 energy_threshold_ = 65.0 * eV; 64 } 65 66 //-------------------------------------------- 67 void G4DNATripleIonisationModel::Initialise( 68 const G4ParticleDefinition* particle, const 69 { 70 if (verbose_level_ > 3) { 71 G4cout << "Calling G4DNATripleIonisationMo 72 } 73 74 proton_def_ = G4Proton::ProtonDefinition(); 75 alpha_def_ = G4DNAGenericIonsManager::Insta 76 carbon_def_ = G4IonTable::GetIonTable()->Get 77 78 constexpr G4double kScaleFactor = 1.0 * m * 79 80 mioni_manager_ = new G4DNAMultipleIonisation 81 82 G4double Z{0.0}, A{0.0}; 83 G4String alpha_param_file{"dna/multipleionis 84 85 if (particle == proton_def_) { 86 87 // *************************************** 88 // for protons 89 const auto& proton = proton_def_->GetParti 90 elow_tab_[proton] = model_elow_tab_[1]; 91 eupp_tab_[proton] = 3.0 * MeV; 92 93 // load cross-section data for single ioni 94 auto xs_proton = new G4DNACrossSectionData 95 new G4LogLogInterpol 96 xs_proton->LoadData("dna/sigma_ionisation_ 97 xs_tab_[proton] = xs_proton; 98 99 // set energy limits 100 SetLowEnergyLimit(elow_tab_[proton]); 101 SetHighEnergyLimit(eupp_tab_[proton]); 102 103 if (!use_champion_param_) { 104 alpha_param_file = "dna/multipleionisati 105 } 106 107 Z = static_cast<G4double>(proton_def_->Get 108 A = static_cast<G4double>(proton_def_->Get 109 110 } else if (particle == alpha_def_) { 111 112 //**************************************** 113 // for alpha particles 114 const auto& alpha = alpha_def_->GetParticl 115 elow_tab_[alpha] = model_elow_tab_[4]; 116 eupp_tab_[alpha] = 23.0 * MeV; 117 118 // load cross-section data for single ioni 119 auto xs_alpha = new G4DNACrossSectionDataS 120 new G4LogLogInterpolat 121 xs_alpha->LoadData("dna/sigma_ionisation_a 122 xs_tab_[alpha] = xs_alpha; 123 124 // set energy limits 125 SetLowEnergyLimit(elow_tab_[alpha]); 126 SetHighEnergyLimit(eupp_tab_[alpha]); 127 128 if (!use_champion_param_) { 129 alpha_param_file = "dna/multipleionisati 130 } 131 132 Z = static_cast<G4double>(alpha_def_->GetA 133 A = static_cast<G4double>(alpha_def_->GetA 134 135 } else if (particle == G4GenericIon::Generic 136 137 // *************************************** 138 // for carbon ions 139 const auto& carbon = carbon_def_->GetParti 140 elow_tab_[carbon] = model_elow_tab_[5] * c 141 eupp_tab_[carbon] = 120.0 * MeV; 142 143 // load cross-section data for single ioni 144 auto xs_carbon = new G4DNACrossSectionData 145 new G4LogLogInterpolat 146 xs_carbon->LoadData("dna/sigma_ionisation_ 147 xs_tab_[carbon] = xs_carbon; 148 149 // set energy limits 150 SetLowEnergyLimit(elow_tab_[carbon]); 151 SetHighEnergyLimit(eupp_tab_[carbon]); 152 153 if (!use_champion_param_) { 154 alpha_param_file = "dna/multipleionisati 155 } 156 157 Z = static_cast<G4double>(carbon_def_->Get 158 A = static_cast<G4double>(carbon_def_->Get 159 160 } 161 162 // load alpha parameter 163 mioni_manager_->LoadAlphaParam(alpha_param_f 164 165 if (verbose_level_ > 0) { 166 G4cout << "G4DNATripleIonisationModel is i 167 << "Energy range: " 168 << LowEnergyLimit() / eV << " eV - 169 << HighEnergyLimit() / keV << " keV 170 << particle->GetParticleName() 171 << G4endl; 172 } 173 174 water_density_ = G4DNAMolecularMaterial::Ins 175 G4Material::GetMaterial( 176 177 atom_deex_ = G4LossTableManager::Instance()- 178 179 if (is_initialized_) { return; } 180 181 particle_change_ = GetParticleChangeForGamma 182 is_initialized_ = true; 183 } 184 185 //-------------------------------------------- 186 G4double G4DNATripleIonisationModel::CrossSect 187 const G4Material* material, const G4Particle 188 G4double ekin, G4double, G4double) 189 { 190 191 if (verbose_level_ > 3) { 192 G4cout << "Calling G4DNATripleIonisationMo 193 << G4endl; 194 } 195 196 // Calculate total cross section for model 197 198 if (pdef != proton_def_ && pdef != alpha_def 199 return 0.0; 200 } 201 202 static G4double water_dens = (*water_density 203 204 const auto& pname = pdef->GetParticleName(); 205 206 const auto low_energy_lim = GetLowEnergyLimi 207 const auto upp_energy_lim = GetUppEnergyLimi 208 209 G4double sigma{0.0}; 210 if (ekin <= upp_energy_lim) { 211 212 if (ekin < low_energy_lim) { ekin = low_en 213 214 CrossSectionDataTable::iterator pos = xs_t 215 if (pos == xs_tab_.end()) { 216 G4Exception("G4DNATripleIonisationModel: 217 "em0002", FatalException, 218 "Model not applicable to par 219 } 220 221 G4DNACrossSectionDataSet* table = pos->sec 222 if (table != nullptr) { 223 const auto a = mioni_manager_->GetAlphaP 224 sigma = table->FindValue(ekin) * a * a; 225 } 226 227 } 228 229 if (verbose_level_ > 2) { 230 231 std::stringstream msg; 232 233 msg << "---------------------------------- 234 msg << " G4DNATripleIonisationModel - XS I 235 msg << " - Kinetic energy(eV): " << ekin/ 236 << pdef->GetParticleName() << "\n"; 237 msg << " - Cross section per water molecu 238 << sigma / cm / cm << "\n"; 239 msg << " - Cross section per water molecu 240 << sigma * water_dens / (1.0 / cm) << 241 msg << " G4DNATripleIonisationModel - XS I 242 msg << "---------------------------------- 243 244 G4cout << msg.str() << G4endl; 245 246 } 247 248 return (sigma * water_dens); 249 250 } 251 252 //-------------------------------------------- 253 void G4DNATripleIonisationModel::SampleSeconda 254 std::vector<G4DynamicParticle*>* vsec, const 255 const G4DynamicParticle* particle, G4double, 256 { 257 258 if (verbose_level_ > 3) { 259 G4cout << "Calling SampleSecondaries() of 260 << G4endl; 261 } 262 263 // get the definition for this parent partic 264 auto pdef = particle->GetDefinition(); 265 266 // get kinetic energy 267 auto ekin = particle->GetKineticEnergy(); 268 269 // get particle name 270 const auto& pname = pdef->GetParticleName(); 271 272 // get energy limits 273 const auto low_energy_lim = GetLowEnergyLimi 274 275 // ***************************************** 276 // stop the transportation process of this p 277 // if its kinetic energy is below the lower 278 if (ekin < low_energy_lim) { 279 particle_change_->SetProposedKineticEnergy 280 particle_change_->ProposeTrackStatus(fStop 281 particle_change_->ProposeLocalEnergyDeposi 282 return; 283 } 284 // ***************************************** 285 286 constexpr G4int kNumSecondaries = 3; 287 constexpr G4double kDeltaTheta = pi * 0.666 288 289 G4int ioni_shell[kNumSecondaries] = {0, 0, 0 290 G4double shell_energy[kNumSecondaries]; 291 292 auto scale_param = mioni_manager_->GetAlphaP 293 scale_param *= scale_param; 294 295 G4bool is_continue{true}; 296 while (1) { 297 ioni_shell[0] = RandomSelect(ekin, scale_p 298 ioni_shell[1] = RandomSelect(ekin, scale_p 299 ioni_shell[2] = RandomSelect(ekin, scale_p 300 is_continue = (ioni_shell[0] == ioni_shell 301 ioni_shell[1] == ioni_shell 302 if (!is_continue) { break; } 303 } 304 305 G4double tot_ioni_energy{0.0}; 306 for (int i = 0; i < kNumSecondaries; i++) { 307 shell_energy[i] = ::water_structure.Ionisa 308 tot_ioni_energy += shell_energy[i]; 309 } 310 311 if (ekin < tot_ioni_energy || tot_ioni_energ 312 return; 313 } 314 315 // generate secondary electrons 316 G4double theta{0.0}, phi{0.0}, tot_ekin2{0.0 317 for (int i = 0; i < kNumSecondaries; i++) { 318 tot_ekin2 += GenerateSecondaries(vsec, cou 319 theta, ph 320 theta += kDeltaTheta; 321 } 322 323 // This should never happen 324 if (mioni_manager_->CheckShellEnergy(eTriple 325 G4Exception("G4DNATripleIonisatioModel::Sa 326 "em2050", FatalException, "Negative lo 327 } 328 329 // ***************************************** 330 // update kinematics for this parent particl 331 const auto primary_dir = particle->GetMoment 332 particle_change_->ProposeMomentumDirection(p 333 334 const auto scattered_energy = ekin - tot_ion 335 336 // update total amount of shell energy 337 tot_ioni_energy = shell_energy[0] + shell_en 338 339 if (stat_code_) { 340 particle_change_->SetProposedKineticEnergy 341 particle_change_->ProposeLocalEnergyDeposi 342 ekin - scatter 343 } else { 344 particle_change_->SetProposedKineticEnergy 345 particle_change_->ProposeLocalEnergyDeposi 346 } 347 348 // ***************************************** 349 // generate triple-ionized water molecules ( 350 const auto the_track = particle_change_->Get 351 mioni_manager_->CreateMultipleIonisedWaterMo 352 eTripleIonisedMolecule, io 353 // ***************************************** 354 355 } 356