| Geant4 Cross Reference |
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>
39 #include "G4DNAOneStepThermalizationModel.hh"
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
52 Meesungnoen2002::gCoeff[13] =
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
61 Meesungnoen2002_amorphous::gCoeff[7] =
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
68 Terrisol1990::gEnergies_T1990[11] =
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
74 Terrisol1990::gStdDev_T1990[11] =
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
95 double Meesungnoen2002::GetRmean(double k){
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
109 double Meesungnoen2002_amorphous::GetRmean(double k){
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
123 void GetGaussianPenetrationFromRmean3D(G4double r_mean,
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
144 void Meesungnoen2002::GetPenetration(G4double k,
145 G4ThreeVector& displacement)
146 {
147 GetGaussianPenetrationFromRmean3D(GetRmean(k), displacement);
148 }
149
150 void Meesungnoen2002_amorphous::GetPenetration(G4double k,
151 G4ThreeVector& displacement)
152 {
153 GetGaussianPenetrationFromRmean3D(GetRmean(k), displacement);
154 }
155
156
157 void Kreipl2009::GetPenetration(G4double k,
158 G4ThreeVector& displacement)
159 {
160 G4double r_mean = Meesungnoen2002::GetRmean(k);
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
176 void Ritchie1994::GetPenetration(G4double k,
177 G4ThreeVector& displacement)
178 {
179 GetGaussianPenetrationFromRmean3D(k/eV * 1.8 * nm, // r_mean
180 displacement);
181 }
182
183 //----------------------------------------------------------------------------
184
185 double Terrisol1990::Get3DStdDeviation(double energy){
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 if(k_eV == 9.){
193 return gStdDev_T1990[10];
194 }
195 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",
200 FatalErrorInArgument,
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
229 double Terrisol1990::GetRmean(double energy){
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
238 void Terrisol1990::GetPenetration(G4double energy,
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
256 G4VEmModel* G4DNASolvationModelFactory::Create(const G4String& penetrationModel)
257 {
258 G4String modelNamePrefix("DNAOneStepThermalizationModel_");
259
260 if(penetrationModel == "Terrisol1990")
261 {
262 return new G4TDNAOneStepThermalizationModel<DNA::Penetration::Terrisol1990>(G4Electron::Definition(), modelNamePrefix + penetrationModel);
263 }
264 if(penetrationModel == "Meesungnoen2002")
265 {
266 return new G4TDNAOneStepThermalizationModel<DNA::Penetration::Meesungnoen2002>(G4Electron::Definition(), modelNamePrefix + penetrationModel);
267 }
268 if(penetrationModel == "Meesungnoen2002_amorphous")
269 {
270 return new G4TDNAOneStepThermalizationModel<DNA::Penetration::Meesungnoen2002_amorphous>(G4Electron::Definition(), modelNamePrefix + penetrationModel);
271 }
272 if(penetrationModel == "Kreipl2009")
273 {
274 return new G4TDNAOneStepThermalizationModel<DNA::Penetration::Kreipl2009>(G4Electron::Definition(), modelNamePrefix + penetrationModel);
275 }
276 if(penetrationModel == "Ritchie1994")
277 {
278 return new G4TDNAOneStepThermalizationModel<DNA::Penetration::Ritchie1994>(G4Electron::Definition(), modelNamePrefix + penetrationModel);
279 }
280
281 G4ExceptionDescription description;
282 description << penetrationModel + " is not a valid model name.";
283 G4Exception("G4DNASolvationModelFactory::Create",
284 "INVALID_ARGUMENT",
285 FatalErrorInArgument,
286 description,
287 "Options are: Terrisol1990, Meesungnoen2002, Ritchie1994.");
288 return nullptr;
289 }
290
291 //------------------------------------------------------------------------------
292 G4VEmModel* G4DNASolvationModelFactory::GetMacroDefinedModel()
293 {
294 auto dnaSubType = G4EmParameters::Instance()->DNAeSolvationSubType();
295
296 switch(dnaSubType)
297 {
298 case fRitchie1994eSolvation:
299 return Create("Ritchie1994");
300 case fTerrisol1990eSolvation:
301 return Create("Terrisol1990");
302 case fKreipl2009eSolvation:
303 return Create("Kreipl2009");
304 case fMeesungnoensolid2002eSolvation:
305 return Create("Meesungnoen2002_amorphous");
306 case fMeesungnoen2002eSolvation:
307 case fDNAUnknownModel:
308 return Create("Meesungnoen2002");
309 default:
310 G4Exception("G4DNASolvationModelFactory::GetMacroDefinedModel",
311 "DnaSubType",
312 FatalErrorInArgument,
313 "The solvation parameter stored in G4EmParameters is unknown. Supported types are: fRitchie1994eSolvation, fTerrisol1990eSolvation, fMeesungnoen2002eSolvation.");
314 }
315
316 return nullptr;
317 }
318