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