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