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