ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNAIonElasticModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNAIonElasticModel.cc
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"
36 #include "G4DNAMolecularMaterial.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 
46 G4DNAIonElasticModel::G4DNAIonElasticModel (const G4ParticleDefinition*,
47  const G4String& nam) :
48  G4VEmModel(nam), isInitialised(false)
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 openings, sampling of atoms
62  // 4 = entering in methods
63 
64  if(verboseLevel > 0)
65  {
66  G4cout << "Ion elastic model is constructed " << G4endl<< "Energy range: "
67  << lowEnergyLimit / eV << " eV - "
68  << highEnergyLimit / MeV << " MeV"
69  << G4endl;
70  }
71 
72  fParticleChangeForGamma = 0;
73  fpMolWaterDensity = 0;
74  fpTableData = 0;
75  fParticle_Mass = -1;
76 
77  // Selection of stationary mode
78 
79  statCode = false;
80 }
81 
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
83 
84 G4DNAIonElasticModel::~G4DNAIonElasticModel ()
85 {
86  // For total cross section
87  if(fpTableData) delete fpTableData;
88 }
89 
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91 
92 void
93 G4DNAIonElasticModel::Initialise (
94  const G4ParticleDefinition* particleDefinition,
95  const G4DataVector& /*cuts*/)
96 {
97 
98  if(verboseLevel > 3)
99  {
100  G4cout << "Calling G4DNAIonElasticModel::Initialise()" << G4endl;
101  }
102 
103  // Energy limits
104 
105  if (LowEnergyLimit() < lowEnergyLimit)
106  {
107  G4cout << "G4DNAIonElasticModel: low energy limit increased from " <<
108  LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
109  SetLowEnergyLimit(lowEnergyLimit);
110  }
111 
112  if (HighEnergyLimit() > highEnergyLimit)
113  {
114  G4cout << "G4DNAIonElasticModel: high energy limit decreased from " <<
115  HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl;
116  SetHighEnergyLimit(highEnergyLimit);
117  }
118 
119  // Reading of data files
120 
121  G4double scaleFactor = 1e-16*cm*cm;
122 
123  char *path = getenv("G4LEDATA");
124 
125  if (!path)
126  {
127  G4Exception("G4IonElasticModel::Initialise","em0006",
128  FatalException,"G4LEDATA environment variable not set.");
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()->FindParticle("proton");
139  G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
140  G4ParticleDefinition* heliumDef = instance->GetIon("helium");
141  G4ParticleDefinition* alphaplusDef = instance->GetIon("alpha+");
142  G4ParticleDefinition* alphaplusplusDef = instance->GetIon("alpha++");
143  G4String proton, hydrogen, helium, alphaplus, alphaplusplus;
144 
145  if (
146  (particleDefinition == protonDef && protonDef != 0)
147  ||
148  (particleDefinition == hydrogenDef && hydrogenDef != 0)
149  )
150  {
151  // For total cross section of p,h
152  fParticle_Mass = 1.;
153  totalXSFile = "dna/sigma_elastic_proton_HTran";
154 
155  // For final state
156  fullFileName << path << "/dna/sigmadiff_cumulated_elastic_proton_HTran.dat";
157  }
158 
159  if (
160  (particleDefinition == instance->GetIon("helium") && heliumDef)
161  ||
162  (particleDefinition == instance->GetIon("alpha+") && alphaplusDef)
163  ||
164  (particleDefinition == instance->GetIon("alpha++") && alphaplusplusDef)
165  )
166  {
167  // For total cross section of he,he+,he++
168  fParticle_Mass = 4.;
169  totalXSFile = "dna/sigma_elastic_alpha_HTran";
170 
171  // For final state
172  fullFileName << path << "/dna/sigmadiff_cumulated_elastic_alpha_HTran.dat";
173  }
174 
175  fpTableData = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
176  fpTableData->LoadData(totalXSFile);
177  std::ifstream diffCrossSection(fullFileName.str().c_str());
178 
179  if (!diffCrossSection)
180  {
181  G4ExceptionDescription description;
182  description << "Missing data file:"
183  <<fullFileName.str().c_str()<< G4endl;
184  G4Exception("G4IonElasticModel::Initialise","em0003",
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[tDummy][eDummy];
214 
215  if (eDummy != eVecm[tDummy].back()) eVecm[tDummy].push_back(eDummy);
216 
217  }
218 
219  // End final state
220  if( verboseLevel>0 )
221  {
222  if (verboseLevel > 2)
223  {
224  G4cout << "Loaded cross section files for ion elastic model" << G4endl;
225  }
226  G4cout << "Ion elastic model is initialized " << G4endl
227  << "Energy range: "
228  << LowEnergyLimit() / eV << " eV - "
229  << HighEnergyLimit() / MeV << " MeV"
230  << G4endl;
231  }
232 
233  // Initialize water density pointer
234  G4DNAMolecularMaterial::Instance()->Initialize();
235  fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->
236  GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
237 
238  if (isInitialised) return;
239  fParticleChangeForGamma = GetParticleChangeForGamma();
240  isInitialised = true;
241 }
242 
243 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
244 
245 G4double
246 G4DNAIonElasticModel::CrossSectionPerVolume (const G4Material* material,
247  const G4ParticleDefinition* p,
248  G4double ekin, G4double, G4double)
249 {
250  if(verboseLevel > 3)
251  {
252  G4cout << "Calling CrossSectionPerVolume() of G4DNAIonElasticModel"
253  << G4endl;
254  }
255 
256  // Calculate total cross section for model
257 
258  G4double sigma=0;
259 
260  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
261 
262  const G4String& particleName = p->GetParticleName();
263 
264  if (ekin <= highEnergyLimit)
265  {
266  //SI : XS must not be zero otherwise sampling of secondaries method ignored
267  if (ekin < killBelowEnergy) return DBL_MAX;
268  //
269 
270  if (fpTableData != 0)
271  {
272  sigma = fpTableData->FindValue(ekin);
273  }
274  else
275  {
276  G4Exception("G4DNAIonElasticModel::ComputeCrossSectionPerVolume","em0002",
277  FatalException,"Model not applicable to particle type.");
278  }
279  }
280 
281  if (verboseLevel > 2)
282  {
283  G4cout << "__________________________________" << G4endl;
284  G4cout << "G4DNAIonElasticModel - XS INFO START" << G4endl;
285  G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
286  G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
287  G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
288  G4cout << "G4DNAIonElasticModel - XS INFO END" << G4endl;
289  }
290 
291  return sigma*waterDensity;
292 }
293 
294 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
295 
296 void
297 G4DNAIonElasticModel::SampleSecondaries (
298  std::vector<G4DynamicParticle*>* /*fvect*/,
299  const G4MaterialCutsCouple* /*couple*/,
300  const G4DynamicParticle* aDynamicParticle, G4double, G4double)
301 {
302 
303  if(verboseLevel > 3)
304  {
305  G4cout << "Calling SampleSecondaries() of G4DNAIonElasticModel" << G4endl;
306  }
307 
308  G4double particleEnergy0 = aDynamicParticle->GetKineticEnergy();
309 
310  if (particleEnergy0 < killBelowEnergy)
311  {
312  fParticleChangeForGamma->SetProposedKineticEnergy(0.);
313  fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
314  fParticleChangeForGamma->ProposeLocalEnergyDeposit(particleEnergy0);
315  return;
316  }
317 
318  if (particleEnergy0>= killBelowEnergy && particleEnergy0 <= highEnergyLimit)
319  {
320  G4double water_mass = 18.;
321 
322  G4double thetaCM = RandomizeThetaCM(particleEnergy0, aDynamicParticle->GetDefinition());
323 
324  //HT:convert to laboratory system
325 
326  G4double theta = std::atan(std::sin(thetaCM*CLHEP::pi/180)
327  /(fParticle_Mass/water_mass+std::cos(thetaCM*CLHEP::pi/180)));
328 
329  G4double cosTheta= std::cos(theta);
330 
331  //
332 
333  G4double phi = 2. * CLHEP::pi * G4UniformRand();
334 
335  G4ThreeVector zVers = aDynamicParticle->GetMomentumDirection();
336  G4ThreeVector xVers = zVers.orthogonal();
337  G4ThreeVector yVers = zVers.cross(xVers);
338 
339  G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
340  G4double yDir = xDir;
341  xDir *= std::cos(phi);
342  yDir *= std::sin(phi);
343 
344  G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
345 
346  fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit());
347 
348  G4double depositEnergyCM = 0;
349 
350  //HT: deposited energy
351  depositEnergyCM = 4. * particleEnergy0 * fParticle_Mass * water_mass *
352  (1-std::cos(thetaCM*CLHEP::pi/180))
353  / (2 * std::pow((fParticle_Mass+water_mass),2));
354 
355  //SI: added protection particleEnergy0 >= depositEnergyCM
356  if (!statCode && (particleEnergy0 >= depositEnergyCM) )
357 
358  fParticleChangeForGamma->SetProposedKineticEnergy(particleEnergy0 - depositEnergyCM);
359 
360  else fParticleChangeForGamma->SetProposedKineticEnergy(particleEnergy0);
361 
362  fParticleChangeForGamma->ProposeLocalEnergyDeposit(depositEnergyCM);
363  }
364 }
365 
366 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
367 
368 G4double
369 G4DNAIonElasticModel::Theta (G4ParticleDefinition * /*particleDefinition*/,
370  G4double k, G4double integrDiff)
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  std::vector<G4double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),
389  eTdummyVec.end(), k);
390  std::vector<G4double>::iterator t1 = t2 - 1;
391 
392  std::vector<G4double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),
393  eVecm[(*t1)].end(),
394  integrDiff);
395  std::vector<G4double>::iterator e11 = e12 - 1;
396 
397  std::vector<G4double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),
398  eVecm[(*t2)].end(),
399  integrDiff);
400  std::vector<G4double>::iterator 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][valueE11];
410  xs12 = fDiffCrossSectionData[valueT1][valueE12];
411  xs21 = fDiffCrossSectionData[valueT2][valueE21];
412  xs22 = fDiffCrossSectionData[valueT2][valueE22];
413 
414  if(xs11 == 0 && xs12 == 0 && xs21 == 0 && xs22 == 0) return (0.);
415 
416  theta = QuadInterpolator(valueE11, valueE12, valueE21, valueE22, xs11, xs12,
417  xs21, xs22, valueT1, valueT2, k, integrDiff);
418 
419  return theta;
420 }
421 
422 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
423 
424 G4double
425 G4DNAIonElasticModel::LinLinInterpolate (G4double e1, G4double e2, G4double e,
426  G4double xs1, G4double xs2)
427 {
428  G4double d1 = xs1;
429  G4double d2 = xs2;
430  G4double value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
431  return value;
432 }
433 
434 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
435 
436 G4double
437 G4DNAIonElasticModel::LinLogInterpolate (G4double e1, G4double e2, G4double e,
438  G4double xs1, G4double xs2)
439 {
440  G4double d1 = std::log(xs1);
441  G4double d2 = std::log(xs2);
442  G4double value = G4Exp(d1 + (d2 - d1) * (e - e1) / (e2 - e1));
443  return value;
444 }
445 
446 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
447 
448 G4double
449 G4DNAIonElasticModel::LogLogInterpolate (G4double e1, G4double e2, G4double e,
450  G4double xs1, G4double xs2)
451 {
452  G4double a = (std::log10(xs2) - std::log10(xs1))
453  / (std::log10(e2) - std::log10(e1));
454  G4double b = std::log10(xs2) - a * std::log10(e2);
455  G4double sigma = a * std::log10(e) + b;
456  G4double value = (std::pow(10., sigma));
457  return value;
458 }
459 
460 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
461 
462 G4double
463 G4DNAIonElasticModel::QuadInterpolator (G4double e11, G4double e12,
464  G4double e21, G4double e22,
465  G4double xs11, G4double xs12,
466  G4double xs21, G4double xs22,
467  G4double t1, G4double t2, G4double t,
468  G4double e)
469 {
470  // Log-Log
471  /*
472  G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
473  G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
474  G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
475  */
476 
477  // Lin-Log
478  /*
479  G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
480  G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
481  G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
482  */
483 
484 // Lin-Lin
485  G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
486  G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
487  G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1,
488  interpolatedvalue2);
489 
490  return value;
491 }
492 
493 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
494 
495 G4double
496 G4DNAIonElasticModel::RandomizeThetaCM (
497  G4double k, G4ParticleDefinition * particleDefinition)
498 {
499  G4double integrdiff = G4UniformRand();
500  return Theta(particleDefinition, k / eV, integrdiff);
501 }
502 
503 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
504 
505 void
506 G4DNAIonElasticModel::SetKillBelowThreshold (G4double threshold)
507 {
508  killBelowEnergy = threshold;
509 
510  if(killBelowEnergy < 100 * eV)
511  {
512  G4cout << "*** WARNING : the G4DNAIonElasticModel class is not "
513  "activated below 100 eV !"
514  << G4endl;
515  }
516 }
517