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 // Created by Z. Francis 29 // Created by Z. Francis 29 30 30 #include "G4DNASancheExcitationModel.hh" 31 #include "G4DNASancheExcitationModel.hh" 31 #include "G4SystemOfUnits.hh" 32 #include "G4SystemOfUnits.hh" 32 #include "G4DNAMolecularMaterial.hh" 33 #include "G4DNAMolecularMaterial.hh" 33 34 34 //....oooOO0OOooo........oooOO0OOooo........oo << 35 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 35 36 36 using namespace std; 37 using namespace std; 37 38 38 //#define SANCHE_VERBOSE << 39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 39 << 40 //....oooOO0OOooo........oooOO0OOooo........oo << 41 40 42 G4DNASancheExcitationModel::G4DNASancheExcitat 41 G4DNASancheExcitationModel::G4DNASancheExcitationModel(const G4ParticleDefinition*, 43 << 42 const G4String& nam) 44 G4VEmModel(nam) << 43 :G4VEmModel(nam),isInitialised(false) 45 { 44 { 46 fpWaterDensity = nullptr; << 45 // nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER"); >> 46 fpWaterDensity = 0; 47 47 48 SetLowEnergyLimit(2.*eV); << 48 lowEnergyLimit = 2 * eV; 49 SetHighEnergyLimit(100*eV); << 49 highEnergyLimit = 100 * eV; 50 nLevels = 9; << 50 SetLowEnergyLimit(lowEnergyLimit); 51 << 51 SetHighEnergyLimit(highEnergyLimit); 52 verboseLevel = 0; << 52 nLevels = 9; 53 // Verbosity scale: << 53 54 // 0 = nothing << 54 verboseLevel= 0; 55 // 1 = warning for energy non-conservation << 55 // Verbosity scale: 56 // 2 = details of energy budget << 56 // 0 = nothing 57 // 3 = calculation of cross sections, file o << 57 // 1 = warning for energy non-conservation 58 // 4 = entering in methods << 58 // 2 = details of energy budget >> 59 // 3 = calculation of cross sections, file openings, sampling of atoms >> 60 // 4 = entering in methods 59 61 60 #ifdef SANCHE_VERBOSE << 62 if (verboseLevel > 0) 61 if (verboseLevel > 0) << 63 { 62 { << 64 G4cout << "Sanche Excitation model is constructed " << G4endl 63 G4cout << "Sanche Excitation model is cons << 65 << "Energy range: " 64 << G4endl << 66 << lowEnergyLimit / eV << " eV - " 65 << "Energy range: " << 67 << highEnergyLimit / eV << " eV" 66 << LowEnergyLimit() / eV << " eV - << 68 << G4endl; 67 << HighEnergyLimit() / eV << " eV" << 69 } 68 << G4endl; << 70 fParticleChangeForGamma = 0; 69 } << 71 fpWaterDensity = 0; 70 #endif << 71 << 72 fParticleChangeForGamma = nullptr; << 73 fpWaterDensity = nullptr; << 74 << 75 // Selection of stationary mode << 76 << 77 statCode = false; << 78 } 72 } 79 73 80 //....oooOO0OOooo........oooOO0OOooo........oo << 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 81 75 82 G4DNASancheExcitationModel::~G4DNASancheExcita 76 G4DNASancheExcitationModel::~G4DNASancheExcitationModel() 83 = default; << 77 {} 84 78 85 //....oooOO0OOooo........oooOO0OOooo........oo << 79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 86 80 87 void << 81 void G4DNASancheExcitationModel::Initialise(const G4ParticleDefinition* /*particle*/, 88 G4DNASancheExcitationModel:: << 82 const G4DataVector& /*cuts*/) 89 Initialise(const G4ParticleDefinition* /*parti << 90 const G4DataVector& /*cuts*/) << 91 { 83 { 92 << 93 #ifdef SANCHE_VERBOSE << 94 if (verboseLevel > 3) << 95 { << 96 G4cout << "Calling G4DNASancheExcitationMo << 97 << G4endl; << 98 } << 99 #endif << 100 << 101 // Energy limits << 102 84 103 if (LowEnergyLimit() < 2.*eV) << 104 { << 105 G4Exception("*** WARNING : the G4DNASanche << 106 "validated below 2 eV !", << 107 "", JustWarning, ""); << 108 } << 109 85 110 if (HighEnergyLimit() > 100.*eV) << 86 if (verboseLevel > 3) 111 { << 87 G4cout << "Calling G4DNASancheExcitationModel::Initialise()" << G4endl; 112 G4cout << "G4DNASancheExcitationModel: hig << 113 HighEnergyLimit()/eV << " eV to " << 100. << 114 SetHighEnergyLimit(100.*eV); << 115 } << 116 88 117 // << 89 // Energy limits 118 #ifdef SANCHE_VERBOSE << 119 if (verboseLevel > 0) << 120 { << 121 G4cout << "Sanche Excitation model is init << 122 << "Energy range: " << 123 << LowEnergyLimit() / eV << " eV - " << 124 << HighEnergyLimit() / eV << " eV" << 125 << G4endl; << 126 } << 127 #endif << 128 << 129 // Initialize water density pointer << 130 fpWaterDensity = G4DNAMolecularMaterial::Ins << 131 GetNumMolPerVolTableFor(G4Material::GetM << 132 << 133 if (isInitialised) {return;} << 134 << 135 fParticleChangeForGamma = GetParticleChangeF << 136 isInitialised = true; << 137 << 138 const char *path = G4FindDataDir("G4LEDATA") << 139 std::ostringstream eFullFileName; << 140 eFullFileName << path << "/dna/sigma_excitat << 141 std::ifstream input(eFullFileName.str().c_st << 142 90 143 if (!input) << 91 if (LowEnergyLimit() < lowEnergyLimit) 144 { << 92 { 145 G4Exception("G4DNASancheExcitationModel::I << 93 G4cout << "G4DNASancheExcitationModel: low energy limit increased from " << 146 FatalException,"Missing data file:/dna << 94 LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl; 147 } << 95 SetLowEnergyLimit(lowEnergyLimit); >> 96 } 148 97 149 // March 25th, 2014 - Vaclav Stepan, Sebasti << 98 if (HighEnergyLimit() > highEnergyLimit) 150 // Added clear for MT << 99 { 151 tdummyVec.clear(); << 100 G4cout << "G4DNASancheExcitationModel: high energy limit decreased from " << 152 // << 101 HighEnergyLimit()/eV << " eV to " << highEnergyLimit/eV << " eV" << G4endl; >> 102 SetHighEnergyLimit(highEnergyLimit); >> 103 } 153 104 154 G4double t; << 105 // 155 G4double xs; << 156 106 157 while(!input.eof()) << 107 if (verboseLevel > 0) 158 { << 108 { 159 input>>t; << 109 G4cout << "Sanche Excitation model is initialized " << G4endl 160 tdummyVec.push_back(t); << 110 << "Energy range: " >> 111 << LowEnergyLimit() / eV << " eV - " >> 112 << HighEnergyLimit() / eV << " eV" >> 113 << G4endl; >> 114 } 161 115 162 fEnergyLevelXS.emplace_back(); << 116 // Initialize water density pointer 163 fEnergyTotalXS.push_back(0); << 117 fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER")); 164 std::vector<G4double>& levelXS = fEnergyLe << 165 levelXS.reserve(9); << 166 118 167 // G4cout<<t; << 119 if (isInitialised) { return; } >> 120 fParticleChangeForGamma = GetParticleChangeForGamma(); >> 121 isInitialised = true; >> 122 >> 123 char *path = getenv("G4LEDATA"); >> 124 std::ostringstream eFullFileName; >> 125 eFullFileName << path << "/dna/sigma_excitationvib_e_sanche.dat"; >> 126 std::ifstream input(eFullFileName.str().c_str()); 168 127 169 for(size_t i = 0 ; i < 9 ;++i) << 128 if (!input) 170 { 129 { 171 input>>xs; << 130 G4Exception("G4DNASancheExcitationModel::Initialise","em0003", 172 levelXS.push_back(xs); << 131 FatalException,"Missing data file:/dna/sigma_excitationvib_e_sanche.dat"); 173 fEnergyTotalXS.back() += xs; << 174 // G4cout <<" " << levelXS[i]; << 175 } 132 } 176 133 177 // G4cout << G4endl; << 134 while(!input.eof()) 178 } << 135 { >> 136 double t; >> 137 input>>t; >> 138 tdummyVec.push_back(t); >> 139 input>>map1[t][0]>>map1[t][1]>>map1[t][2]>>map1[t][3]>>map1[t][4]>>map1[t][5]>>map1[t][6]>>map1[t][7]>>map1[t][8]; >> 140 //G4cout<<t<<" "<<map1[t][0]<<map1[t][1]<<map1[t][2]<<map1[t][3]<<map1[t][4]<<map1[t][5]<<map1[t][6]<<map1[t][7]<<map1[t][8]<<G4endl; >> 141 } 179 } 142 } 180 143 181 //....oooOO0OOooo........oooOO0OOooo........oo << 144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 182 145 183 G4double G4DNASancheExcitationModel::CrossSect << 146 G4double G4DNASancheExcitationModel::CrossSectionPerVolume(const G4Material* material, 184 << 147 const G4ParticleDefinition* particleDefinition, 185 #ifdef SANCHE_VERBOSE << 186 << 187 #endif << 188 << 189 148 G4double ekin, 190 149 G4double, 191 150 G4double) 192 { 151 { 193 #ifdef SANCHE_VERBOSE << 152 if (verboseLevel > 3) 194 if (verboseLevel > 3) << 153 G4cout << "Calling CrossSectionPerVolume() of G4DNASancheExcitationModel" << G4endl; 195 { << 196 G4cout << "Calling CrossSectionPerVolume() << 197 << G4endl; << 198 } << 199 #endif << 200 154 201 // Calculate total cross section for model << 155 // Calculate total cross section for model 202 156 203 G4double sigma = 0.; << 157 G4double sigma=0; 204 158 205 G4double waterDensity = (*fpWaterDensity)[ma << 159 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()]; 206 160 207 if (ekin >= LowEnergyLimit() && ekin <= High << 161 if(waterDensity!= 0.0) 208 sigma = TotalCrossSection(ekin); << 162 // if (material == nistwater || material->GetBaseMaterial() == nistwater) >> 163 { 209 164 210 #ifdef SANCHE_VERBOSE << 165 if (particleDefinition == G4Electron::ElectronDefinition()) 211 if (verboseLevel > 2) << 166 { 212 { << 167 if (ekin >= lowEnergyLimit && ekin < highEnergyLimit) 213 G4cout << "_______________________________ << 168 { 214 G4cout << "=== G4DNASancheExcitationModel << 169 sigma = Sum(ekin); 215 G4cout << "=== Kinetic energy(eV)=" << eki << 170 } 216 G4cout << "=== Cross section per water mol << 171 } 217 G4cout << "=== Cross section per water mol << 172 218 G4cout << "=== G4DNASancheExcitationModel << 173 if (verboseLevel > 2) 219 } << 174 { 220 #endif << 175 G4cout << "__________________________________" << G4endl; >> 176 G4cout << "°°° G4DNASancheExcitationModel - XS INFO START" << G4endl; >> 177 G4cout << "°°° Kinetic energy(eV)=" << ekin/eV << " particle : " << particleDefinition->GetParticleName() << G4endl; >> 178 G4cout << "°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl; >> 179 G4cout << "°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl; >> 180 // G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl; >> 181 G4cout << "°°° G4DNASancheExcitationModel - XS INFO END" << G4endl; >> 182 } >> 183 >> 184 } // if water >> 185 >> 186 >> 187 // return sigma*2*material->GetAtomicNumDensityVector()[1]; >> 188 return sigma*2*waterDensity; >> 189 // see papers for factor 2 description 221 190 222 return sigma*2.*waterDensity; << 223 // see papers for factor 2 description (liqu << 224 } 191 } 225 192 226 //....oooOO0OOooo........oooOO0OOooo........oo << 193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 227 194 228 void G4DNASancheExcitationModel::SampleSeconda << 195 void G4DNASancheExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>*, 229 << 230 196 const G4MaterialCutsCouple*, 231 197 const G4DynamicParticle* aDynamicElectron, 232 198 G4double, 233 199 G4double) 234 { 200 { 235 #ifdef SANCHE_VERBOSE << 236 if (verboseLevel > 3) << 237 { << 238 G4cout << "Calling SampleSecondaries() of << 239 << G4endl; << 240 } << 241 #endif << 242 201 243 G4double electronEnergy0 = aDynamicElectron- << 202 if (verboseLevel > 3) 244 G4int level = RandomSelect(electronEnergy0); << 203 G4cout << "Calling SampleSecondaries() of G4DNASancheExcitationModel" << G4endl; 245 G4double excitationEnergy = VibrationEnergy( << 246 G4double newEnergy = electronEnergy0 - excit << 247 << 248 /* << 249 if (electronEnergy0 < highEnergyLimit) << 250 { << 251 if (newEnergy >= lowEnergyLimit) << 252 { << 253 fParticleChangeForGamma->ProposeMomentu << 254 fParticleChangeForGamma->SetProposedKin << 255 fParticleChangeForGamma->ProposeLocalEn << 256 } << 257 << 258 else << 259 { << 260 fParticleChangeForGamma->ProposeTrackSt << 261 fParticleChangeForGamma->ProposeLocalEn << 262 } << 263 } << 264 */ << 265 204 266 if (electronEnergy0 <= HighEnergyLimit() && << 205 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy(); 267 { << 206 G4int level = RandomSelect(electronEnergy0); >> 207 G4double excitationEnergy = VibrationEnergy(level); // levels go from 0 to 8 >> 208 G4double newEnergy = electronEnergy0 - excitationEnergy; 268 209 269 if (!statCode) << 210 /* 270 { << 211 if (electronEnergy0 < highEnergyLimit) >> 212 { >> 213 if (newEnergy >= lowEnergyLimit) >> 214 { 271 fParticleChangeForGamma->ProposeMomentum 215 fParticleChangeForGamma->ProposeMomentumDirection(aDynamicElectron->GetMomentumDirection()); 272 fParticleChangeForGamma->SetProposedKine 216 fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy); 273 fParticleChangeForGamma->ProposeLocalEne 217 fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy); 274 } 218 } 275 219 276 else << 220 else 277 { 221 { 278 fParticleChangeForGamma->ProposeMomentum << 222 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill); 279 fParticleChangeForGamma->SetProposedKine << 223 fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0); 280 fParticleChangeForGamma->ProposeLocalEne << 281 } 224 } 282 << 283 } 225 } >> 226 */ 284 227 285 // << 228 if (electronEnergy0 < highEnergyLimit && newEnergy>0.) 286 } << 229 { 287 << 230 fParticleChangeForGamma->ProposeMomentumDirection(aDynamicElectron->GetMomentumDirection()); 288 //....oooOO0OOooo........oooOO0OOooo........oo << 231 fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy); >> 232 fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy); >> 233 } 289 234 290 G4double G4DNASancheExcitationModel::PartialCr << 235 // 291 << 292 { << 293 // Protection against out of boundary access << 294 if (t/eV==tdummyVec.back()) t=t*(1.-1e-12); << 295 // << 296 << 297 auto t2 = std::upper_bound(tdummyVec.begin() << 298 << 299 auto t1 = t2 - 1; << 300 << 301 size_t i1 = t1 - tdummyVec.begin(); << 302 size_t i2 = t2 - tdummyVec.begin(); << 303 << 304 G4double sigma = LinInterpolate((*t1), (*t2) << 305 t / eV, << 306 fEnergyLevelXS << 307 fEnergyLevelXS << 308 << 309 static const G4double conv_factor = 1e-16 * << 310 << 311 sigma *= conv_factor; << 312 if (sigma == 0.) sigma = 1e-30; << 313 return (sigma); << 314 } 236 } 315 237 316 //....oooOO0OOooo........oooOO0OOooo........oo 238 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 317 239 318 G4double G4DNASancheExcitationModel::TotalCros << 240 G4double G4DNASancheExcitationModel::PartialCrossSection(G4double t, G4int level) 319 { 241 { 320 // Protection against out of boundary access << 242 std::vector<double>::iterator t2 = std::upper_bound(tdummyVec.begin(),tdummyVec.end(), t/eV); 321 if (t/eV==tdummyVec.back()) t=t*(1.-1e-12); << 243 std::vector<double>::iterator t1 = t2-1; 322 // << 244 323 << 245 double sigma = LinInterpolate((*t1), (*t2), t/eV, map1[*t1][level], map1[*t2][level]); 324 auto t2 = std::upper_bound(tdummyVec.begin() << 246 sigma*=1e-16*cm*cm; 325 << 247 if(sigma==0.)sigma=1e-30; 326 auto t1 = t2 - 1; << 248 return (sigma); 327 << 328 size_t i1 = t1 - tdummyVec.begin(); << 329 size_t i2 = t2 - tdummyVec.begin(); << 330 << 331 G4double sigma = LinInterpolate((*t1), (*t2) << 332 t / eV, << 333 fEnergyTotalXS << 334 fEnergyTotalXS << 335 << 336 static const G4double conv_factor = 1e-16 * << 337 << 338 sigma *= conv_factor; << 339 if (sigma == 0.) sigma = 1e-30; << 340 return (sigma); << 341 } 249 } 342 250 343 //....oooOO0OOooo........oooOO0OOooo........oo 251 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 344 252 345 G4double G4DNASancheExcitationModel::Vibration 253 G4double G4DNASancheExcitationModel::VibrationEnergy(G4int level) 346 { 254 { 347 static G4double energies[9] = { 0.01, 0.024, << 255 G4double energies[9] = {0.01, 0.024, 0.061, 0.092, 0.204, 0.417, 0.460, 0.500, 0.835}; 348 0.500, 0.835 }; << 256 return(energies[level]*eV); 349 return (energies[level] * eV); << 350 } 257 } 351 258 352 //....oooOO0OOooo........oooOO0OOooo........oo 259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 353 260 354 G4int G4DNASancheExcitationModel::RandomSelect 261 G4int G4DNASancheExcitationModel::RandomSelect(G4double k) 355 { 262 { 356 263 357 // Level Selection Counting can be done here << 264 // Level Selection Counting can be done here ! 358 265 359 G4int i = nLevels; << 266 G4int i = nLevels; 360 G4double value = 0.; << 267 G4double value = 0.; 361 std::deque<G4double> values; << 268 std::deque<double> values; 362 269 363 while (i > 0) << 270 while (i > 0) 364 { << 271 { 365 i--; << 272 i--; 366 G4double partial = PartialCrossSection(k, << 273 G4double partial = PartialCrossSection(k,i); 367 values.push_front(partial); << 274 values.push_front(partial); 368 value += partial; << 275 value += partial; 369 } << 276 } 370 << 371 value *= G4UniformRand(); << 372 277 373 i = nLevels; << 278 value *= G4UniformRand(); >> 279 >> 280 i = nLevels; 374 281 375 while (i > 0) << 282 while (i > 0) 376 { << 377 i--; << 378 if (values[i] > value) << 379 { 283 { 380 //outcount<<i<<" "<<VibrationEnergy(i)< << 284 i--; 381 return i; << 285 if (values[i] > value) >> 286 { >> 287 //outcount<<i<<" "<<VibrationEnergy(i)<<G4endl; >> 288 return i; >> 289 } >> 290 value -= values[i]; 382 } 291 } 383 value -= values[i]; << 384 } << 385 292 386 //outcount<<0<<" "<<VibrationEnergy(0)<<G4e << 293 //outcount<<0<<" "<<VibrationEnergy(0)<<G4endl; 387 294 388 return 0; << 295 return 0; 389 } 296 } 390 297 391 //....oooOO0OOooo........oooOO0OOooo........oo 298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 392 299 393 G4double G4DNASancheExcitationModel::Sum(G4dou 300 G4double G4DNASancheExcitationModel::Sum(G4double k) 394 { 301 { 395 G4double totalCrossSection = 0.; << 302 G4double totalCrossSection = 0.; 396 << 397 for (G4int i = 0; i < nLevels; i++) << 398 { << 399 totalCrossSection += PartialCrossSection(k << 400 } << 401 303 402 return totalCrossSection; << 304 for (G4int i=0; i<nLevels; i++) >> 305 { >> 306 totalCrossSection += PartialCrossSection(k,i); >> 307 } >> 308 return totalCrossSection; 403 } 309 } 404 310 405 //....oooOO0OOooo........oooOO0OOooo........oo 311 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 406 312 407 G4double G4DNASancheExcitationModel::LinInterp << 313 G4double G4DNASancheExcitationModel::LinInterpolate(G4double e1, 408 314 G4double e2, 409 315 G4double e, 410 316 G4double xs1, 411 317 G4double xs2) 412 { 318 { 413 G4double a = (xs2 - xs1) / (e2 - e1); << 319 G4double a = (xs2 - xs1) / (e2 - e1); 414 G4double b = xs2 - a * e2; << 320 G4double b = xs2 - a*e2; 415 G4double value = a * e + b; << 321 G4double value = a*e + b; 416 // G4cout<<"interP >> "<<e1<<" "<<e2<<" " << 322 // G4cout<<"interP >> "<<e1<<" "<<e2<<" "<<e<<" "<<xs1<<" "<<xs2<<" "<<a<<" "<<b<<" "<<value<<G4endl; 417 // <<xs1<<" "<<xs2<<" "<<a<<" "<<b<<" "< << 418 323 419 return value; << 324 return value; 420 } 325 } 421 326 422 327