Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNAIonElasticModel.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Author: H. N. Tran (Ton Duc Thang University)
27 // p, H, He, He+ and He++ models are assumed identical
28 // NIMB 343, 132-137 (2015)
29 //
30 // The Geant4-DNA web site is available at http://geant4-dna.org
31 //
32 
33 #include "G4DNAIonElasticModel.hh"
34 #include "G4PhysicalConstants.hh"
35 #include "G4SystemOfUnits.hh"
37 #include "G4ParticleTable.hh"
38 #include "G4Exp.hh"
39 
40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
41 
42 using namespace std;
43 
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45 
47  const G4String& nam) :
48  G4VEmModel(nam), isInitialised(false)
49 {
50  //nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
51 
52  killBelowEnergy = 100 * eV;
53  lowEnergyLimit = 0 * eV;
54  highEnergyLimit = 1 * MeV;
55  SetLowEnergyLimit(lowEnergyLimit);
56  SetHighEnergyLimit(highEnergyLimit);
57 
58  verboseLevel = 0;
59  // Verbosity scale:
60  // 0 = nothing
61  // 1 = warning for energy non-conservation
62  // 2 = details of energy budget
63  // 3 = calculation of cross sections, file openings, sampling of atoms
64  // 4 = entering in methods
65 
66  if(verboseLevel > 0)
67  {
68  G4cout << "Ion elastic model is constructed " << G4endl<< "Energy range: "
69  << lowEnergyLimit / eV << " eV - "
70  << highEnergyLimit / MeV << " MeV"
71  << G4endl;
72  }
74  fpMolWaterDensity = 0;
75  fpTableData = 0;
76  fParticle_Mass = -1;
77 
78  // Selection of stationary mode
79 
80  statCode = false;
81 }
82 
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84 
86 {
87  // For total cross section
88  if(fpTableData) delete fpTableData;
89 }
90 
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92 
93 void
95  const G4ParticleDefinition* particleDefinition,
96  const G4DataVector& /*cuts*/)
97 {
98 
99  if(verboseLevel > 3)
100  {
101  G4cout << "Calling G4DNAIonElasticModel::Initialise()" << G4endl;
102  }
103 
104  // Energy limits
105 
106  if (LowEnergyLimit() < lowEnergyLimit)
107  {
108  G4cout << "G4DNAIonElasticModel: low energy limit increased from " <<
109  LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
110  SetLowEnergyLimit(lowEnergyLimit);
111  }
112 
113  if (HighEnergyLimit() > highEnergyLimit)
114  {
115  G4cout << "G4DNAIonElasticModel: high energy limit decreased from " <<
116  HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl;
117  SetHighEnergyLimit(highEnergyLimit);
118  }
119 
120  // Reading of data files
121 
122  G4double scaleFactor = 1e-16*cm*cm;
123 
124  char *path = getenv("G4LEDATA");
125 
126  if (!path)
127  {
128  G4Exception("G4IonElasticModel::Initialise","em0006",
129  FatalException,"G4LEDATA environment variable not set.");
130  return;
131  }
132 
133  G4String totalXSFile;
134  std::ostringstream fullFileName;
135 
138  G4ParticleDefinition* protonDef =
140  G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
141  G4ParticleDefinition* heliumDef = instance->GetIon("helium");
142  G4ParticleDefinition* alphaplusDef = instance->GetIon("alpha+");
143  G4ParticleDefinition* alphaplusplusDef = instance->GetIon("alpha++");
144  G4String proton, hydrogen, helium, alphaplus, alphaplusplus;
145 
146  if (
147  (particleDefinition == protonDef && protonDef != 0)
148  ||
149  (particleDefinition == hydrogenDef && hydrogenDef != 0)
150  )
151  {
152  // For total cross section of p,h
153  fParticle_Mass = 1.;
154  totalXSFile = "dna/sigma_elastic_proton_HTran";
155 
156  // For final state
157  fullFileName << path << "/dna/sigmadiff_cumulated_elastic_proton_HTran.dat";
158  }
159 
160  if (
161  (particleDefinition == instance->GetIon("helium") && heliumDef)
162  ||
163  (particleDefinition == instance->GetIon("alpha+") && alphaplusDef)
164  ||
165  (particleDefinition == instance->GetIon("alpha++") && alphaplusplusDef)
166  )
167  {
168  fParticle_Mass = 4.;
169  // For total cross section of he,he+,he++
170  totalXSFile = "dna/sigma_elastic_alpha_HTran";
171 
172  // For final state
173  fullFileName << path << "/dna/sigmadiff_cumulated_elastic_alpha_HTran.dat";
174  }
175 
176  fpTableData = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
177  fpTableData->LoadData(totalXSFile);
178  std::ifstream diffCrossSection(fullFileName.str().c_str());
179 
180  if (!diffCrossSection)
181  {
182  G4ExceptionDescription description;
183  description << "Missing data file:"
184  <<fullFileName.str().c_str()<< G4endl;
185  G4Exception("G4IonElasticModel::Initialise","em0003",
187  description);
188  }
189 
190  // Added clear for MT
191 
192  eTdummyVec.clear();
193  eVecm.clear();
194  fDiffCrossSectionData.clear();
195 
196  //
197 
198  eTdummyVec.push_back(0.);
199 
200  while(!diffCrossSection.eof())
201  {
202  double tDummy;
203  double eDummy;
204  diffCrossSection>>tDummy>>eDummy;
205 
206  // SI : mandatory eVecm initialization
207 
208  if (tDummy != eTdummyVec.back())
209  {
210  eTdummyVec.push_back(tDummy);
211  eVecm[tDummy].push_back(0.);
212  }
213 
214  diffCrossSection>>fDiffCrossSectionData[tDummy][eDummy];
215 
216  if (eDummy != eVecm[tDummy].back()) eVecm[tDummy].push_back(eDummy);
217 
218  }
219 
220  // End final state
221  if( verboseLevel>0 )
222  {
223  if (verboseLevel > 2)
224  {
225  G4cout << "Loaded cross section files for ion elastic model" << G4endl;
226  }
227  G4cout << "Ion elastic model is initialized " << G4endl
228  << "Energy range: "
229  << LowEnergyLimit() / eV << " eV - "
230  << HighEnergyLimit() / MeV << " MeV"
231  << G4endl;
232  }
233 
234  // Initialize water density pointer
236  fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->
237  GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
238 
239  if (isInitialised)
240  { return;}
242  isInitialised = true;
243 }
244 
245 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
246 
247 G4double
249  const G4ParticleDefinition* p,
250  G4double ekin, G4double, G4double)
251 {
252  if(verboseLevel > 3)
253  {
254  G4cout << "Calling CrossSectionPerVolume() of G4DNAIonElasticModel"
255  << G4endl;
256  }
257 
258  // Calculate total cross section for model
259 
260  G4double sigma=0;
261 
262  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
263 
264  if(waterDensity!= 0.0)
265  {
266  const G4String& particleName = p->GetParticleName();
267 
268  if (ekin < highEnergyLimit)
269  {
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  }
283  }
284 
285  if (verboseLevel > 2)
286  {
287  G4cout << "__________________________________" << G4endl;
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;
293  }
294 
295  }
296 
297  return sigma*waterDensity;
298 }
299 
300 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
301 
302 void
304  std::vector<G4DynamicParticle*>* /*fvect*/,
305  const G4MaterialCutsCouple* /*couple*/,
306  const G4DynamicParticle* aDynamicParticle, G4double, G4double)
307 {
308 
309  if(verboseLevel > 3)
310  {
311  G4cout << "Calling SampleSecondaries() of G4DNAIonElasticModel" << G4endl;
312  }
313 
314  G4double particleEnergy0 = aDynamicParticle->GetKineticEnergy();
315 
316  if (particleEnergy0 < killBelowEnergy)
317  {
321  return;
322  }
323 
324  if (particleEnergy0>= killBelowEnergy && particleEnergy0 < highEnergyLimit)
325  {
326  G4double water_mass = 18.;
327 
328  G4double thetaCM = RandomizeThetaCM(particleEnergy0, aDynamicParticle->GetDefinition());
329 
330  //HT:convert to laboratory system
331 
332  G4double theta = std::atan(std::sin(thetaCM*CLHEP::pi/180)
333  /(fParticle_Mass/water_mass+std::cos(thetaCM*CLHEP::pi/180)));
334 
335  G4double cosTheta= std::cos(theta);
336 
337  //
338 
339  G4double phi = 2. * CLHEP::pi * G4UniformRand();
340 
341  G4ThreeVector zVers = aDynamicParticle->GetMomentumDirection();
342  G4ThreeVector xVers = zVers.orthogonal();
343  G4ThreeVector yVers = zVers.cross(xVers);
344 
345  G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
346  G4double yDir = xDir;
347  xDir *= std::cos(phi);
348  yDir *= std::sin(phi);
349 
350  G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
351 
353 
354  G4double depositEnergyCM = 0;
355 
356  //HT: deposited energy
357  depositEnergyCM = 4. * particleEnergy0 * fParticle_Mass * water_mass *
358  (1-std::cos(thetaCM*CLHEP::pi/180))
359  / (2 * std::pow((fParticle_Mass+water_mass),2));
360 
361  //SI: added protection particleEnergy0 >= depositEnergyCM
362  if (!statCode && (particleEnergy0 >= depositEnergyCM) )
363 
364  fParticleChangeForGamma->SetProposedKineticEnergy(particleEnergy0 - depositEnergyCM);
365 
366  else fParticleChangeForGamma->SetProposedKineticEnergy(particleEnergy0);
367 
369  }
370 }
371 
372 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
373 
374 G4double
375 G4DNAIonElasticModel::Theta (G4ParticleDefinition * /*particleDefinition*/,
376  G4double k, G4double integrDiff)
377 {
378  G4double theta = 0.;
379  G4double valueT1 = 0;
380  G4double valueT2 = 0;
381  G4double valueE21 = 0;
382  G4double valueE22 = 0;
383  G4double valueE12 = 0;
384  G4double valueE11 = 0;
385  G4double xs11 = 0;
386  G4double xs12 = 0;
387  G4double xs21 = 0;
388  G4double xs22 = 0;
389 
390  std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),
391  eTdummyVec.end(), k);
392  std::vector<double>::iterator t1 = t2 - 1;
393 
394  std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),
395  eVecm[(*t1)].end(),
396  integrDiff);
397  std::vector<double>::iterator e11 = e12 - 1;
398 
399  std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),
400  eVecm[(*t2)].end(),
401  integrDiff);
402  std::vector<double>::iterator e21 = e22 - 1;
403 
404  valueT1 = *t1;
405  valueT2 = *t2;
406  valueE21 = *e21;
407  valueE22 = *e22;
408  valueE12 = *e12;
409  valueE11 = *e11;
410 
411  xs11 = fDiffCrossSectionData[valueT1][valueE11];
412  xs12 = fDiffCrossSectionData[valueT1][valueE12];
413  xs21 = fDiffCrossSectionData[valueT2][valueE21];
414  xs22 = fDiffCrossSectionData[valueT2][valueE22];
415 
416  if(xs11 == 0 && xs12 == 0 && xs21 == 0 && xs22 == 0) return (0.);
417 
418  theta = QuadInterpolator(valueE11, valueE12, valueE21, valueE22, xs11, xs12,
419  xs21, xs22, valueT1, valueT2, k, integrDiff);
420 
421  return theta;
422 }
423 
424 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
425 
426 G4double
427 G4DNAIonElasticModel::LinLinInterpolate (G4double e1, G4double e2, G4double e,
428  G4double xs1, G4double xs2)
429 {
430  G4double d1 = xs1;
431  G4double d2 = xs2;
432  G4double value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
433  return value;
434 }
435 
436 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
437 
438 G4double
439 G4DNAIonElasticModel::LinLogInterpolate (G4double e1, G4double e2, G4double e,
440  G4double xs1, G4double xs2)
441 {
442  G4double d1 = std::log(xs1);
443  G4double d2 = std::log(xs2);
444  G4double value = G4Exp(d1 + (d2 - d1) * (e - e1) / (e2 - e1));
445  return value;
446 }
447 
448 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
449 
450 G4double
451 G4DNAIonElasticModel::LogLogInterpolate (G4double e1, G4double e2, G4double e,
452  G4double xs1, G4double xs2)
453 {
454  G4double a = (std::log10(xs2) - std::log10(xs1))
455  / (std::log10(e2) - std::log10(e1));
456  G4double b = std::log10(xs2) - a * std::log10(e2);
457  G4double sigma = a * std::log10(e) + b;
458  G4double value = (std::pow(10., sigma));
459  return value;
460 }
461 
462 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
463 
464 G4double
465 G4DNAIonElasticModel::QuadInterpolator (G4double e11, G4double e12,
466  G4double e21, G4double e22,
467  G4double xs11, G4double xs12,
468  G4double xs21, G4double xs22,
469  G4double t1, G4double t2, G4double t,
470  G4double e)
471 {
472  // Log-Log
473  /*
474  G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
475  G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
476  G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
477  */
478 
479  // Lin-Log
480  /*
481  G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
482  G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
483  G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
484  */
485 
486 // Lin-Lin
487  G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
488  G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
489  G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1,
490  interpolatedvalue2);
491 
492  return value;
493 }
494 
495 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
496 
497 G4double
498 G4DNAIonElasticModel::RandomizeThetaCM (
499  G4double k, G4ParticleDefinition * particleDefinition)
500 {
501  G4double integrdiff = G4UniformRand();
502  return Theta(particleDefinition, k / eV, integrdiff);
503 }
504 
505 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
506 
507 void
509 {
510  killBelowEnergy = threshold;
511 
512  if(killBelowEnergy < 100 * eV)
513  {
514  G4cout << "*** WARNING : the G4DNAIonElasticModel class is not "
515  "activated below 100 eV !"
516  << G4endl;
517  }
518 }
519 
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:636
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:602
size_t GetIndex() const
Definition: G4Material.hh:262
const char * p
Definition: xmltok.h:285
static const G4double d2
virtual G4bool LoadData(const G4String &argFileName)
G4ParticleDefinition * GetDefinition() const
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4String & GetParticleName() const
void SetKillBelowThreshold(G4double threshold)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:732
string material
Definition: eplot.py:19
tuple b
Definition: test.py:12
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const XML_Char int const XML_Char * value
Definition: expat.h:331
G4ParticleChangeForGamma * fParticleChangeForGamma
const G4ThreeVector & GetMomentumDirection() const
virtual void Initialise(const G4ParticleDefinition *particuleDefinition, const G4DataVector &)
static constexpr double cm
Definition: G4SIunits.hh:119
static constexpr double eV
Definition: G4SIunits.hh:215
static G4DNAGenericIonsManager * Instance(void)
virtual G4double FindValue(G4double e, G4int componentId=0) const
static const G4double d1
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static G4DNAMolecularMaterial * Instance()
static G4ParticleTable * GetParticleTable()
Hep3Vector unit() const
tuple t1
Definition: plottest35.py:33
Hep3Vector orthogonal() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
static MCTruthManager * instance
Hep3Vector cross(const Hep3Vector &) const
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:739
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
#define DBL_MAX
Definition: templates.hh:83
static constexpr double pi
Definition: SystemOfUnits.h:54
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:132
G4DNAIonElasticModel(const G4ParticleDefinition *p=0, const G4String &nam="DNAIonElasticModel")
G4ParticleDefinition * GetIon(const G4String &name)