Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // Author: H. N. Tran (Ton Duc Thang University) 27 // p, H, He, He+ and He++ models are assumed identical 28 // NIMB 343, 132-137 (2015) 29 // 30 // The Geant4-DNA web site is available at http://geant4-dna.org 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........oooOO0OOooo........oooOO0OOooo.... 41 42 using namespace std; 43 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 45 46 G4DNAIonElasticModel::G4DNAIonElasticModel (const G4ParticleDefinition*, 47 const G4String& nam) : 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 openings, sampling of atoms 62 // 4 = entering in methods 63 64 if(verboseLevel > 0) 65 { 66 G4cout << "Ion elastic model is constructed " << G4endl<< "Energy range: " 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........oooOO0OOooo........oooOO0OOooo.... 83 84 G4DNAIonElasticModel::~G4DNAIonElasticModel () 85 { 86 // For total cross section 87 delete fpTableData; 88 } 89 90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 91 92 void 93 G4DNAIonElasticModel::Initialise ( 94 const G4ParticleDefinition* particleDefinition, 95 const G4DataVector& /*cuts*/) 96 { 97 98 if(verboseLevel > 3) 99 { 100 G4cout << "Calling G4DNAIonElasticModel::Initialise()" << G4endl; 101 } 102 103 // Energy limits 104 105 if (LowEnergyLimit() < lowEnergyLimit) 106 { 107 G4cout << "G4DNAIonElasticModel: low energy limit increased from " << 108 LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl; 109 SetLowEnergyLimit(lowEnergyLimit); 110 } 111 112 if (HighEnergyLimit() > highEnergyLimit) 113 { 114 G4cout << "G4DNAIonElasticModel: high energy limit decreased from " << 115 HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl; 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","em0006", 128 FatalException,"G4LEDATA environment variable not set."); 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()->FindParticle("proton"); 139 G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen"); 140 G4ParticleDefinition* heliumDef = instance->GetIon("helium"); 141 G4ParticleDefinition* alphaplusDef = instance->GetIon("alpha+"); 142 G4ParticleDefinition* alphaplusplusDef = instance->GetIon("alpha++"); 143 G4String proton, hydrogen, helium, alphaplus, alphaplusplus; 144 145 if ( 146 (particleDefinition == protonDef && protonDef != nullptr) 147 || 148 (particleDefinition == hydrogenDef && hydrogenDef != nullptr) 149 ) 150 { 151 // For total cross section of p,h 152 fParticle_Mass = 1.; 153 totalXSFile = "dna/sigma_elastic_proton_HTran"; 154 155 // For final state 156 fullFileName << path << "/dna/sigmadiff_cumulated_elastic_proton_HTran.dat"; 157 } 158 159 if ( 160 (particleDefinition == instance->GetIon("helium") && (heliumDef != nullptr)) 161 || 162 (particleDefinition == instance->GetIon("alpha+") && (alphaplusDef != nullptr)) 163 || 164 (particleDefinition == instance->GetIon("alpha++") && (alphaplusplusDef != nullptr)) 165 ) 166 { 167 // For total cross section of he,he+,he++ 168 fParticle_Mass = 4.; 169 totalXSFile = "dna/sigma_elastic_alpha_HTran"; 170 171 // For final state 172 fullFileName << path << "/dna/sigmadiff_cumulated_elastic_alpha_HTran.dat"; 173 } 174 175 fpTableData = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor ); 176 fpTableData->LoadData(totalXSFile); 177 std::ifstream diffCrossSection(fullFileName.str().c_str()); 178 179 if (!diffCrossSection) 180 { 181 G4ExceptionDescription description; 182 description << "Missing data file:" 183 <<fullFileName.str().c_str()<< G4endl; 184 G4Exception("G4IonElasticModel::Initialise","em0003", 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[tDummy][eDummy]; 214 215 if (eDummy != eVecm[tDummy].back()) eVecm[tDummy].push_back(eDummy); 216 217 } 218 219 // End final state 220 if( verboseLevel>0 ) 221 { 222 if (verboseLevel > 2) 223 { 224 G4cout << "Loaded cross section files for ion elastic model" << G4endl; 225 } 226 G4cout << "Ion elastic model is initialized " << G4endl 227 << "Energy range: " 228 << LowEnergyLimit() / eV << " eV - " 229 << HighEnergyLimit() / MeV << " MeV" 230 << G4endl; 231 } 232 233 // Initialize water density pointer 234 G4DNAMolecularMaterial::Instance()->Initialize(); 235 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()-> 236 GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER")); 237 238 if (isInitialised) return; 239 fParticleChangeForGamma = GetParticleChangeForGamma(); 240 isInitialised = true; 241 } 242 243 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 244 245 G4double 246 G4DNAIonElasticModel::CrossSectionPerVolume (const G4Material* material, 247 const G4ParticleDefinition* p, 248 G4double ekin, G4double, G4double) 249 { 250 if(verboseLevel > 3) 251 { 252 G4cout << "Calling CrossSectionPerVolume() of G4DNAIonElasticModel" 253 << G4endl; 254 } 255 256 // Calculate total cross section for model 257 258 G4double sigma=0; 259 260 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()]; 261 262 const G4String& particleName = p->GetParticleName(); 263 264 if (ekin <= highEnergyLimit) 265 { 266 //SI : XS must not be zero otherwise sampling of secondaries method ignored 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::ComputeCrossSectionPerVolume","em0002", 277 FatalException,"Model not applicable to particle type."); 278 } 279 } 280 281 if (verboseLevel > 2) 282 { 283 G4cout << "__________________________________" << G4endl; 284 G4cout << "G4DNAIonElasticModel - XS INFO START" << G4endl; 285 G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl; 286 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl; 287 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl; 288 G4cout << "G4DNAIonElasticModel - XS INFO END" << G4endl; 289 } 290 291 return sigma*waterDensity; 292 } 293 294 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 295 296 void 297 G4DNAIonElasticModel::SampleSecondaries ( 298 std::vector<G4DynamicParticle*>* /*fvect*/, 299 const G4MaterialCutsCouple* /*couple*/, 300 const G4DynamicParticle* aDynamicParticle, G4double, G4double) 301 { 302 303 if(verboseLevel > 3) 304 { 305 G4cout << "Calling SampleSecondaries() of G4DNAIonElasticModel" << G4endl; 306 } 307 308 G4double particleEnergy0 = aDynamicParticle->GetKineticEnergy(); 309 310 if (particleEnergy0 < killBelowEnergy) 311 { 312 fParticleChangeForGamma->SetProposedKineticEnergy(0.); 313 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill); 314 fParticleChangeForGamma->ProposeLocalEnergyDeposit(particleEnergy0); 315 return; 316 } 317 318 if (particleEnergy0>= killBelowEnergy && particleEnergy0 <= highEnergyLimit) 319 { 320 G4double water_mass = 18.; 321 322 G4double thetaCM = RandomizeThetaCM(particleEnergy0, aDynamicParticle->GetDefinition()); 323 324 //HT:convert to laboratory system 325 326 G4double theta = std::atan(std::sin(thetaCM*CLHEP::pi/180) 327 /(fParticle_Mass/water_mass+std::cos(thetaCM*CLHEP::pi/180))); 328 329 G4double cosTheta= std::cos(theta); 330 331 // 332 333 G4double phi = 2. * CLHEP::pi * G4UniformRand(); 334 335 G4ThreeVector zVers = aDynamicParticle->GetMomentumDirection(); 336 G4ThreeVector xVers = zVers.orthogonal(); 337 G4ThreeVector yVers = zVers.cross(xVers); 338 339 G4double xDir = std::sqrt(1. - cosTheta*cosTheta); 340 G4double yDir = xDir; 341 xDir *= std::cos(phi); 342 yDir *= std::sin(phi); 343 344 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers)); 345 346 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit()); 347 348 G4double depositEnergyCM = 0; 349 350 //HT: deposited energy 351 depositEnergyCM = 4. * particleEnergy0 * fParticle_Mass * water_mass * 352 (1-std::cos(thetaCM*CLHEP::pi/180)) 353 / (2 * std::pow((fParticle_Mass+water_mass),2)); 354 355 //SI: added protection particleEnergy0 >= depositEnergyCM 356 if (!statCode && (particleEnergy0 >= depositEnergyCM) ) 357 358 fParticleChangeForGamma->SetProposedKineticEnergy(particleEnergy0 - depositEnergyCM); 359 360 else fParticleChangeForGamma->SetProposedKineticEnergy(particleEnergy0); 361 362 fParticleChangeForGamma->ProposeLocalEnergyDeposit(depositEnergyCM); 363 } 364 } 365 366 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 367 368 G4double 369 G4DNAIonElasticModel::Theta (G4ParticleDefinition * /*particleDefinition*/, 370 G4double k, G4double integrDiff) 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 eTdummyVec.end(), k); 390 auto t1 = t2 - 1; 391 392 auto e12 = std::upper_bound(eVecm[(*t1)].begin(), 393 eVecm[(*t1)].end(), 394 integrDiff); 395 auto e11 = e12 - 1; 396 397 auto e22 = std::upper_bound(eVecm[(*t2)].begin(), 398 eVecm[(*t2)].end(), 399 integrDiff); 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][valueE11]; 410 xs12 = fDiffCrossSectionData[valueT1][valueE12]; 411 xs21 = fDiffCrossSectionData[valueT2][valueE21]; 412 xs22 = fDiffCrossSectionData[valueT2][valueE22]; 413 414 if(xs11 == 0 && xs12 == 0 && xs21 == 0 && xs22 == 0) return (0.); 415 416 theta = QuadInterpolator(valueE11, valueE12, valueE21, valueE22, xs11, xs12, 417 xs21, xs22, valueT1, valueT2, k, integrDiff); 418 419 return theta; 420 } 421 422 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 423 424 G4double 425 G4DNAIonElasticModel::LinLinInterpolate (G4double e1, G4double e2, G4double e, 426 G4double xs1, G4double xs2) 427 { 428 G4double d1 = xs1; 429 G4double d2 = xs2; 430 G4double value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1)); 431 return value; 432 } 433 434 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 435 436 G4double 437 G4DNAIonElasticModel::LinLogInterpolate (G4double e1, G4double e2, G4double e, 438 G4double xs1, G4double xs2) 439 { 440 G4double d1 = std::log(xs1); 441 G4double d2 = std::log(xs2); 442 G4double value = G4Exp(d1 + (d2 - d1) * (e - e1) / (e2 - e1)); 443 return value; 444 } 445 446 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 447 448 G4double 449 G4DNAIonElasticModel::LogLogInterpolate (G4double e1, G4double e2, G4double e, 450 G4double xs1, G4double xs2) 451 { 452 G4double a = (std::log10(xs2) - std::log10(xs1)) 453 / (std::log10(e2) - std::log10(e1)); 454 G4double b = std::log10(xs2) - a * std::log10(e2); 455 G4double sigma = a * std::log10(e) + b; 456 G4double value = (std::pow(10., sigma)); 457 return value; 458 } 459 460 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 461 462 G4double 463 G4DNAIonElasticModel::QuadInterpolator (G4double e11, G4double e12, 464 G4double e21, G4double e22, 465 G4double xs11, G4double xs12, 466 G4double xs21, G4double xs22, 467 G4double t1, G4double t2, G4double t, 468 G4double e) 469 { 470 // Log-Log 471 /* 472 G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12); 473 G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22); 474 G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2); 475 */ 476 477 // Lin-Log 478 /* 479 G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12); 480 G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22); 481 G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2); 482 */ 483 484 // Lin-Lin 485 G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12); 486 G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22); 487 G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1, 488 interpolatedvalue2); 489 490 return value; 491 } 492 493 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 494 495 G4double 496 G4DNAIonElasticModel::RandomizeThetaCM ( 497 G4double k, G4ParticleDefinition * particleDefinition) 498 { 499 G4double integrdiff = G4UniformRand(); 500 return Theta(particleDefinition, k / eV, integrdiff); 501 } 502 503 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 504 505 void 506 G4DNAIonElasticModel::SetKillBelowThreshold (G4double threshold) 507 { 508 killBelowEnergy = threshold; 509 510 if(killBelowEnergy < 100 * eV) 511 { 512 G4cout << "*** WARNING : the G4DNAIonElasticModel class is not " 513 "activated below 100 eV !" 514 << G4endl; 515 } 516 } 517 518