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 // Created on 2016/04/08 27 // 28 // Authors: D. Sakata, S. Incerti 29 // 30 // This class perform transmission term of vol 31 // based on Quinn Model, see Phys. Rev. vol 12 32 33 #include "G4DNAQuinnPlasmonExcitationModel.hh" 34 #include "G4SystemOfUnits.hh" 35 #include "G4RandomDirection.hh" 36 37 //....oooOO0OOooo........oooOO0OOooo........oo 38 39 using namespace std; 40 41 //....oooOO0OOooo........oooOO0OOooo........oo 42 43 G4DNAQuinnPlasmonExcitationModel::G4DNAQuinnPl 44 (const G4ParticleDefin 45 const G4String& nam): 46 G4VEmModel(nam) 47 { 48 fpMaterialDensity = nullptr; 49 fLowEnergyLimit = 10 * eV; 50 fHighEnergyLimit = 1.0 * GeV; 51 52 for(int & i : nValenceElectron) i=0; 53 54 verboseLevel = 0; 55 56 if (verboseLevel > 0) 57 { 58 G4cout << "Quinn plasmon excitation model 59 } 60 fParticleChangeForGamma = nullptr; 61 statCode = false; 62 } 63 64 //....oooOO0OOooo........oooOO0OOooo........oo 65 66 G4DNAQuinnPlasmonExcitationModel::~G4DNAQuinnP 67 = default; 68 69 //....oooOO0OOooo........oooOO0OOooo........oo 70 71 void G4DNAQuinnPlasmonExcitationModel::Initial 72 (const G4ParticleDefin 73 const G4DataVector& / 74 { 75 for(int & i : nValenceElectron) i=0; 76 77 if (verboseLevel > 3) 78 { 79 G4cout << 80 "Calling G4DNAQuinnPlasmonExcitatio 81 << G4endl; 82 } 83 84 if(particle == G4Electron::ElectronDefinitio 85 { 86 fLowEnergyLimit = 10 * eV; 87 fHighEnergyLimit = 1.0 * GeV; 88 } 89 else 90 { 91 G4Exception("G4DNAQuinnPlasmonExcitationMo 92 FatalException,"Not defined for other part 93 return; 94 } 95 96 // Get Number of valence electrons 97 G4ProductionCutsTable* theCoupleTable = 98 G4ProductionCutsTable::GetProductionCuts 99 100 auto numOfCouples = (G4int)theCoupleTable-> 101 102 for(G4int i=0;i<numOfCouples;i++){ 103 104 const G4MaterialCutsCouple* couple = 105 theCoupleTable->GetMaterialCutsCoupl 106 107 const G4Material* material = couple->GetMa 108 109 const G4ElementVector* theElementVector =m 110 111 std::size_t nelm = material->GetNumberOfEl 112 if (nelm==1) // Protection: only for singl 113 { 114 G4int z = G4lrint((*theElementVector)[0] 115 if(z<=100) 116 { 117 nValenceElectron[z] = GetNValenceElec 118 } 119 else 120 { 121 G4Exception("G4DNAQuinnPlasmonExcitati 122 FatalException,"The model is not app 123 } 124 } 125 //for(G4int j=0;j<nelm;j++){ 126 // G4int z=G4lrint((*theElementVector)[ 127 // if(z<=100){nValenceElectron[z] = Ge 128 //} 129 } 130 131 if( verboseLevel>0 ) 132 { 133 G4cout << "Quinn plasmon excitation model 134 << "Energy range: " 135 << LowEnergyLimit() / eV << " eV - " 136 << HighEnergyLimit() / keV << " keV for " 137 << particle->GetParticleName() 138 << G4endl; 139 } 140 141 if (isInitialised){return;} 142 fParticleChangeForGamma = GetParticleChangeF 143 isInitialised = true; 144 } 145 146 //....oooOO0OOooo........oooOO0OOooo........oo 147 148 G4double G4DNAQuinnPlasmonExcitationModel::Cro 149 (const G4Material* ma 150 const G4ParticleDefi 151 G4double ekin, 152 G4double, 153 G4double) 154 { 155 if (verboseLevel > 3) 156 { 157 G4cout << 158 "Calling CrossSectionPerVolume() of G 159 << G4endl; 160 } 161 162 // Protection: only for single element 163 if(material->GetNumberOfElements()>1) return 164 G4double z = material->GetZ(); 165 166 // Protection: only for Gold 167 if (z!=79){return 0.;} 168 169 170 G4double sigma = 0; 171 G4double atomicNDensity = material->GetAtomi 172 173 if(atomicNDensity!= 0.0) 174 { 175 if (ekin >= fLowEnergyLimit && ekin < fHig 176 { 177 sigma = GetCrossSection(material,particl 178 } 179 180 if (verboseLevel > 2) 181 { 182 G4cout<<"_______________________________ 183 G4cout<<"=== G4DNAQuinnPlasmonExcitation 184 G4cout<<"=== Kinetic energy (eV)=" << ek 185 <<particleDefinition->GetParticleN 186 G4cout<<"=== Cross section per atom for 187 <<sigma/cm/cm << G4endl; 188 G4cout<<"=== Cross section per atom for 189 <<sigma*atomicNDensity/(1./cm) << 190 G4cout<<"=== G4DNAQuinnPlasmonExcitation 191 } 192 } 193 194 return sigma*atomicNDensity; 195 } 196 197 //....oooOO0OOooo........oooOO0OOooo........oo 198 199 void G4DNAQuinnPlasmonExcitationModel::SampleS 200 (std::vector<G4Dynami 201 const G4MaterialCuts 202 const G4DynamicParti 203 G4double,G4double) 204 { 205 206 if (verboseLevel > 3) 207 { 208 G4cout << 209 "Calling SampleSecondaries() of G4D 210 << G4endl; 211 } 212 213 const G4Material *material = couple->GetMate 214 215 G4ParticleDefinition* particle = aDynamicPar 216 217 G4double k = aDynamicPar 218 219 if(particle == G4Electron::ElectronDefinitio 220 { 221 G4double e = 1.; 222 G4int z = material->GetZ(); 223 G4int Nve = 0; 224 225 //TODO: have to be change to realistic!! 226 if(z<100) Nve = nValenceElectron[z]; 227 228 G4double A = material->GetA()/g/mole 229 G4double Dens = material->GetDensity()/ 230 G4double veDens = Dens*CLHEP::Avogadro*Nv 231 232 G4double omega_p = std::sqrt(veDens*std::p 233 (CLHEP::epsilon0/(1./cm) 234 /(CLHEP::c_squared/cm/cm 235 236 G4double excitationEnergy = CLHEP::hbar_ 237 G4double newEnergy = k - excitati 238 239 240 if (newEnergy > 0) 241 { 242 fParticleChangeForGamma-> 243 ProposeMomentumDirection(aDynamicP 244 245 fParticleChangeForGamma->ProposeLocalEne 246 247 if(!statCode) 248 { 249 fParticleChangeForGamma->SetPropose 250 } 251 else 252 { 253 fParticleChangeForGamma->SetPropose 254 255 } 256 } 257 } 258 } 259 260 //....oooOO0OOooo........oooOO0OOooo........oo 261 262 G4double G4DNAQuinnPlasmonExcitationModel::Get 263 (const G4Material* ma 264 const G4ParticleDefi 265 G4double kineticEner 266 { 267 G4double value=0; 268 269 if(particle == G4Electron::ElectronDefinitio 270 { 271 G4double e = 1.; 272 G4int z = material->GetZ(); 273 G4int Nve = 0; 274 if(z<100) Nve = nValenceElectron[z]; 275 G4double A = material->GetA()/g/mol 276 G4double Dens = material->GetDensity() 277 G4double veDens = Dens*CLHEP::Avogadro*N 278 279 G4double omega_p = std::sqrt(veDens*std:: 280 /(CLHEP::epsilon0/(1./ 281 (CLHEP::c_squared/cm/c 282 283 G4double fEnergy = std::pow(CLHEP::h_Plan 284 std::pow(3*veDens/CLHE 285 *(CLHEP::c_squared/cm/ 286 287 G4double p0 = sqrt(2*CLHEP::electron 288 /(CLHEP::c_squared/cm/ 289 290 G4double p = sqrt(2*CLHEP::electron 291 /(CLHEP::c_squared/cm/ 292 293 G4double mfp = 2*CLHEP::Bohr_radius/c 294 /(CLHEP::hbar_Planck*o 295 (G4Log((std::pow(std:: 296 +2*CLHEP::electron_mas 297 (CLHEP::c_squared/cm/c 298 *CLHEP::hbar_Planck,1. 299 /(p-std::pow(std::pow( 300 (CLHEP::c_squared/cm/c 301 *CLHEP::hbar_Planck,1. 302 303 G4double excitationEnergy = CLHEP::hbar 304 305 if((0<mfp)&&(0<veDens)&&(excitationEnergy 306 value = 1./(veDens*mfp); 307 } 308 } 309 return value*cm*cm; 310 } 311 312 //....oooOO0OOooo........oooOO0OOooo........oo 313 314 G4int G4DNAQuinnPlasmonExcitationModel::GetNVa 315 { 316 317 G4int Nve=0; 318 319 // Current limitation to gold 320 if (z!=79){return 0.;} 321 322 if (verboseLevel > 3) 323 { 324 G4cout << 325 "Calling GetNValenceElectron() of G 326 << G4endl; 327 } 328 329 const char *datadir=nullptr; 330 331 if(datadir == nullptr) 332 { 333 datadir = G4FindDataDir("G4LEDATA"); 334 if(datadir == nullptr) 335 { 336 G4Exception("G4DNAQuinnPlasmonExcitatio 337 ,"em0002",FatalException, 338 "Enviroment variable G4LEDA 339 return 0; 340 } 341 } 342 343 std::ostringstream targetfile; 344 targetfile.str(""); 345 targetfile.clear(stringstream::goodbit); 346 targetfile << datadir <<"/dna/atomicstate_Z" 347 std::ifstream fin(targetfile.str().c_str()); 348 349 if(!fin) 350 { 351 G4cout<< " Error : "<< targetfile.str() << 352 G4Exception("G4DNAQuinnPlasmonExcitationMo 353 ,"em0003",FatalException, 354 "There is no target file"); 355 return 0; 356 } 357 358 string buff0,buff1,buff2,buff3,buff4,buff5,b 359 fin >> buff0 >>buff1>>buff2>>buff3>>buff4>>b 360 361 while(true){ 362 fin >> buff0 >>buff1>>buff2>>buff3>>buff4> 363 if(!fin.eof()) 364 { 365 Nve = stoi(buff3); 366 } 367 else 368 { 369 break; 370 } 371 } 372 return Nve; 373 } 374 375