| 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 char *path = getenv("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