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