| 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 // Author: H. N. Tran (Ton Duc Thang Universit 26 // Author: H. N. Tran (Ton Duc Thang University)
27 // p, H, He, He+ and He++ models are assumed i 27 // p, H, He, He+ and He++ models are assumed identical
28 // NIMB 343, 132-137 (2015) 28 // NIMB 343, 132-137 (2015)
29 // 29 //
30 // The Geant4-DNA web site is available at htt 30 // The Geant4-DNA web site is available at http://geant4-dna.org
31 // 31 //
32 32
33 #include "G4DNAIonElasticModel.hh" 33 #include "G4DNAIonElasticModel.hh"
34 #include "G4PhysicalConstants.hh" 34 #include "G4PhysicalConstants.hh"
35 #include "G4SystemOfUnits.hh" 35 #include "G4SystemOfUnits.hh"
36 #include "G4DNAMolecularMaterial.hh" 36 #include "G4DNAMolecularMaterial.hh"
37 #include "G4ParticleTable.hh" 37 #include "G4ParticleTable.hh"
38 #include "G4Exp.hh" <<
39 38
40 //....oooOO0OOooo........oooOO0OOooo........oo 39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
41 40
42 using namespace std; 41 using namespace std;
43 42
44 //....oooOO0OOooo........oooOO0OOooo........oo 43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45 44
46 G4DNAIonElasticModel::G4DNAIonElasticModel (co 45 G4DNAIonElasticModel::G4DNAIonElasticModel (const G4ParticleDefinition*,
47 co 46 const G4String& nam) :
48 G4VEmModel(nam) << 47 G4VEmModel(nam), isInitialised(false)
49 { 48 {
>> 49 //nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
>> 50
50 killBelowEnergy = 100 * eV; 51 killBelowEnergy = 100 * eV;
51 lowEnergyLimit = 0 * eV; 52 lowEnergyLimit = 0 * eV;
52 highEnergyLimit = 1 * MeV; 53 highEnergyLimit = 1 * MeV;
53 SetLowEnergyLimit(lowEnergyLimit); 54 SetLowEnergyLimit(lowEnergyLimit);
54 SetHighEnergyLimit(highEnergyLimit); 55 SetHighEnergyLimit(highEnergyLimit);
55 56
56 verboseLevel = 0; 57 verboseLevel = 0;
57 // Verbosity scale: 58 // Verbosity scale:
58 // 0 = nothing 59 // 0 = nothing
59 // 1 = warning for energy non-conservation 60 // 1 = warning for energy non-conservation
60 // 2 = details of energy budget 61 // 2 = details of energy budget
61 // 3 = calculation of cross sections, file o 62 // 3 = calculation of cross sections, file openings, sampling of atoms
62 // 4 = entering in methods 63 // 4 = entering in methods
63 64
64 if(verboseLevel > 0) 65 if(verboseLevel > 0)
65 { 66 {
66 G4cout << "Ion elastic model is constructe 67 G4cout << "Ion elastic model is constructed " << G4endl<< "Energy range: "
67 << lowEnergyLimit / eV << " eV - " 68 << lowEnergyLimit / eV << " eV - "
68 << highEnergyLimit / MeV << " MeV" 69 << highEnergyLimit / MeV << " MeV"
69 << G4endl; 70 << G4endl;
70 } 71 }
71 << 72 fParticleChangeForGamma = 0;
72 fParticleChangeForGamma = nullptr; << 73 fpMolWaterDensity = 0;
73 fpMolWaterDensity = nullptr; << 74 fpTableData = 0;
74 fpTableData = nullptr; <<
75 fParticle_Mass = -1; 75 fParticle_Mass = -1;
76 76
77 // Selection of stationary mode 77 // Selection of stationary mode
78 78
79 statCode = false; 79 statCode = false;
80 } 80 }
81 81
82 //....oooOO0OOooo........oooOO0OOooo........oo 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
83 83
84 G4DNAIonElasticModel::~G4DNAIonElasticModel () 84 G4DNAIonElasticModel::~G4DNAIonElasticModel ()
85 { 85 {
86 // For total cross section 86 // For total cross section
87 delete fpTableData; << 87 if(fpTableData) delete fpTableData;
88 } 88 }
89 89
90 //....oooOO0OOooo........oooOO0OOooo........oo 90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91 91
92 void 92 void
93 G4DNAIonElasticModel::Initialise ( 93 G4DNAIonElasticModel::Initialise (
94 const G4ParticleDefinition* particleDefini 94 const G4ParticleDefinition* particleDefinition,
95 const G4DataVector& /*cuts*/) 95 const G4DataVector& /*cuts*/)
96 { 96 {
97 97
98 if(verboseLevel > 3) 98 if(verboseLevel > 3)
99 { 99 {
100 G4cout << "Calling G4DNAIonElasticModel::I 100 G4cout << "Calling G4DNAIonElasticModel::Initialise()" << G4endl;
101 } 101 }
102 102
103 // Energy limits 103 // Energy limits
104 104
105 if (LowEnergyLimit() < lowEnergyLimit) 105 if (LowEnergyLimit() < lowEnergyLimit)
106 { 106 {
107 G4cout << "G4DNAIonElasticModel: low energ 107 G4cout << "G4DNAIonElasticModel: low energy limit increased from " <<
108 LowEnergyLimit()/eV << " eV to " << lowEne 108 LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
109 SetLowEnergyLimit(lowEnergyLimit); 109 SetLowEnergyLimit(lowEnergyLimit);
110 } 110 }
111 111
112 if (HighEnergyLimit() > highEnergyLimit) 112 if (HighEnergyLimit() > highEnergyLimit)
113 { 113 {
114 G4cout << "G4DNAIonElasticModel: high ener 114 G4cout << "G4DNAIonElasticModel: high energy limit decreased from " <<
115 HighEnergyLimit()/MeV << " MeV to " << hig 115 HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl;
116 SetHighEnergyLimit(highEnergyLimit); 116 SetHighEnergyLimit(highEnergyLimit);
117 } 117 }
118 118
119 // Reading of data files 119 // Reading of data files
120 120
121 G4double scaleFactor = 1e-16*cm*cm; 121 G4double scaleFactor = 1e-16*cm*cm;
122 122
123 const char *path = G4FindDataDir("G4LEDATA") << 123 char *path = getenv("G4LEDATA");
124 124
125 if (path == nullptr) << 125 if (!path)
126 { 126 {
127 G4Exception("G4IonElasticModel::Initialise 127 G4Exception("G4IonElasticModel::Initialise","em0006",
128 FatalException,"G4LEDATA environment v 128 FatalException,"G4LEDATA environment variable not set.");
129 return; 129 return;
130 } 130 }
131 131
132 G4String totalXSFile; 132 G4String totalXSFile;
133 std::ostringstream fullFileName; 133 std::ostringstream fullFileName;
134 134
135 G4DNAGenericIonsManager *instance; 135 G4DNAGenericIonsManager *instance;
136 instance = G4DNAGenericIonsManager::Instance 136 instance = G4DNAGenericIonsManager::Instance();
137 G4ParticleDefinition* protonDef = 137 G4ParticleDefinition* protonDef =
138 G4ParticleTable::GetParticleTable()->FindPar 138 G4ParticleTable::GetParticleTable()->FindParticle("proton");
139 G4ParticleDefinition* hydrogenDef = instance 139 G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
140 G4ParticleDefinition* heliumDef = instance-> 140 G4ParticleDefinition* heliumDef = instance->GetIon("helium");
141 G4ParticleDefinition* alphaplusDef = instanc 141 G4ParticleDefinition* alphaplusDef = instance->GetIon("alpha+");
142 G4ParticleDefinition* alphaplusplusDef = ins 142 G4ParticleDefinition* alphaplusplusDef = instance->GetIon("alpha++");
143 G4String proton, hydrogen, helium, alphaplus 143 G4String proton, hydrogen, helium, alphaplus, alphaplusplus;
144 144
145 if ( 145 if (
146 (particleDefinition == protonDef && prot << 146 (particleDefinition == protonDef && protonDef != 0)
147 || 147 ||
148 (particleDefinition == hydrogenDef && hy << 148 (particleDefinition == hydrogenDef && hydrogenDef != 0)
149 ) 149 )
150 { 150 {
151 // For total cross section of p,h 151 // For total cross section of p,h
152 fParticle_Mass = 1.; << 152 fParticle_Mass = 1.;
153 totalXSFile = "dna/sigma_elastic_proton_HT 153 totalXSFile = "dna/sigma_elastic_proton_HTran";
154 154
155 // For final state 155 // For final state
156 fullFileName << path << "/dna/sigmadiff_cu 156 fullFileName << path << "/dna/sigmadiff_cumulated_elastic_proton_HTran.dat";
157 } 157 }
158 158
159 if ( 159 if (
160 (particleDefinition == instance->GetIon( << 160 (particleDefinition == instance->GetIon("helium") && heliumDef)
161 || 161 ||
162 (particleDefinition == instance->GetIon( << 162 (particleDefinition == instance->GetIon("alpha+") && alphaplusDef)
163 || 163 ||
164 (particleDefinition == instance->GetIon( << 164 (particleDefinition == instance->GetIon("alpha++") && alphaplusplusDef)
165 ) 165 )
166 { 166 {
>> 167 fParticle_Mass = 4.;
167 // For total cross section of he,he+,he++ 168 // For total cross section of he,he+,he++
168 fParticle_Mass = 4.; <<
169 totalXSFile = "dna/sigma_elastic_alpha_HTr 169 totalXSFile = "dna/sigma_elastic_alpha_HTran";
170 170
171 // For final state 171 // For final state
172 fullFileName << path << "/dna/sigmadiff_cu 172 fullFileName << path << "/dna/sigmadiff_cumulated_elastic_alpha_HTran.dat";
173 } 173 }
174 174
175 fpTableData = new G4DNACrossSectionDataSet(n 175 fpTableData = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
176 fpTableData->LoadData(totalXSFile); 176 fpTableData->LoadData(totalXSFile);
177 std::ifstream diffCrossSection(fullFileName. 177 std::ifstream diffCrossSection(fullFileName.str().c_str());
178 178
179 if (!diffCrossSection) 179 if (!diffCrossSection)
180 { 180 {
181 G4ExceptionDescription description; 181 G4ExceptionDescription description;
182 description << "Missing data file:" 182 description << "Missing data file:"
183 <<fullFileName.str().c_str()<< G4endl; 183 <<fullFileName.str().c_str()<< G4endl;
184 G4Exception("G4IonElasticModel::Initialise 184 G4Exception("G4IonElasticModel::Initialise","em0003",
185 FatalException, 185 FatalException,
186 description); 186 description);
187 } 187 }
188 188
189 // Added clear for MT 189 // Added clear for MT
190 190
191 eTdummyVec.clear(); 191 eTdummyVec.clear();
192 eVecm.clear(); 192 eVecm.clear();
193 fDiffCrossSectionData.clear(); 193 fDiffCrossSectionData.clear();
194 194
195 // 195 //
196 196
197 eTdummyVec.push_back(0.); 197 eTdummyVec.push_back(0.);
198 198
199 while(!diffCrossSection.eof()) 199 while(!diffCrossSection.eof())
200 { 200 {
201 G4double tDummy; << 201 double tDummy;
202 G4double eDummy; << 202 double eDummy;
203 diffCrossSection>>tDummy>>eDummy; 203 diffCrossSection>>tDummy>>eDummy;
204 204
205 // SI : mandatory eVecm initialization 205 // SI : mandatory eVecm initialization
206 206
207 if (tDummy != eTdummyVec.back()) 207 if (tDummy != eTdummyVec.back())
208 { 208 {
209 eTdummyVec.push_back(tDummy); 209 eTdummyVec.push_back(tDummy);
210 eVecm[tDummy].push_back(0.); 210 eVecm[tDummy].push_back(0.);
211 } 211 }
212 212
213 diffCrossSection>>fDiffCrossSectionData[tD 213 diffCrossSection>>fDiffCrossSectionData[tDummy][eDummy];
214 214
215 if (eDummy != eVecm[tDummy].back()) eVecm[ 215 if (eDummy != eVecm[tDummy].back()) eVecm[tDummy].push_back(eDummy);
216 216
217 } 217 }
218 218
219 // End final state 219 // End final state
220 if( verboseLevel>0 ) 220 if( verboseLevel>0 )
221 { 221 {
222 if (verboseLevel > 2) 222 if (verboseLevel > 2)
223 { 223 {
224 G4cout << "Loaded cross section files fo 224 G4cout << "Loaded cross section files for ion elastic model" << G4endl;
225 } 225 }
226 G4cout << "Ion elastic model is initialize 226 G4cout << "Ion elastic model is initialized " << G4endl
227 << "Energy range: " 227 << "Energy range: "
228 << LowEnergyLimit() / eV << " eV - " 228 << LowEnergyLimit() / eV << " eV - "
229 << HighEnergyLimit() / MeV << " MeV" 229 << HighEnergyLimit() / MeV << " MeV"
230 << G4endl; 230 << G4endl;
231 } 231 }
232 232
233 // Initialize water density pointer 233 // Initialize water density pointer
234 G4DNAMolecularMaterial::Instance()->Initiali 234 G4DNAMolecularMaterial::Instance()->Initialize();
235 fpMolWaterDensity = G4DNAMolecularMaterial:: 235 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->
236 GetNumMolPerVolTableFor(G4Material::GetMater 236 GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
237 237
238 if (isInitialised) return; << 238 if (isInitialised)
>> 239 { return;}
239 fParticleChangeForGamma = GetParticleChangeF 240 fParticleChangeForGamma = GetParticleChangeForGamma();
240 isInitialised = true; 241 isInitialised = true;
241 } 242 }
242 243
243 //....oooOO0OOooo........oooOO0OOooo........oo 244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
244 245
245 G4double 246 G4double
246 G4DNAIonElasticModel::CrossSectionPerVolume (c 247 G4DNAIonElasticModel::CrossSectionPerVolume (const G4Material* material,
247 c 248 const G4ParticleDefinition* p,
248 G 249 G4double ekin, G4double, G4double)
249 { 250 {
250 if(verboseLevel > 3) 251 if(verboseLevel > 3)
251 { 252 {
252 G4cout << "Calling CrossSectionPerVolume() 253 G4cout << "Calling CrossSectionPerVolume() of G4DNAIonElasticModel"
253 << G4endl; 254 << G4endl;
254 } 255 }
255 256
256 // Calculate total cross section for model 257 // Calculate total cross section for model
257 258
258 G4double sigma=0; 259 G4double sigma=0;
259 260
260 G4double waterDensity = (*fpMolWaterDensity) 261 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
261 262
262 const G4String& particleName = p->GetParticl << 263 if(waterDensity!= 0.0)
263 <<
264 if (ekin <= highEnergyLimit) <<
265 { 264 {
266 //SI : XS must not be zero otherwise sampl << 265 const G4String& particleName = p->GetParticleName();
267 if (ekin < killBelowEnergy) return DBL_MAX <<
268 // <<
269 266
270 if (fpTableData != nullptr) << 267 if (ekin < highEnergyLimit)
271 { 268 {
272 sigma = fpTableData->FindValue(ekin); << 269 //SI : XS must not be zero otherwise sampling of secondaries method ignored
>> 270 if (ekin < killBelowEnergy) return DBL_MAX;
>> 271 //
>> 272
>> 273 if (fpTableData != 0) // MK: is this check necessary?
>> 274 {
>> 275 sigma = fpTableData->FindValue(ekin);
>> 276 }
>> 277 else
>> 278 {
>> 279 G4Exception("G4DNAIonElasticModel::ComputeCrossSectionPerVolume","em0002",
>> 280 FatalException,"Model not applicable to particle type.");
>> 281 }
273 } 282 }
274 else << 283
>> 284 if (verboseLevel > 2)
275 { 285 {
276 G4Exception("G4DNAIonElasticModel::Compu << 286 G4cout << "__________________________________" << G4endl;
277 FatalException,"Model not applicable << 287 G4cout << "G4DNAIonElasticModel - XS INFO START" << G4endl;
>> 288 G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
>> 289 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
>> 290 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
>> 291 G4cout << "G4DNAIonElasticModel - XS INFO END" << G4endl;
278 } 292 }
279 } <<
280 293
281 if (verboseLevel > 2) <<
282 { <<
283 G4cout << "_______________________________ <<
284 G4cout << "G4DNAIonElasticModel - XS INFO <<
285 G4cout << "Kinetic energy(eV)=" << ekin/eV <<
286 G4cout << "Cross section per water molecul <<
287 G4cout << "Cross section per water molecul <<
288 G4cout << "G4DNAIonElasticModel - XS INFO <<
289 } 294 }
290 295
291 return sigma*waterDensity; 296 return sigma*waterDensity;
292 } 297 }
293 298
294 //....oooOO0OOooo........oooOO0OOooo........oo 299 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
295 300
296 void 301 void
297 G4DNAIonElasticModel::SampleSecondaries ( 302 G4DNAIonElasticModel::SampleSecondaries (
298 std::vector<G4DynamicParticle*>* /*fvect*/ 303 std::vector<G4DynamicParticle*>* /*fvect*/,
299 const G4MaterialCutsCouple* /*couple*/, 304 const G4MaterialCutsCouple* /*couple*/,
300 const G4DynamicParticle* aDynamicParticle, 305 const G4DynamicParticle* aDynamicParticle, G4double, G4double)
301 { 306 {
302 307
303 if(verboseLevel > 3) 308 if(verboseLevel > 3)
304 { 309 {
305 G4cout << "Calling SampleSecondaries() of 310 G4cout << "Calling SampleSecondaries() of G4DNAIonElasticModel" << G4endl;
306 } 311 }
307 312
308 G4double particleEnergy0 = aDynamicParticle- 313 G4double particleEnergy0 = aDynamicParticle->GetKineticEnergy();
309 314
310 if (particleEnergy0 < killBelowEnergy) 315 if (particleEnergy0 < killBelowEnergy)
311 { 316 {
312 fParticleChangeForGamma->SetProposedKineti 317 fParticleChangeForGamma->SetProposedKineticEnergy(0.);
313 fParticleChangeForGamma->ProposeTrackStatu 318 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
314 fParticleChangeForGamma->ProposeLocalEnerg 319 fParticleChangeForGamma->ProposeLocalEnergyDeposit(particleEnergy0);
315 return; 320 return;
316 } 321 }
317 322
318 if (particleEnergy0>= killBelowEnergy && par << 323 if (particleEnergy0>= killBelowEnergy && particleEnergy0 < highEnergyLimit)
319 { 324 {
320 G4double water_mass = 18.; 325 G4double water_mass = 18.;
321 326
322 G4double thetaCM = RandomizeThetaCM(partic 327 G4double thetaCM = RandomizeThetaCM(particleEnergy0, aDynamicParticle->GetDefinition());
323 328
324 //HT:convert to laboratory system 329 //HT:convert to laboratory system
325 330
326 G4double theta = std::atan(std::sin(thetaC 331 G4double theta = std::atan(std::sin(thetaCM*CLHEP::pi/180)
327 /(fParticle_Mass/water_mass+std::cos(t 332 /(fParticle_Mass/water_mass+std::cos(thetaCM*CLHEP::pi/180)));
328 333
329 G4double cosTheta= std::cos(theta); 334 G4double cosTheta= std::cos(theta);
330 335
331 // 336 //
332 337
333 G4double phi = 2. * CLHEP::pi * G4UniformR 338 G4double phi = 2. * CLHEP::pi * G4UniformRand();
334 339
335 G4ThreeVector zVers = aDynamicParticle->Ge 340 G4ThreeVector zVers = aDynamicParticle->GetMomentumDirection();
336 G4ThreeVector xVers = zVers.orthogonal(); 341 G4ThreeVector xVers = zVers.orthogonal();
337 G4ThreeVector yVers = zVers.cross(xVers); 342 G4ThreeVector yVers = zVers.cross(xVers);
338 343
339 G4double xDir = std::sqrt(1. - cosTheta*co 344 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
340 G4double yDir = xDir; 345 G4double yDir = xDir;
341 xDir *= std::cos(phi); 346 xDir *= std::cos(phi);
342 yDir *= std::sin(phi); 347 yDir *= std::sin(phi);
343 348
344 G4ThreeVector zPrimeVers((xDir*xVers + yDi 349 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
345 350
346 fParticleChangeForGamma->ProposeMomentumDi 351 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit());
347 352
348 G4double depositEnergyCM = 0; 353 G4double depositEnergyCM = 0;
349 354
350 //HT: deposited energy 355 //HT: deposited energy
351 depositEnergyCM = 4. * particleEnergy0 * f 356 depositEnergyCM = 4. * particleEnergy0 * fParticle_Mass * water_mass *
352 (1-std::cos(thetaCM*CLHEP::pi/180)) 357 (1-std::cos(thetaCM*CLHEP::pi/180))
353 / (2 * std::pow((fParticle_Mass+water_mass 358 / (2 * std::pow((fParticle_Mass+water_mass),2));
354 359
355 //SI: added protection particleEnergy0 >= << 360 if (!statCode) fParticleChangeForGamma->SetProposedKineticEnergy(particleEnergy0 - depositEnergyCM);
356 if (!statCode && (particleEnergy0 >= depos <<
357 <<
358 fParticleChangeForGamma->SetProposedKine <<
359 361
360 else fParticleChangeForGamma->SetProposedK 362 else fParticleChangeForGamma->SetProposedKineticEnergy(particleEnergy0);
361 363
362 fParticleChangeForGamma->ProposeLocalEnerg 364 fParticleChangeForGamma->ProposeLocalEnergyDeposit(depositEnergyCM);
363 } 365 }
364 } 366 }
365 367
366 //....oooOO0OOooo........oooOO0OOooo........oo 368 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
367 369
368 G4double 370 G4double
369 G4DNAIonElasticModel::Theta (G4ParticleDefinit 371 G4DNAIonElasticModel::Theta (G4ParticleDefinition * /*particleDefinition*/,
370 G4double k, G4dou 372 G4double k, G4double integrDiff)
371 { 373 {
372 G4double theta = 0.; 374 G4double theta = 0.;
373 G4double valueT1 = 0; 375 G4double valueT1 = 0;
374 G4double valueT2 = 0; 376 G4double valueT2 = 0;
375 G4double valueE21 = 0; 377 G4double valueE21 = 0;
376 G4double valueE22 = 0; 378 G4double valueE22 = 0;
377 G4double valueE12 = 0; 379 G4double valueE12 = 0;
378 G4double valueE11 = 0; 380 G4double valueE11 = 0;
379 G4double xs11 = 0; 381 G4double xs11 = 0;
380 G4double xs12 = 0; 382 G4double xs12 = 0;
381 G4double xs21 = 0; 383 G4double xs21 = 0;
382 G4double xs22 = 0; 384 G4double xs22 = 0;
383 385
384 // Protection against out of boundary access << 386 std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),
385 if (k==eTdummyVec.back()) k=k*(1.-1e-12); <<
386 // <<
387 <<
388 auto t2 = std::upper_bound(eTdummyVec.begin( <<
389 387 eTdummyVec.end(), k);
390 auto t1 = t2 - 1; << 388 std::vector<double>::iterator t1 = t2 - 1;
391 389
392 auto e12 = std::upper_bound(eVecm[(*t1)].beg << 390 std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),
393 391 eVecm[(*t1)].end(),
394 392 integrDiff);
395 auto e11 = e12 - 1; << 393 std::vector<double>::iterator e11 = e12 - 1;
396 394
397 auto e22 = std::upper_bound(eVecm[(*t2)].beg << 395 std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),
398 396 eVecm[(*t2)].end(),
399 397 integrDiff);
400 auto e21 = e22 - 1; << 398 std::vector<double>::iterator e21 = e22 - 1;
401 399
402 valueT1 = *t1; 400 valueT1 = *t1;
403 valueT2 = *t2; 401 valueT2 = *t2;
404 valueE21 = *e21; 402 valueE21 = *e21;
405 valueE22 = *e22; 403 valueE22 = *e22;
406 valueE12 = *e12; 404 valueE12 = *e12;
407 valueE11 = *e11; 405 valueE11 = *e11;
408 406
409 xs11 = fDiffCrossSectionData[valueT1][valueE 407 xs11 = fDiffCrossSectionData[valueT1][valueE11];
410 xs12 = fDiffCrossSectionData[valueT1][valueE 408 xs12 = fDiffCrossSectionData[valueT1][valueE12];
411 xs21 = fDiffCrossSectionData[valueT2][valueE 409 xs21 = fDiffCrossSectionData[valueT2][valueE21];
412 xs22 = fDiffCrossSectionData[valueT2][valueE 410 xs22 = fDiffCrossSectionData[valueT2][valueE22];
413 411
414 if(xs11 == 0 && xs12 == 0 && xs21 == 0 && xs 412 if(xs11 == 0 && xs12 == 0 && xs21 == 0 && xs22 == 0) return (0.);
415 413
416 theta = QuadInterpolator(valueE11, valueE12, 414 theta = QuadInterpolator(valueE11, valueE12, valueE21, valueE22, xs11, xs12,
417 xs21, xs22, valueT1 415 xs21, xs22, valueT1, valueT2, k, integrDiff);
418 416
419 return theta; 417 return theta;
420 } 418 }
421 419
422 //....oooOO0OOooo........oooOO0OOooo........oo 420 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
423 421
424 G4double 422 G4double
425 G4DNAIonElasticModel::LinLinInterpolate (G4dou 423 G4DNAIonElasticModel::LinLinInterpolate (G4double e1, G4double e2, G4double e,
426 G4dou 424 G4double xs1, G4double xs2)
427 { 425 {
428 G4double d1 = xs1; 426 G4double d1 = xs1;
429 G4double d2 = xs2; 427 G4double d2 = xs2;
430 G4double value = (d1 + (d2 - d1) * (e - e1) 428 G4double value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
431 return value; 429 return value;
432 } 430 }
433 431
434 //....oooOO0OOooo........oooOO0OOooo........oo 432 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
435 433
436 G4double 434 G4double
437 G4DNAIonElasticModel::LinLogInterpolate (G4dou 435 G4DNAIonElasticModel::LinLogInterpolate (G4double e1, G4double e2, G4double e,
438 G4dou 436 G4double xs1, G4double xs2)
439 { 437 {
440 G4double d1 = std::log(xs1); 438 G4double d1 = std::log(xs1);
441 G4double d2 = std::log(xs2); 439 G4double d2 = std::log(xs2);
442 G4double value = G4Exp(d1 + (d2 - d1) * (e - << 440 G4double value = std::exp(d1 + (d2 - d1) * (e - e1) / (e2 - e1));
443 return value; 441 return value;
444 } 442 }
445 443
446 //....oooOO0OOooo........oooOO0OOooo........oo 444 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
447 445
448 G4double 446 G4double
449 G4DNAIonElasticModel::LogLogInterpolate (G4dou 447 G4DNAIonElasticModel::LogLogInterpolate (G4double e1, G4double e2, G4double e,
450 G4dou 448 G4double xs1, G4double xs2)
451 { 449 {
452 G4double a = (std::log10(xs2) - std::log10(x 450 G4double a = (std::log10(xs2) - std::log10(xs1))
453 / (std::log10(e2) - std::log10(e1)); 451 / (std::log10(e2) - std::log10(e1));
454 G4double b = std::log10(xs2) - a * std::log1 452 G4double b = std::log10(xs2) - a * std::log10(e2);
455 G4double sigma = a * std::log10(e) + b; 453 G4double sigma = a * std::log10(e) + b;
456 G4double value = (std::pow(10., sigma)); 454 G4double value = (std::pow(10., sigma));
457 return value; 455 return value;
458 } 456 }
459 457
460 //....oooOO0OOooo........oooOO0OOooo........oo 458 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
461 459
462 G4double 460 G4double
463 G4DNAIonElasticModel::QuadInterpolator (G4doub 461 G4DNAIonElasticModel::QuadInterpolator (G4double e11, G4double e12,
464 G4doub 462 G4double e21, G4double e22,
465 G4doub 463 G4double xs11, G4double xs12,
466 G4doub 464 G4double xs21, G4double xs22,
467 G4doub 465 G4double t1, G4double t2, G4double t,
468 G4doub 466 G4double e)
469 { 467 {
470 // Log-Log 468 // Log-Log
471 /* 469 /*
472 G4double interpolatedvalue1 = LogLogInterpo 470 G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
473 G4double interpolatedvalue2 = LogLogInterpo 471 G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
474 G4double value = LogLogInterpolate(t1, t2, 472 G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
475 */ 473 */
476 474
477 // Lin-Log 475 // Lin-Log
478 /* 476 /*
479 G4double interpolatedvalue1 = LinLogInterpo 477 G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
480 G4double interpolatedvalue2 = LinLogInterpo 478 G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
481 G4double value = LinLogInterpolate(t1, t2, 479 G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
482 */ 480 */
483 481
484 // Lin-Lin 482 // Lin-Lin
485 G4double interpolatedvalue1 = LinLinInterpol 483 G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
486 G4double interpolatedvalue2 = LinLinInterpol 484 G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
487 G4double value = LinLinInterpolate(t1, t2, t 485 G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1,
488 interpola 486 interpolatedvalue2);
489 487
490 return value; 488 return value;
491 } 489 }
492 490
493 //....oooOO0OOooo........oooOO0OOooo........oo 491 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
494 492
495 G4double 493 G4double
496 G4DNAIonElasticModel::RandomizeThetaCM ( 494 G4DNAIonElasticModel::RandomizeThetaCM (
497 G4double k, G4ParticleDefinition * particl 495 G4double k, G4ParticleDefinition * particleDefinition)
498 { 496 {
499 G4double integrdiff = G4UniformRand(); 497 G4double integrdiff = G4UniformRand();
500 return Theta(particleDefinition, k / eV, int 498 return Theta(particleDefinition, k / eV, integrdiff);
501 } 499 }
502 500
503 //....oooOO0OOooo........oooOO0OOooo........oo 501 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
504 502
505 void 503 void
506 G4DNAIonElasticModel::SetKillBelowThreshold (G 504 G4DNAIonElasticModel::SetKillBelowThreshold (G4double threshold)
507 { 505 {
508 killBelowEnergy = threshold; 506 killBelowEnergy = threshold;
509 507
510 if(killBelowEnergy < 100 * eV) 508 if(killBelowEnergy < 100 * eV)
511 { 509 {
512 G4cout << "*** WARNING : the G4DNAIonElast 510 G4cout << "*** WARNING : the G4DNAIonElasticModel class is not "
513 "activated below 100 eV !" 511 "activated below 100 eV !"
514 << G4endl; 512 << G4endl;
515 } 513 }
516 } 514 }
517 515
518 516