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 // 26 // 27 27 28 #include "G4DNAChampionElasticModel.hh" 28 #include "G4DNAChampionElasticModel.hh" 29 #include "G4PhysicalConstants.hh" 29 #include "G4PhysicalConstants.hh" 30 #include "G4SystemOfUnits.hh" 30 #include "G4SystemOfUnits.hh" 31 #include "G4DNAMolecularMaterial.hh" 31 #include "G4DNAMolecularMaterial.hh" 32 #include "G4Exp.hh" 32 #include "G4Exp.hh" 33 33 34 //....oooOO0OOooo........oooOO0OOooo........oo 34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 35 35 36 using namespace std; 36 using namespace std; 37 37 38 #define CHAMPION_VERBOSE 38 #define CHAMPION_VERBOSE 39 39 40 //....oooOO0OOooo........oooOO0OOooo........oo 40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 41 41 42 G4DNAChampionElasticModel:: 42 G4DNAChampionElasticModel:: 43 G4DNAChampionElasticModel(const G4ParticleDefi 43 G4DNAChampionElasticModel(const G4ParticleDefinition*, const G4String& nam) : 44 G4VEmModel(nam) << 44 G4VEmModel(nam), isInitialised(false) 45 { 45 { 46 SetLowEnergyLimit(7.4 * eV); 46 SetLowEnergyLimit(7.4 * eV); 47 SetHighEnergyLimit(1. * MeV); 47 SetHighEnergyLimit(1. * MeV); 48 48 49 verboseLevel = 0; 49 verboseLevel = 0; 50 // Verbosity scale: 50 // Verbosity scale: 51 // 0 = nothing 51 // 0 = nothing 52 // 1 = warning for energy non-conservation 52 // 1 = warning for energy non-conservation 53 // 2 = details of energy budget 53 // 2 = details of energy budget 54 // 3 = calculation of cross sections, file o 54 // 3 = calculation of cross sections, file openings, sampling of atoms 55 // 4 = entering in methods 55 // 4 = entering in methods 56 56 57 #ifdef CHAMPION_VERBOSE 57 #ifdef CHAMPION_VERBOSE 58 if (verboseLevel > 0) 58 if (verboseLevel > 0) 59 { 59 { 60 G4cout << "Champion Elastic model is const 60 G4cout << "Champion Elastic model is constructed " 61 << G4endl 61 << G4endl 62 << "Energy range: " 62 << "Energy range: " 63 << LowEnergyLimit() / eV << " eV - 63 << LowEnergyLimit() / eV << " eV - " 64 << HighEnergyLimit() / MeV << " MeV 64 << HighEnergyLimit() / MeV << " MeV" 65 << G4endl; 65 << G4endl; 66 } 66 } 67 #endif 67 #endif 68 68 69 fParticleChangeForGamma = nullptr; << 69 fParticleChangeForGamma = 0; 70 fpMolWaterDensity = nullptr; << 70 fpMolWaterDensity = 0; 71 fpData = nullptr; << 71 fpData = 0; 72 } 72 } 73 73 74 //....oooOO0OOooo........oooOO0OOooo........oo 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 75 75 76 G4DNAChampionElasticModel::~G4DNAChampionElast 76 G4DNAChampionElasticModel::~G4DNAChampionElasticModel() 77 { 77 { 78 // For total cross section 78 // For total cross section 79 delete fpData; << 79 if(fpData) delete fpData; 80 80 81 // For final state 81 // For final state 82 eVecm.clear(); 82 eVecm.clear(); 83 } 83 } 84 84 85 //....oooOO0OOooo........oooOO0OOooo........oo 85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 86 86 87 void G4DNAChampionElasticModel::Initialise(con 87 void G4DNAChampionElasticModel::Initialise(const G4ParticleDefinition* particle, 88 con 88 const G4DataVector& /*cuts*/) 89 { 89 { 90 #ifdef CHAMPION_VERBOSE 90 #ifdef CHAMPION_VERBOSE 91 if (verboseLevel > 3) 91 if (verboseLevel > 3) 92 { 92 { 93 G4cout << "Calling G4DNAChampionElasticMod 93 G4cout << "Calling G4DNAChampionElasticModel::Initialise()" << G4endl; 94 } 94 } 95 #endif 95 #endif 96 96 97 if(particle->GetParticleName() != "e-") 97 if(particle->GetParticleName() != "e-") 98 { 98 { 99 G4Exception("G4DNAChampionElasticModel::In 99 G4Exception("G4DNAChampionElasticModel::Initialise", 100 "em0002", 100 "em0002", 101 FatalException, 101 FatalException, 102 "Model not applicable to parti 102 "Model not applicable to particle type."); 103 } 103 } 104 104 105 // Energy limits 105 // Energy limits 106 106 107 if (LowEnergyLimit() < 7.4*eV) 107 if (LowEnergyLimit() < 7.4*eV) 108 { 108 { 109 G4cout << "G4DNAChampionElasticModel: low 109 G4cout << "G4DNAChampionElasticModel: low energy limit increased from " 110 << LowEnergyLimit()/eV << " eV to " 110 << LowEnergyLimit()/eV << " eV to " << 7.4 << " eV" 111 << G4endl; 111 << G4endl; 112 SetLowEnergyLimit(7.4*eV); 112 SetLowEnergyLimit(7.4*eV); 113 } 113 } 114 114 115 if (HighEnergyLimit() > 1.*MeV) 115 if (HighEnergyLimit() > 1.*MeV) 116 { 116 { 117 G4cout << "G4DNAChampionElasticModel: high 117 G4cout << "G4DNAChampionElasticModel: high energy limit decreased from " 118 << HighEnergyLimit()/MeV << " MeV t 118 << HighEnergyLimit()/MeV << " MeV to " << 1. << " MeV" 119 << G4endl; 119 << G4endl; 120 SetHighEnergyLimit(1.*MeV); 120 SetHighEnergyLimit(1.*MeV); 121 } 121 } 122 122 123 if (isInitialised) { return; } 123 if (isInitialised) { return; } 124 124 125 // *** ELECTRON 125 // *** ELECTRON 126 // For total cross section 126 // For total cross section 127 // Reading of data files 127 // Reading of data files 128 128 129 G4double scaleFactor = 1e-16*cm*cm; 129 G4double scaleFactor = 1e-16*cm*cm; 130 130 131 G4String fileElectron("dna/sigma_elastic_e_c 131 G4String fileElectron("dna/sigma_elastic_e_champion"); 132 132 133 fpData = new G4DNACrossSectionDataSet(new G4 133 fpData = new G4DNACrossSectionDataSet(new G4LogLogInterpolation(), 134 eV, 134 eV, 135 scaleF 135 scaleFactor ); 136 fpData->LoadData(fileElectron); 136 fpData->LoadData(fileElectron); 137 137 138 // For final state 138 // For final state 139 139 140 const char *path = G4FindDataDir("G4LEDATA") 140 const char *path = G4FindDataDir("G4LEDATA"); 141 141 142 if (path == nullptr) << 142 if (!path) 143 { 143 { 144 G4Exception("G4ChampionElasticModel::Initi 144 G4Exception("G4ChampionElasticModel::Initialise", 145 "em0006", 145 "em0006", 146 FatalException, 146 FatalException, 147 "G4LEDATA environment variable 147 "G4LEDATA environment variable not set."); 148 return; 148 return; 149 } 149 } 150 150 151 std::ostringstream eFullFileName; 151 std::ostringstream eFullFileName; 152 eFullFileName << path << "/dna/sigmadiff_cum 152 eFullFileName << path << "/dna/sigmadiff_cumulated_elastic_e_champion_hp.dat"; 153 std::ifstream eDiffCrossSection(eFullFileNam 153 std::ifstream eDiffCrossSection(eFullFileName.str().c_str()); 154 154 155 if (!eDiffCrossSection) 155 if (!eDiffCrossSection) 156 { 156 { 157 G4ExceptionDescription errMsg; 157 G4ExceptionDescription errMsg; 158 errMsg << "Missing data file:/dna/sigmadif 158 errMsg << "Missing data file:/dna/sigmadiff_cumulated_elastic_e_champion_hp.dat; " 159 << "please use G4EMLOW6.36 and abov 159 << "please use G4EMLOW6.36 and above."; 160 160 161 G4Exception("G4DNAChampionElasticModel::In 161 G4Exception("G4DNAChampionElasticModel::Initialise", 162 "em0003", 162 "em0003", 163 FatalException, 163 FatalException, 164 errMsg); 164 errMsg); 165 } 165 } 166 166 167 // March 25th, 2014 - Vaclav Stepan, Sebasti 167 // March 25th, 2014 - Vaclav Stepan, Sebastien Incerti 168 // Added clear for MT 168 // Added clear for MT 169 169 170 eTdummyVec.clear(); 170 eTdummyVec.clear(); 171 eVecm.clear(); 171 eVecm.clear(); 172 eDiffCrossSectionData.clear(); 172 eDiffCrossSectionData.clear(); 173 173 174 // 174 // 175 175 176 eTdummyVec.push_back(0.); 176 eTdummyVec.push_back(0.); 177 177 178 while(!eDiffCrossSection.eof()) 178 while(!eDiffCrossSection.eof()) 179 { 179 { 180 G4double tDummy; 180 G4double tDummy; 181 G4double eDummy; 181 G4double eDummy; 182 eDiffCrossSection >> tDummy >> eDummy; 182 eDiffCrossSection >> tDummy >> eDummy; 183 183 184 // SI : mandatory eVecm initialization 184 // SI : mandatory eVecm initialization 185 185 186 if (tDummy != eTdummyVec.back()) 186 if (tDummy != eTdummyVec.back()) 187 { 187 { 188 eTdummyVec.push_back(tDummy); 188 eTdummyVec.push_back(tDummy); 189 eVecm[tDummy].push_back(0.); 189 eVecm[tDummy].push_back(0.); 190 } 190 } 191 191 192 eDiffCrossSection >> eDiffCrossSectionData 192 eDiffCrossSection >> eDiffCrossSectionData[tDummy][eDummy]; 193 193 194 if (eDummy != eVecm[tDummy].back()) eVecm[ 194 if (eDummy != eVecm[tDummy].back()) eVecm[tDummy].push_back(eDummy); 195 } 195 } 196 196 197 // End final state 197 // End final state 198 198 199 #ifdef CHAMPION_VERBOSE 199 #ifdef CHAMPION_VERBOSE 200 if (verboseLevel>0) 200 if (verboseLevel>0) 201 { 201 { 202 if (verboseLevel > 2) 202 if (verboseLevel > 2) 203 { 203 { 204 G4cout << "Loaded cross section files fo 204 G4cout << "Loaded cross section files for Champion Elastic model" << G4endl; 205 } 205 } 206 206 207 G4cout << "Champion Elastic model is initi 207 G4cout << "Champion Elastic model is initialized " << G4endl 208 << "Energy range: " 208 << "Energy range: " 209 << LowEnergyLimit() / eV << " eV - 209 << LowEnergyLimit() / eV << " eV - " 210 << HighEnergyLimit() / MeV << " MeV 210 << HighEnergyLimit() / MeV << " MeV" 211 << G4endl; 211 << G4endl; 212 } 212 } 213 #endif 213 #endif 214 214 215 // Initialize water density pointer 215 // Initialize water density pointer 216 G4DNAMolecularMaterial::Instance()->Initiali 216 G4DNAMolecularMaterial::Instance()->Initialize(); 217 217 218 fpMolWaterDensity = G4DNAMolecularMaterial:: 218 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()-> 219 GetNumMolPerVolTableFor(G4Material::GetMat 219 GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER")); 220 220 221 fParticleChangeForGamma = GetParticleChangeF 221 fParticleChangeForGamma = GetParticleChangeForGamma(); 222 isInitialised = true; 222 isInitialised = true; 223 223 224 } 224 } 225 225 226 //....oooOO0OOooo........oooOO0OOooo........oo 226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 227 227 228 G4double 228 G4double 229 G4DNAChampionElasticModel:: 229 G4DNAChampionElasticModel:: 230 CrossSectionPerVolume(const G4Material* materi 230 CrossSectionPerVolume(const G4Material* material, 231 #ifdef CHAMPION_VERBOSE 231 #ifdef CHAMPION_VERBOSE 232 const G4ParticleDefiniti 232 const G4ParticleDefinition* p, 233 #else 233 #else 234 const G4ParticleDefiniti 234 const G4ParticleDefinition*, 235 #endif 235 #endif 236 G4double ekin, 236 G4double ekin, 237 G4double, 237 G4double, 238 G4double) 238 G4double) 239 { 239 { 240 #ifdef CHAMPION_VERBOSE 240 #ifdef CHAMPION_VERBOSE 241 if (verboseLevel > 3) 241 if (verboseLevel > 3) 242 { 242 { 243 G4cout << "Calling CrossSectionPerVolume() 243 G4cout << "Calling CrossSectionPerVolume() of G4DNAChampionElasticModel" 244 << G4endl; 244 << G4endl; 245 } 245 } 246 #endif 246 #endif 247 247 248 // Calculate total cross section for model 248 // Calculate total cross section for model 249 249 250 G4double sigma = 0.; 250 G4double sigma = 0.; 251 G4double waterDensity = (*fpMolWaterDensity) 251 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()]; 252 252 253 if (ekin <= HighEnergyLimit() && ekin >= Low 253 if (ekin <= HighEnergyLimit() && ekin >= LowEnergyLimit()) 254 { 254 { 255 //SI : XS must not be zero otherwise sam 255 //SI : XS must not be zero otherwise sampling of secondaries method ignored 256 // 256 // 257 sigma = fpData->FindValue(ekin); 257 sigma = fpData->FindValue(ekin); 258 } 258 } 259 259 260 #ifdef CHAMPION_VERBOSE 260 #ifdef CHAMPION_VERBOSE 261 if (verboseLevel > 2) 261 if (verboseLevel > 2) 262 { 262 { 263 G4cout << "_______________________________ 263 G4cout << "__________________________________" << G4endl; 264 G4cout << "=== G4DNAChampionElasticModel - 264 G4cout << "=== G4DNAChampionElasticModel - XS INFO START" << G4endl; 265 G4cout << "=== Kinetic energy(eV)=" << eki 265 G4cout << "=== Kinetic energy(eV)=" << ekin/eV << " particle : " << p->GetParticleName() << G4endl; 266 G4cout << "=== Cross section per water mol 266 G4cout << "=== Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl; 267 G4cout << "=== Cross section per water mol 267 G4cout << "=== Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl; 268 G4cout << "=== G4DNAChampionElasticModel - 268 G4cout << "=== G4DNAChampionElasticModel - XS INFO END" << G4endl; 269 } 269 } 270 #endif 270 #endif 271 271 272 return sigma*waterDensity; 272 return sigma*waterDensity; 273 } 273 } 274 274 275 //....oooOO0OOooo........oooOO0OOooo........oo 275 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 276 276 277 void G4DNAChampionElasticModel::SampleSecondar 277 void G4DNAChampionElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/, 278 278 const G4MaterialCutsCouple* /*couple*/, 279 279 const G4DynamicParticle* aDynamicElectron, 280 280 G4double, 281 281 G4double) 282 { 282 { 283 283 284 #ifdef CHAMPION_VERBOSE 284 #ifdef CHAMPION_VERBOSE 285 if (verboseLevel > 3) 285 if (verboseLevel > 3) 286 { 286 { 287 G4cout << "Calling SampleSecondaries() of 287 G4cout << "Calling SampleSecondaries() of G4DNAChampionElasticModel" << G4endl; 288 } 288 } 289 #endif 289 #endif 290 290 291 G4double electronEnergy0 = aDynamicElectron- 291 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy(); 292 292 293 G4double cosTheta = RandomizeCosTheta(electr 293 G4double cosTheta = RandomizeCosTheta(electronEnergy0); 294 294 295 G4double phi = 2. * pi * G4UniformRand(); 295 G4double phi = 2. * pi * G4UniformRand(); 296 296 297 G4ThreeVector zVers = aDynamicElectron->GetM 297 G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection(); 298 G4ThreeVector xVers = zVers.orthogonal(); 298 G4ThreeVector xVers = zVers.orthogonal(); 299 G4ThreeVector yVers = zVers.cross(xVers); 299 G4ThreeVector yVers = zVers.cross(xVers); 300 300 301 G4double xDir = std::sqrt(1. - cosTheta*cosT 301 G4double xDir = std::sqrt(1. - cosTheta*cosTheta); 302 G4double yDir = xDir; 302 G4double yDir = xDir; 303 xDir *= std::cos(phi); 303 xDir *= std::cos(phi); 304 yDir *= std::sin(phi); 304 yDir *= std::sin(phi); 305 305 306 G4ThreeVector zPrimeVers((xDir*xVers + yDir* 306 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers)); 307 307 308 fParticleChangeForGamma->ProposeMomentumDire 308 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit()); 309 309 310 fParticleChangeForGamma->SetProposedKineticE 310 fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0); 311 311 312 } 312 } 313 313 314 //....oooOO0OOooo........oooOO0OOooo........oo 314 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 315 315 316 G4double G4DNAChampionElasticModel::Theta(G4do 316 G4double G4DNAChampionElasticModel::Theta(G4double k, 317 G4do 317 G4double integrDiff) 318 { 318 { 319 G4double theta = 0.; 319 G4double theta = 0.; 320 G4double valueT1 = 0; 320 G4double valueT1 = 0; 321 G4double valueT2 = 0; 321 G4double valueT2 = 0; 322 G4double valueE21 = 0; 322 G4double valueE21 = 0; 323 G4double valueE22 = 0; 323 G4double valueE22 = 0; 324 G4double valueE12 = 0; 324 G4double valueE12 = 0; 325 G4double valueE11 = 0; 325 G4double valueE11 = 0; 326 G4double xs11 = 0; 326 G4double xs11 = 0; 327 G4double xs12 = 0; 327 G4double xs12 = 0; 328 G4double xs21 = 0; 328 G4double xs21 = 0; 329 G4double xs22 = 0; 329 G4double xs22 = 0; 330 330 331 auto t2 = std::upper_bound(eTdummyVec.begin( << 331 std::vector<G4double>::iterator t2 = std::upper_bound(eTdummyVec.begin(), 332 332 eTdummyVec.end(), k); 333 auto t1 = t2 - 1; << 333 std::vector<G4double>::iterator t1 = t2 - 1; 334 334 335 auto e12 = std::upper_bound(eVecm[(*t1)].beg << 335 std::vector<G4double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(), 336 336 eVecm[(*t1)].end(), 337 337 integrDiff); 338 auto e11 = e12 - 1; << 338 std::vector<G4double>::iterator e11 = e12 - 1; 339 339 340 auto e22 = std::upper_bound(eVecm[(*t2)].beg << 340 std::vector<G4double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(), 341 341 eVecm[(*t2)].end(), 342 342 integrDiff); 343 auto e21 = e22 - 1; << 343 std::vector<G4double>::iterator e21 = e22 - 1; 344 344 345 valueT1 = *t1; 345 valueT1 = *t1; 346 valueT2 = *t2; 346 valueT2 = *t2; 347 valueE21 = *e21; 347 valueE21 = *e21; 348 valueE22 = *e22; 348 valueE22 = *e22; 349 valueE12 = *e12; 349 valueE12 = *e12; 350 valueE11 = *e11; 350 valueE11 = *e11; 351 351 352 xs11 = eDiffCrossSectionData[valueT1][valueE 352 xs11 = eDiffCrossSectionData[valueT1][valueE11]; 353 xs12 = eDiffCrossSectionData[valueT1][valueE 353 xs12 = eDiffCrossSectionData[valueT1][valueE12]; 354 xs21 = eDiffCrossSectionData[valueT2][valueE 354 xs21 = eDiffCrossSectionData[valueT2][valueE21]; 355 xs22 = eDiffCrossSectionData[valueT2][valueE 355 xs22 = eDiffCrossSectionData[valueT2][valueE22]; 356 356 357 if (xs11 == 0 && xs12 == 0 && xs21 == 0 && x 357 if (xs11 == 0 && xs12 == 0 && xs21 == 0 && xs22 == 0) return (0.); 358 358 359 theta = QuadInterpolator(valueE11, valueE12, 359 theta = QuadInterpolator(valueE11, valueE12, valueE21, valueE22, xs11, xs12, 360 xs21, xs22, valueT1 360 xs21, xs22, valueT1, valueT2, k, integrDiff); 361 361 362 return theta; 362 return theta; 363 } 363 } 364 364 365 //....oooOO0OOooo........oooOO0OOooo........oo 365 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 366 366 367 G4double G4DNAChampionElasticModel::LinLogInte 367 G4double G4DNAChampionElasticModel::LinLogInterpolate(G4double e1, 368 368 G4double e2, 369 369 G4double e, 370 370 G4double xs1, 371 371 G4double xs2) 372 { 372 { 373 G4double d1 = std::log(xs1); 373 G4double d1 = std::log(xs1); 374 G4double d2 = std::log(xs2); 374 G4double d2 = std::log(xs2); 375 G4double value = G4Exp(d1 + (d2 - d1) * (e - 375 G4double value = G4Exp(d1 + (d2 - d1) * (e - e1) / (e2 - e1)); 376 return value; 376 return value; 377 } 377 } 378 378 379 //....oooOO0OOooo........oooOO0OOooo........oo 379 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 380 380 381 G4double G4DNAChampionElasticModel::LinLinInte 381 G4double G4DNAChampionElasticModel::LinLinInterpolate(G4double e1, 382 382 G4double e2, 383 383 G4double e, 384 384 G4double xs1, 385 385 G4double xs2) 386 { 386 { 387 G4double d1 = xs1; 387 G4double d1 = xs1; 388 G4double d2 = xs2; 388 G4double d2 = xs2; 389 G4double value = (d1 + (d2 - d1) * (e - e1) 389 G4double value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1)); 390 return value; 390 return value; 391 } 391 } 392 392 393 //....oooOO0OOooo........oooOO0OOooo........oo 393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 394 394 395 G4double G4DNAChampionElasticModel::LogLogInte 395 G4double G4DNAChampionElasticModel::LogLogInterpolate(G4double e1, 396 396 G4double e2, 397 397 G4double e, 398 398 G4double xs1, 399 399 G4double xs2) 400 { 400 { 401 G4double a = (std::log10(xs2) - std::log10(x 401 G4double a = (std::log10(xs2) - std::log10(xs1)) 402 / (std::log10(e2) - std::log10(e1)); 402 / (std::log10(e2) - std::log10(e1)); 403 G4double b = std::log10(xs2) - a * std::log1 403 G4double b = std::log10(xs2) - a * std::log10(e2); 404 G4double sigma = a * std::log10(e) + b; 404 G4double sigma = a * std::log10(e) + b; 405 G4double value = (std::pow(10., sigma)); 405 G4double value = (std::pow(10., sigma)); 406 return value; 406 return value; 407 } 407 } 408 408 409 //....oooOO0OOooo........oooOO0OOooo........oo 409 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 410 410 411 G4double G4DNAChampionElasticModel::QuadInterp 411 G4double G4DNAChampionElasticModel::QuadInterpolator(G4double e11, 412 412 G4double e12, 413 413 G4double e21, 414 414 G4double e22, 415 415 G4double xs11, 416 416 G4double xs12, 417 417 G4double xs21, 418 418 G4double xs22, 419 419 G4double t1, 420 420 G4double t2, 421 421 G4double t, 422 422 G4double e) 423 { 423 { 424 // Log-Log 424 // Log-Log 425 /* 425 /* 426 G4double interpolatedvalue1 = LogLogInterpo 426 G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12); 427 G4double interpolatedvalue2 = LogLogInterpo 427 G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22); 428 G4double value = LogLogInterpolate(t1, t2, 428 G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2); 429 429 430 430 431 // Lin-Log 431 // Lin-Log 432 G4double interpolatedvalue1 = LinLogInterpo 432 G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12); 433 G4double interpolatedvalue2 = LinLogInterpo 433 G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22); 434 G4double value = LinLogInterpolate(t1, t2, 434 G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2); 435 */ 435 */ 436 436 437 // Lin-Lin 437 // Lin-Lin 438 G4double interpolatedvalue1 = LinLinInterpol 438 G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12); 439 G4double interpolatedvalue2 = LinLinInterpol 439 G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22); 440 G4double value = LinLinInterpolate(t1, t2, t 440 G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1, 441 interpola 441 interpolatedvalue2); 442 442 443 return value; 443 return value; 444 } 444 } 445 445 446 //....oooOO0OOooo........oooOO0OOooo........oo 446 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 447 447 448 G4double G4DNAChampionElasticModel::RandomizeC 448 G4double G4DNAChampionElasticModel::RandomizeCosTheta(G4double k) 449 { 449 { 450 450 451 G4double integrdiff = 0; 451 G4double integrdiff = 0; 452 G4double uniformRand = G4UniformRand(); 452 G4double uniformRand = G4UniformRand(); 453 integrdiff = uniformRand; 453 integrdiff = uniformRand; 454 454 455 G4double theta = 0.; 455 G4double theta = 0.; 456 G4double cosTheta = 0.; 456 G4double cosTheta = 0.; 457 theta = Theta(k / eV, integrdiff); 457 theta = Theta(k / eV, integrdiff); 458 458 459 cosTheta = std::cos(theta * pi / 180); 459 cosTheta = std::cos(theta * pi / 180); 460 460 461 return cosTheta; 461 return cosTheta; 462 } 462 } 463 463 464 //....oooOO0OOooo........oooOO0OOooo........oo 464 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 465 465 466 void G4DNAChampionElasticModel::SetKillBelowTh 466 void G4DNAChampionElasticModel::SetKillBelowThreshold(G4double) 467 { 467 { 468 G4ExceptionDescription errMsg; 468 G4ExceptionDescription errMsg; 469 errMsg << "The method G4DNAChampionElasticMo 469 errMsg << "The method G4DNAChampionElasticModel::SetKillBelowThreshold is deprecated"; 470 470 471 G4Exception("G4DNAChampionElasticModel::SetK 471 G4Exception("G4DNAChampionElasticModel::SetKillBelowThreshold", 472 "deprecated", 472 "deprecated", 473 JustWarning, 473 JustWarning, 474 errMsg); 474 errMsg); 475 } 475 } 476 476