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