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