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 // Author: H. N. Tran (Ton Duc Thang Universit 27 // p, H, He, He+ and He++ models are assumed i 28 // NIMB 343, 132-137 (2015) 29 // 30 // The Geant4-DNA web site is available at htt 31 // 32 33 #include "G4DNAIonElasticModel.hh" 34 #include "G4PhysicalConstants.hh" 35 #include "G4SystemOfUnits.hh" 36 #include "G4DNAMolecularMaterial.hh" 37 #include "G4ParticleTable.hh" 38 #include "G4Exp.hh" 39 40 //....oooOO0OOooo........oooOO0OOooo........oo 41 42 using namespace std; 43 44 //....oooOO0OOooo........oooOO0OOooo........oo 45 46 G4DNAIonElasticModel::G4DNAIonElasticModel (co 47 co 48 G4VEmModel(nam) 49 { 50 killBelowEnergy = 100 * eV; 51 lowEnergyLimit = 0 * eV; 52 highEnergyLimit = 1 * MeV; 53 SetLowEnergyLimit(lowEnergyLimit); 54 SetHighEnergyLimit(highEnergyLimit); 55 56 verboseLevel = 0; 57 // Verbosity scale: 58 // 0 = nothing 59 // 1 = warning for energy non-conservation 60 // 2 = details of energy budget 61 // 3 = calculation of cross sections, file o 62 // 4 = entering in methods 63 64 if(verboseLevel > 0) 65 { 66 G4cout << "Ion elastic model is constructe 67 << lowEnergyLimit / eV << " eV - " 68 << highEnergyLimit / MeV << " MeV" 69 << G4endl; 70 } 71 72 fParticleChangeForGamma = nullptr; 73 fpMolWaterDensity = nullptr; 74 fpTableData = nullptr; 75 fParticle_Mass = -1; 76 77 // Selection of stationary mode 78 79 statCode = false; 80 } 81 82 //....oooOO0OOooo........oooOO0OOooo........oo 83 84 G4DNAIonElasticModel::~G4DNAIonElasticModel () 85 { 86 // For total cross section 87 delete fpTableData; 88 } 89 90 //....oooOO0OOooo........oooOO0OOooo........oo 91 92 void 93 G4DNAIonElasticModel::Initialise ( 94 const G4ParticleDefinition* particleDefini 95 const G4DataVector& /*cuts*/) 96 { 97 98 if(verboseLevel > 3) 99 { 100 G4cout << "Calling G4DNAIonElasticModel::I 101 } 102 103 // Energy limits 104 105 if (LowEnergyLimit() < lowEnergyLimit) 106 { 107 G4cout << "G4DNAIonElasticModel: low energ 108 LowEnergyLimit()/eV << " eV to " << lowEne 109 SetLowEnergyLimit(lowEnergyLimit); 110 } 111 112 if (HighEnergyLimit() > highEnergyLimit) 113 { 114 G4cout << "G4DNAIonElasticModel: high ener 115 HighEnergyLimit()/MeV << " MeV to " << hig 116 SetHighEnergyLimit(highEnergyLimit); 117 } 118 119 // Reading of data files 120 121 G4double scaleFactor = 1e-16*cm*cm; 122 123 const char *path = G4FindDataDir("G4LEDATA") 124 125 if (path == nullptr) 126 { 127 G4Exception("G4IonElasticModel::Initialise 128 FatalException,"G4LEDATA environment v 129 return; 130 } 131 132 G4String totalXSFile; 133 std::ostringstream fullFileName; 134 135 G4DNAGenericIonsManager *instance; 136 instance = G4DNAGenericIonsManager::Instance 137 G4ParticleDefinition* protonDef = 138 G4ParticleTable::GetParticleTable()->FindPar 139 G4ParticleDefinition* hydrogenDef = instance 140 G4ParticleDefinition* heliumDef = instance-> 141 G4ParticleDefinition* alphaplusDef = instanc 142 G4ParticleDefinition* alphaplusplusDef = ins 143 G4String proton, hydrogen, helium, alphaplus 144 145 if ( 146 (particleDefinition == protonDef && prot 147 || 148 (particleDefinition == hydrogenDef && hy 149 ) 150 { 151 // For total cross section of p,h 152 fParticle_Mass = 1.; 153 totalXSFile = "dna/sigma_elastic_proton_HT 154 155 // For final state 156 fullFileName << path << "/dna/sigmadiff_cu 157 } 158 159 if ( 160 (particleDefinition == instance->GetIon( 161 || 162 (particleDefinition == instance->GetIon( 163 || 164 (particleDefinition == instance->GetIon( 165 ) 166 { 167 // For total cross section of he,he+,he++ 168 fParticle_Mass = 4.; 169 totalXSFile = "dna/sigma_elastic_alpha_HTr 170 171 // For final state 172 fullFileName << path << "/dna/sigmadiff_cu 173 } 174 175 fpTableData = new G4DNACrossSectionDataSet(n 176 fpTableData->LoadData(totalXSFile); 177 std::ifstream diffCrossSection(fullFileName. 178 179 if (!diffCrossSection) 180 { 181 G4ExceptionDescription description; 182 description << "Missing data file:" 183 <<fullFileName.str().c_str()<< G4endl; 184 G4Exception("G4IonElasticModel::Initialise 185 FatalException, 186 description); 187 } 188 189 // Added clear for MT 190 191 eTdummyVec.clear(); 192 eVecm.clear(); 193 fDiffCrossSectionData.clear(); 194 195 // 196 197 eTdummyVec.push_back(0.); 198 199 while(!diffCrossSection.eof()) 200 { 201 G4double tDummy; 202 G4double eDummy; 203 diffCrossSection>>tDummy>>eDummy; 204 205 // SI : mandatory eVecm initialization 206 207 if (tDummy != eTdummyVec.back()) 208 { 209 eTdummyVec.push_back(tDummy); 210 eVecm[tDummy].push_back(0.); 211 } 212 213 diffCrossSection>>fDiffCrossSectionData[tD 214 215 if (eDummy != eVecm[tDummy].back()) eVecm[ 216 217 } 218 219 // End final state 220 if( verboseLevel>0 ) 221 { 222 if (verboseLevel > 2) 223 { 224 G4cout << "Loaded cross section files fo 225 } 226 G4cout << "Ion elastic model is initialize 227 << "Energy range: " 228 << LowEnergyLimit() / eV << " eV - " 229 << HighEnergyLimit() / MeV << " MeV" 230 << G4endl; 231 } 232 233 // Initialize water density pointer 234 G4DNAMolecularMaterial::Instance()->Initiali 235 fpMolWaterDensity = G4DNAMolecularMaterial:: 236 GetNumMolPerVolTableFor(G4Material::GetMater 237 238 if (isInitialised) return; 239 fParticleChangeForGamma = GetParticleChangeF 240 isInitialised = true; 241 } 242 243 //....oooOO0OOooo........oooOO0OOooo........oo 244 245 G4double 246 G4DNAIonElasticModel::CrossSectionPerVolume (c 247 c 248 G 249 { 250 if(verboseLevel > 3) 251 { 252 G4cout << "Calling CrossSectionPerVolume() 253 << G4endl; 254 } 255 256 // Calculate total cross section for model 257 258 G4double sigma=0; 259 260 G4double waterDensity = (*fpMolWaterDensity) 261 262 const G4String& particleName = p->GetParticl 263 264 if (ekin <= highEnergyLimit) 265 { 266 //SI : XS must not be zero otherwise sampl 267 if (ekin < killBelowEnergy) return DBL_MAX 268 // 269 270 if (fpTableData != nullptr) 271 { 272 sigma = fpTableData->FindValue(ekin); 273 } 274 else 275 { 276 G4Exception("G4DNAIonElasticModel::Compu 277 FatalException,"Model not applicable 278 } 279 } 280 281 if (verboseLevel > 2) 282 { 283 G4cout << "_______________________________ 284 G4cout << "G4DNAIonElasticModel - XS INFO 285 G4cout << "Kinetic energy(eV)=" << ekin/eV 286 G4cout << "Cross section per water molecul 287 G4cout << "Cross section per water molecul 288 G4cout << "G4DNAIonElasticModel - XS INFO 289 } 290 291 return sigma*waterDensity; 292 } 293 294 //....oooOO0OOooo........oooOO0OOooo........oo 295 296 void 297 G4DNAIonElasticModel::SampleSecondaries ( 298 std::vector<G4DynamicParticle*>* /*fvect*/ 299 const G4MaterialCutsCouple* /*couple*/, 300 const G4DynamicParticle* aDynamicParticle, 301 { 302 303 if(verboseLevel > 3) 304 { 305 G4cout << "Calling SampleSecondaries() of 306 } 307 308 G4double particleEnergy0 = aDynamicParticle- 309 310 if (particleEnergy0 < killBelowEnergy) 311 { 312 fParticleChangeForGamma->SetProposedKineti 313 fParticleChangeForGamma->ProposeTrackStatu 314 fParticleChangeForGamma->ProposeLocalEnerg 315 return; 316 } 317 318 if (particleEnergy0>= killBelowEnergy && par 319 { 320 G4double water_mass = 18.; 321 322 G4double thetaCM = RandomizeThetaCM(partic 323 324 //HT:convert to laboratory system 325 326 G4double theta = std::atan(std::sin(thetaC 327 /(fParticle_Mass/water_mass+std::cos(t 328 329 G4double cosTheta= std::cos(theta); 330 331 // 332 333 G4double phi = 2. * CLHEP::pi * G4UniformR 334 335 G4ThreeVector zVers = aDynamicParticle->Ge 336 G4ThreeVector xVers = zVers.orthogonal(); 337 G4ThreeVector yVers = zVers.cross(xVers); 338 339 G4double xDir = std::sqrt(1. - cosTheta*co 340 G4double yDir = xDir; 341 xDir *= std::cos(phi); 342 yDir *= std::sin(phi); 343 344 G4ThreeVector zPrimeVers((xDir*xVers + yDi 345 346 fParticleChangeForGamma->ProposeMomentumDi 347 348 G4double depositEnergyCM = 0; 349 350 //HT: deposited energy 351 depositEnergyCM = 4. * particleEnergy0 * f 352 (1-std::cos(thetaCM*CLHEP::pi/180)) 353 / (2 * std::pow((fParticle_Mass+water_mass 354 355 //SI: added protection particleEnergy0 >= 356 if (!statCode && (particleEnergy0 >= depos 357 358 fParticleChangeForGamma->SetProposedKine 359 360 else fParticleChangeForGamma->SetProposedK 361 362 fParticleChangeForGamma->ProposeLocalEnerg 363 } 364 } 365 366 //....oooOO0OOooo........oooOO0OOooo........oo 367 368 G4double 369 G4DNAIonElasticModel::Theta (G4ParticleDefinit 370 G4double k, G4dou 371 { 372 G4double theta = 0.; 373 G4double valueT1 = 0; 374 G4double valueT2 = 0; 375 G4double valueE21 = 0; 376 G4double valueE22 = 0; 377 G4double valueE12 = 0; 378 G4double valueE11 = 0; 379 G4double xs11 = 0; 380 G4double xs12 = 0; 381 G4double xs21 = 0; 382 G4double xs22 = 0; 383 384 // Protection against out of boundary access 385 if (k==eTdummyVec.back()) k=k*(1.-1e-12); 386 // 387 388 auto t2 = std::upper_bound(eTdummyVec.begin( 389 390 auto t1 = t2 - 1; 391 392 auto e12 = std::upper_bound(eVecm[(*t1)].beg 393 394 395 auto e11 = e12 - 1; 396 397 auto e22 = std::upper_bound(eVecm[(*t2)].beg 398 399 400 auto e21 = e22 - 1; 401 402 valueT1 = *t1; 403 valueT2 = *t2; 404 valueE21 = *e21; 405 valueE22 = *e22; 406 valueE12 = *e12; 407 valueE11 = *e11; 408 409 xs11 = fDiffCrossSectionData[valueT1][valueE 410 xs12 = fDiffCrossSectionData[valueT1][valueE 411 xs21 = fDiffCrossSectionData[valueT2][valueE 412 xs22 = fDiffCrossSectionData[valueT2][valueE 413 414 if(xs11 == 0 && xs12 == 0 && xs21 == 0 && xs 415 416 theta = QuadInterpolator(valueE11, valueE12, 417 xs21, xs22, valueT1 418 419 return theta; 420 } 421 422 //....oooOO0OOooo........oooOO0OOooo........oo 423 424 G4double 425 G4DNAIonElasticModel::LinLinInterpolate (G4dou 426 G4dou 427 { 428 G4double d1 = xs1; 429 G4double d2 = xs2; 430 G4double value = (d1 + (d2 - d1) * (e - e1) 431 return value; 432 } 433 434 //....oooOO0OOooo........oooOO0OOooo........oo 435 436 G4double 437 G4DNAIonElasticModel::LinLogInterpolate (G4dou 438 G4dou 439 { 440 G4double d1 = std::log(xs1); 441 G4double d2 = std::log(xs2); 442 G4double value = G4Exp(d1 + (d2 - d1) * (e - 443 return value; 444 } 445 446 //....oooOO0OOooo........oooOO0OOooo........oo 447 448 G4double 449 G4DNAIonElasticModel::LogLogInterpolate (G4dou 450 G4dou 451 { 452 G4double a = (std::log10(xs2) - std::log10(x 453 / (std::log10(e2) - std::log10(e1)); 454 G4double b = std::log10(xs2) - a * std::log1 455 G4double sigma = a * std::log10(e) + b; 456 G4double value = (std::pow(10., sigma)); 457 return value; 458 } 459 460 //....oooOO0OOooo........oooOO0OOooo........oo 461 462 G4double 463 G4DNAIonElasticModel::QuadInterpolator (G4doub 464 G4doub 465 G4doub 466 G4doub 467 G4doub 468 G4doub 469 { 470 // Log-Log 471 /* 472 G4double interpolatedvalue1 = LogLogInterpo 473 G4double interpolatedvalue2 = LogLogInterpo 474 G4double value = LogLogInterpolate(t1, t2, 475 */ 476 477 // Lin-Log 478 /* 479 G4double interpolatedvalue1 = LinLogInterpo 480 G4double interpolatedvalue2 = LinLogInterpo 481 G4double value = LinLogInterpolate(t1, t2, 482 */ 483 484 // Lin-Lin 485 G4double interpolatedvalue1 = LinLinInterpol 486 G4double interpolatedvalue2 = LinLinInterpol 487 G4double value = LinLinInterpolate(t1, t2, t 488 interpola 489 490 return value; 491 } 492 493 //....oooOO0OOooo........oooOO0OOooo........oo 494 495 G4double 496 G4DNAIonElasticModel::RandomizeThetaCM ( 497 G4double k, G4ParticleDefinition * particl 498 { 499 G4double integrdiff = G4UniformRand(); 500 return Theta(particleDefinition, k / eV, int 501 } 502 503 //....oooOO0OOooo........oooOO0OOooo........oo 504 505 void 506 G4DNAIonElasticModel::SetKillBelowThreshold (G 507 { 508 killBelowEnergy = threshold; 509 510 if(killBelowEnergy < 100 * eV) 511 { 512 G4cout << "*** WARNING : the G4DNAIonElast 513 "activated below 100 eV !" 514 << G4endl; 515 } 516 } 517 518