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