ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNAOneStepThermalizationModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNAOneStepThermalizationModel.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 //
27 // Author: Mathieu Karamitros
28 //
29 // WARNING : This class is released as a prototype.
30 // It might strongly evolve or even disapear in the next releases.
31 //
32 // History:
33 // -----------
34 // 10 Oct 2011 M.Karamitros created
35 //
36 // -------------------------------------------------------------------
37 
38 #include <algorithm>
40 #include "globals.hh"
41 #include "G4Exp.hh"
42 #include "G4RandomDirection.hh"
43 #include "G4Electron.hh"
44 #include "G4EmParameters.hh"
45 
46 //------------------------------------------------------------------------------
47 
48 namespace DNA {
49 namespace Penetration {
50 
51 const double
53 { -4.06217193e-08, 3.06848412e-06, -9.93217814e-05,
54  1.80172797e-03, -2.01135480e-02, 1.42939448e-01,
55  -6.48348714e-01, 1.85227848e+00, -3.36450378e+00,
56  4.37785068e+00, -4.20557339e+00, 3.81679083e+00,
57  -2.34069784e-01 };
58 // fit from Meesungnoen, 2002
59 
60 const double
62 { 7.3144e-05, -2.2474e-03, 3.4555e-02,
63  -4.3574e-01, 2.8954e+00, -1.0381e+00,
64  1.4300e+00 };
65 // fit from Meesungnoen, 2002
66 
67 const double
69 { 0.2, 0.5, 1, 2, 3, 4, 5, 6, 7,
70  // The two last are not in the dataset
71  8, 9}; // eV
72 
73 const double
75 { 17.68*CLHEP::angstrom,
76  22.3*CLHEP::angstrom,
77  28.49*CLHEP::angstrom,
78  45.35*CLHEP::angstrom,
79  70.03*CLHEP::angstrom,
80  98.05*CLHEP::angstrom,
81  120.56*CLHEP::angstrom,
82  132.73*CLHEP::angstrom,
83  142.60*CLHEP::angstrom,
84  // the above value as given in the paper's table does not match
85  // b=27.22 nm nor the mean value. 129.62*CLHEP::angstrom could be
86  // a better fit.
87  //
88  // The two last are made up
89  137.9*CLHEP::angstrom,
90  120.7*CLHEP::angstrom
91 }; // angstrom
92 
93 //----------------------------------------------------------------------------
94 
96  G4double k_eV = k/eV;
97 
98  if(k_eV>0.1){ // data until 0.2 eV
99  G4double r_mean = 0;
100  for(int8_t i=12; i!=-1 ; --i){
101  r_mean+=gCoeff[12-i]*std::pow(k_eV,i);
102  }
103  r_mean*=CLHEP::nanometer;
104  return r_mean;
105  }
106  return 0;
107 }
108 
110  G4double k_eV = k/eV;
111 
112  if(k_eV>0.1){ // data until 0.2 eV
113  G4double r_mean = 0;
114  for(int8_t i=6; i!=-1 ; --i){
115  r_mean+=gCoeff[6-i]*std::pow(k_eV,i);
116  }
117  r_mean*=CLHEP::nanometer;
118  return r_mean;
119  }
120  return 0;
121 }
122 
124  G4ThreeVector& displacement)
125 {
126  if(r_mean == 0)
127  {
128  // rare events:
129  // prevent H2O and secondary electron from being placed at the same position
130  displacement = G4RandomDirection() * (1e-3*CLHEP::nanometer);
131  return;
132  }
133 
134  static constexpr double convertRmean3DToSigma1D = 0.62665706865775006;
135  // = sqrt(CLHEP::pi)/pow(2,3./2.)
136 
137  // Use r_mean to build a 3D gaussian
138  const double sigma1D = r_mean * convertRmean3DToSigma1D;
139  displacement = G4ThreeVector(G4RandGauss::shoot(0, sigma1D),
140  G4RandGauss::shoot(0, sigma1D),
141  G4RandGauss::shoot(0, sigma1D));
142 }
143 
145  G4ThreeVector& displacement)
146 {
148 }
149 
151  G4ThreeVector& displacement)
152 {
154 }
155 
156 
158  G4ThreeVector& displacement)
159 {
161 
162  if(r_mean == 0)
163  {
164  // rare events:
165  // prevent H2O and secondary electron from being placed at the same position
166  displacement = G4RandomDirection() * (1e-3*CLHEP::nanometer);
167  return;
168  }
169 
170  double r = G4RandGamma::shoot(2,2);
171 
172  displacement = G4RandomDirection() * r * r_mean;
173 }
174 //----------------------------------------------------------------------------
175 
177  G4ThreeVector& displacement)
178 {
179  GetGaussianPenetrationFromRmean3D(k/eV * 1.8 * nm, // r_mean
180  displacement);
181 }
182 
183 //----------------------------------------------------------------------------
184 
186  G4double k_eV = energy/eV;
187  if(k_eV < 0.2){
188  // rare events:
189  // prevent H2O and secondary electron to be at the spot
190  return 1e-3*CLHEP::nanometer;
191  }
192  else if(k_eV == 9.){
193  return gStdDev_T1990[10];
194  }
195  else if(k_eV > 9.){
196  G4ExceptionDescription description;
197  description << "Terrisol1990 is not tabulated for energies greater than 9eV";
198  G4Exception("Terrisol1990::Get3DStdDeviation",
199  "INVALID_ARGUMENT",
201  description);
202  }
203 
204  size_t lowBin, upBin;
205 
206  if(k_eV >= 1.){
207  lowBin=std::floor(k_eV)+1;
208  upBin=std::min(lowBin+1, size_t(10));
209  }
210  else{
211  auto it=std::lower_bound(&gEnergies_T1990[0],
212  &gEnergies_T1990[2],
213  k_eV);
214  lowBin = it-&gEnergies_T1990[0];
215  upBin = lowBin+1;
216  }
217 
218  double lowE = gEnergies_T1990[lowBin];
219  double upE = gEnergies_T1990[upBin];
220 
221  double lowS = gStdDev_T1990[lowBin];
222  double upS = gStdDev_T1990[upBin];
223 
224  double tanA = (lowS-upS)/(lowE-upE);
225  double sigma3D = lowS + (k_eV-lowE)*tanA;
226  return sigma3D;
227 }
228 
230  double sigma3D=Get3DStdDeviation(energy);
231 
232  static constexpr double s2r=1.595769121605731; // = pow(2,3./2.)/sqrt(CLHEP::pi)
233 
234  double r_mean=sigma3D*s2r;
235  return r_mean;
236 }
237 
239  G4ThreeVector& displacement){
240  double sigma3D = Get3DStdDeviation(energy);
241 
242  static constexpr double factor = 2.20496999539; // = 1./(3. - 8./CLHEP::pi);
243 
244  double sigma1D = std::sqrt(std::pow(sigma3D, 2.)*factor);
245 
246  displacement = G4ThreeVector(G4RandGauss::shoot(0, sigma1D),
247  G4RandGauss::shoot(0, sigma1D),
248  G4RandGauss::shoot(0, sigma1D));
249 }
250 
251 } // Penetration
252 } // DNA
253 
254 //------------------------------------------------------------------------------
255 
257 {
258  G4String modelNamePrefix("DNAOneStepThermalizationModel_");
259 
260  if(penetrationModel == "Terrisol1990")
261  {
262  return new G4TDNAOneStepThermalizationModel<DNA::Penetration::Terrisol1990>(G4Electron::Definition(), modelNamePrefix + penetrationModel);
263  }
264  else if(penetrationModel == "Meesungnoen2002")
265  {
267  }
268  else if(penetrationModel == "Meesungnoen2002_amorphous")
269  {
271  }
272  else if(penetrationModel == "Kreipl2009")
273  {
274  return new G4TDNAOneStepThermalizationModel<DNA::Penetration::Kreipl2009>(G4Electron::Definition(), modelNamePrefix + penetrationModel);
275  }
276  else if(penetrationModel == "Ritchie1994")
277  {
278  return new G4TDNAOneStepThermalizationModel<DNA::Penetration::Ritchie1994>(G4Electron::Definition(), modelNamePrefix + penetrationModel);
279  }
280  else
281  {
282  G4ExceptionDescription description;
283  description << penetrationModel + " is not a valid model name.";
284  G4Exception("G4DNASolvationModelFactory::Create",
285  "INVALID_ARGUMENT",
287  description,
288  "Options are: Terrisol1990, Meesungnoen2002, Ritchie1994.");
289  }
290  return nullptr;
291 }
292 
293 //------------------------------------------------------------------------------
295 {
296  auto dnaSubType = G4EmParameters::Instance()->DNAeSolvationSubType();
297 
298  switch(dnaSubType)
299  {
301  return Create("Ritchie1994");
303  return Create("Terrisol1990");
305  return Create("Kreipl2009");
307  return Create("Meesungnoen2002_amorphous");
309  case fDNAUnknownModel:
310  return Create("Meesungnoen2002");
311  default:
312  G4Exception("G4DNASolvationModelFactory::GetMacroDefinedModel",
313  "DnaSubType",
315  "The solvation parameter stored in G4EmParameters is unknown. Supported types are: fRitchie1994eSolvation, fTerrisol1990eSolvation, fMeesungnoen2002eSolvation.");
316  }
317 
318  return nullptr;
319 }