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