| 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 //
27 // Author: Mathieu Karamitros
28 //
29 // WARNING : This class is released as a proto
30 // It might strongly evolve or even disapear i
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.932178
54 1.80172797e-03, -2.01135480e-02, 1.429394
55 -6.48348714e-01, 1.85227848e+00, -3.364503
56 4.37785068e+00, -4.20557339e+00, 3.816790
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 t
85 // b=27.22 nm nor the mean value. 129.62*CLH
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(dou
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(G4doubl
124 G4Three
125 {
126 if(r_mean == 0)
127 {
128 // rare events:
129 // prevent H2O and secondary electron from
130 displacement = G4RandomDirection() * (1e-3
131 return;
132 }
133
134 static constexpr double convertRmean3DToSigm
135 // = sqrt(CLHEP::pi)/pow(2,3./2.)
136
137 // Use r_mean to build a 3D gaussian
138 const double sigma1D = r_mean * convertRmean
139 displacement = G4ThreeVector(G4RandGauss::sh
140 G4RandGauss::sh
141 G4RandGauss::sh
142 }
143
144 void Meesungnoen2002::GetPenetration(G4double
145 G4ThreeVe
146 {
147 GetGaussianPenetrationFromRmean3D(GetRmean(k
148 }
149
150 void Meesungnoen2002_amorphous::GetPenetration
151 G4ThreeVe
152 {
153 GetGaussianPenetrationFromRmean3D(GetRmean(k
154 }
155
156
157 void Kreipl2009::GetPenetration(G4double k,
158 G4ThreeVe
159 {
160 G4double r_mean = Meesungnoen2002::GetRmean(
161
162 if(r_mean == 0)
163 {
164 // rare events:
165 // prevent H2O and secondary electron from b
166 displacement = G4RandomDirection() * (1e-3*C
167 return;
168 }
169
170 double r = G4RandGamma::shoot(2,2);
171
172 displacement = G4RandomDirection() * r * r_m
173 }
174 //--------------------------------------------
175
176 void Ritchie1994::GetPenetration(G4double k,
177 G4ThreeVector
178 {
179 GetGaussianPenetrationFromRmean3D(k/eV * 1.8
180 displaceme
181 }
182
183 //--------------------------------------------
184
185 double Terrisol1990::Get3DStdDeviation(double
186 G4double k_eV = energy/eV;
187 if(k_eV < 0.2){
188 // rare events:
189 // prevent H2O and secondary electron to
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 tabula
198 G4Exception("Terrisol1990::Get3DStdDeviati
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[
212 &gEnergies_T1990[
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.59576912160573
233
234 double r_mean=sigma3D*s2r;
235 return r_mean;
236 }
237
238 void Terrisol1990::GetPenetration(G4double ene
239 G4ThreeVecto
240 double sigma3D = Get3DStdDeviation(energy);
241
242 static constexpr double factor = 2.204969995
243
244 double sigma1D = std::sqrt(std::pow(sigma3D,
245
246 displacement = G4ThreeVector(G4RandGauss::sh
247 G4RandGauss::sh
248 G4RandGauss::sh
249 }
250
251 } // Penetration
252 } // DNA
253
254 //--------------------------------------------
255
256 G4VEmModel* G4DNASolvationModelFactory::Create
257 {
258 G4String modelNamePrefix("DNAOneStepThermali
259
260 if(penetrationModel == "Terrisol1990")
261 {
262 return new G4TDNAOneStepThermalizationMode
263 }
264 if(penetrationModel == "Meesungnoen2002")
265 {
266 return new G4TDNAOneStepThermalizationMode
267 }
268 if(penetrationModel == "Meesungnoen2002_amor
269 {
270 return new G4TDNAOneStepThermalizationModel<
271 }
272 if(penetrationModel == "Kreipl2009")
273 {
274 return new G4TDNAOneStepThermalizationModel<
275 }
276 if(penetrationModel == "Ritchie1994")
277 {
278 return new G4TDNAOneStepThermalizationMode
279 }
280
281 G4ExceptionDescription description;
282 description << penetrationModel + " is not a
283 G4Exception("G4DNASolvationModelFactory::Cre
284 "INVALID_ARGUMENT",
285 FatalErrorInArgument,
286 description,
287 "Options are: Terrisol1990, Mees
288 return nullptr;
289 }
290
291 //--------------------------------------------
292 G4VEmModel* G4DNASolvationModelFactory::GetMac
293 {
294 auto dnaSubType = G4EmParameters::Instance()
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::G
311 "DnaSubType",
312 FatalErrorInArgument,
313 "The solvation parameter store
314 }
315
316 return nullptr;
317 }
318