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 // $Id: G4DNARuddIonisationExtendedModel.cc 96060 2016-03-11 12:58:04Z gcosmo $ >> 27 // GEANT4 tag $Name: $ 26 // 28 // 27 // Modified by Z. Francis, S. Incerti to handl 29 // Modified by Z. Francis, S. Incerti to handle HZE 28 // && inverse rudd function sampling 26-10-201 30 // && inverse rudd function sampling 26-10-2010 29 // << 30 // Rewitten by V.Ivanchenko 21.05.2023 << 31 // << 32 31 33 #include "G4EmCorrections.hh" << 34 #include "G4DNARuddIonisationExtendedModel.hh" 32 #include "G4DNARuddIonisationExtendedModel.hh" 35 #include "G4PhysicalConstants.hh" 33 #include "G4PhysicalConstants.hh" 36 #include "G4SystemOfUnits.hh" 34 #include "G4SystemOfUnits.hh" 37 #include "G4UAtomicDeexcitation.hh" 35 #include "G4UAtomicDeexcitation.hh" 38 #include "G4LossTableManager.hh" 36 #include "G4LossTableManager.hh" 39 #include "G4NistManager.hh" << 40 #include "G4DNAChemistryManager.hh" 37 #include "G4DNAChemistryManager.hh" 41 #include "G4DNAMolecularMaterial.hh" 38 #include "G4DNAMolecularMaterial.hh" 42 39 >> 40 //SEB 43 #include "G4IonTable.hh" 41 #include "G4IonTable.hh" 44 #include "G4DNARuddAngle.hh" 42 #include "G4DNARuddAngle.hh" 45 #include "G4DeltaAngle.hh" 43 #include "G4DeltaAngle.hh" 46 #include "G4Exp.hh" 44 #include "G4Exp.hh" 47 #include "G4Log.hh" << 45 // 48 #include "G4Pow.hh" << 49 #include "G4Alpha.hh" << 50 #include "G4Proton.hh" << 51 46 52 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 53 48 54 G4DNACrossSectionDataSet* G4DNARuddIonisationE << 49 using namespace std; 55 G4DNACrossSectionDataSet* G4DNARuddIonisationE << 56 G4DNACrossSectionDataSet* G4DNARuddIonisationE << 57 const std::vector<G4double>* G4DNARuddIonisati << 58 << 59 namespace << 60 { << 61 const G4double scaleFactor = CLHEP::m*CLHEP: << 62 const G4double tolerance = 1*CLHEP::eV; << 63 const G4double Ry = 13.6*CLHEP::eV; << 64 << 65 // Following values provided by M. Dingfelde << 66 const G4double Bj[5] = {12.60*CLHEP::eV, 14. << 67 32.20*CLHEP::eV, 539 << 68 } << 69 50 70 //....oooOO0OOooo........oooOO0OOooo........oo 51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 71 52 72 G4DNARuddIonisationExtendedModel::G4DNARuddIon 53 G4DNARuddIonisationExtendedModel::G4DNARuddIonisationExtendedModel(const G4ParticleDefinition*, 73 54 const G4String& nam) 74 : G4VEmModel(nam) << 55 :G4VEmModel(nam),isInitialised(false) 75 { 56 { 76 fEmCorrections = G4LossTableManager::Instanc << 57 // nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER"); 77 fGpow = G4Pow::GetInstance(); << 58 fpWaterDensity = 0; 78 fLowestEnergy = 100*CLHEP::eV; << 59 79 fLimitEnergy = 1*CLHEP::keV; << 60 slaterEffectiveCharge[0]=0.; >> 61 slaterEffectiveCharge[1]=0.; >> 62 slaterEffectiveCharge[2]=0.; >> 63 sCoefficient[0]=0.; >> 64 sCoefficient[1]=0.; >> 65 sCoefficient[2]=0.; >> 66 >> 67 lowEnergyLimitForA[1] = 0 * eV; >> 68 lowEnergyLimitForA[2] = 0 * eV; >> 69 lowEnergyLimitForA[3] = 0 * eV; >> 70 lowEnergyLimitOfModelForA[1] = 100 * eV; >> 71 lowEnergyLimitOfModelForA[4] = 1 * keV; >> 72 lowEnergyLimitOfModelForA[5] = 0.5 * MeV; // For A = 3 or above, limit is MeV/uma >> 73 killBelowEnergyForA[1] = lowEnergyLimitOfModelForA[1]; >> 74 killBelowEnergyForA[4] = lowEnergyLimitOfModelForA[4]; >> 75 killBelowEnergyForA[5] = lowEnergyLimitOfModelForA[5]; >> 76 >> 77 verboseLevel= 0; >> 78 // Verbosity scale: >> 79 // 0 = nothing >> 80 // 1 = warning for energy non-conservation >> 81 // 2 = details of energy budget >> 82 // 3 = calculation of cross sections, file openings, sampling of atoms >> 83 // 4 = entering in methods >> 84 >> 85 if( verboseLevel>0 ) >> 86 { >> 87 G4cout << "Rudd ionisation model is constructed " << G4endl; >> 88 } 80 89 81 // Mark this model as "applicable" for atomi << 90 // Define default angular generator 82 SetDeexcitationFlag(true); << 91 SetAngularDistribution(new G4DNARuddAngle()); 83 92 84 // Define default angular generator << 93 // Mark this model as "applicable" for atomic deexcitation 85 SetAngularDistribution(new G4DNARuddAngle()) << 94 SetDeexcitationFlag(true); >> 95 fAtomDeexcitation = 0; >> 96 fParticleChangeForGamma = 0; 86 97 87 if (nullptr == xshelium) { LoadData(); } << 98 // Selection of stationary mode >> 99 >> 100 statCode = false; 88 } 101 } 89 102 90 //....oooOO0OOooo........oooOO0OOooo........oo 103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 91 104 92 G4DNARuddIonisationExtendedModel::~G4DNARuddIo 105 G4DNARuddIonisationExtendedModel::~G4DNARuddIonisationExtendedModel() 93 { 106 { 94 if(isFirst) { << 107 // Cross section 95 for(auto & i : xsdata) { delete i; } << 108 >> 109 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos; >> 110 for (pos = tableData.begin(); pos != tableData.end(); ++pos) >> 111 { >> 112 G4DNACrossSectionDataSet* table = pos->second; >> 113 delete table; >> 114 } >> 115 >> 116 // The following removal is forbidden G4VEnergyLossModel takes care of deletion >> 117 // however coverity will signal this as an error >> 118 // if (fAtomDeexcitation) {delete fAtomDeexcitation;} >> 119 >> 120 } >> 121 >> 122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 123 >> 124 void G4DNARuddIonisationExtendedModel::Initialise(const G4ParticleDefinition* particle, >> 125 const G4DataVector& /*cuts*/) >> 126 { >> 127 if (verboseLevel > 3) >> 128 G4cout << "Calling G4DNARuddIonisationExtendedModel::Initialise()" << G4endl; >> 129 >> 130 // Energy limits >> 131 >> 132 G4String fileProton("dna/sigma_ionisation_p_rudd"); >> 133 G4String fileHydrogen("dna/sigma_ionisation_h_rudd"); >> 134 G4String fileAlphaPlusPlus("dna/sigma_ionisation_alphaplusplus_rudd"); >> 135 G4String fileAlphaPlus("dna/sigma_ionisation_alphaplus_rudd"); >> 136 G4String fileHelium("dna/sigma_ionisation_he_rudd"); >> 137 G4String fileLithium("dna/sigma_ionisation_li_rudd"); >> 138 G4String fileBeryllium("dna/sigma_ionisation_be_rudd"); >> 139 G4String fileBoron("dna/sigma_ionisation_b_rudd"); >> 140 G4String fileCarbon("dna/sigma_ionisation_c_rudd"); >> 141 G4String fileNitrogen("dna/sigma_ionisation_n_rudd"); >> 142 G4String fileOxygen("dna/sigma_ionisation_o_rudd"); >> 143 G4String fileSilicon("dna/sigma_ionisation_si_rudd"); >> 144 G4String fileIron("dna/sigma_ionisation_fe_rudd"); >> 145 >> 146 G4DNAGenericIonsManager *instance; >> 147 instance = G4DNAGenericIonsManager::Instance(); >> 148 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition(); >> 149 G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen"); >> 150 G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++"); >> 151 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+"); >> 152 G4ParticleDefinition* heliumDef = instance->GetIon("helium"); >> 153 >> 154 //SEB >> 155 //G4ParticleDefinition* carbonDef = instance->GetIon("carbon"); >> 156 //G4ParticleDefinition* nitrogenDef = instance->GetIon("nitrogen"); >> 157 //G4ParticleDefinition* oxygenDef = instance->GetIon("oxygen"); >> 158 //G4ParticleDefinition* siliconDef = instance->GetIon("silicon"); >> 159 //G4ParticleDefinition* ironDef = instance->GetIon("iron"); >> 160 G4ParticleDefinition* lithiumDef = G4IonTable::GetIonTable()->GetIon(3,7); >> 161 G4ParticleDefinition* berylliumDef = G4IonTable::GetIonTable()->GetIon(4,9); >> 162 G4ParticleDefinition* boronDef = G4IonTable::GetIonTable()->GetIon(5,11); >> 163 G4ParticleDefinition* carbonDef = G4IonTable::GetIonTable()->GetIon(6,12); >> 164 G4ParticleDefinition* nitrogenDef = G4IonTable::GetIonTable()->GetIon(7,14); >> 165 G4ParticleDefinition* oxygenDef = G4IonTable::GetIonTable()->GetIon(8,16); >> 166 G4ParticleDefinition* siliconDef = G4IonTable::GetIonTable()->GetIon(14,28); >> 167 G4ParticleDefinition* ironDef = G4IonTable::GetIonTable()->GetIon(26,56); >> 168 // >> 169 >> 170 G4String proton; >> 171 G4String hydrogen; >> 172 G4String alphaPlusPlus; >> 173 G4String alphaPlus; >> 174 G4String helium; >> 175 G4String lithium; >> 176 G4String beryllium; >> 177 G4String boron; >> 178 G4String carbon; >> 179 G4String nitrogen; >> 180 G4String oxygen; >> 181 G4String silicon; >> 182 G4String iron; >> 183 >> 184 G4double scaleFactor = 1 * m*m; >> 185 >> 186 // LIMITS AND DATA >> 187 >> 188 // ********************************************************************************************** >> 189 >> 190 proton = protonDef->GetParticleName(); >> 191 tableFile[proton] = fileProton; >> 192 lowEnergyLimit[proton] = lowEnergyLimitForA[1]; >> 193 highEnergyLimit[proton] = 500. * keV; >> 194 >> 195 // Cross section >> 196 >> 197 G4DNACrossSectionDataSet* tableProton = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, >> 198 eV, >> 199 scaleFactor ); >> 200 tableProton->LoadData(fileProton); >> 201 tableData[proton] = tableProton; >> 202 >> 203 // ********************************************************************************************** >> 204 >> 205 hydrogen = hydrogenDef->GetParticleName(); >> 206 tableFile[hydrogen] = fileHydrogen; >> 207 >> 208 lowEnergyLimit[hydrogen] = lowEnergyLimitForA[1]; >> 209 highEnergyLimit[hydrogen] = 100. * MeV; >> 210 >> 211 // Cross section >> 212 >> 213 G4DNACrossSectionDataSet* tableHydrogen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, >> 214 eV, >> 215 scaleFactor ); >> 216 tableHydrogen->LoadData(fileHydrogen); >> 217 >> 218 tableData[hydrogen] = tableHydrogen; >> 219 >> 220 // ********************************************************************************************** >> 221 >> 222 alphaPlusPlus = alphaPlusPlusDef->GetParticleName(); >> 223 tableFile[alphaPlusPlus] = fileAlphaPlusPlus; >> 224 >> 225 lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForA[4]; >> 226 highEnergyLimit[alphaPlusPlus] = 400. * MeV; >> 227 >> 228 // Cross section >> 229 >> 230 G4DNACrossSectionDataSet* tableAlphaPlusPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, >> 231 eV, >> 232 scaleFactor ); >> 233 tableAlphaPlusPlus->LoadData(fileAlphaPlusPlus); >> 234 >> 235 tableData[alphaPlusPlus] = tableAlphaPlusPlus; >> 236 >> 237 // ********************************************************************************************** >> 238 >> 239 alphaPlus = alphaPlusDef->GetParticleName(); >> 240 tableFile[alphaPlus] = fileAlphaPlus; >> 241 >> 242 lowEnergyLimit[alphaPlus] = lowEnergyLimitForA[4]; >> 243 highEnergyLimit[alphaPlus] = 400. * MeV; >> 244 >> 245 // Cross section >> 246 >> 247 G4DNACrossSectionDataSet* tableAlphaPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, >> 248 eV, >> 249 scaleFactor ); >> 250 tableAlphaPlus->LoadData(fileAlphaPlus); >> 251 tableData[alphaPlus] = tableAlphaPlus; >> 252 >> 253 // ********************************************************************************************** >> 254 >> 255 helium = heliumDef->GetParticleName(); >> 256 tableFile[helium] = fileHelium; >> 257 >> 258 lowEnergyLimit[helium] = lowEnergyLimitForA[4]; >> 259 highEnergyLimit[helium] = 400. * MeV; >> 260 >> 261 // Cross section >> 262 >> 263 G4DNACrossSectionDataSet* tableHelium = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, >> 264 eV, >> 265 scaleFactor ); >> 266 tableHelium->LoadData(fileHelium); >> 267 tableData[helium] = tableHelium; >> 268 >> 269 // ********************************************************************************************** >> 270 >> 271 lithium = lithiumDef->GetParticleName(); >> 272 tableFile[lithium] = fileLithium; >> 273 >> 274 //SEB >> 275 //lowEnergyLimit[carbon] = lowEnergyLimitForA[5] * particle->GetAtomicMass(); >> 276 //highEnergyLimit[carbon] = 1e6* particle->GetAtomicMass() * MeV; >> 277 lowEnergyLimit[lithium] = 0.5*7*MeV; >> 278 highEnergyLimit[lithium] = 1e6*7*MeV; >> 279 // >> 280 >> 281 // Cross section >> 282 >> 283 G4DNACrossSectionDataSet* tableLithium = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, >> 284 eV, >> 285 scaleFactor ); >> 286 tableLithium->LoadData(fileLithium); >> 287 tableData[lithium] = tableLithium; >> 288 >> 289 // ********************************************************************************************** >> 290 >> 291 beryllium = berylliumDef->GetParticleName(); >> 292 tableFile[beryllium] = fileBeryllium; >> 293 >> 294 //SEB >> 295 //lowEnergyLimit[carbon] = lowEnergyLimitForA[5] * particle->GetAtomicMass(); >> 296 //highEnergyLimit[carbon] = 1e6* particle->GetAtomicMass() * MeV; >> 297 lowEnergyLimit[beryllium] = 0.5*9*MeV; >> 298 highEnergyLimit[beryllium] = 1e6*9*MeV; >> 299 // >> 300 >> 301 // Cross section >> 302 >> 303 G4DNACrossSectionDataSet* tableBeryllium = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, >> 304 eV, >> 305 scaleFactor ); >> 306 tableBeryllium->LoadData(fileBeryllium); >> 307 tableData[beryllium] = tableBeryllium; >> 308 >> 309 // ********************************************************************************************** >> 310 >> 311 boron = boronDef->GetParticleName(); >> 312 tableFile[boron] = fileBoron; >> 313 >> 314 //SEB >> 315 //lowEnergyLimit[carbon] = lowEnergyLimitForA[5] * particle->GetAtomicMass(); >> 316 //highEnergyLimit[carbon] = 1e6* particle->GetAtomicMass() * MeV; >> 317 lowEnergyLimit[boron] = 0.5*11*MeV; >> 318 highEnergyLimit[boron] = 1e6*11*MeV; >> 319 // >> 320 >> 321 // Cross section >> 322 >> 323 G4DNACrossSectionDataSet* tableBoron = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, >> 324 eV, >> 325 scaleFactor ); >> 326 tableBoron->LoadData(fileBoron); >> 327 tableData[boron] = tableBoron; >> 328 >> 329 // ********************************************************************************************** >> 330 >> 331 carbon = carbonDef->GetParticleName(); >> 332 tableFile[carbon] = fileCarbon; >> 333 >> 334 //SEB >> 335 //lowEnergyLimit[carbon] = lowEnergyLimitForA[5] * particle->GetAtomicMass(); >> 336 //highEnergyLimit[carbon] = 1e6* particle->GetAtomicMass() * MeV; >> 337 lowEnergyLimit[carbon] = 0.5*12*MeV; >> 338 highEnergyLimit[carbon] = 1e6*12*MeV; >> 339 // >> 340 >> 341 // Cross section >> 342 >> 343 G4DNACrossSectionDataSet* tableCarbon = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, >> 344 eV, >> 345 scaleFactor ); >> 346 tableCarbon->LoadData(fileCarbon); >> 347 tableData[carbon] = tableCarbon; >> 348 >> 349 // ********************************************************************************************** >> 350 >> 351 oxygen = oxygenDef->GetParticleName(); >> 352 tableFile[oxygen] = fileOxygen; >> 353 >> 354 //SEB >> 355 //lowEnergyLimit[oxygen] = lowEnergyLimitForA[5]* particle->GetAtomicMass(); >> 356 //highEnergyLimit[oxygen] = 1e6* particle->GetAtomicMass()* MeV; >> 357 lowEnergyLimit[oxygen] = 0.5*16*MeV; >> 358 highEnergyLimit[oxygen] = 1e6*16*MeV; >> 359 // >> 360 >> 361 // Cross section >> 362 >> 363 G4DNACrossSectionDataSet* tableOxygen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, >> 364 eV, >> 365 scaleFactor ); >> 366 tableOxygen->LoadData(fileOxygen); >> 367 tableData[oxygen] = tableOxygen; >> 368 >> 369 // ********************************************************************************************** >> 370 >> 371 nitrogen = nitrogenDef->GetParticleName(); >> 372 tableFile[nitrogen] = fileNitrogen; >> 373 >> 374 //SEB >> 375 //lowEnergyLimit[nitrogen] = lowEnergyLimitForA[5]* particle->GetAtomicMass(); >> 376 //highEnergyLimit[nitrogen] = 1e6* particle->GetAtomicMass()* MeV; >> 377 lowEnergyLimit[nitrogen] = 0.5*14*MeV; >> 378 highEnergyLimit[nitrogen] = 1e6*14*MeV; >> 379 // >> 380 >> 381 // Cross section >> 382 >> 383 G4DNACrossSectionDataSet* tableNitrogen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, >> 384 eV, >> 385 scaleFactor ); >> 386 tableNitrogen->LoadData(fileNitrogen); >> 387 tableData[nitrogen] = tableNitrogen; >> 388 >> 389 // ********************************************************************************************** >> 390 >> 391 silicon = siliconDef->GetParticleName(); >> 392 tableFile[silicon] = fileSilicon; >> 393 >> 394 //lowEnergyLimit[silicon] = lowEnergyLimitForA[5]* particle->GetAtomicMass(); >> 395 //highEnergyLimit[silicon] = 1e6* particle->GetAtomicMass()* MeV; >> 396 lowEnergyLimit[silicon] = 0.5*28*MeV; >> 397 highEnergyLimit[silicon] = 1e6*28*MeV; >> 398 // >> 399 >> 400 // Cross section >> 401 >> 402 G4DNACrossSectionDataSet* tableSilicon = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, >> 403 eV, >> 404 scaleFactor ); >> 405 tableSilicon->LoadData(fileSilicon); >> 406 tableData[silicon] = tableSilicon; >> 407 >> 408 // ********************************************************************************************** >> 409 >> 410 iron = ironDef->GetParticleName(); >> 411 tableFile[iron] = fileIron; >> 412 >> 413 //SEB >> 414 //lowEnergyLimit[iron] = lowEnergyLimitForA[5]* particle->GetAtomicMass(); >> 415 //highEnergyLimit[iron] = 1e6* particle->GetAtomicMass()* MeV; >> 416 lowEnergyLimit[iron] = 0.5*56*MeV; >> 417 highEnergyLimit[iron] = 1e6*56*MeV; >> 418 // >> 419 >> 420 // Cross section >> 421 >> 422 G4DNACrossSectionDataSet* tableIron = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, >> 423 eV, >> 424 scaleFactor ); >> 425 tableIron->LoadData(fileIron); >> 426 tableData[iron] = tableIron; >> 427 >> 428 // ********************************************************************************************** >> 429 >> 430 //SEB: not anymore >> 431 // ZF Following lines can be replaced by: >> 432 //SetLowEnergyLimit(lowEnergyLimit[particle->GetParticleName()]); >> 433 //SetHighEnergyLimit(highEnergyLimit[particle->GetParticleName()]); >> 434 // at least for HZE >> 435 >> 436 if (particle==protonDef) >> 437 { >> 438 SetLowEnergyLimit(lowEnergyLimit[proton]); >> 439 SetHighEnergyLimit(highEnergyLimit[proton]); 96 } 440 } >> 441 >> 442 if (particle==hydrogenDef) >> 443 { >> 444 SetLowEnergyLimit(lowEnergyLimit[hydrogen]); >> 445 SetHighEnergyLimit(highEnergyLimit[hydrogen]); >> 446 } >> 447 >> 448 if (particle==heliumDef) >> 449 { >> 450 SetLowEnergyLimit(lowEnergyLimit[helium]); >> 451 SetHighEnergyLimit(highEnergyLimit[helium]); >> 452 } >> 453 >> 454 if (particle==alphaPlusDef) >> 455 { >> 456 SetLowEnergyLimit(lowEnergyLimit[alphaPlus]); >> 457 SetHighEnergyLimit(highEnergyLimit[alphaPlus]); >> 458 } >> 459 >> 460 if (particle==alphaPlusPlusDef) >> 461 { >> 462 SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]); >> 463 SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]); >> 464 } >> 465 >> 466 if (particle==lithiumDef) >> 467 { >> 468 SetLowEnergyLimit(lowEnergyLimit[lithium]); >> 469 SetHighEnergyLimit(highEnergyLimit[lithium]); >> 470 } >> 471 >> 472 if (particle==berylliumDef) >> 473 { >> 474 SetLowEnergyLimit(lowEnergyLimit[beryllium]); >> 475 SetHighEnergyLimit(highEnergyLimit[beryllium]); >> 476 } >> 477 >> 478 if (particle==boronDef) >> 479 { >> 480 SetLowEnergyLimit(lowEnergyLimit[boron]); >> 481 SetHighEnergyLimit(highEnergyLimit[boron]); >> 482 } >> 483 >> 484 if (particle==carbonDef) >> 485 { >> 486 SetLowEnergyLimit(lowEnergyLimit[carbon]); >> 487 SetHighEnergyLimit(highEnergyLimit[carbon]); >> 488 } >> 489 >> 490 if (particle==nitrogenDef) >> 491 { >> 492 SetLowEnergyLimit(lowEnergyLimit[nitrogen]); >> 493 SetHighEnergyLimit(highEnergyLimit[nitrogen]); >> 494 } >> 495 >> 496 if (particle==oxygenDef) >> 497 { >> 498 SetLowEnergyLimit(lowEnergyLimit[oxygen]); >> 499 SetHighEnergyLimit(highEnergyLimit[oxygen]); >> 500 } >> 501 >> 502 if (particle==siliconDef) >> 503 { >> 504 SetLowEnergyLimit(lowEnergyLimit[silicon]); >> 505 SetHighEnergyLimit(highEnergyLimit[silicon]); >> 506 } >> 507 >> 508 if (particle==ironDef) >> 509 { >> 510 SetLowEnergyLimit(lowEnergyLimit[iron]); >> 511 SetHighEnergyLimit(highEnergyLimit[iron]); >> 512 } >> 513 >> 514 //---------------------------------------------------------------------- >> 515 >> 516 if( verboseLevel>0 ) >> 517 { >> 518 G4cout << "Rudd ionisation model is initialized " << G4endl >> 519 << "Energy range: " >> 520 << LowEnergyLimit() / eV << " eV - " >> 521 << HighEnergyLimit() / keV << " keV for " >> 522 << particle->GetParticleName() >> 523 << G4endl; >> 524 } >> 525 >> 526 // Initialize water density pointer >> 527 fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER")); >> 528 >> 529 // >> 530 >> 531 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation(); >> 532 >> 533 if (isInitialised) { return; } >> 534 fParticleChangeForGamma = GetParticleChangeForGamma(); >> 535 isInitialised = true; 97 } 536 } 98 537 99 //....oooOO0OOooo........oooOO0OOooo........oo 538 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 100 539 101 void G4DNARuddIonisationExtendedModel::LoadDat << 540 G4double G4DNARuddIonisationExtendedModel::CrossSectionPerVolume(const G4Material* material, >> 541 const G4ParticleDefinition* particleDefinition, >> 542 G4double k, >> 543 G4double, >> 544 G4double) 102 { 545 { 103 // initialisation of static data once << 546 //SEB: particleDefinition->GetParticleName() is for eg. Fe56 104 isFirst = true; << 547 // particleDefinition->GetPDGMass() is correct 105 G4String filename("dna/sigma_ionisation_h_ru << 548 // particleDefinition->GetAtomicNumber() is correct 106 xsdata[0] = new G4DNACrossSectionDataSet(new << 549 107 xsdata[0]->LoadData(filename); << 550 if (verboseLevel > 3) 108 << 551 G4cout << "Calling CrossSectionPerVolume() of G4DNARuddIonisationExtendedModel" << G4endl; 109 filename = "dna/sigma_ionisation_p_rudd"; << 552 110 xsdata[1] = new G4DNACrossSectionDataSet(new << 553 // Calculate total cross section for model 111 xsdata[1]->LoadData(filename); << 554 112 << 555 G4DNAGenericIonsManager *instance; 113 filename = "dna/sigma_ionisation_alphapluspl << 556 instance = G4DNAGenericIonsManager::Instance(); 114 xsdata[2] = new G4DNACrossSectionDataSet(new << 557 115 xsdata[2]->LoadData(filename); << 558 if ( 116 << 559 particleDefinition != G4Proton::ProtonDefinition() 117 filename = "dna/sigma_ionisation_li_rudd"; << 560 && 118 xsdata[3] = new G4DNACrossSectionDataSet(new << 561 particleDefinition != instance->GetIon("hydrogen") 119 xsdata[3]->LoadData(filename); << 562 && 120 << 563 particleDefinition != instance->GetIon("alpha++") 121 filename = "dna/sigma_ionisation_be_rudd"; << 564 && 122 xsdata[4] = new G4DNACrossSectionDataSet(new << 565 particleDefinition != instance->GetIon("alpha+") 123 xsdata[4]->LoadData(filename); << 566 && 124 << 567 particleDefinition != instance->GetIon("helium") 125 filename = "dna/sigma_ionisation_b_rudd"; << 568 && 126 xsdata[5] = new G4DNACrossSectionDataSet(new << 569 //SEB 127 xsdata[5]->LoadData(filename); << 570 //particleDefinition != instance->GetIon("carbon") 128 << 571 //&& 129 filename = "dna/sigma_ionisation_c_rudd"; << 572 //particleDefinition != instance->GetIon("nitrogen") 130 xsdata[6] = new G4DNACrossSectionDataSet(new << 573 //&& 131 xsdata[6]->LoadData(filename); << 574 //particleDefinition != instance->GetIon("oxygen") 132 << 575 //&& 133 filename = "dna/sigma_ionisation_n_rudd"; << 576 //particleDefinition != instance->GetIon("iron") 134 xsdata[7] = new G4DNACrossSectionDataSet(new << 577 particleDefinition != G4IonTable::GetIonTable()->GetIon(3,7) 135 xsdata[7]->LoadData(filename); << 578 && 136 << 579 particleDefinition != G4IonTable::GetIonTable()->GetIon(4,9) 137 filename = "dna/sigma_ionisation_o_rudd"; << 580 && 138 xsdata[8] = new G4DNACrossSectionDataSet(new << 581 particleDefinition != G4IonTable::GetIonTable()->GetIon(5,11) 139 xsdata[8]->LoadData(filename); << 582 && 140 << 583 particleDefinition != G4IonTable::GetIonTable()->GetIon(6,12) 141 filename = "dna/sigma_ionisation_si_rudd"; << 584 && 142 xsdata[14] = new G4DNACrossSectionDataSet(ne << 585 particleDefinition != G4IonTable::GetIonTable()->GetIon(7,14) 143 xsdata[14]->LoadData(filename); << 586 && 144 << 587 particleDefinition != G4IonTable::GetIonTable()->GetIon(8,16) 145 filename = "dna/sigma_ionisation_fe_rudd"; << 588 && 146 xsdata[26] = new G4DNACrossSectionDataSet(ne << 589 particleDefinition != G4IonTable::GetIonTable()->GetIon(14,28) 147 xsdata[26]->LoadData(filename); << 590 && 148 << 591 particleDefinition != G4IonTable::GetIonTable()->GetIon(26,56) 149 filename = "dna/sigma_ionisation_alphaplus_r << 592 // 150 xsalphaplus = new G4DNACrossSectionDataSet(n << 593 ) 151 xsalphaplus->LoadData(filename); << 594 152 << 595 return 0; 153 filename = "dna/sigma_ionisation_he_rudd"; << 596 154 xshelium = new G4DNACrossSectionDataSet(new << 597 G4double lowLim = 0; 155 xshelium->LoadData(filename); << 598 156 << 599 if ( particleDefinition == G4Proton::ProtonDefinition() 157 // to avoid possible threading problem fill << 600 || particleDefinition == instance->GetIon("hydrogen") 158 auto water = G4NistManager::Instance()->Find << 601 ) 159 fpWaterDensity = << 602 160 G4DNAMolecularMaterial::Instance()->GetNum << 603 lowLim = lowEnergyLimitOfModelForA[1]; >> 604 >> 605 else if ( particleDefinition == instance->GetIon("alpha++") >> 606 || particleDefinition == instance->GetIon("alpha+") >> 607 || particleDefinition == instance->GetIon("helium") >> 608 ) >> 609 >> 610 lowLim = lowEnergyLimitOfModelForA[4]; >> 611 >> 612 else lowLim = lowEnergyLimitOfModelForA[5]; >> 613 >> 614 G4double highLim = 0; >> 615 G4double sigma=0; >> 616 >> 617 >> 618 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()]; >> 619 >> 620 if(waterDensity!= 0.0) >> 621 // if (material == nistwater || material->GetBaseMaterial() == nistwater) >> 622 { >> 623 const G4String& particleName = particleDefinition->GetParticleName(); >> 624 >> 625 std::map< G4String,G4double,std::less<G4String> >::iterator pos2; >> 626 pos2 = highEnergyLimit.find(particleName); >> 627 >> 628 if (pos2 != highEnergyLimit.end()) >> 629 { >> 630 highLim = pos2->second; >> 631 } >> 632 >> 633 if (k <= highLim) >> 634 { >> 635 >> 636 //SI : XS must not be zero otherwise sampling of secondaries method ignored >> 637 >> 638 if (k < lowLim) k = lowLim; >> 639 >> 640 // >> 641 >> 642 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos; >> 643 pos = tableData.find(particleName); >> 644 >> 645 if (pos != tableData.end()) >> 646 { >> 647 G4DNACrossSectionDataSet* table = pos->second; >> 648 if (table != 0) >> 649 { >> 650 sigma = table->FindValue(k); >> 651 } >> 652 } >> 653 else >> 654 { >> 655 G4Exception("G4DNARuddIonisationExtendedModel::CrossSectionPerVolume","em0002", >> 656 FatalException,"Model not applicable to particle type."); >> 657 } >> 658 >> 659 } // if (k >= lowLim && k < highLim) >> 660 >> 661 if (verboseLevel > 2) >> 662 { >> 663 G4cout << "__________________________________" << G4endl; >> 664 G4cout << "G4DNARuddIonisationExtendedModel - XS INFO START" << G4endl; >> 665 G4cout << "Kinetic energy(eV)=" << k/eV << " particle : " << particleDefinition->GetParticleName() << G4endl; >> 666 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl; >> 667 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl; >> 668 //G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl; >> 669 G4cout << "G4DNARuddIonisationExtendedModel - XS INFO END" << G4endl; >> 670 >> 671 } >> 672 >> 673 } // if (waterMaterial) >> 674 >> 675 return sigma*waterDensity; >> 676 // return sigma*material->GetAtomicNumDensityVector()[1]; >> 677 161 } 678 } 162 679 163 //....oooOO0OOooo........oooOO0OOooo........oo 680 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 164 681 165 void G4DNARuddIonisationExtendedModel::Initial << 682 void G4DNARuddIonisationExtendedModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect, 166 << 683 const G4MaterialCutsCouple* couple, >> 684 const G4DynamicParticle* particle, >> 685 G4double, >> 686 G4double) 167 { 687 { 168 if (p != fParticle) { SetParticle(p); } << 688 //SEB: particle->GetDefinition()->GetParticleName() is for eg. Fe56 >> 689 // particle->GetDefinition()->GetPDGMass() is correct >> 690 // particle->GetDefinition()->GetAtomicNumber() is correct >> 691 // particle->GetDefinition()->GetAtomicMass() is correct >> 692 >> 693 if (verboseLevel > 3) >> 694 G4cout << "Calling SampleSecondaries() of G4DNARuddIonisationExtendedModel" << G4endl; >> 695 >> 696 G4double lowLim = 0; >> 697 G4double highLim = 0; >> 698 >> 699 // ZF: the following line summarizes the commented part >> 700 >> 701 if(particle->GetDefinition()->GetAtomicMass() <= 4) lowLim = killBelowEnergyForA[particle->GetDefinition()->GetAtomicMass()]; >> 702 >> 703 else lowLim = killBelowEnergyForA[5]*particle->GetDefinition()->GetAtomicMass(); >> 704 >> 705 /* >> 706 >> 707 if(particle->GetDefinition()->GetAtomicMass() >= 5) lowLim = killBelowEnergyForA[5]*particle->GetDefinition()->GetAtomicMass(); >> 708 >> 709 if ( particle->GetDefinition() == G4Proton::ProtonDefinition() >> 710 || particle->GetDefinition() == instance->GetIon("hydrogen") >> 711 ) >> 712 >> 713 lowLim = killBelowEnergyForA[1]; >> 714 >> 715 if ( particle->GetDefinition() == instance->GetIon("alpha++") >> 716 || particle->GetDefinition() == instance->GetIon("alpha+") >> 717 || particle->GetDefinition() == instance->GetIon("helium") >> 718 ) >> 719 >> 720 lowLim = killBelowEnergyForA[4]; >> 721 >> 722 */ >> 723 // >> 724 >> 725 G4double k = particle->GetKineticEnergy(); >> 726 >> 727 const G4String& particleName = particle->GetDefinition()->GetParticleName(); >> 728 >> 729 // SI - the following is useless since lowLim is already defined >> 730 /* >> 731 std::map< G4String,G4double,std::less<G4String> >::iterator pos1; >> 732 pos1 = lowEnergyLimit.find(particleName); >> 733 >> 734 if (pos1 != lowEnergyLimit.end()) >> 735 { >> 736 lowLim = pos1->second; >> 737 } >> 738 */ 169 739 170 // particle change object may be externally << 740 std::map< G4String,G4double,std::less<G4String> >::iterator pos2; 171 if (nullptr == fParticleChangeForGamma) { << 741 pos2 = highEnergyLimit.find(particleName); 172 fParticleChangeForGamma = GetParticleChang << 173 } << 174 742 175 // initialisation once in each thread << 743 if (pos2 != highEnergyLimit.end())highLim = pos2->second; 176 if (!isInitialised) { << 744 177 isInitialised = true; << 745 if (k >= lowLim && k < highLim) 178 const G4String& pname = fParticle->GetPart << 746 { 179 if (pname == "proton") { << 747 G4ParticleDefinition* definition = particle->GetDefinition(); 180 idx = 1; << 748 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection(); 181 xscurrent = xsdata[1]; << 749 /* 182 fElow = fLowestEnergy; << 750 G4double particleMass = definition->GetPDGMass(); 183 } else if (pname == "hydrogen") { << 751 G4double totalEnergy = k + particleMass; 184 idx = 0; << 752 G4double pSquare = k*(totalEnergy+particleMass); 185 xscurrent = xsdata[0]; << 753 G4double totalMomentum = std::sqrt(pSquare); 186 fElow = fLowestEnergy; << 754 */ 187 } else if (pname == "alpha") { << 755 188 idx = 1; << 756 G4int ionizationShell = RandomSelect(k,particleName); 189 xscurrent = xsdata[2]; << 757 190 isHelium = true; << 758 // sample deexcitation 191 fElow = fLimitEnergy; << 759 // here we assume that H_{2}O electronic levels are the same as Oxygen. 192 } else if (pname == "alpha+") { << 760 // this can be considered true with a rough 10% error in energy on K-shell, 193 idx = 1; << 761 194 isHelium = true; << 762 G4int secNumberInit = 0; // need to know at a certain point the energy of secondaries 195 xscurrent = xsalphaplus; << 763 G4int secNumberFinal = 0; // So I'll make the diference and then sum the energies 196 fElow = fLimitEnergy; << 764 G4double bindingEnergy = 0; 197 // The following values are provided by << 765 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell); 198 slaterEffectiveCharge[0]=2.0; << 766 199 slaterEffectiveCharge[1]=2.0; << 767 //SI: additional protection if tcs interpolation method is modified 200 slaterEffectiveCharge[2]=2.0; << 768 if (k<bindingEnergy) return; 201 sCoefficient[0]=0.7; << 769 // 202 sCoefficient[1]=0.15; << 770 203 sCoefficient[2]=0.15; << 771 G4int Z = 8; 204 } else if (pname == "helium") { << 772 205 idx = 0; << 773 if(fAtomDeexcitation) { 206 isHelium = true; << 774 G4AtomicShellEnumerator as = fKShell; 207 fElow = fLimitEnergy; << 775 208 xscurrent = xshelium; << 776 if (ionizationShell <5 && ionizationShell >1) 209 slaterEffectiveCharge[0]=1.7; << 777 { 210 slaterEffectiveCharge[1]=1.15; << 778 as = G4AtomicShellEnumerator(4-ionizationShell); 211 slaterEffectiveCharge[2]=1.15; << 779 } 212 sCoefficient[0]=0.5; << 780 else if (ionizationShell <2) 213 sCoefficient[1]=0.25; << 781 { 214 sCoefficient[2]=0.25; << 782 as = G4AtomicShellEnumerator(3); 215 } else { << 783 } 216 isIon = true; << 784 217 idx = -1; << 785 // DEBUG 218 xscurrent = xsdata[1]; << 786 // if (ionizationShell == 4) { 219 fElow = fLowestEnergy; << 787 // 220 } << 788 // G4cout << "Z: " << Z << " as: " << as 221 // defined stationary mode << 789 // << " ionizationShell: " << ionizationShell << " bindingEnergy: "<< bindingEnergy/eV << G4endl; 222 statCode = G4EmParameters::Instance()->DNA << 790 // G4cout << "Press <Enter> key to continue..." << G4endl; 223 << 791 // G4cin.ignore(); 224 // initialise atomic de-excitation << 792 // } 225 fAtomDeexcitation = G4LossTableManager::In << 793 226 << 794 227 if (verbose > 0) { << 795 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as); 228 G4cout << "### G4DNARuddIonisationExtend << 796 secNumberInit = fvect->size(); 229 << "/n idx=" << idx << " isIon=" << << 797 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0); 230 << " isHelium=" << isHelium << G4endl; << 798 secNumberFinal = fvect->size(); >> 799 } >> 800 >> 801 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell); >> 802 >> 803 G4ThreeVector deltaDirection = >> 804 GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic, >> 805 Z, ionizationShell, >> 806 couple->GetMaterial()); >> 807 >> 808 // SI: the following lines are not needed anymore >> 809 /* >> 810 G4double cosTheta = 0.; >> 811 G4double phi = 0.; >> 812 RandomizeEjectedElectronDirection(definition, k,secondaryKinetic, cosTheta, phi, ionizationShell); >> 813 >> 814 G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta); >> 815 G4double dirX = sinTheta*std::cos(phi); >> 816 G4double dirY = sinTheta*std::sin(phi); >> 817 G4double dirZ = cosTheta; >> 818 G4ThreeVector deltaDirection(dirX,dirY,dirZ); >> 819 deltaDirection.rotateUz(primaryDirection); >> 820 */ >> 821 >> 822 // Ignored for ions on electrons >> 823 /* >> 824 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 )); >> 825 >> 826 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x(); >> 827 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y(); >> 828 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z(); >> 829 G4double finalMomentum = std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz); >> 830 finalPx /= finalMomentum; >> 831 finalPy /= finalMomentum; >> 832 finalPz /= finalMomentum; >> 833 >> 834 G4ThreeVector direction; >> 835 direction.set(finalPx,finalPy,finalPz); >> 836 >> 837 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ; >> 838 */ >> 839 >> 840 fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection); >> 841 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic; >> 842 G4double deexSecEnergy = 0; >> 843 for (G4int j=secNumberInit; j < secNumberFinal; j++) { >> 844 >> 845 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy(); >> 846 >> 847 } >> 848 >> 849 if (!statCode) >> 850 { >> 851 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy); >> 852 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy-secondaryKinetic-deexSecEnergy); >> 853 } >> 854 else >> 855 { >> 856 fParticleChangeForGamma->SetProposedKineticEnergy(k); >> 857 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy); >> 858 } >> 859 >> 860 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ; >> 861 fvect->push_back(dp); >> 862 >> 863 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack(); >> 864 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule, >> 865 ionizationShell, >> 866 theIncomingTrack); 231 } 867 } 232 } << 868 >> 869 // SI - not useful since low energy of model is 0 eV >> 870 >> 871 if (k < lowLim) >> 872 { >> 873 fParticleChangeForGamma->SetProposedKineticEnergy(0.); >> 874 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill); >> 875 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k); >> 876 } >> 877 233 } 878 } 234 879 235 //....oooOO0OOooo........oooOO0OOooo........oo << 880 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 236 881 237 void G4DNARuddIonisationExtendedModel::SetPart << 882 G4double G4DNARuddIonisationExtendedModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition, >> 883 G4double k, >> 884 G4int shell) 238 { 885 { 239 fParticle = p; << 886 //-- Fast sampling method ----- 240 fMass = p->GetPDGMass(); << 887 G4double proposed_energy; 241 fMassRate = (isIon) ? CLHEP::proton_mass_c2/ << 888 G4double random1; >> 889 G4double value_sampling; >> 890 G4double max1; >> 891 >> 892 do >> 893 { >> 894 proposed_energy = ProposedSampledEnergy(particleDefinition, k, shell); // Proposed energy by inverse function sampling >> 895 >> 896 max1=0.; >> 897 >> 898 for(G4double en=0.; en<20.; en+=1.) if(RejectionFunction(particleDefinition, k, en, shell) > max1) >> 899 max1=RejectionFunction(particleDefinition, k, en, shell); >> 900 >> 901 random1 = G4UniformRand()*max1; >> 902 >> 903 value_sampling = RejectionFunction(particleDefinition, k, proposed_energy, shell); >> 904 >> 905 } while(random1 > value_sampling); >> 906 >> 907 return(proposed_energy); 242 } 908 } 243 909 244 //....oooOO0OOooo........oooOO0OOooo........oo << 910 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 245 911 246 G4double << 912 // The following section is not used anymore but is kept for memory 247 G4DNARuddIonisationExtendedModel::CrossSection << 913 // GetAngularDistribution()->SampleDirectionForShell is used instead 248 << 914 249 << 915 /* 250 << 916 void G4DNARuddIonisationExtendedModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition, >> 917 G4double k, >> 918 G4double secKinetic, >> 919 G4double & cosTheta, >> 920 G4double & phi, >> 921 G4int shell ) 251 { 922 { 252 // check if model is applicable for given ma << 923 G4double maxSecKinetic = 0.; 253 G4double density = (material->GetIndex() < f << 924 G4double maximumEnergyTransfer = 0.; 254 ? (*fpWaterDensity)[material->GetIndex()] << 925 255 if (0.0 == density) { return 0.0; } << 926 // ZF. generalized & relativistic version 256 927 257 // ion may be different << 928 if( (k/MeV)/(particleDefinition->GetPDGMass()/MeV) <= 0.1 ) 258 if (fParticle != part) { SetParticle(part); << 929 { >> 930 maximumEnergyTransfer= 4.* (electron_mass_c2 / particleDefinition->GetPDGMass()) * k; >> 931 maximumEnergyTransfer+=waterStructure.IonisationEnergy(shell); >> 932 } >> 933 else >> 934 { >> 935 G4double approx_nuc_number = particleDefinition->GetPDGMass() / proton_mass_c2; >> 936 G4double en_per_nucleon = k/approx_nuc_number; >> 937 G4double beta2 = 1. - 1./pow( (1.+(en_per_nucleon/electron_mass_c2)*(electron_mass_c2/proton_mass_c2)), 2.); >> 938 G4double gamma = 1./sqrt(1.-beta2); >> 939 maximumEnergyTransfer = 2.*electron_mass_c2*(gamma*gamma-1.)/(1.+2.*gamma*(electron_mass_c2/particleDefinition->GetPDGMass())+pow(electron_mass_c2/particleDefinition->GetPDGMass(), 2.) ); >> 940 maximumEnergyTransfer+=waterStructure.IonisationEnergy(shell); >> 941 } 259 942 260 // ion shoud be stopped - check on kinetic e << 943 maxSecKinetic = maximumEnergyTransfer-waterStructure.IonisationEnergy(shell); 261 if (kinE < fLowestEnergy) { return DBL_MAX; << 262 944 263 G4double e = kinE*fMassRate; << 945 phi = twopi * G4UniformRand(); 264 946 265 G4double sigma = (e > fElow) ? xscurrent->Fi << 947 if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic); 266 : xscurrent->FindValue(fElow) * e / fElow; << 948 else cosTheta = (2.*G4UniformRand())-1.; 267 949 268 if (idx == -1) { << 950 } 269 sigma *= fEmCorrections->EffectiveChargeSq << 951 */ 270 } << 271 952 272 sigma *= density; << 953 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 273 954 274 if (verbose > 1) { << 955 G4double G4DNARuddIonisationExtendedModel::RejectionFunction(G4ParticleDefinition* particleDefinition, 275 G4cout << "G4DNARuddIonisationExtendedMode << 956 G4double k, 276 << " Ekin(keV)=" << kinE/CLHEP::keV << 957 G4double proposed_ws, 277 << " sigma(cm^2)=" << sigma/CLHEP:: << 958 G4int ionizationLevelIndex) 278 } << 959 { 279 return sigma; << 960 const G4int j=ionizationLevelIndex; >> 961 G4double Bj_energy, alphaConst; >> 962 G4double Ry = 13.6*eV; >> 963 const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.}; >> 964 >> 965 // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV}; //Ding Paper >> 966 >> 967 // Following values provided by M. Dingfelder (priv. comm) >> 968 const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV}; >> 969 >> 970 if (j == 4) >> 971 { >> 972 alphaConst = 0.66; >> 973 //---Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm) >> 974 Bj_energy = waterStructure.IonisationEnergy(ionizationLevelIndex); >> 975 //--- >> 976 } >> 977 else >> 978 { >> 979 alphaConst = 0.64; >> 980 Bj_energy = Bj[ionizationLevelIndex]; >> 981 } >> 982 >> 983 G4double energyTransfer = proposed_ws + Bj_energy; >> 984 proposed_ws/=Bj_energy; >> 985 G4DNAGenericIonsManager *instance; >> 986 instance = G4DNAGenericIonsManager::Instance(); >> 987 G4double tau = 0.; >> 988 G4double A_ion = 0.; >> 989 tau = (electron_mass_c2 / particleDefinition->GetPDGMass()) * k; >> 990 A_ion = particleDefinition->GetAtomicMass(); >> 991 >> 992 G4double v2; >> 993 G4double beta2; >> 994 >> 995 if((tau/MeV)<5.447761194e-2) >> 996 { >> 997 v2 = tau / Bj_energy; >> 998 beta2 = 2.*tau / electron_mass_c2; >> 999 } >> 1000 // Relativistic >> 1001 else >> 1002 { >> 1003 v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) )); >> 1004 beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion)); >> 1005 } >> 1006 >> 1007 G4double v = std::sqrt(v2); >> 1008 G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj_energy)); >> 1009 G4double rejection_term = 1.+G4Exp(alphaConst*(proposed_ws - wc) / v); >> 1010 rejection_term = (1./rejection_term)*CorrectionFactor(particleDefinition,k,ionizationLevelIndex) * Gj[j]; >> 1011 //* (S/Bj_energy) ; Not needed anymore >> 1012 >> 1013 G4bool isHelium = false; >> 1014 >> 1015 if ( particleDefinition == G4Proton::ProtonDefinition() >> 1016 || particleDefinition == instance->GetIon("hydrogen") >> 1017 ) >> 1018 { >> 1019 return(rejection_term); >> 1020 } >> 1021 >> 1022 else if(particleDefinition->GetAtomicMass() > 4) // anything above Helium >> 1023 { >> 1024 G4double Z = particleDefinition->GetAtomicNumber(); >> 1025 >> 1026 G4double x = 100.*std::sqrt(beta2)/std::pow(Z,(2./3.)); >> 1027 G4double Zeffion = Z*(1.-G4Exp(-1.316*x+0.112*x*x-0.0650*x*x*x)); >> 1028 rejection_term*=Zeffion*Zeffion; >> 1029 } >> 1030 >> 1031 else if (particleDefinition == instance->GetIon("alpha++") ) >> 1032 { >> 1033 isHelium = true; >> 1034 slaterEffectiveCharge[0]=0.; >> 1035 slaterEffectiveCharge[1]=0.; >> 1036 slaterEffectiveCharge[2]=0.; >> 1037 sCoefficient[0]=0.; >> 1038 sCoefficient[1]=0.; >> 1039 sCoefficient[2]=0.; >> 1040 } >> 1041 >> 1042 else if (particleDefinition == instance->GetIon("alpha+") ) >> 1043 { >> 1044 isHelium = true; >> 1045 slaterEffectiveCharge[0]=2.0; >> 1046 // The following values are provided by M. Dingfelder (priv. comm) >> 1047 slaterEffectiveCharge[1]=2.0; >> 1048 slaterEffectiveCharge[2]=2.0; >> 1049 // >> 1050 sCoefficient[0]=0.7; >> 1051 sCoefficient[1]=0.15; >> 1052 sCoefficient[2]=0.15; >> 1053 } >> 1054 >> 1055 else if (particleDefinition == instance->GetIon("helium") ) >> 1056 { >> 1057 isHelium = true; >> 1058 slaterEffectiveCharge[0]=1.7; >> 1059 slaterEffectiveCharge[1]=1.15; >> 1060 slaterEffectiveCharge[2]=1.15; >> 1061 sCoefficient[0]=0.5; >> 1062 sCoefficient[1]=0.25; >> 1063 sCoefficient[2]=0.25; >> 1064 } >> 1065 >> 1066 // if ( particleDefinition == instance->GetIon("helium") >> 1067 // || particleDefinition == instance->GetIon("alpha+") >> 1068 // || particleDefinition == instance->GetIon("alpha++") >> 1069 // ) >> 1070 >> 1071 if (isHelium) >> 1072 { >> 1073 >> 1074 G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber(); >> 1075 >> 1076 zEff -= ( sCoefficient[0] * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.) + >> 1077 sCoefficient[1] * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.) + >> 1078 sCoefficient[2] * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.) ); >> 1079 >> 1080 rejection_term*= zEff * zEff; >> 1081 } >> 1082 >> 1083 return (rejection_term); 280 } 1084 } 281 1085 282 //....oooOO0OOooo........oooOO0OOooo........oo << 283 1086 284 void << 1087 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 285 G4DNARuddIonisationExtendedModel::SampleSecond << 286 << 287 << 288 << 289 { << 290 const G4ParticleDefinition* pd = dpart->GetD << 291 if (fParticle != pd) { SetParticle(pd); } << 292 << 293 // stop ion with energy below low energy lim << 294 G4double kinE = dpart->GetKineticEnergy(); << 295 // ion shoud be stopped - check on kinetic e << 296 if (kinE <= fLowestEnergy) { << 297 fParticleChangeForGamma->SetProposedKineti << 298 fParticleChangeForGamma->ProposeTrackStatu << 299 fParticleChangeForGamma->ProposeLocalEnerg << 300 return; << 301 } << 302 << 303 G4int shell = SelectShell(kinE*fMassRate); << 304 G4double bindingEnergy = (useDNAWaterStructu << 305 ? waterStructure.IonisationEnergy(shell) : << 306 1088 307 //Si: additional protection if tcs interpola << 308 if (kinE < bindingEnergy) { return; } << 309 << 310 G4double esec = SampleElectronEnergy(kinE, s << 311 G4double esum = 0.0; << 312 1089 313 // sample deexcitation << 1090 G4double G4DNARuddIonisationExtendedModel::ProposedSampledEnergy(G4ParticleDefinition* particle, 314 // here we assume that H2O electronic levels << 1091 G4double k, 315 // this can be considered true with a rough << 1092 G4int ionizationLevelIndex) 316 G4int Z = 8; << 317 G4ThreeVector deltaDir = << 318 GetAngularDistribution()->SampleDirectionF << 319 << 320 // SI: only atomic deexcitation from K shell << 321 if(fAtomDeexcitation != nullptr && shell == << 322 auto as = G4AtomicShellEnumerator(0); << 323 auto ashell = fAtomDeexcitation->GetAtomic << 324 fAtomDeexcitation->GenerateParticles(fvect << 325 << 326 // compute energy sum from de-excitation << 327 for (auto const & ptr : *fvect) { << 328 esum += ptr->GetKineticEnergy(); << 329 } << 330 } << 331 // check energy balance << 332 // remaining excitation energy of water mole << 333 G4double exc = bindingEnergy - esum; << 334 << 335 // remaining projectile energy << 336 G4double scatteredEnergy = kinE - bindingEne << 337 if(scatteredEnergy < -tolerance || exc < -to << 338 G4cout << "G4DNARuddIonisationExtendedMode << 339 << "negative final E(keV)=" << scat << 340 << kinE/CLHEP::keV << " " << pd->G << 341 << " Edelta(keV)=" << esec/CLHEP::k << 342 << G4endl; << 343 } << 344 << 345 // projectile << 346 if (!statCode) { << 347 fParticleChangeForGamma->SetProposedKineti << 348 fParticleChangeForGamma->ProposeLocalEnerg << 349 } else { << 350 fParticleChangeForGamma->SetProposedKineti << 351 fParticleChangeForGamma->ProposeLocalEnerg << 352 } << 353 << 354 // delta-electron << 355 auto dp = new G4DynamicParticle(G4Electron: << 356 fvect->push_back(dp); << 357 << 358 // create radical << 359 const G4Track* theIncomingTrack = fParticleC << 360 G4DNAChemistryManager::Instance()->CreateWat << 361 theIncomingTrack); << 362 } << 363 << 364 //....oooOO0OOooo........oooOO0OOooo........oo << 365 << 366 G4int G4DNARuddIonisationExtendedModel::Select << 367 { << 368 G4double sum = 0.0; << 369 G4double xs; << 370 for (G4int i=0; i<5; ++i) { << 371 auto ptr = xscurrent->GetComponent(i); << 372 xs = (e > fElow) ? ptr->FindValue(e) : ptr << 373 sum += xs; << 374 fTemp[i] = sum; << 375 } << 376 sum *= G4UniformRand(); << 377 for (G4int i=0; i<5; ++i) { << 378 if (sum <= fTemp[i]) { return i; } << 379 } << 380 return 0; << 381 } << 382 << 383 //....oooOO0OOooo........oooOO0OOooo........oo << 384 << 385 G4double G4DNARuddIonisationExtendedModel::Max << 386 { << 387 // kinematic limit << 388 G4double tau = kine/fMass; << 389 G4double gam = 1.0 + tau; << 390 G4double emax = 2.0*CLHEP::electron_mass_c2* << 391 << 392 // Initialisation of sampling << 393 G4double A1, B1, C1, D1, E1, A2, B2, C2, D2; << 394 if (shell == 4) { << 395 //Data For Liquid Water K SHELL from Dingf << 396 A1 = 1.25; << 397 B1 = 0.5; << 398 C1 = 1.00; << 399 D1 = 1.00; << 400 E1 = 3.00; << 401 A2 = 1.10; << 402 B2 = 1.30; << 403 C2 = 1.00; << 404 D2 = 0.00; << 405 alphaConst = 0.66; << 406 } else { << 407 //Data For Liquid Water from Dingfelder (P << 408 A1 = 1.02; << 409 B1 = 82.0; << 410 C1 = 0.45; << 411 D1 = -0.80; << 412 E1 = 0.38; << 413 A2 = 1.07; << 414 // Value provided by M. Dingfelder (priv. << 415 B2 = 11.6; << 416 C2 = 0.60; << 417 D2 = 0.04; << 418 alphaConst = 0.64; << 419 } << 420 bEnergy = Bj[shell]; << 421 G4double v2 = 0.25*emax/(bEnergy*gam*gam); << 422 v = std::sqrt(v2); << 423 u = Ry/bEnergy; << 424 wc = 4.*v2 - 2.*v - 0.25*u; << 425 << 426 G4double L1 = (C1 * fGpow->powA(v, D1)) / (1 << 427 G4double L2 = C2 * fGpow->powA(v, D2); << 428 G4double H1 = (A1 * G4Log(1. + v2)) / (v2 + << 429 G4double H2 = (A2 / v2) + (B2 / (v2 * v2)); << 430 << 431 F1 = L1 + H1; << 432 F2 = (L2 * H2) / (L2 + H2); << 433 return emax; << 434 } << 435 << 436 //....oooOO0OOooo........oooOO0OOooo........oo << 437 << 438 G4double G4DNARuddIonisationExtendedModel::Sam << 439 << 440 { << 441 G4double emax = MaxEnergy(kine, shell); << 442 // compute cumulative probability function << 443 G4double step = 1*CLHEP::eV; << 444 auto nn = (G4int)(emax/step); << 445 nn = std::min(std::max(nn, 10), 100); << 446 step = emax/(G4double)nn; << 447 << 448 // find max probability << 449 G4double pmax = ProbabilityFunction(kine, 0. << 450 //G4cout << "## E(keV)=" << kine/keV << " em << 451 // << " pmax(0)=" << pmax << " shell=" << 452 << 453 G4double e0 = 0.0; // energy with max probab << 454 // 2 areas after point with max probability << 455 G4double e1 = emax; << 456 G4double e2 = emax; << 457 G4double p1 = 0.0; << 458 G4double p2 = 0.0; << 459 const G4double f = 0.25; << 460 << 461 // find max probability << 462 G4double e = 0.0; << 463 G4double p = 0.0; << 464 for (G4int i=0; i<nn; ++i) { << 465 e += step; << 466 p = ProbabilityFunction(kine, e, shell); << 467 if (p > pmax) { << 468 pmax = p; << 469 e0 = e; << 470 } else { << 471 break; << 472 } << 473 } << 474 // increase step to be more effective << 475 step *= 2.0; << 476 // 2-nd area << 477 for (G4int i=0; i<nn; ++i) { << 478 e += step; << 479 if (std::abs(e - emax) < step) { << 480 e1 = emax; << 481 break; << 482 } << 483 p = ProbabilityFunction(kine, e, shell); << 484 if (p < f*pmax) { << 485 p1 = p; << 486 e1 = e; << 487 break; << 488 } << 489 } << 490 // 3-d area << 491 if (e < emax) { << 492 for (G4int i=0; i<nn; ++i) { << 493 e += step; << 494 if (std::abs(e - emax) < step) { << 495 e2 = emax; << 496 break; << 497 } << 498 p = ProbabilityFunction(kine, e, shell); << 499 if (p < f*p1) { << 500 p2 = p; << 501 e2 = e; << 502 break; << 503 } << 504 } << 505 } << 506 pmax *= 1.05; << 507 // regression method with 3 regions << 508 G4double s0 = pmax*e1; << 509 G4double s1 = s0 + p1 * (e2 - e1); << 510 G4double s2 = s1 + p2 * (emax - e2); << 511 s0 = (s0 == s1) ? 1.0 : s0 / s2; << 512 s1 = (s1 == s2) ? 1.0 : s1 / s2; << 513 << 514 //G4cout << "pmax=" << pmax << " e1(keV)=" < << 515 // << " p2=" << p2 << " s0=" << s0 << " s1 << 516 << 517 // sampling << 518 G4int count = 0; << 519 G4double ymax, y, deltae; << 520 for (G4int i = 0; i<100000; ++i) { << 521 G4double q = G4UniformRand(); << 522 if (q <= s0) { << 523 ymax = pmax; << 524 deltae = e1 * q / s0; << 525 } else if (q <= s1) { << 526 ymax = p1; << 527 deltae = e1 + (e2 - e1) * (q - s0) / (s1 << 528 } else { << 529 ymax = p2; << 530 deltae = e2 + (emax - e2) * (q - s1) / ( << 531 } << 532 y = ProbabilityFunction(kine, deltae, shel << 533 //G4cout << " " << i << ". deltae=" << << 534 // << " y=" << y << " ymax=" << ymax << 535 if (y > ymax && count < 10) { << 536 ++count; << 537 G4cout << "G4DNARuddIonisationExtendedMo << 538 << fParticle->GetParticleName() << " E( << 539 << " Edelta(keV)=" << deltae/CLHEP::keV << 540 << " y=" << y << " ymax=" << ymax << " << 541 } << 542 if (ymax * G4UniformRand() <= y) { << 543 return deltae; << 544 } << 545 } << 546 deltae = std::min(e0 + step, 0.5*emax); << 547 return deltae; << 548 } << 549 << 550 //....oooOO0OOooo........oooOO0OOooo........oo << 551 << 552 G4double G4DNARuddIonisationExtendedModel::Pro << 553 << 554 << 555 { << 556 // Shells ids are 0 1 2 3 4 (4 is k shell) << 557 // !!Attention, "energyTransfer" here is the << 558 // that the secondary kinetic en << 559 // << 560 // ds S F1(nu) + << 561 // ---- = G(k) * ---- ----------------- << 562 // dw Bj (1+w)^3 * [1 + e << 563 // << 564 // w is the secondary electron kinetic Energ << 565 // << 566 // All the other parameters can be found in << 567 // << 568 // M.Eugene Rudd, 1988, User-Friendly model << 569 // electrons from protons or electron collis << 570 // << 571 G4double w = deltae/bEnergy; << 572 G4double x = alphaConst*(w - wc)/v; << 573 G4double y = (x > -15.) ? 1.0 + G4Exp(x) : 1 << 574 << 575 G4double res = CorrectionFactor(kine, shell) << 576 (fGpow->powN((1. + w)/u, 3) * y); << 577 << 578 if (isHelium) { << 579 G4double energyTransfer = deltae + bEnergy << 580 G4double Zeff = 2.0 - << 581 (sCoefficient[0] * S_1s(kine, energyTran << 582 sCoefficient[1] * S_2s(kine, energyTran << 583 sCoefficient[2] * S_2p(kine, energyTran << 584 << 585 res *= Zeff * Zeff; << 586 } << 587 return std::max(res, 0.0); << 588 } << 589 << 590 //....oooOO0OOooo........oooOO0OOooo........oo << 591 << 592 G4double G4DNARuddIonisationExtendedModel::Com << 593 const G4ParticleDefinition* p, G4doub << 594 { << 595 if (fParticle != p) { SetParticle(p); } << 596 MaxEnergy(e, shell); << 597 return ProbabilityFunction(e, deltae, shell) << 598 } << 599 << 600 //....oooOO0OOooo........oooOO0OOooo........oo << 601 << 602 G4double G4DNARuddIonisationExtendedModel::S_1 << 603 << 604 << 605 << 606 { 1093 { 607 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2) << 608 // Dingfelder, in Chattanooga 2005 proceedin << 609 1094 610 G4double r = Rh(kine, energyTransfer, slater << 1095 const G4int j=ionizationLevelIndex; 611 G4double value = 1. - G4Exp(-2 * r) * ( ( 2. << 1096 612 return value; << 1097 G4double A1, B1, C1, D1, E1, A2, B2, C2, D2; >> 1098 //G4double alphaConst ; >> 1099 G4double Bj_energy; >> 1100 >> 1101 // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV}; //Ding Paper >> 1102 // Following values provided by M. Dingfelder (priv. comm) >> 1103 >> 1104 const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV}; >> 1105 >> 1106 if (j == 4) >> 1107 { >> 1108 //Data For Liquid Water K SHELL from Dingfelder (Protons in Water) >> 1109 A1 = 1.25; >> 1110 B1 = 0.5; >> 1111 C1 = 1.00; >> 1112 D1 = 1.00; >> 1113 E1 = 3.00; >> 1114 A2 = 1.10; >> 1115 B2 = 1.30; >> 1116 C2 = 1.00; >> 1117 D2 = 0.00; >> 1118 //alphaConst = 0.66; >> 1119 //---Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm) >> 1120 Bj_energy = waterStructure.IonisationEnergy(ionizationLevelIndex); >> 1121 //--- >> 1122 } >> 1123 else >> 1124 { >> 1125 //Data For Liquid Water from Dingfelder (Protons in Water) >> 1126 A1 = 1.02; >> 1127 B1 = 82.0; >> 1128 C1 = 0.45; >> 1129 D1 = -0.80; >> 1130 E1 = 0.38; >> 1131 A2 = 1.07; >> 1132 //B2 = 14.6; From Ding Paper >> 1133 // Value provided by M. Dingfelder (priv. comm) >> 1134 B2 = 11.6; >> 1135 // >> 1136 C2 = 0.60; >> 1137 D2 = 0.04; >> 1138 //alphaConst = 0.64; >> 1139 >> 1140 Bj_energy = Bj[ionizationLevelIndex]; >> 1141 } >> 1142 >> 1143 G4double tau = 0.; >> 1144 G4double A_ion = 0.; >> 1145 tau = (electron_mass_c2 / particle->GetPDGMass()) * k; >> 1146 >> 1147 A_ion = particle->GetAtomicMass(); >> 1148 >> 1149 G4double v2; >> 1150 G4double beta2; >> 1151 if((tau/MeV)<5.447761194e-2) >> 1152 { >> 1153 v2 = tau / Bj_energy; >> 1154 beta2 = 2.*tau / electron_mass_c2; >> 1155 } >> 1156 // Relativistic >> 1157 else >> 1158 { >> 1159 v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) )); >> 1160 beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion)); >> 1161 } >> 1162 >> 1163 G4double v = std::sqrt(v2); >> 1164 //G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj_energy)); >> 1165 G4double L1 = (C1* std::pow(v,(D1))) / (1.+ E1*std::pow(v, (D1+4.))); >> 1166 G4double L2 = C2*std::pow(v,(D2)); >> 1167 G4double H1 = (A1*std::log(1.+v2)) / (v2+(B1/v2)); >> 1168 G4double H2 = (A2/v2) + (B2/(v2*v2)); >> 1169 G4double F1 = L1+H1; >> 1170 G4double F2 = (L2*H2)/(L2+H2); >> 1171 >> 1172 // ZF. generalized & relativistic version >> 1173 G4double maximumEnergy; >> 1174 >> 1175 //---- maximum kinetic energy , non relativistic ------ >> 1176 if( (k/MeV)/(particle->GetPDGMass()/MeV) <= 0.1 ) >> 1177 { >> 1178 maximumEnergy = 4.* (electron_mass_c2 / particle->GetPDGMass()) * k; >> 1179 } >> 1180 //---- relativistic ----------------------------------- >> 1181 else >> 1182 { >> 1183 G4double gamma = 1./sqrt(1.-beta2); >> 1184 maximumEnergy = 2.*electron_mass_c2*(gamma*gamma-1.)/ >> 1185 (1.+2.*gamma*(electron_mass_c2/particle->GetPDGMass())+pow(electron_mass_c2/particle->GetPDGMass(), 2.) ); >> 1186 } >> 1187 >> 1188 //either it is transfered energy or secondary electron energy ... >> 1189 //maximumEnergy-=Bj_energy; >> 1190 >> 1191 //----------------------------------------------------- >> 1192 G4double wmax = maximumEnergy/Bj_energy; >> 1193 G4double c = wmax*(F2*wmax+F1*(2.+wmax))/(2.*(1.+wmax)*(1.+wmax)); >> 1194 c=1./c; //!!!!!!!!!!! manual calculus leads to c=1/c >> 1195 G4double randVal = G4UniformRand(); >> 1196 G4double proposed_ws = F1*F1*c*c + 2.*F2*c*randVal - 2.*F1*c*randVal; >> 1197 proposed_ws = -F1*c+2.*randVal+std::sqrt(proposed_ws); >> 1198 // proposed_ws = -F1*c+2.*randVal-std::sqrt(proposed_ws); >> 1199 proposed_ws/= ( F1*c + F2*c - 2.*randVal ); >> 1200 proposed_ws*=Bj_energy; >> 1201 >> 1202 return(proposed_ws); >> 1203 613 } 1204 } 614 1205 615 //....oooOO0OOooo........oooOO0OOooo........oo 1206 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 616 1207 617 G4double G4DNARuddIonisationExtendedModel::S_2 << 1208 G4double G4DNARuddIonisationExtendedModel::S_1s(G4double t, 618 << 1209 G4double energyTransferred, 619 << 1210 G4double slaterEffectiveChg, 620 1211 G4double shellNumber) 621 { 1212 { 622 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4) << 1213 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2) 623 // Dingfelder, in Chattanooga 2005 proceedin << 1214 // Dingfelder, in Chattanooga 2005 proceedings, formula (7) 624 1215 625 G4double r = Rh(kine, energyTransfer, slater << 1216 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber); 626 G4double value = << 1217 G4double value = 1. - G4Exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. ); 627 1. - G4Exp(-2 * r) * (((2. * r * r + 2.) * << 628 1218 629 return value; << 1219 return value; 630 } 1220 } 631 1221 632 //....oooOO0OOooo........oooOO0OOooo........oo 1222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 633 1223 634 G4double G4DNARuddIonisationExtendedModel::S_2 << 1224 G4double G4DNARuddIonisationExtendedModel::S_2s(G4double t, 635 << 1225 G4double energyTransferred, 636 << 1226 G4double slaterEffectiveChg, 637 1227 G4double shellNumber) 638 { 1228 { 639 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^ << 1229 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4) 640 // Dingfelder, in Chattanooga 2005 proceedin << 1230 // Dingfelder, in Chattanooga 2005 proceedings, formula (8) >> 1231 >> 1232 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber); >> 1233 G4double value = 1. - G4Exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.); 641 1234 642 G4double r = Rh(kine, energyTransfer, slater << 1235 return value; 643 G4double value = << 644 1. - G4Exp(-2 * r) * (((( 2./3. * r + 4./3 << 645 1236 646 return value; << 647 } 1237 } 648 1238 649 //....oooOO0OOooo........oooOO0OOooo........oo 1239 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 650 1240 651 G4double G4DNARuddIonisationExtendedModel::Rh( << 1241 G4double G4DNARuddIonisationExtendedModel::S_2p(G4double t, 652 << 1242 G4double energyTransferred, >> 1243 G4double slaterEffectiveChg, >> 1244 G4double shellNumber) 653 { 1245 { 654 // The following values are provided by M. D << 1246 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4) 655 // Dingfelder, in Chattanooga 2005 proceedin << 1247 // Dingfelder, in Chattanooga 2005 proceedings, formula (9) 656 1248 657 G4double escaled = CLHEP::electron_mass_c2/f << 1249 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber); 658 const G4double H = 13.60569172 * CLHEP::eV; << 1250 G4double value = 1. - G4Exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.); 659 G4double value = 2.0*std::sqrt(escaled / H)* << 660 1251 661 return value; << 1252 return value; 662 } 1253 } 663 1254 664 //....oooOO0OOooo........oooOO0OOooo........oo 1255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 665 1256 666 G4double << 1257 G4double G4DNARuddIonisationExtendedModel::R(G4double t, 667 G4DNARuddIonisationExtendedModel::CorrectionFa << 1258 G4double energyTransferred, >> 1259 G4double slaterEffectiveChg, >> 1260 G4double shellNumber) 668 { 1261 { 669 // ZF Shortened << 1262 // tElectron = m_electron / m_alpha * t 670 G4double res = 1.0; << 1263 // Dingfelder, in Chattanooga 2005 proceedings, p 4 671 if (shell < 4 && 0 != idx) { << 1264 672 const G4double ln10 = fGpow->logZ(10); << 1265 G4double tElectron = 0.511/3728. * t; 673 G4double x = 2.0*((G4Log(kine/CLHEP::eV)/l << 674 // The following values are provided by M. 1266 // The following values are provided by M. Dingfelder (priv. comm) 675 res = 0.6/(1.0 + G4Exp(x)) + 0.9; << 1267 G4double H = 2.*13.60569172 * eV; 676 } << 1268 G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) * (slaterEffectiveChg/shellNumber); 677 return res; << 1269 >> 1270 return value; >> 1271 } >> 1272 >> 1273 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 1274 >> 1275 G4double G4DNARuddIonisationExtendedModel::CorrectionFactor(G4ParticleDefinition* particleDefinition, G4double k, G4int shell) >> 1276 { >> 1277 // ZF Shortened >> 1278 G4DNAGenericIonsManager *instance; >> 1279 instance = G4DNAGenericIonsManager::Instance(); >> 1280 >> 1281 if (particleDefinition == instance->GetIon("hydrogen") && shell < 4) >> 1282 { >> 1283 G4double value = (std::log10(k/eV)-4.2)/0.5; >> 1284 // The following values are provided by M. Dingfelder (priv. comm) >> 1285 return((0.6/(1+G4Exp(value))) + 0.9); >> 1286 } >> 1287 else >> 1288 { >> 1289 return(1.); >> 1290 } 678 } 1291 } 679 1292 680 //....oooOO0OOooo........oooOO0OOooo........oo 1293 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 1294 >> 1295 G4int G4DNARuddIonisationExtendedModel::RandomSelect(G4double k, const G4String& particle ) >> 1296 { >> 1297 >> 1298 G4int level = 0; >> 1299 >> 1300 // Retrieve data table corresponding to the current particle type >> 1301 >> 1302 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos; >> 1303 pos = tableData.find(particle); >> 1304 >> 1305 if (pos != tableData.end()) >> 1306 { >> 1307 G4DNACrossSectionDataSet* table = pos->second; >> 1308 >> 1309 if (table != 0) >> 1310 { >> 1311 G4double* valuesBuffer = new G4double[table->NumberOfComponents()]; >> 1312 >> 1313 const size_t n(table->NumberOfComponents()); >> 1314 size_t i(n); >> 1315 G4double value = 0.; >> 1316 >> 1317 while (i>0) >> 1318 { >> 1319 i--; >> 1320 valuesBuffer[i] = table->GetComponent(i)->FindValue(k); >> 1321 >> 1322 value += valuesBuffer[i]; >> 1323 } >> 1324 >> 1325 value *= G4UniformRand(); >> 1326 >> 1327 i = n; >> 1328 >> 1329 while (i > 0) >> 1330 { >> 1331 i--; >> 1332 >> 1333 if (valuesBuffer[i] > value) >> 1334 { >> 1335 delete[] valuesBuffer; >> 1336 return i; >> 1337 } >> 1338 value -= valuesBuffer[i]; >> 1339 } >> 1340 >> 1341 if (valuesBuffer) delete[] valuesBuffer; >> 1342 >> 1343 } >> 1344 } >> 1345 else >> 1346 { >> 1347 G4Exception("G4DNARuddIonisationExtendedModel::RandomSelect","em0002", >> 1348 FatalException,"Model not applicable to particle type."); >> 1349 } >> 1350 >> 1351 return level; >> 1352 } >> 1353 >> 1354 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 1355 >> 1356 G4double G4DNARuddIonisationExtendedModel::PartialCrossSection(const G4Track& track ) >> 1357 { >> 1358 G4double sigma = 0.; >> 1359 >> 1360 const G4DynamicParticle* particle = track.GetDynamicParticle(); >> 1361 G4double k = particle->GetKineticEnergy(); >> 1362 >> 1363 G4double lowLim = 0; >> 1364 G4double highLim = 0; >> 1365 >> 1366 const G4String& particleName = particle->GetDefinition()->GetParticleName(); >> 1367 >> 1368 std::map< G4String,G4double,std::less<G4String> >::iterator pos1; >> 1369 pos1 = lowEnergyLimit.find(particleName); >> 1370 >> 1371 if (pos1 != lowEnergyLimit.end()) >> 1372 { >> 1373 lowLim = pos1->second; >> 1374 } >> 1375 >> 1376 std::map< G4String,G4double,std::less<G4String> >::iterator pos2; >> 1377 pos2 = highEnergyLimit.find(particleName); >> 1378 >> 1379 if (pos2 != highEnergyLimit.end()) >> 1380 { >> 1381 highLim = pos2->second; >> 1382 } >> 1383 >> 1384 if (k >= lowLim && k <= highLim) >> 1385 { >> 1386 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos; >> 1387 pos = tableData.find(particleName); >> 1388 >> 1389 if (pos != tableData.end()) >> 1390 { >> 1391 G4DNACrossSectionDataSet* table = pos->second; >> 1392 if (table != 0) >> 1393 { >> 1394 sigma = table->FindValue(k); >> 1395 } >> 1396 } >> 1397 else >> 1398 { >> 1399 G4Exception("G4DNARuddIonisationExtendedModel::PartialCrossSection","em0002", >> 1400 FatalException,"Model not applicable to particle type."); >> 1401 } >> 1402 } >> 1403 >> 1404 return sigma; >> 1405 } >> 1406 >> 1407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 1408 >> 1409 G4double G4DNARuddIonisationExtendedModel::Sum(G4double /* energy */, const G4String& /* particle */) >> 1410 { >> 1411 return 0; >> 1412 } >> 1413 681 1414