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 // CPA100 elastic model class for electrons 27 // 28 // Based on the work of M. Terrissol and M. C. 29 // 30 // Users are requested to cite the following p 31 // - M. Terrissol, A. Baudre, Radiat. Prot. Do 32 // - M.C. Bordage, J. Bordes, S. Edel, M. Terr 33 // M. Bardies, N. Lampe, S. Incerti, Phys. M 34 // 35 // Authors of this class: 36 // M.C. Bordage, M. Terrissol, S. Edel, J. Bor 37 // 38 // 15.01.2014: creation 39 // 40 // Based on the study by S. Zein et. al. Nucl. 41 // 1/2/2023 : Hoang added modification 42 43 #include "G4DNACPA100ElasticModel.hh" 44 45 #include "G4DNAMaterialManager.hh" 46 #include "G4DNAMolecularMaterial.hh" 47 #include "G4EnvironmentUtils.hh" 48 #include "G4SystemOfUnits.hh" 49 using namespace std; 50 G4DNACPA100ElasticModel::G4DNACPA100ElasticMod 51 : G4VDNAModel(nam, "all") 52 { 53 fpGuanine = G4Material::GetMaterial("G4_GUAN 54 fpG4_WATER = G4Material::GetMaterial("G4_WAT 55 fpDeoxyribose = G4Material::GetMaterial("G4_ 56 fpCytosine = G4Material::GetMaterial("G4_CYT 57 fpThymine = G4Material::GetMaterial("G4_THYM 58 fpAdenine = G4Material::GetMaterial("G4_ADEN 59 fpPhosphate = G4Material::GetMaterial("G4_PH 60 fpParticle = G4Electron::ElectronDefinition( 61 } 62 63 void G4DNACPA100ElasticModel::Initialise(const 64 const 65 { 66 if (isInitialised) { 67 return; 68 } 69 if (verboseLevel > 3) { 70 G4cout << "Calling G4DNACPA100ExcitationMo 71 } 72 73 if (!G4DNAMaterialManager::Instance()->IsLoc 74 if (p != fpParticle) { 75 std::ostringstream oss; 76 oss << " Model is not applied for this p 77 G4Exception("G4DNACPA100ElasticModel::G4 78 oss.str().c_str()); 79 } 80 81 const char* path = G4FindDataDir("G4LEDATA 82 83 if (path == nullptr) { 84 G4Exception("G4DNACPA100ElasticModel::In 85 "G4LEDATA environment variab 86 return; 87 } 88 89 std::size_t index; 90 if (fpG4_WATER != nullptr) { 91 index = fpG4_WATER->GetIndex(); 92 fLevels[index] = 1.214e-4; 93 AddCrossSectionData(index, p, "dna/sigma 94 "dna/sigmadiff_cumul 95 SetLowELimit(index, p, 11. * eV); 96 SetHighELimit(index, p, 255955. * eV); 97 } 98 if (fpGuanine != nullptr) { 99 index = fpGuanine->GetIndex(); 100 fLevels[index] = 1.4504480e-05; 101 AddCrossSectionData(index, p, "dna/sigma 102 "dna/sigmadiff_cumul 103 SetLowELimit(index, p, 11 * eV); 104 SetHighELimit(index, p, 1 * MeV); 105 } 106 if (fpDeoxyribose != nullptr) { 107 index = fpDeoxyribose->GetIndex(); 108 fLevels[index] = 1.6343100e-05; 109 AddCrossSectionData(index, p, "dna/sigma 110 "dna/sigmadiff_cumul 111 SetLowELimit(index, p, 11 * eV); 112 SetHighELimit(index, p, 1 * MeV); 113 } 114 if (fpCytosine != nullptr) { 115 index = fpCytosine->GetIndex(); 116 fLevels[index] = 1.9729660e-05; 117 AddCrossSectionData(index, p, "dna/sigma 118 "dna/sigmadiff_cumul 119 SetLowELimit(index, p, 11 * eV); 120 SetHighELimit(index, p, 1 * MeV); 121 } 122 if (fpThymine != nullptr) { 123 index = fpThymine->GetIndex(); 124 fLevels[index] = 1.7381300e-05; 125 AddCrossSectionData(index, p, "dna/sigma 126 "dna/sigmadiff_cumul 127 SetLowELimit(index, p, 11 * eV); 128 SetHighELimit(index, p, 1 * MeV); 129 } 130 if (fpAdenine != nullptr) { 131 index = fpAdenine->GetIndex(); 132 fLevels[index] = 1.6221800e-05; 133 AddCrossSectionData(index, p, "dna/sigma 134 "dna/sigmadiff_cumul 135 SetLowELimit(index, p, 11 * eV); 136 SetHighELimit(index, p, 1 * MeV); 137 } 138 if (fpPhosphate != nullptr) { 139 index = fpPhosphate->GetIndex(); 140 fLevels[index] = 2.2369600e-05; 141 AddCrossSectionData(index, p, "dna/sigma 142 "dna/sigmadiff_cumul 143 SetLowELimit(index, p, 11 * eV); 144 SetHighELimit(index, p, 1 * MeV); 145 } 146 147 // Load data 148 LoadCrossSectionData(p); 149 G4DNAMaterialManager::Instance()->SetMaste 150 fpModelData = this; 151 } 152 else { 153 auto dataModel = dynamic_cast<G4DNACPA100E 154 G4DNAMaterialManager::Instance()->GetMod 155 if (dataModel == nullptr) { 156 G4cout << "G4DNACPA100ElasticModel::Cros 157 G4Exception("G4DNACPA100ElasticModel::Cr 158 "no modelData is registered" 159 } 160 else { 161 fpModelData = dataModel; 162 } 163 } 164 165 fParticleChangeForGamma = GetParticleChangeF 166 isInitialised = true; 167 } 168 169 //....oooOO0OOooo........oooOO0OOooo........oo 170 171 G4double G4DNACPA100ElasticModel::CrossSection 172 173 174 { 175 // Get the name of the current particle 176 const G4String& particleName = p->GetParticl 177 auto materialID = pMaterial->GetIndex(); 178 179 // set killBelowEnergy value for current mat 180 fKillBelowEnergy = fpModelData->GetLowELimit 181 182 G4double sigma = 0.; 183 184 if (ekin < fpModelData->GetHighELimit(materi 185 if (ekin < fKillBelowEnergy) { 186 return DBL_MAX; 187 } 188 189 auto tableData = fpModelData->GetData(); 190 191 if ((*tableData)[materialID][p] == nullptr 192 G4Exception("G4DNACPA100ElasticModel::Cr 193 "No model is registered"); 194 } 195 sigma = (*tableData)[materialID][p]->FindV 196 } 197 198 if (verboseLevel > 2) { 199 auto MolDensity = 200 (*G4DNAMolecularMaterial::Instance()->Ge 201 G4cout << "_______________________________ 202 G4cout << "°°° G4DNACPA100ElasticModel 203 G4cout << "°°° Kinetic energy(eV)=" << 204 G4cout << "°°° lowLim (eV) = " << GetLo 205 << " highLim (eV) : " << GetHighELi 206 G4cout << "°°° Materials = " << (*G4Mat 207 << G4endl; 208 G4cout << "°°° Cross section per molecu 209 G4cout << "°°° Cross section per Phosph 210 << G4endl; 211 G4cout << "°°° G4DNACPA100ElasticModel 212 } 213 214 auto MolDensity = 215 (*G4DNAMolecularMaterial::Instance()->GetN 216 return sigma * MolDensity; 217 } 218 219 //....oooOO0OOooo........oooOO0OOooo........oo 220 221 void G4DNACPA100ElasticModel::SampleSecondarie 222 223 224 225 { 226 G4double electronEnergy0 = aDynamicElectron- 227 auto materialID = couple->GetMaterial()->Get 228 auto p = aDynamicElectron->GetParticleDefini 229 230 if (p != fpParticle) { 231 G4Exception("G4DNACPA100ElasticModel::Samp 232 "This particle is not applied 233 } 234 if (electronEnergy0 < fKillBelowEnergy) { 235 return; 236 } 237 G4double cosTheta = fpModelData->RandomizeCo 238 G4double phi = 2. * CLHEP::pi * G4UniformRan 239 240 const G4ThreeVector& zVers = aDynamicElectro 241 242 G4double CT1, ST1, CF1, SF1, CT2, ST2, CF2, 243 G4double sinTheta = std::sqrt(1 - cosTheta * 244 245 CT1 = zVers.z(); 246 ST1 = std::sqrt(1. - CT1 * CT1); 247 248 if (ST1 != 0) 249 CF1 = zVers.x() / ST1; 250 else 251 CF1 = std::cos(2. * CLHEP::pi * G4UniformR 252 if (ST1 != 0) 253 SF1 = zVers.y() / ST1; 254 else 255 SF1 = std::sqrt(1. - CF1 * CF1); 256 257 G4double A3, A4, A5, A2, A1; 258 259 A3 = sinTheta * std::cos(phi); 260 A4 = A3 * CT1 + ST1 * cosTheta; 261 A5 = sinTheta * std::sin(phi); 262 A2 = A4 * SF1 + A5 * CF1; 263 A1 = A4 * CF1 - A5 * SF1; 264 265 CT2 = CT1 * cosTheta - ST1 * A3; 266 ST2 = std::sqrt(1. - CT2 * CT2); 267 268 if (ST2 == 0) ST2 = 1E-6; 269 CF2 = A1 / ST2; 270 SF2 = A2 / ST2; 271 G4ThreeVector zPrimeVers(ST2 * CF2, ST2 * SF 272 273 fParticleChangeForGamma->ProposeMomentumDire 274 275 auto EnergyDeposit = fpModelData->GetElastic 276 fParticleChangeForGamma->ProposeLocalEnergyD 277 if (statCode) { 278 fParticleChangeForGamma->SetProposedKineti 279 } 280 else { 281 auto newEnergy = electronEnergy0 - EnergyD 282 fParticleChangeForGamma->SetProposedKineti 283 } 284 } 285 286 G4double G4DNACPA100ElasticModel::Theta(const 287 G4doub 288 { 289 G4double theta, valueT1, valueT2, valueE21, 290 G4double xs11 = 0; 291 G4double xs12 = 0; 292 G4double xs21 = 0; 293 G4double xs22 = 0; 294 if (p == G4Electron::ElectronDefinition()) { 295 if (k == tValuesVec[materialID][p].back()) 296 k = k * (1. - 1e-12); 297 } 298 auto t2 = 299 std::upper_bound(tValuesVec[materialID][ 300 auto t1 = t2 - 1; 301 302 auto e12 = std::upper_bound(eValuesVect[ma 303 eValuesVect[ma 304 auto e11 = e12 - 1; 305 306 auto e22 = std::upper_bound(eValuesVect[ma 307 eValuesVect[ma 308 auto e21 = e22 - 1; 309 310 valueT1 = *t1; 311 valueT2 = *t2; 312 valueE21 = *e21; 313 valueE22 = *e22; 314 valueE12 = *e12; 315 valueE11 = *e11; 316 317 xs11 = diffCrossSectionData[materialID][p] 318 xs12 = diffCrossSectionData[materialID][p] 319 xs21 = diffCrossSectionData[materialID][p] 320 xs22 = diffCrossSectionData[materialID][p] 321 } 322 323 if (xs11 == 0 && xs12 == 0 && xs21 == 0 && x 324 return (0.); 325 } 326 327 theta = QuadInterpolator(valueE11, valueE12, 328 valueT2, k, integrD 329 330 return theta; 331 } 332 333 //....oooOO0OOooo........oooOO0OOooo........oo 334 335 G4double G4DNACPA100ElasticModel::LinLogInterp 336 337 { 338 G4double d1 = std::log(xs1); 339 G4double d2 = std::log(xs2); 340 G4double value = std::exp(d1 + (d2 - d1) * ( 341 return value; 342 } 343 344 //....oooOO0OOooo........oooOO0OOooo........oo 345 346 G4double G4DNACPA100ElasticModel::LinLinInterp 347 348 { 349 G4double d1 = xs1; 350 G4double d2 = xs2; 351 G4double value = (d1 + (d2 - d1) * (e - e1) 352 return value; 353 } 354 355 //....oooOO0OOooo........oooOO0OOooo........oo 356 357 G4double G4DNACPA100ElasticModel::LogLogInterp 358 359 { 360 G4double a = (std::log10(xs2) - std::log10(x 361 G4double b = std::log10(xs2) - a * std::log1 362 G4double sigma = a * std::log10(e) + b; 363 G4double value = (std::pow(10., sigma)); 364 return value; 365 } 366 367 //....oooOO0OOooo........oooOO0OOooo........oo 368 369 G4double G4DNACPA100ElasticModel::QuadInterpol 370 371 372 373 { 374 // Log-Log 375 /* 376 G4double interpolatedvalue1 = LogLogInterp 377 G4double interpolatedvalue2 = LogLogInterp 378 G4double value = LogLogInterpolate(t1, t2, 379 380 381 // Lin-Log 382 G4double interpolatedvalue1 = LinLogInterp 383 G4double interpolatedvalue2 = LinLogInterp 384 G4double value = LinLogInterpolate(t1, t2, 385 */ 386 387 // Lin-Lin 388 G4double interpolatedvalue1 = LinLinInterpol 389 G4double interpolatedvalue2 = LinLinInterpol 390 G4double value = LinLinInterpolate(t1, t2, t 391 392 return value; 393 } 394 395 //....oooOO0OOooo........oooOO0OOooo........oo 396 397 G4double G4DNACPA100ElasticModel::RandomizeCos 398 { 399 G4double integrdiff = 0; // PROBABILITY bet 400 G4double uniformRand = G4UniformRand(); 401 integrdiff = uniformRand; 402 G4double cosTheta = 0.; 403 cosTheta = 1 - Theta(G4Electron::ElectronDef 404 // cosTheta = std::cos(theta * CLHEP::pi / 1 405 return cosTheta; 406 } 407 //....oooOO0OOooo........oooOO0OOooo........oo 408 409 void G4DNACPA100ElasticModel::ReadDiffCSFile(c 410 c 411 c 412 { 413 const char* path = G4FindDataDir("G4LEDATA") 414 if (path == nullptr) { 415 G4Exception("G4DNACPA100ElasticModel::Read 416 "G4LEDATA environment variable 417 return; 418 } 419 420 std::ostringstream fullFileName; 421 fullFileName << path << "/" << file << ".dat 422 423 std::ifstream diffCrossSection(fullFileName. 424 // error if file is not there 425 std::stringstream endPath; 426 if (!diffCrossSection) { 427 endPath << "Missing data file: " << file; 428 G4Exception("G4DNACPA100ElasticModel::Init 429 endPath.str().c_str()); 430 } 431 432 tValuesVec[materialName][particleName].push_ 433 434 G4String line; 435 while (std::getline(diffCrossSection, line)) 436 // 437 std::istringstream testIss(line); 438 G4String test; 439 testIss >> test; 440 if (test == "#") { 441 continue; 442 } 443 // check if line is empty 444 if (line.empty()) { 445 continue; 446 } 447 std::istringstream iss(line); 448 449 G4double tDummy; 450 G4double eDummy; 451 452 iss >> tDummy >> eDummy; 453 454 if (tDummy != tValuesVec[materialName][par 455 // Add the current T value 456 tValuesVec[materialName][particleName].p 457 // Make it correspond to a default zero 458 eValuesVect[materialName][particleName][ 459 } 460 iss >> diffCrossSectionData[materialName][ 461 462 if (eDummy != eValuesVect[materialName][pa 463 eValuesVect[materialName][particleName][ 464 } 465 } 466 } 467