| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // Author: Mathieu Karamitros 27 // Author: Mathieu Karamitros
28 // 28 //
29 // WARNING : This class is released as a proto 29 // WARNING : This class is released as a prototype.
30 // It might strongly evolve or even disapear i 30 // It might strongly evolve or even disapear in the next releases.
31 // 31 //
32 // History: 32 // History:
33 // ----------- 33 // -----------
34 // 10 Oct 2011 M.Karamitros created 34 // 10 Oct 2011 M.Karamitros created
35 // 35 //
36 // ------------------------------------------- 36 // -------------------------------------------------------------------
37 37
38 #include <algorithm> 38 #include <algorithm>
39 #include "G4DNAOneStepThermalizationModel.hh" 39 #include "G4DNAOneStepThermalizationModel.hh"
40 #include "globals.hh" 40 #include "globals.hh"
41 #include "G4Exp.hh" 41 #include "G4Exp.hh"
42 #include "G4RandomDirection.hh" 42 #include "G4RandomDirection.hh"
43 #include "G4Electron.hh" 43 #include "G4Electron.hh"
44 #include "G4EmParameters.hh" 44 #include "G4EmParameters.hh"
45 45
46 //-------------------------------------------- 46 //------------------------------------------------------------------------------
47 47
48 namespace DNA { 48 namespace DNA {
49 namespace Penetration { 49 namespace Penetration {
50 50
51 const double 51 const double
52 Meesungnoen2002::gCoeff[13] = 52 Meesungnoen2002::gCoeff[13] =
53 { -4.06217193e-08, 3.06848412e-06, -9.932178 53 { -4.06217193e-08, 3.06848412e-06, -9.93217814e-05,
54 1.80172797e-03, -2.01135480e-02, 1.429394 54 1.80172797e-03, -2.01135480e-02, 1.42939448e-01,
55 -6.48348714e-01, 1.85227848e+00, -3.364503 55 -6.48348714e-01, 1.85227848e+00, -3.36450378e+00,
56 4.37785068e+00, -4.20557339e+00, 3.816790 56 4.37785068e+00, -4.20557339e+00, 3.81679083e+00,
57 -2.34069784e-01 }; 57 -2.34069784e-01 };
58 // fit from Meesungnoen, 2002 58 // fit from Meesungnoen, 2002
59 59
60 const double 60 const double
61 Meesungnoen2002_amorphous::gCoeff[7] = 61 Meesungnoen2002_amorphous::gCoeff[7] =
62 { 7.3144e-05, -2.2474e-03, 3.4555e-02, 62 { 7.3144e-05, -2.2474e-03, 3.4555e-02,
63 -4.3574e-01, 2.8954e+00, -1.0381e+00, 63 -4.3574e-01, 2.8954e+00, -1.0381e+00,
64 1.4300e+00 }; 64 1.4300e+00 };
65 // fit from Meesungnoen, 2002 65 // fit from Meesungnoen, 2002
66 66
67 const double 67 const double
68 Terrisol1990::gEnergies_T1990[11] = 68 Terrisol1990::gEnergies_T1990[11] =
69 { 0.2, 0.5, 1, 2, 3, 4, 5, 6, 7, 69 { 0.2, 0.5, 1, 2, 3, 4, 5, 6, 7,
70 // The two last are not in the dataset 70 // The two last are not in the dataset
71 8, 9}; // eV 71 8, 9}; // eV
72 72
73 const double 73 const double
74 Terrisol1990::gStdDev_T1990[11] = 74 Terrisol1990::gStdDev_T1990[11] =
75 { 17.68*CLHEP::angstrom, 75 { 17.68*CLHEP::angstrom,
76 22.3*CLHEP::angstrom, 76 22.3*CLHEP::angstrom,
77 28.49*CLHEP::angstrom, 77 28.49*CLHEP::angstrom,
78 45.35*CLHEP::angstrom, 78 45.35*CLHEP::angstrom,
79 70.03*CLHEP::angstrom, 79 70.03*CLHEP::angstrom,
80 98.05*CLHEP::angstrom, 80 98.05*CLHEP::angstrom,
81 120.56*CLHEP::angstrom, 81 120.56*CLHEP::angstrom,
82 132.73*CLHEP::angstrom, 82 132.73*CLHEP::angstrom,
83 142.60*CLHEP::angstrom, 83 142.60*CLHEP::angstrom,
84 // the above value as given in the paper's t 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*CLH 85 // b=27.22 nm nor the mean value. 129.62*CLHEP::angstrom could be
86 // a better fit. 86 // a better fit.
87 // 87 //
88 // The two last are made up 88 // The two last are made up
89 137.9*CLHEP::angstrom, 89 137.9*CLHEP::angstrom,
90 120.7*CLHEP::angstrom 90 120.7*CLHEP::angstrom
91 }; // angstrom 91 }; // angstrom
92 92
93 //-------------------------------------------- 93 //----------------------------------------------------------------------------
94 94
95 double Meesungnoen2002::GetRmean(double k){ 95 double Meesungnoen2002::GetRmean(double k){
96 G4double k_eV = k/eV; 96 G4double k_eV = k/eV;
97 97
98 if(k_eV>0.1){ // data until 0.2 eV 98 if(k_eV>0.1){ // data until 0.2 eV
99 G4double r_mean = 0; 99 G4double r_mean = 0;
100 for(int8_t i=12; i!=-1 ; --i){ 100 for(int8_t i=12; i!=-1 ; --i){
101 r_mean+=gCoeff[12-i]*std::pow(k_eV,i); 101 r_mean+=gCoeff[12-i]*std::pow(k_eV,i);
102 } 102 }
103 r_mean*=CLHEP::nanometer; 103 r_mean*=CLHEP::nanometer;
104 return r_mean; 104 return r_mean;
105 } 105 }
106 return 0; 106 return 0;
107 } 107 }
108 108
109 double Meesungnoen2002_amorphous::GetRmean(dou 109 double Meesungnoen2002_amorphous::GetRmean(double k){
110 G4double k_eV = k/eV; 110 G4double k_eV = k/eV;
111 111
112 if(k_eV>0.1){ // data until 0.2 eV 112 if(k_eV>0.1){ // data until 0.2 eV
113 G4double r_mean = 0; 113 G4double r_mean = 0;
114 for(int8_t i=6; i!=-1 ; --i){ 114 for(int8_t i=6; i!=-1 ; --i){
115 r_mean+=gCoeff[6-i]*std::pow(k_eV,i); 115 r_mean+=gCoeff[6-i]*std::pow(k_eV,i);
116 } 116 }
117 r_mean*=CLHEP::nanometer; 117 r_mean*=CLHEP::nanometer;
118 return r_mean; 118 return r_mean;
119 } 119 }
120 return 0; 120 return 0;
121 } 121 }
122 122
123 void GetGaussianPenetrationFromRmean3D(G4doubl 123 void GetGaussianPenetrationFromRmean3D(G4double r_mean,
124 G4Three 124 G4ThreeVector& displacement)
125 { 125 {
126 if(r_mean == 0) 126 if(r_mean == 0)
127 { 127 {
128 // rare events: 128 // rare events:
129 // prevent H2O and secondary electron from 129 // prevent H2O and secondary electron from being placed at the same position
130 displacement = G4RandomDirection() * (1e-3 130 displacement = G4RandomDirection() * (1e-3*CLHEP::nanometer);
131 return; 131 return;
132 } 132 }
133 133
134 static constexpr double convertRmean3DToSigm 134 static constexpr double convertRmean3DToSigma1D = 0.62665706865775006;
135 // = sqrt(CLHEP::pi)/pow(2,3./2.) 135 // = sqrt(CLHEP::pi)/pow(2,3./2.)
136 136
137 // Use r_mean to build a 3D gaussian 137 // Use r_mean to build a 3D gaussian
138 const double sigma1D = r_mean * convertRmean 138 const double sigma1D = r_mean * convertRmean3DToSigma1D;
139 displacement = G4ThreeVector(G4RandGauss::sh 139 displacement = G4ThreeVector(G4RandGauss::shoot(0, sigma1D),
140 G4RandGauss::sh 140 G4RandGauss::shoot(0, sigma1D),
141 G4RandGauss::sh 141 G4RandGauss::shoot(0, sigma1D));
142 } 142 }
143 143
144 void Meesungnoen2002::GetPenetration(G4double 144 void Meesungnoen2002::GetPenetration(G4double k,
145 G4ThreeVe 145 G4ThreeVector& displacement)
146 { 146 {
147 GetGaussianPenetrationFromRmean3D(GetRmean(k 147 GetGaussianPenetrationFromRmean3D(GetRmean(k), displacement);
148 } 148 }
149 149
150 void Meesungnoen2002_amorphous::GetPenetration 150 void Meesungnoen2002_amorphous::GetPenetration(G4double k,
151 G4ThreeVe 151 G4ThreeVector& displacement)
152 { 152 {
153 GetGaussianPenetrationFromRmean3D(GetRmean(k 153 GetGaussianPenetrationFromRmean3D(GetRmean(k), displacement);
154 } 154 }
155 155
156 156
157 void Kreipl2009::GetPenetration(G4double k, 157 void Kreipl2009::GetPenetration(G4double k,
158 G4ThreeVe 158 G4ThreeVector& displacement)
159 { 159 {
160 G4double r_mean = Meesungnoen2002::GetRmean( 160 G4double r_mean = Meesungnoen2002::GetRmean(k);
161 161
162 if(r_mean == 0) 162 if(r_mean == 0)
163 { 163 {
164 // rare events: 164 // rare events:
165 // prevent H2O and secondary electron from b 165 // prevent H2O and secondary electron from being placed at the same position
166 displacement = G4RandomDirection() * (1e-3*C 166 displacement = G4RandomDirection() * (1e-3*CLHEP::nanometer);
167 return; 167 return;
168 } 168 }
169 169
170 double r = G4RandGamma::shoot(2,2); 170 double r = G4RandGamma::shoot(2,2);
171 171
172 displacement = G4RandomDirection() * r * r_m 172 displacement = G4RandomDirection() * r * r_mean;
173 } 173 }
174 //-------------------------------------------- 174 //----------------------------------------------------------------------------
175 175
176 void Ritchie1994::GetPenetration(G4double k, 176 void Ritchie1994::GetPenetration(G4double k,
177 G4ThreeVector 177 G4ThreeVector& displacement)
178 { 178 {
179 GetGaussianPenetrationFromRmean3D(k/eV * 1.8 179 GetGaussianPenetrationFromRmean3D(k/eV * 1.8 * nm, // r_mean
180 displaceme 180 displacement);
181 } 181 }
182 182
183 //-------------------------------------------- 183 //----------------------------------------------------------------------------
184 184
185 double Terrisol1990::Get3DStdDeviation(double 185 double Terrisol1990::Get3DStdDeviation(double energy){
186 G4double k_eV = energy/eV; 186 G4double k_eV = energy/eV;
187 if(k_eV < 0.2){ 187 if(k_eV < 0.2){
188 // rare events: 188 // rare events:
189 // prevent H2O and secondary electron to 189 // prevent H2O and secondary electron to be at the spot
190 return 1e-3*CLHEP::nanometer; 190 return 1e-3*CLHEP::nanometer;
191 } 191 }
192 if(k_eV == 9.){ << 192 else if(k_eV == 9.){
193 return gStdDev_T1990[10]; 193 return gStdDev_T1990[10];
194 } 194 }
195 if(k_eV > 9.){ << 195 else if(k_eV > 9.){
196 G4ExceptionDescription description; 196 G4ExceptionDescription description;
197 description << "Terrisol1990 is not tabula 197 description << "Terrisol1990 is not tabulated for energies greater than 9eV";
198 G4Exception("Terrisol1990::Get3DStdDeviati 198 G4Exception("Terrisol1990::Get3DStdDeviation",
199 "INVALID_ARGUMENT", 199 "INVALID_ARGUMENT",
200 FatalErrorInArgument, 200 FatalErrorInArgument,
201 description); 201 description);
202 } 202 }
203 203
204 size_t lowBin, upBin; 204 size_t lowBin, upBin;
205 205
206 if(k_eV >= 1.){ 206 if(k_eV >= 1.){
207 lowBin=std::floor(k_eV)+1; 207 lowBin=std::floor(k_eV)+1;
208 upBin=std::min(lowBin+1, size_t(10)); 208 upBin=std::min(lowBin+1, size_t(10));
209 } 209 }
210 else{ 210 else{
211 auto it=std::lower_bound(&gEnergies_T1990[ 211 auto it=std::lower_bound(&gEnergies_T1990[0],
212 &gEnergies_T1990[ 212 &gEnergies_T1990[2],
213 k_eV); 213 k_eV);
214 lowBin = it-&gEnergies_T1990[0]; 214 lowBin = it-&gEnergies_T1990[0];
215 upBin = lowBin+1; 215 upBin = lowBin+1;
216 } 216 }
217 217
218 double lowE = gEnergies_T1990[lowBin]; 218 double lowE = gEnergies_T1990[lowBin];
219 double upE = gEnergies_T1990[upBin]; 219 double upE = gEnergies_T1990[upBin];
220 220
221 double lowS = gStdDev_T1990[lowBin]; 221 double lowS = gStdDev_T1990[lowBin];
222 double upS = gStdDev_T1990[upBin]; 222 double upS = gStdDev_T1990[upBin];
223 223
224 double tanA = (lowS-upS)/(lowE-upE); 224 double tanA = (lowS-upS)/(lowE-upE);
225 double sigma3D = lowS + (k_eV-lowE)*tanA; 225 double sigma3D = lowS + (k_eV-lowE)*tanA;
226 return sigma3D; 226 return sigma3D;
227 } 227 }
228 228
229 double Terrisol1990::GetRmean(double energy){ 229 double Terrisol1990::GetRmean(double energy){
230 double sigma3D=Get3DStdDeviation(energy); 230 double sigma3D=Get3DStdDeviation(energy);
231 231
232 static constexpr double s2r=1.59576912160573 232 static constexpr double s2r=1.595769121605731; // = pow(2,3./2.)/sqrt(CLHEP::pi)
233 233
234 double r_mean=sigma3D*s2r; 234 double r_mean=sigma3D*s2r;
235 return r_mean; 235 return r_mean;
236 } 236 }
237 237
238 void Terrisol1990::GetPenetration(G4double ene 238 void Terrisol1990::GetPenetration(G4double energy,
239 G4ThreeVecto 239 G4ThreeVector& displacement){
240 double sigma3D = Get3DStdDeviation(energy); 240 double sigma3D = Get3DStdDeviation(energy);
241 241
242 static constexpr double factor = 2.204969995 242 static constexpr double factor = 2.20496999539; // = 1./(3. - 8./CLHEP::pi);
243 243
244 double sigma1D = std::sqrt(std::pow(sigma3D, 244 double sigma1D = std::sqrt(std::pow(sigma3D, 2.)*factor);
245 245
246 displacement = G4ThreeVector(G4RandGauss::sh 246 displacement = G4ThreeVector(G4RandGauss::shoot(0, sigma1D),
247 G4RandGauss::sh 247 G4RandGauss::shoot(0, sigma1D),
248 G4RandGauss::sh 248 G4RandGauss::shoot(0, sigma1D));
249 } 249 }
250 250
251 } // Penetration 251 } // Penetration
252 } // DNA 252 } // DNA
253 253
254 //-------------------------------------------- 254 //------------------------------------------------------------------------------
255 255
256 G4VEmModel* G4DNASolvationModelFactory::Create 256 G4VEmModel* G4DNASolvationModelFactory::Create(const G4String& penetrationModel)
257 { 257 {
258 G4String modelNamePrefix("DNAOneStepThermali 258 G4String modelNamePrefix("DNAOneStepThermalizationModel_");
259 259
260 if(penetrationModel == "Terrisol1990") 260 if(penetrationModel == "Terrisol1990")
261 { 261 {
262 return new G4TDNAOneStepThermalizationMode 262 return new G4TDNAOneStepThermalizationModel<DNA::Penetration::Terrisol1990>(G4Electron::Definition(), modelNamePrefix + penetrationModel);
263 } 263 }
264 if(penetrationModel == "Meesungnoen2002") << 264 else if(penetrationModel == "Meesungnoen2002")
265 { 265 {
266 return new G4TDNAOneStepThermalizationMode 266 return new G4TDNAOneStepThermalizationModel<DNA::Penetration::Meesungnoen2002>(G4Electron::Definition(), modelNamePrefix + penetrationModel);
267 } 267 }
268 if(penetrationModel == "Meesungnoen2002_amor << 268 else if(penetrationModel == "Meesungnoen2002_amorphous")
269 { 269 {
270 return new G4TDNAOneStepThermalizationModel< 270 return new G4TDNAOneStepThermalizationModel<DNA::Penetration::Meesungnoen2002_amorphous>(G4Electron::Definition(), modelNamePrefix + penetrationModel);
271 } 271 }
272 if(penetrationModel == "Kreipl2009") << 272 else if(penetrationModel == "Kreipl2009")
273 { 273 {
274 return new G4TDNAOneStepThermalizationModel< 274 return new G4TDNAOneStepThermalizationModel<DNA::Penetration::Kreipl2009>(G4Electron::Definition(), modelNamePrefix + penetrationModel);
275 } 275 }
276 if(penetrationModel == "Ritchie1994") << 276 else if(penetrationModel == "Ritchie1994")
277 { 277 {
278 return new G4TDNAOneStepThermalizationMode 278 return new G4TDNAOneStepThermalizationModel<DNA::Penetration::Ritchie1994>(G4Electron::Definition(), modelNamePrefix + penetrationModel);
279 } 279 }
280 << 280 else
281 G4ExceptionDescription description; << 281 {
282 description << penetrationModel + " is not a << 282 G4ExceptionDescription description;
283 G4Exception("G4DNASolvationModelFactory::Cre << 283 description << penetrationModel + " is not a valid model name.";
284 "INVALID_ARGUMENT", << 284 G4Exception("G4DNASolvationModelFactory::Create",
285 FatalErrorInArgument, << 285 "INVALID_ARGUMENT",
286 description, << 286 FatalErrorInArgument,
287 "Options are: Terrisol1990, Mees << 287 description,
>> 288 "Options are: Terrisol1990, Meesungnoen2002, Ritchie1994.");
>> 289 }
288 return nullptr; 290 return nullptr;
289 } 291 }
290 292
291 //-------------------------------------------- 293 //------------------------------------------------------------------------------
292 G4VEmModel* G4DNASolvationModelFactory::GetMac 294 G4VEmModel* G4DNASolvationModelFactory::GetMacroDefinedModel()
293 { 295 {
294 auto dnaSubType = G4EmParameters::Instance() 296 auto dnaSubType = G4EmParameters::Instance()->DNAeSolvationSubType();
295 297
296 switch(dnaSubType) 298 switch(dnaSubType)
297 { 299 {
298 case fRitchie1994eSolvation: 300 case fRitchie1994eSolvation:
299 return Create("Ritchie1994"); 301 return Create("Ritchie1994");
300 case fTerrisol1990eSolvation: 302 case fTerrisol1990eSolvation:
301 return Create("Terrisol1990"); 303 return Create("Terrisol1990");
302 case fKreipl2009eSolvation: 304 case fKreipl2009eSolvation:
303 return Create("Kreipl2009"); 305 return Create("Kreipl2009");
304 case fMeesungnoensolid2002eSolvation: 306 case fMeesungnoensolid2002eSolvation:
305 return Create("Meesungnoen2002_amorphous"); 307 return Create("Meesungnoen2002_amorphous");
306 case fMeesungnoen2002eSolvation: 308 case fMeesungnoen2002eSolvation:
307 case fDNAUnknownModel: 309 case fDNAUnknownModel:
308 return Create("Meesungnoen2002"); 310 return Create("Meesungnoen2002");
309 default: 311 default:
310 G4Exception("G4DNASolvationModelFactory::G 312 G4Exception("G4DNASolvationModelFactory::GetMacroDefinedModel",
311 "DnaSubType", 313 "DnaSubType",
312 FatalErrorInArgument, 314 FatalErrorInArgument,
313 "The solvation parameter store 315 "The solvation parameter stored in G4EmParameters is unknown. Supported types are: fRitchie1994eSolvation, fTerrisol1990eSolvation, fMeesungnoen2002eSolvation.");
314 } 316 }
315 317
316 return nullptr; 318 return nullptr;
317 } 319 }
318 320