Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // Authors: S. Meylan and C. Villagrasa (IRSN, 27 // Models come from 28 // M. Bug et al, Rad. Phys and Chem. 130, 459- 29 // 30 31 #include "G4DNAPTBElasticModel.hh" 32 33 #include "G4DNAChampionElasticModel.hh" 34 #include "G4DNAMaterialManager.hh" 35 #include "G4DNAMolecularMaterial.hh" 36 #include "G4Proton.hh" 37 #include "G4SystemOfUnits.hh" 38 39 G4DNAPTBElasticModel::G4DNAPTBElasticModel(con 40 con 41 : G4VDNAModel(nam, applyToMaterial) 42 { 43 if (verboseLevel > 0) { 44 G4cout << "PTB Elastic model is constructe 45 } 46 fpTHF = G4Material::GetMaterial("THF", false 47 fpPY = G4Material::GetMaterial("PY", false); 48 fpPU = G4Material::GetMaterial("PU", false); 49 fpTMP = G4Material::GetMaterial("TMP", false 50 fpG4_WATER = G4Material::GetMaterial("G4_WAT 51 fpBackbone_THF = G4Material::GetMaterial("ba 52 fpCytosine_PY = G4Material::GetMaterial("cyt 53 fpThymine_PY = G4Material::GetMaterial("thym 54 fpAdenine_PU = G4Material::GetMaterial("aden 55 fpBackbone_TMP = G4Material::GetMaterial("ba 56 fpGuanine_PU = G4Material::GetMaterial("guan 57 fpN2 = G4Material::GetMaterial("N2", false); 58 } 59 60 //....oooOO0OOooo........oooOO0OOooo........oo 61 62 void G4DNAPTBElasticModel::Initialise(const G4 63 const G4 64 { 65 if (isInitialised) { 66 return; 67 } 68 if (verboseLevel > 3) 69 { 70 G4cout << "Calling G4DNAPTBElasticModel::I 71 } 72 if (particle != G4Electron::ElectronDefiniti 73 std::ostringstream oss; 74 oss << " Model is not applied for this par 75 G4Exception("G4DNAPTBElasticModel::G4DNAPT 76 oss.str().c_str()); 77 } 78 G4double scaleFactor = 1e-16 * cm * cm; 79 //****************************************** 80 // Cross section data 81 //****************************************** 82 83 std::size_t index; 84 // MPietrzak, adding paths for N2 85 if (fpN2 != nullptr) { 86 index = fpN2->GetIndex(); 87 AddCrossSectionData(index, particle, "dna/ 88 "dna/sigmadiff_cumulat 89 SetLowELimit(index, particle, 10 * eV); 90 SetHighELimit(index, particle, 1.02 * MeV) 91 } 92 // MPietrzak 93 94 if (fpTHF != nullptr) { 95 index = fpTHF->GetIndex(); 96 AddCrossSectionData(index, particle, "dna/ 97 "dna/sigmadiff_cumulat 98 SetLowELimit(index, particle, 10 * eV); 99 SetHighELimit(index, particle, 1 * keV); 100 } 101 102 if (fpPY != nullptr) { 103 index = fpPY->GetIndex(); 104 AddCrossSectionData(index, particle, "dna/ 105 "dna/sigmadiff_cumulat 106 SetLowELimit(index, particle, 10 * eV); 107 SetHighELimit(index, particle, 1 * keV); 108 } 109 110 if (fpPU != nullptr) { 111 index = fpPU->GetIndex(); 112 AddCrossSectionData(index, particle, "dna/ 113 "dna/sigmadiff_cumulat 114 SetLowELimit(index, particle, 10 * eV); 115 SetHighELimit(index, particle, 1 * keV); 116 } 117 118 if (fpTMP != nullptr) { 119 index = fpTMP->GetIndex(); 120 AddCrossSectionData(index, particle, "dna/ 121 "dna/sigmadiff_cumulat 122 SetLowELimit(index, particle, 10 * eV); 123 SetHighELimit(index, particle, 1 * keV); 124 } 125 //???? 126 if (fpG4_WATER != nullptr) { 127 index = fpG4_WATER->GetIndex(); 128 AddCrossSectionData(index, particle, "dna/ 129 "dna/sigmadiff_cumulat 130 SetLowELimit(index, particle, 10 * eV); 131 SetHighELimit(index, particle, 1 * keV); 132 } 133 // DNA materials 134 // 135 if (fpBackbone_THF != nullptr) { 136 index = fpBackbone_THF->GetIndex(); 137 AddCrossSectionData(index, particle, "dna/ 138 "dna/sigmadiff_cumulat 139 SetLowELimit(index, particle, 10 * eV); 140 SetHighELimit(index, particle, 1 * keV); 141 } 142 143 if (fpCytosine_PY != nullptr) { 144 index = fpCytosine_PY->GetIndex(); 145 AddCrossSectionData(index, particle, "dna/ 146 "dna/sigmadiff_cumulat 147 SetLowELimit(index, particle, 10 * eV); 148 SetHighELimit(index, particle, 1 * keV); 149 } 150 151 if (fpThymine_PY != nullptr) { 152 index = fpThymine_PY->GetIndex(); 153 AddCrossSectionData(index, particle, "dna/ 154 "dna/sigmadiff_cumulat 155 SetLowELimit(index, particle, 10 * eV); 156 SetHighELimit(index, particle, 1 * keV); 157 } 158 159 if (fpAdenine_PU != nullptr) { 160 index = fpAdenine_PU->GetIndex(); 161 AddCrossSectionData(index, particle, "dna/ 162 "dna/sigmadiff_cumulat 163 SetLowELimit(index, particle, 10 * eV); 164 SetHighELimit(index, particle, 1 * keV); 165 } 166 if (fpGuanine_PU != nullptr) { 167 index = fpGuanine_PU->GetIndex(); 168 AddCrossSectionData(index, particle, "dna/ 169 "dna/sigmadiff_cumulat 170 SetLowELimit(index, particle, 10 * eV); 171 SetHighELimit(index, particle, 1 * keV); 172 } 173 174 if (fpBackbone_TMP != nullptr) { 175 index = fpBackbone_TMP->GetIndex(); 176 AddCrossSectionData(index, particle, "dna/ 177 "dna/sigmadiff_cumulat 178 SetLowELimit(index, particle, 10 * eV); 179 SetHighELimit(index, particle, 1 * keV); 180 } 181 182 if (!G4DNAMaterialManager::Instance()->IsLoc 183 // Load the data 184 LoadCrossSectionData(particle); 185 G4DNAMaterialManager::Instance()->SetMaste 186 fpModelData = this; 187 } 188 else { 189 auto dataModel = dynamic_cast<G4DNAPTBElas 190 G4DNAMaterialManager::Instance()->GetMod 191 if (dataModel == nullptr) { 192 G4cout << "G4DNAPTBElasticModel::Initial 193 G4Exception("G4DNAPTBElasticModel::Initi 194 "not good modelData"); 195 } 196 else { 197 fpModelData = dataModel; 198 } 199 } 200 201 if (verboseLevel > 2) { 202 G4cout << "Loaded cross section files for 203 } 204 205 fParticleChangeForGamma = GetParticleChangeF 206 isInitialised = true; 207 } 208 209 //....oooOO0OOooo........oooOO0OOooo........oo 210 211 void G4DNAPTBElasticModel::ReadDiffCSFile(cons 212 cons 213 cons 214 { 215 // Method to read and save the information c 216 // This method is not yet standard. 217 218 // get the path of the G4LEDATA data folder 219 const char* path = G4FindDataDir("G4LEDATA") 220 // if it is not found then quit and print er 221 if (path == nullptr) { 222 G4Exception("G4DNAPTBElasticModel::ReadAll 223 "G4LEDATA environment variable 224 return; 225 } 226 227 // build the fullFileName path of the data f 228 std::ostringstream fullFileName; 229 fullFileName << path << "/" << file << ".dat 230 231 // open the data file 232 std::ifstream diffCrossSection(fullFileName. 233 // error if file is not there 234 std::stringstream endPath; 235 if (!diffCrossSection) { 236 endPath << "Missing data file: " << file; 237 G4Exception("G4DNAPTBElasticModel::Initial 238 endPath.str().c_str()); 239 } 240 241 tValuesVec[materialName][particleName].push_ 242 243 G4String line; 244 245 // read the file line by line until we reach 246 while (std::getline(diffCrossSection, line)) 247 // check if the line is comment or empty 248 // 249 std::istringstream testIss(line); 250 G4String test; 251 testIss >> test; 252 // check first caracter to determine if fo 253 if (test == "#") { 254 // skip the line by beginning a new whil 255 continue; 256 } 257 // check if line is empty 258 if (line.empty()) { 259 // skip the line by beginning a new whil 260 continue; 261 } 262 // 263 // end of the check 264 265 // transform the line into a iss 266 std::istringstream iss(line); 267 268 // Variables to be filled by the input fil 269 G4double tDummy; 270 G4double eDummy; 271 272 // fill the variables with the content of 273 iss >> tDummy >> eDummy; 274 275 // SI : mandatory Vecm initialization 276 277 // Fill two vectors contained in maps of t 278 // [materialName][particleName]=vector 279 // [materialName][particleName][T]=vector 280 // to list all the incident energies (tVal 281 // file 282 // 283 // Check if we already have the current T 284 // If not then add it 285 if (tDummy != tValuesVec[materialName][par 286 // Add the current T value 287 tValuesVec[materialName][particleName].p 288 // Make it correspond to a default zero 289 eValuesVect[materialName][particleName][ 290 } 291 292 // Put the differential cross section valu 293 // map 294 iss >> diffCrossSectionData[materialName][ 295 296 // If the current E value (eDummy) is diff 297 // then add it to the vector 298 if (eDummy != eValuesVect[materialName][pa 299 eValuesVect[materialName][particleName][ 300 } 301 } 302 } 303 304 //....oooOO0OOooo........oooOO0OOooo........oo 305 306 G4double G4DNAPTBElasticModel::CrossSectionPer 307 308 309 { 310 if (verboseLevel > 3){ 311 G4cout << "Calling CrossSectionPerVolume() 312 } 313 314 // Get the name of the current particle 315 const G4String& particleName = p->GetParticl 316 const std::size_t& materialID = pMaterial->G 317 318 // set killBelowEnergy value for current mat 319 fKillBelowEnergy = fpModelData->GetLowELimit 320 // initialise the return value (cross sectio 321 G4double sigma = 0.; 322 323 // check if we are below the high energy lim 324 if (ekin < fpModelData->GetHighELimit(materi 325 // This is used to kill the particle if it 326 // If the energy is lower then we return a 327 // method will be called for sure. SampleS 328 // simulation. 329 // 330 // SI : XS must not be zero otherwise samp 331 if (ekin < fKillBelowEnergy) { 332 return DBL_MAX; 333 } 334 335 // Get the tables with the cross section d 336 auto tableData = fpModelData->GetData(); 337 if ((*tableData)[materialID][p] == nullptr 338 G4Exception("G4DNAPTBElasticModel::Cross 339 "No model is registered"); 340 } 341 // Retrieve the cross section value 342 sigma = (*tableData)[materialID][p]->FindV 343 } 344 345 if (verboseLevel > 2) { 346 G4cout << "_______________________________ 347 G4cout << "°°° G4DNAPTBElasticModel - X 348 G4cout << "°°° Kinetic energy(eV)=" << 349 G4cout << "°°° Cross section per molecu 350 G4cout << "°°° G4DNAPTBElasticModel - X 351 } 352 353 // Return the cross section 354 auto MolDensity = 355 (*G4DNAMolecularMaterial::Instance()->GetN 356 return sigma * MolDensity; 357 } 358 359 //....oooOO0OOooo........oooOO0OOooo........oo 360 361 void G4DNAPTBElasticModel::SampleSecondaries(s 362 c 363 c 364 G 365 { 366 if (verboseLevel > 3) { 367 G4cout << "Calling SampleSecondaries() of 368 } 369 370 G4double electronEnergy0 = aDynamicElectron- 371 const std::size_t& materialID = couple->GetM 372 auto p = aDynamicElectron->GetParticleDefini 373 374 // set killBelowEnergy value for material 375 fKillBelowEnergy = fpModelData->GetLowELimit 376 377 // If the particle (electron here) energy is 378 // simulation 379 if (electronEnergy0 < fKillBelowEnergy) { 380 fParticleChangeForGamma->SetProposedKineti 381 fParticleChangeForGamma->ProposeTrackStatu 382 fParticleChangeForGamma->ProposeLocalEnerg 383 } 384 // If we are above the kill limite and below 385 else if (electronEnergy0 >= fKillBelowEnergy 386 // Random sampling of the cosTheta 387 G4double cosTheta = fpModelData->Randomize 388 389 // Random sampling of phi 390 G4double phi = 2. * CLHEP::pi * G4UniformR 391 392 auto zVers = aDynamicElectron->GetMomentum 393 auto xVers = zVers.orthogonal(); 394 auto yVers = zVers.cross(xVers); 395 396 G4double xDir = std::sqrt(1. - cosTheta * 397 G4double yDir = xDir; 398 xDir *= std::cos(phi); 399 yDir *= std::sin(phi); 400 401 // Particle direction after ModelInterface 402 G4ThreeVector zPrikeVers((xDir * xVers + y 403 404 // Give the new direction 405 fParticleChangeForGamma->ProposeMomentumDi 406 407 // Update the energy which does not change 408 fParticleChangeForGamma->SetProposedKineti 409 } 410 } 411 412 //....oooOO0OOooo........oooOO0OOooo........oo 413 414 G4double G4DNAPTBElasticModel::Theta(const G4P 415 const std 416 { 417 G4double theta = 0.; 418 G4double valueT1 = 0; 419 G4double valueT2 = 0; 420 G4double valueE21 = 0; 421 G4double valueE22 = 0; 422 G4double valueE12 = 0; 423 G4double valueE11 = 0; 424 G4double xs11 = 0; 425 G4double xs12 = 0; 426 G4double xs21 = 0; 427 G4double xs22 = 0; 428 if (p == G4Electron::ElectronDefinition()) { 429 auto t2 = 430 std::upper_bound(tValuesVec[materialID][ 431 auto t1 = t2 - 1; 432 433 auto e12 = std::upper_bound(eValuesVect[ma 434 eValuesVect[ma 435 auto e11 = e12 - 1; 436 437 auto e22 = std::upper_bound(eValuesVect[ma 438 eValuesVect[ma 439 auto e21 = e22 - 1; 440 441 valueT1 = *t1; 442 valueT2 = *t2; 443 valueE21 = *e21; 444 valueE22 = *e22; 445 valueE12 = *e12; 446 valueE11 = *e11; 447 448 xs11 = diffCrossSectionData[materialID][p] 449 xs12 = diffCrossSectionData[materialID][p] 450 xs21 = diffCrossSectionData[materialID][p] 451 xs22 = diffCrossSectionData[materialID][p] 452 } 453 454 if (xs11 == 0 && xs12 == 0 && xs21 == 0 && x 455 return (0.); 456 } 457 458 theta = QuadInterpolator(valueE11, valueE12, 459 valueT2, k, integrD 460 461 return theta; 462 } 463 464 //....oooOO0OOooo........oooOO0OOooo........oo 465 466 G4double G4DNAPTBElasticModel::LinLogInterpola 467 468 { 469 G4double d1 = std::log(xs1); 470 G4double d2 = std::log(xs2); 471 G4double value = std::exp(d1 + (d2 - d1) * ( 472 return value; 473 } 474 475 //....oooOO0OOooo........oooOO0OOooo........oo 476 477 G4double G4DNAPTBElasticModel::LinLinInterpola 478 479 { 480 G4double d1 = xs1; 481 G4double d2 = xs2; 482 G4double value = (d1 + (d2 - d1) * (e - e1) 483 return value; 484 } 485 486 //....oooOO0OOooo........oooOO0OOooo........oo 487 488 G4double G4DNAPTBElasticModel::LogLogInterpola 489 490 { 491 G4double a = (std::log10(xs2) - std::log10(x 492 G4double b = std::log10(xs2) - a * std::log1 493 G4double sigma = a * std::log10(e) + b; 494 G4double value = (std::pow(10., sigma)); 495 return value; 496 } 497 498 //....oooOO0OOooo........oooOO0OOooo........oo 499 500 G4double G4DNAPTBElasticModel::QuadInterpolato 501 502 503 504 { 505 // Log-Log 506 /* 507 G4double interpolatedvalue1 = LogLogInterpo 508 G4double interpolatedvalue2 = LogLogInterpo 509 G4double value = LogLogInterpolate(t1, t2, 510 511 512 // Lin-Log 513 G4double interpolatedvalue1 = LinLogInterpo 514 G4double interpolatedvalue2 = LinLogInterpo 515 G4double value = LinLogInterpolate(t1, t2, 516 */ 517 518 // Lin-Lin 519 G4double interpolatedvalue1 = LinLinInterpol 520 G4double interpolatedvalue2 = LinLinInterpol 521 G4double value = LinLinInterpolate(t1, t2, t 522 523 return value; 524 } 525 526 //....oooOO0OOooo........oooOO0OOooo........oo 527 528 G4double G4DNAPTBElasticModel::RandomizeCosThe 529 { 530 G4double integrdiff = 0; 531 G4double uniformRand = G4UniformRand(); 532 integrdiff = uniformRand; 533 534 G4double theta = 0.; 535 G4double cosTheta = 0.; 536 theta = Theta(G4Electron::ElectronDefinition 537 538 cosTheta = std::cos(theta * CLHEP::pi / 180) 539 540 return cosTheta; 541 } 542