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 // G4MicroElecElasticModel.cc, 2011/08/29 A.Va 27 // G4MicroElecElasticModel.cc, 2011/08/29 A.Valentin, M. Raine 28 // 28 // 29 // Based on the following publications 29 // Based on the following publications 30 // - Geant4 physics processes for microdo 30 // - Geant4 physics processes for microdosimetry simulation: 31 // very low energy electromagnetic models 31 // very low energy electromagnetic models for electrons in Si, 32 // NIM B, vol. 288, pp. 66 - 73, 2012. 32 // NIM B, vol. 288, pp. 66 - 73, 2012. 33 // 33 // 34 // 34 // 35 //....oooOO0OOooo........oooOO0OOooo........oo 35 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 36 36 37 37 38 #include "G4MicroElecElasticModel.hh" 38 #include "G4MicroElecElasticModel.hh" 39 #include "G4PhysicalConstants.hh" 39 #include "G4PhysicalConstants.hh" 40 #include "G4SystemOfUnits.hh" 40 #include "G4SystemOfUnits.hh" 41 #include "G4Exp.hh" 41 #include "G4Exp.hh" 42 42 43 //....oooOO0OOooo........oooOO0OOooo........oo 43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 44 44 45 using namespace std; 45 using namespace std; 46 46 47 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 48 48 49 G4MicroElecElasticModel::G4MicroElecElasticMod 49 G4MicroElecElasticModel::G4MicroElecElasticModel(const G4ParticleDefinition*, 50 const G4String& nam) << 50 const G4String& nam) 51 :G4VEmModel(nam),isInitialised(false) << 51 :G4VEmModel(nam),isInitialised(false) 52 { 52 { 53 nistSi = G4NistManager::Instance()->FindOrBu 53 nistSi = G4NistManager::Instance()->FindOrBuildMaterial("G4_Si"); 54 54 55 killBelowEnergy = 16.7 * eV; // Minimum e- e 55 killBelowEnergy = 16.7 * eV; // Minimum e- energy for energy loss by excitation 56 lowEnergyLimit = 0 * eV; 56 lowEnergyLimit = 0 * eV; 57 lowEnergyLimitOfModel = 5 * eV; // The model 57 lowEnergyLimitOfModel = 5 * eV; // The model lower energy is 5 eV 58 highEnergyLimit = 100. * MeV; 58 highEnergyLimit = 100. * MeV; 59 SetLowEnergyLimit(lowEnergyLimit); 59 SetLowEnergyLimit(lowEnergyLimit); 60 SetHighEnergyLimit(highEnergyLimit); 60 SetHighEnergyLimit(highEnergyLimit); 61 61 62 verboseLevel= 0; 62 verboseLevel= 0; 63 // Verbosity scale: 63 // Verbosity scale: 64 // 0 = nothing 64 // 0 = nothing 65 // 1 = warning for energy non-conservation 65 // 1 = warning for energy non-conservation 66 // 2 = details of energy budget 66 // 2 = details of energy budget 67 // 3 = calculation of cross sections, file o 67 // 3 = calculation of cross sections, file openings, sampling of atoms 68 // 4 = entering in methods 68 // 4 = entering in methods 69 69 70 if( verboseLevel>0 ) 70 if( verboseLevel>0 ) 71 { << 71 { 72 G4cout << "MicroElec Elastic model is co << 72 G4cout << "MicroElec Elastic model is constructed " << G4endl 73 << "Energy range: " << 73 << "Energy range: " 74 << lowEnergyLimit / eV << " eV - " << 74 << lowEnergyLimit / eV << " eV - " 75 << highEnergyLimit / MeV << " MeV" << 75 << highEnergyLimit / MeV << " MeV" 76 << G4endl; << 76 << G4endl; 77 } << 77 } 78 fParticleChangeForGamma = 0; 78 fParticleChangeForGamma = 0; 79 } 79 } 80 80 81 //....oooOO0OOooo........oooOO0OOooo........oo 81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 82 82 83 G4MicroElecElasticModel::~G4MicroElecElasticMo 83 G4MicroElecElasticModel::~G4MicroElecElasticModel() 84 { 84 { 85 // For total cross section 85 // For total cross section 86 for (auto & pos : tableData) << 87 { << 88 G4MicroElecCrossSectionDataSet* table = << 89 delete table; << 90 } << 91 86 92 // For final state << 87 std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator pos; 93 eVecm.clear(); << 88 for (pos = tableData.begin(); pos != tableData.end(); ++pos) >> 89 { >> 90 G4MicroElecCrossSectionDataSet* table = pos->second; >> 91 delete table; >> 92 } >> 93 >> 94 // For final state >> 95 >> 96 eVecm.clear(); >> 97 94 } 98 } 95 99 96 //....oooOO0OOooo........oooOO0OOooo........oo 100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 97 101 98 void G4MicroElecElasticModel::Initialise(const 102 void G4MicroElecElasticModel::Initialise(const G4ParticleDefinition* /*particle*/, 99 const G4DataVector& /*cuts*/) << 103 const G4DataVector& /*cuts*/) 100 { 104 { >> 105 101 if (verboseLevel > 3) 106 if (verboseLevel > 3) 102 G4cout << "Calling G4MicroElecElasticModel 107 G4cout << "Calling G4MicroElecElasticModel::Initialise()" << G4endl; 103 108 104 // Energy limits 109 // Energy limits >> 110 105 if (LowEnergyLimit() < lowEnergyLimit) 111 if (LowEnergyLimit() < lowEnergyLimit) 106 { << 112 { 107 G4cout << "G4MicroElecElasticModel: low << 113 G4cout << "G4MicroElecElasticModel: low energy limit increased from " << 108 LowEnergyLimit()/eV << " eV to " << lowEnerg 114 LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl; 109 SetLowEnergyLimit(lowEnergyLimit); << 115 SetLowEnergyLimit(lowEnergyLimit); 110 } 116 } 111 117 112 if (HighEnergyLimit() > highEnergyLimit) 118 if (HighEnergyLimit() > highEnergyLimit) 113 { << 119 { 114 G4cout << "G4MicroElecElasticModel: high << 120 G4cout << "G4MicroElecElasticModel: high energy limit decreased from " << 115 HighEnergyLimit()/MeV << " MeV to " << 121 HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl; 116 SetHighEnergyLimit(highEnergyLimit); << 122 SetHighEnergyLimit(highEnergyLimit); 117 } << 123 } 118 124 119 // Reading of data files 125 // Reading of data files 120 126 121 G4double scaleFactor = 1e-18 * cm * cm; 127 G4double scaleFactor = 1e-18 * cm * cm; >> 128 122 G4String fileElectron("microelec/sigma_elast 129 G4String fileElectron("microelec/sigma_elastic_e_Si"); 123 130 124 G4ParticleDefinition* electronDef = G4Electr 131 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition(); 125 G4String electron; 132 G4String electron; 126 133 127 // For total cross section << 134 // For total cross section 128 electron = electronDef->GetParticleName(); << 129 tableFile[electron] = fileElectron; << 130 135 131 G4MicroElecCrossSectionDataSet* tableE = new << 136 electron = electronDef->GetParticleName(); 132 tableE->LoadData(fileElectron); << 133 tableData[electron] = tableE; << 134 137 135 // For final state << 138 tableFile[electron] = fileElectron; 136 const char* path = G4FindDataDir("G4LEDATA") << 137 139 138 if (!path) << 140 G4MicroElecCrossSectionDataSet* tableE = new G4MicroElecCrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor ); >> 141 tableE->LoadData(fileElectron); >> 142 tableData[electron] = tableE; >> 143 >> 144 // For final state >> 145 >> 146 char *path = std::getenv("G4LEDATA"); >> 147 >> 148 if (!path) 139 { 149 { 140 G4Exception("G4MicroElecElasticModel::In 150 G4Exception("G4MicroElecElasticModel::Initialise","em0006",FatalException,"G4LEDATA environment variable not set."); 141 return; 151 return; 142 } 152 } 143 153 144 std::ostringstream eFullFileName; << 154 std::ostringstream eFullFileName; 145 eFullFileName << path << "/microelec/sigmadi << 155 eFullFileName << path << "/microelec/sigmadiff_cumulated_elastic_e_Si.dat"; 146 std::ifstream eDiffCrossSection(eFullFileNam << 156 std::ifstream eDiffCrossSection(eFullFileName.str().c_str()); 147 << 157 148 if (!eDiffCrossSection) << 158 if (!eDiffCrossSection) 149 G4Exception("G4MicroElecElasticModel::Init << 159 G4Exception("G4MicroElecElasticModel::Initialise","em0003",FatalException,"Missing data file: /microelec/sigmadiff_cumulated_elastic_e_Si.dat"); 150 "Missing data file: /microelec/sigmadiff_c << 160 151 << 161 152 // Added clear for MT << 162 // October 21th, 2014 - Melanie Raine 153 eTdummyVec.clear(); << 163 // Added clear for MT 154 eVecm.clear(); << 164 155 eDiffCrossSectionData.clear(); << 165 eTdummyVec.clear(); >> 166 eVecm.clear(); >> 167 eDiffCrossSectionData.clear(); 156 168 157 // << 169 // 158 eTdummyVec.push_back(0.); << 159 170 160 while(!eDiffCrossSection.eof()) << 171 >> 172 eTdummyVec.push_back(0.); >> 173 >> 174 while(!eDiffCrossSection.eof()) 161 { 175 { 162 double tDummy; << 176 double tDummy; 163 double eDummy; << 177 double eDummy; 164 eDiffCrossSection>>tDummy>>eDummy; << 178 eDiffCrossSection>>tDummy>>eDummy; 165 179 166 if (tDummy != eTdummyVec.back()) << 180 // SI : mandatory eVecm initialization >> 181 >> 182 if (tDummy != eTdummyVec.back()) 167 { 183 { 168 eTdummyVec.push_back(tDummy); 184 eTdummyVec.push_back(tDummy); 169 eVecm[tDummy].push_back(0.); 185 eVecm[tDummy].push_back(0.); 170 } 186 } 171 187 172 eDiffCrossSection>>eDiffCrossSectionData << 188 eDiffCrossSection>>eDiffCrossSectionData[tDummy][eDummy]; >> 189 >> 190 if (eDummy != eVecm[tDummy].back()) eVecm[tDummy].push_back(eDummy); 173 191 174 if (eDummy != eVecm[tDummy].back()) eVec << 175 } 192 } 176 // End final state << 193 >> 194 // End final state 177 195 178 if (verboseLevel > 2) 196 if (verboseLevel > 2) 179 G4cout << "Loaded cross section files for 197 G4cout << "Loaded cross section files for MicroElec Elastic model" << G4endl; 180 198 181 if( verboseLevel>0 ) 199 if( verboseLevel>0 ) 182 { << 200 { 183 G4cout << "MicroElec Elastic model is in << 201 G4cout << "MicroElec Elastic model is initialized " << G4endl 184 << "Energy range: " << 202 << "Energy range: " 185 << LowEnergyLimit() / eV << " eV - " << 203 << LowEnergyLimit() / eV << " eV - " 186 << HighEnergyLimit() / MeV << " MeV" << 204 << HighEnergyLimit() / MeV << " MeV" 187 << G4endl; << 205 << G4endl; 188 } << 206 } 189 207 190 if (isInitialised) { return; } 208 if (isInitialised) { return; } 191 fParticleChangeForGamma = GetParticleChangeF 209 fParticleChangeForGamma = GetParticleChangeForGamma(); 192 isInitialised = true; 210 isInitialised = true; >> 211 193 } 212 } 194 213 195 //....oooOO0OOooo........oooOO0OOooo........oo 214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 196 215 197 G4double G4MicroElecElasticModel::CrossSection 216 G4double G4MicroElecElasticModel::CrossSectionPerVolume(const G4Material* material, 198 const G4ParticleDefinition* p, << 217 const G4ParticleDefinition* p, 199 G4double ekin, << 218 G4double ekin, 200 G4double, << 219 G4double, 201 G4double) << 220 G4double) 202 { 221 { 203 if (verboseLevel > 3) 222 if (verboseLevel > 3) 204 G4cout << "Calling CrossSectionPerVolume() 223 G4cout << "Calling CrossSectionPerVolume() of G4MicroElecElasticModel" << G4endl; 205 224 206 // Calculate total cross section for model << 225 // Calculate total cross section for model 207 G4double sigma=0; << 208 G4double density = material->GetTotNbOfAtoms << 209 226 210 if (material == nistSi || material->GetBaseM << 227 G4double sigma=0; 211 { << 228 212 const G4String& particleName = p->GetPar << 229 G4double density = material->GetTotNbOfAtomsPerVolume(); >> 230 >> 231 if (material == nistSi || material->GetBaseMaterial() == nistSi) >> 232 { >> 233 const G4String& particleName = p->GetParticleName(); 213 234 214 if (ekin < highEnergyLimit) << 235 if (ekin < highEnergyLimit) >> 236 { >> 237 //SI : XS must not be zero otherwise sampling of secondaries method ignored >> 238 if (ekin < killBelowEnergy) return DBL_MAX; >> 239 // >> 240 >> 241 std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator pos; >> 242 pos = tableData.find(particleName); >> 243 >> 244 if (pos != tableData.end()) 215 { 245 { 216 //SI : XS must not be zero otherwise sampl << 246 G4MicroElecCrossSectionDataSet* table = pos->second; 217 if (ekin < killBelowEnergy) return DBL_MAX << 247 if (table != 0) 218 // << 248 { 219 << 249 sigma = table->FindValue(ekin); 220 auto pos = tableData.find(particleName); << 250 } 221 if (pos != tableData.end()) << 222 { << 223 G4MicroElecCrossSectionDataSet* table << 224 if (table != nullptr) << 225 { << 226 sigma = table->FindValue(ekin); << 227 } << 228 } << 229 else << 230 { << 231 G4Exception("G4MicroElecElasticModel:: << 232 FatalException,"Model not applicable t << 233 } << 234 } 251 } 235 << 252 else 236 if (verboseLevel > 3) << 237 { 253 { 238 G4cout << "---> Kinetic energy(eV)=" << ek << 254 G4Exception("G4MicroElecElasticModel::ComputeCrossSectionPerVolume","em0002",FatalException,"Model not applicable to particle type."); 239 G4cout << " - Cross section per Si atom (c << 240 G4cout << " - Cross section per Si atom (c << 241 } 255 } 242 } << 256 } 243 return sigma*density; << 257 >> 258 if (verboseLevel > 3) >> 259 { >> 260 G4cout << "---> Kinetic energy(eV)=" << ekin/eV << G4endl; >> 261 G4cout << " - Cross section per Si atom (cm^2)=" << sigma/cm/cm << G4endl; >> 262 G4cout << " - Cross section per Si atom (cm^-1)=" << sigma*density/(1./cm) << G4endl; >> 263 } >> 264 >> 265 } >> 266 >> 267 return sigma*density; 244 } 268 } 245 269 246 //....oooOO0OOooo........oooOO0OOooo........oo 270 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 247 271 248 void G4MicroElecElasticModel::SampleSecondarie 272 void G4MicroElecElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/, 249 const G4MaterialCutsCouple* /*coup << 273 const G4MaterialCutsCouple* /*couple*/, 250 const G4DynamicParticle* aDynamicE << 274 const G4DynamicParticle* aDynamicElectron, 251 G4double, << 275 G4double, 252 G4double) << 276 G4double) 253 { 277 { 254 278 255 if (verboseLevel > 3) 279 if (verboseLevel > 3) 256 G4cout << "Calling SampleSecondaries() of 280 G4cout << "Calling SampleSecondaries() of G4MicroElecElasticModel" << G4endl; 257 281 258 G4double electronEnergy0 = aDynamicElectron- 282 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy(); 259 283 260 if (electronEnergy0 < killBelowEnergy) 284 if (electronEnergy0 < killBelowEnergy) 261 { << 285 { 262 fParticleChangeForGamma->SetProposedKine << 286 fParticleChangeForGamma->SetProposedKineticEnergy(0.); 263 fParticleChangeForGamma->ProposeTrackSta << 287 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill); 264 fParticleChangeForGamma->ProposeLocalEne << 288 fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0); 265 return ; << 289 return ; 266 } << 290 } 267 291 268 if (electronEnergy0>= killBelowEnergy && ele 292 if (electronEnergy0>= killBelowEnergy && electronEnergy0 < highEnergyLimit) 269 { << 293 { 270 G4double cosTheta = RandomizeCosTheta(el << 294 G4double cosTheta = RandomizeCosTheta(electronEnergy0); 271 G4double phi = 2. * pi * G4UniformRand() << 272 G4ThreeVector zVers = aDynamicElectron-> << 273 G4ThreeVector xVers = zVers.orthogonal() << 274 G4ThreeVector yVers = zVers.cross(xVers) << 275 << 276 G4double xDir = std::sqrt(1. - cosTheta* << 277 G4double yDir = xDir; << 278 xDir *= std::cos(phi); << 279 yDir *= std::sin(phi); << 280 295 281 G4ThreeVector zPrimeVers((xDir*xVers + y << 296 G4double phi = 2. * pi * G4UniformRand(); >> 297 >> 298 G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection(); >> 299 G4ThreeVector xVers = zVers.orthogonal(); >> 300 G4ThreeVector yVers = zVers.cross(xVers); >> 301 >> 302 G4double xDir = std::sqrt(1. - cosTheta*cosTheta); >> 303 G4double yDir = xDir; >> 304 xDir *= std::cos(phi); >> 305 yDir *= std::sin(phi); >> 306 >> 307 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers)); >> 308 >> 309 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit()) ; >> 310 >> 311 fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0); >> 312 } 282 313 283 fParticleChangeForGamma->ProposeMomentum << 284 fParticleChangeForGamma->SetProposedKine << 285 } << 286 } 314 } 287 315 288 //....oooOO0OOooo........oooOO0OOooo........oo 316 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 289 317 290 G4double G4MicroElecElasticModel::Theta 318 G4double G4MicroElecElasticModel::Theta 291 (G4ParticleDefinition * particleDefinition, G4 << 319 (G4ParticleDefinition * particleDefinition, G4double k, G4double integrDiff) 292 { 320 { >> 321 293 G4double theta = 0.; 322 G4double theta = 0.; 294 G4double valueT1 = 0; 323 G4double valueT1 = 0; 295 G4double valueT2 = 0; 324 G4double valueT2 = 0; 296 G4double valueE21 = 0; 325 G4double valueE21 = 0; 297 G4double valueE22 = 0; 326 G4double valueE22 = 0; 298 G4double valueE12 = 0; 327 G4double valueE12 = 0; 299 G4double valueE11 = 0; 328 G4double valueE11 = 0; 300 G4double xs11 = 0; 329 G4double xs11 = 0; 301 G4double xs12 = 0; 330 G4double xs12 = 0; 302 G4double xs21 = 0; 331 G4double xs21 = 0; 303 G4double xs22 = 0; 332 G4double xs22 = 0; 304 333 >> 334 305 if (particleDefinition == G4Electron::Electr 335 if (particleDefinition == G4Electron::ElectronDefinition()) 306 { << 336 { 307 auto t2 = std::upper_bound(eTdummyVec.be << 337 std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k); 308 auto t1 = t2-1; << 338 std::vector<double>::iterator t1 = t2-1; 309 auto e12 = std::upper_bound(eVecm[(*t1)] << 339 310 auto e11 = e12-1; << 340 std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), integrDiff); 311 auto e22 = std::upper_bound(eVecm[(*t2)] << 341 std::vector<double>::iterator e11 = e12-1; 312 auto e21 = e22-1; << 342 313 << 343 std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), integrDiff); 314 valueT1 =*t1; << 344 std::vector<double>::iterator e21 = e22-1; 315 valueT2 =*t2; << 345 316 valueE21 =*e21; << 346 valueT1 =*t1; 317 valueE22 =*e22; << 347 valueT2 =*t2; 318 valueE12 =*e12; << 348 valueE21 =*e21; 319 valueE11 =*e11; << 349 valueE22 =*e22; 320 << 350 valueE12 =*e12; 321 xs11 = eDiffCrossSectionData[valueT1][va << 351 valueE11 =*e11; 322 xs12 = eDiffCrossSectionData[valueT1][va << 352 323 xs21 = eDiffCrossSectionData[valueT2][va << 353 xs11 = eDiffCrossSectionData[valueT1][valueE11]; 324 xs22 = eDiffCrossSectionData[valueT2][va << 354 xs12 = eDiffCrossSectionData[valueT1][valueE12]; 325 } << 355 xs21 = eDiffCrossSectionData[valueT2][valueE21]; >> 356 xs22 = eDiffCrossSectionData[valueT2][valueE22]; >> 357 >> 358 } >> 359 326 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) 360 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.); 327 361 328 theta = QuadInterpolator( valueE11, valueE1 362 theta = QuadInterpolator( valueE11, valueE12, 329 valueE21, valueE22, 363 valueE21, valueE22, 330 xs11, xs12, 364 xs11, xs12, 331 xs21, xs22, 365 xs21, xs22, 332 valueT1, valueT2, 366 valueT1, valueT2, 333 k, integrDiff ); 367 k, integrDiff ); 334 368 335 return theta; 369 return theta; 336 } 370 } 337 371 338 //....oooOO0OOooo........oooOO0OOooo........oo 372 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 339 373 340 G4double G4MicroElecElasticModel::LinLogInterp 374 G4double G4MicroElecElasticModel::LinLogInterpolate(G4double e1, 341 G4double e2, << 375 G4double e2, 342 G4double e, << 376 G4double e, 343 G4double xs1, << 377 G4double xs1, 344 G4double xs2) << 378 G4double xs2) 345 { 379 { 346 G4double d1 = std::log(xs1); 380 G4double d1 = std::log(xs1); 347 G4double d2 = std::log(xs2); 381 G4double d2 = std::log(xs2); 348 G4double value = G4Exp(d1 + (d2 - d1)*(e - e 382 G4double value = G4Exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1)); 349 return value; 383 return value; 350 } 384 } 351 385 352 //....oooOO0OOooo........oooOO0OOooo........oo 386 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 353 387 354 G4double G4MicroElecElasticModel::LinLinInterp 388 G4double G4MicroElecElasticModel::LinLinInterpolate(G4double e1, 355 G4double e2, << 389 G4double e2, 356 G4double e, << 390 G4double e, 357 G4double xs1, << 391 G4double xs1, 358 G4double xs2) << 392 G4double xs2) 359 { 393 { 360 G4double d1 = xs1; 394 G4double d1 = xs1; 361 G4double d2 = xs2; 395 G4double d2 = xs2; 362 G4double value = (d1 + (d2 - d1)*(e - e1)/ ( 396 G4double value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1)); 363 return value; 397 return value; 364 } 398 } 365 399 366 //....oooOO0OOooo........oooOO0OOooo........oo 400 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 367 401 368 G4double G4MicroElecElasticModel::LogLogInterp 402 G4double G4MicroElecElasticModel::LogLogInterpolate(G4double e1, 369 G4double e2, << 403 G4double e2, 370 G4double e, << 404 G4double e, 371 G4double xs1, << 405 G4double xs1, 372 G4double xs2) << 406 G4double xs2) 373 { 407 { 374 G4double a = (std::log10(xs2)-std::log10(xs1 408 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1)); 375 G4double b = std::log10(xs2) - a*std::log10( 409 G4double b = std::log10(xs2) - a*std::log10(e2); 376 G4double sigma = a*std::log10(e) + b; 410 G4double sigma = a*std::log10(e) + b; 377 G4double value = (std::pow(10.,sigma)); 411 G4double value = (std::pow(10.,sigma)); 378 return value; 412 return value; 379 } 413 } 380 414 381 //....oooOO0OOooo........oooOO0OOooo........oo 415 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 382 416 383 G4double G4MicroElecElasticModel::QuadInterpol 417 G4double G4MicroElecElasticModel::QuadInterpolator(G4double e11, G4double e12, 384 G4double e21, G4double e22, << 418 G4double e21, G4double e22, 385 G4double xs11, G4double xs12, << 419 G4double xs11, G4double xs12, 386 G4double xs21, G4double xs22, << 420 G4double xs21, G4double xs22, 387 G4double t1, G4double t2, << 421 G4double t1, G4double t2, 388 G4double t, G4double e) << 422 G4double t, G4double e) 389 { << 423 { 390 // Log-Log << 424 // Log-Log 391 /* << 425 /* 392 G4double interpolatedvalue1 = LogLogInterp << 426 G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12); 393 G4double interpolatedvalue2 = LogLogInterp << 427 G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22); 394 G4double value = LogLogInterpolate(t1, t2, << 428 G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2); 395 << 429 396 << 430 397 // Lin-Log << 431 // Lin-Log 398 G4double interpolatedvalue1 = LinLogInterp << 432 G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12); 399 G4double interpolatedvalue2 = LinLogInterp << 433 G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22); 400 G4double value = LinLogInterpolate(t1, t2, << 434 G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2); 401 */ << 435 */ 402 436 403 // Lin-Lin 437 // Lin-Lin 404 G4double interpolatedvalue1 = LinLinInterpol 438 G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12); 405 G4double interpolatedvalue2 = LinLinInterpol 439 G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22); 406 G4double value = LinLinInterpolate(t1, t2, t 440 G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2); 407 441 408 return value; 442 return value; 409 } 443 } 410 444 411 //....oooOO0OOooo........oooOO0OOooo........oo 445 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 412 446 413 G4double G4MicroElecElasticModel::RandomizeCos 447 G4double G4MicroElecElasticModel::RandomizeCosTheta(G4double k) 414 { 448 { 415 G4double integrdiff=0; 449 G4double integrdiff=0; 416 G4double uniformRand=G4UniformRand(); << 450 G4double uniformRand=G4UniformRand(); 417 integrdiff = uniformRand; << 451 integrdiff = uniformRand; 418 452 419 G4double theta=0.; << 453 G4double theta=0.; 420 G4double cosTheta=0.; << 454 G4double cosTheta=0.; 421 theta = Theta(G4Electron::ElectronDefinition << 455 theta = Theta(G4Electron::ElectronDefinition(),k/eV,integrdiff); 422 456 423 cosTheta= std::cos(theta*pi/180); << 457 cosTheta= std::cos(theta*pi/180); 424 458 425 return cosTheta; << 459 return cosTheta; 426 } 460 } 427 461