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 27 28 #include "G4DNAMillerGreenExcitationModel.hh" 28 #include "G4DNAMillerGreenExcitationModel.hh" 29 #include "G4SystemOfUnits.hh" 29 #include "G4SystemOfUnits.hh" 30 #include "G4DNAChemistryManager.hh" 30 #include "G4DNAChemistryManager.hh" 31 #include "G4DNAMolecularMaterial.hh" 31 #include "G4DNAMolecularMaterial.hh" 32 #include "G4Exp.hh" 32 #include "G4Exp.hh" 33 #include "G4Pow.hh" << 34 #include "G4Alpha.hh" << 35 33 36 static G4Pow * gpow = G4Pow::GetInstance(); << 37 //....oooOO0OOooo........oooOO0OOooo........oo 34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 38 35 39 using namespace std; 36 using namespace std; 40 37 41 //....oooOO0OOooo........oooOO0OOooo........oo 38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 42 39 43 G4DNAMillerGreenExcitationModel::G4DNAMillerGr 40 G4DNAMillerGreenExcitationModel::G4DNAMillerGreenExcitationModel(const G4ParticleDefinition*, 44 41 const G4String& nam) 45 :G4VEmModel(nam) << 42 :G4VEmModel(nam),isInitialised(false) 46 { 43 { 47 fpMolWaterDensity = nullptr; << 44 fpMolWaterDensity = 0; 48 45 49 nLevels=0; 46 nLevels=0; 50 kineticEnergyCorrection[0]=0.; 47 kineticEnergyCorrection[0]=0.; 51 kineticEnergyCorrection[1]=0.; 48 kineticEnergyCorrection[1]=0.; 52 kineticEnergyCorrection[2]=0.; 49 kineticEnergyCorrection[2]=0.; 53 kineticEnergyCorrection[3]=0.; 50 kineticEnergyCorrection[3]=0.; 54 51 55 verboseLevel= 0; 52 verboseLevel= 0; 56 // Verbosity scale: 53 // Verbosity scale: 57 // 0 = nothing 54 // 0 = nothing 58 // 1 = warning for energy non-conservation 55 // 1 = warning for energy non-conservation 59 // 2 = details of energy budget 56 // 2 = details of energy budget 60 // 3 = calculation of cross sections, file o 57 // 3 = calculation of cross sections, file openings, sampling of atoms 61 // 4 = entering in methods 58 // 4 = entering in methods 62 59 63 if( verboseLevel>0 ) 60 if( verboseLevel>0 ) 64 { 61 { 65 G4cout << "Miller & Green excitation model 62 G4cout << "Miller & Green excitation model is constructed " << G4endl; 66 } 63 } 67 fParticleChangeForGamma = nullptr; << 64 fParticleChangeForGamma = 0; 68 65 69 // Selection of stationary mode 66 // Selection of stationary mode 70 67 71 statCode = false; 68 statCode = false; 72 } 69 } 73 70 74 //....oooOO0OOooo........oooOO0OOooo........oo 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 75 72 76 G4DNAMillerGreenExcitationModel::~G4DNAMillerG 73 G4DNAMillerGreenExcitationModel::~G4DNAMillerGreenExcitationModel() 77 = default; << 74 {} 78 75 79 //....oooOO0OOooo........oooOO0OOooo........oo 76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 80 77 81 void G4DNAMillerGreenExcitationModel::Initiali 78 void G4DNAMillerGreenExcitationModel::Initialise(const G4ParticleDefinition* particle, 82 79 const G4DataVector& /*cuts*/) 83 { 80 { 84 81 85 if (verboseLevel > 3) 82 if (verboseLevel > 3) 86 G4cout << "Calling G4DNAMillerGreenExcitat 83 G4cout << "Calling G4DNAMillerGreenExcitationModel::Initialise()" << G4endl; 87 84 88 // Energy limits 85 // Energy limits 89 86 90 G4DNAGenericIonsManager *instance; 87 G4DNAGenericIonsManager *instance; 91 instance = G4DNAGenericIonsManager::Instance 88 instance = G4DNAGenericIonsManager::Instance(); 92 protonDef = G4Proton::ProtonDefinition(); << 89 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition(); 93 hydrogenDef = instance->GetIon("hydrogen"); << 90 G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen"); 94 alphaPlusPlusDef = G4Alpha::Alpha(); << 91 G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++"); 95 alphaPlusDef = instance->GetIon("alpha+"); << 92 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+"); 96 heliumDef = instance->GetIon("helium"); << 93 G4ParticleDefinition* heliumDef = instance->GetIon("helium"); 97 94 98 G4String proton; 95 G4String proton; 99 G4String hydrogen; 96 G4String hydrogen; 100 G4String alphaPlusPlus; 97 G4String alphaPlusPlus; 101 G4String alphaPlus; 98 G4String alphaPlus; 102 G4String helium; 99 G4String helium; 103 100 104 // LIMITS AND CONSTANTS 101 // LIMITS AND CONSTANTS 105 102 106 proton = protonDef->GetParticleName(); 103 proton = protonDef->GetParticleName(); 107 lowEnergyLimit[proton] = 10. * eV; 104 lowEnergyLimit[proton] = 10. * eV; 108 highEnergyLimit[proton] = 500. * keV; 105 highEnergyLimit[proton] = 500. * keV; 109 106 110 kineticEnergyCorrection[0] = 1.; 107 kineticEnergyCorrection[0] = 1.; 111 slaterEffectiveCharge[0][0] = 0.; 108 slaterEffectiveCharge[0][0] = 0.; 112 slaterEffectiveCharge[1][0] = 0.; 109 slaterEffectiveCharge[1][0] = 0.; 113 slaterEffectiveCharge[2][0] = 0.; 110 slaterEffectiveCharge[2][0] = 0.; 114 sCoefficient[0][0] = 0.; 111 sCoefficient[0][0] = 0.; 115 sCoefficient[1][0] = 0.; 112 sCoefficient[1][0] = 0.; 116 sCoefficient[2][0] = 0.; 113 sCoefficient[2][0] = 0.; 117 114 118 hydrogen = hydrogenDef->GetParticleName(); 115 hydrogen = hydrogenDef->GetParticleName(); 119 lowEnergyLimit[hydrogen] = 10. * eV; 116 lowEnergyLimit[hydrogen] = 10. * eV; 120 highEnergyLimit[hydrogen] = 500. * keV; 117 highEnergyLimit[hydrogen] = 500. * keV; 121 118 122 kineticEnergyCorrection[0] = 1.; 119 kineticEnergyCorrection[0] = 1.; 123 slaterEffectiveCharge[0][0] = 0.; 120 slaterEffectiveCharge[0][0] = 0.; 124 slaterEffectiveCharge[1][0] = 0.; 121 slaterEffectiveCharge[1][0] = 0.; 125 slaterEffectiveCharge[2][0] = 0.; 122 slaterEffectiveCharge[2][0] = 0.; 126 sCoefficient[0][0] = 0.; 123 sCoefficient[0][0] = 0.; 127 sCoefficient[1][0] = 0.; 124 sCoefficient[1][0] = 0.; 128 sCoefficient[2][0] = 0.; 125 sCoefficient[2][0] = 0.; 129 126 130 alphaPlusPlus = alphaPlusPlusDef->GetParticl 127 alphaPlusPlus = alphaPlusPlusDef->GetParticleName(); 131 lowEnergyLimit[alphaPlusPlus] = 1. * keV; 128 lowEnergyLimit[alphaPlusPlus] = 1. * keV; 132 highEnergyLimit[alphaPlusPlus] = 400. * MeV; 129 highEnergyLimit[alphaPlusPlus] = 400. * MeV; 133 130 134 kineticEnergyCorrection[1] = 0.9382723/3.727 131 kineticEnergyCorrection[1] = 0.9382723/3.727417; 135 slaterEffectiveCharge[0][1]=0.; 132 slaterEffectiveCharge[0][1]=0.; 136 slaterEffectiveCharge[1][1]=0.; 133 slaterEffectiveCharge[1][1]=0.; 137 slaterEffectiveCharge[2][1]=0.; 134 slaterEffectiveCharge[2][1]=0.; 138 sCoefficient[0][1]=0.; 135 sCoefficient[0][1]=0.; 139 sCoefficient[1][1]=0.; 136 sCoefficient[1][1]=0.; 140 sCoefficient[2][1]=0.; 137 sCoefficient[2][1]=0.; 141 138 142 alphaPlus = alphaPlusDef->GetParticleName(); 139 alphaPlus = alphaPlusDef->GetParticleName(); 143 lowEnergyLimit[alphaPlus] = 1. * keV; 140 lowEnergyLimit[alphaPlus] = 1. * keV; 144 highEnergyLimit[alphaPlus] = 400. * MeV; 141 highEnergyLimit[alphaPlus] = 400. * MeV; 145 142 146 kineticEnergyCorrection[2] = 0.9382723/3.727 143 kineticEnergyCorrection[2] = 0.9382723/3.727417; 147 slaterEffectiveCharge[0][2]=2.0; 144 slaterEffectiveCharge[0][2]=2.0; 148 145 149 // Following values provided by M. Dingfelde 146 // Following values provided by M. Dingfelder 150 slaterEffectiveCharge[1][2]=2.00; 147 slaterEffectiveCharge[1][2]=2.00; 151 slaterEffectiveCharge[2][2]=2.00; 148 slaterEffectiveCharge[2][2]=2.00; 152 // 149 // 153 sCoefficient[0][2]=0.7; 150 sCoefficient[0][2]=0.7; 154 sCoefficient[1][2]=0.15; 151 sCoefficient[1][2]=0.15; 155 sCoefficient[2][2]=0.15; 152 sCoefficient[2][2]=0.15; 156 153 157 helium = heliumDef->GetParticleName(); 154 helium = heliumDef->GetParticleName(); 158 lowEnergyLimit[helium] = 1. * keV; 155 lowEnergyLimit[helium] = 1. * keV; 159 highEnergyLimit[helium] = 400. * MeV; 156 highEnergyLimit[helium] = 400. * MeV; 160 157 161 kineticEnergyCorrection[3] = 0.9382723/3.727 158 kineticEnergyCorrection[3] = 0.9382723/3.727417; 162 slaterEffectiveCharge[0][3]=1.7; 159 slaterEffectiveCharge[0][3]=1.7; 163 slaterEffectiveCharge[1][3]=1.15; 160 slaterEffectiveCharge[1][3]=1.15; 164 slaterEffectiveCharge[2][3]=1.15; 161 slaterEffectiveCharge[2][3]=1.15; 165 sCoefficient[0][3]=0.5; 162 sCoefficient[0][3]=0.5; 166 sCoefficient[1][3]=0.25; 163 sCoefficient[1][3]=0.25; 167 sCoefficient[2][3]=0.25; 164 sCoefficient[2][3]=0.25; 168 165 169 // 166 // 170 167 171 if (particle==protonDef) 168 if (particle==protonDef) 172 { 169 { 173 SetLowEnergyLimit(lowEnergyLimit[proton]); 170 SetLowEnergyLimit(lowEnergyLimit[proton]); 174 SetHighEnergyLimit(highEnergyLimit[proton] 171 SetHighEnergyLimit(highEnergyLimit[proton]); 175 } 172 } 176 173 177 if (particle==hydrogenDef) 174 if (particle==hydrogenDef) 178 { 175 { 179 SetLowEnergyLimit(lowEnergyLimit[hydrogen] 176 SetLowEnergyLimit(lowEnergyLimit[hydrogen]); 180 SetHighEnergyLimit(highEnergyLimit[hydroge 177 SetHighEnergyLimit(highEnergyLimit[hydrogen]); 181 } 178 } 182 179 183 if (particle==alphaPlusPlusDef) 180 if (particle==alphaPlusPlusDef) 184 { 181 { 185 SetLowEnergyLimit(lowEnergyLimit[alphaPlus 182 SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]); 186 SetHighEnergyLimit(highEnergyLimit[alphaPl 183 SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]); 187 } 184 } 188 185 189 if (particle==alphaPlusDef) 186 if (particle==alphaPlusDef) 190 { 187 { 191 SetLowEnergyLimit(lowEnergyLimit[alphaPlus 188 SetLowEnergyLimit(lowEnergyLimit[alphaPlus]); 192 SetHighEnergyLimit(highEnergyLimit[alphaPl 189 SetHighEnergyLimit(highEnergyLimit[alphaPlus]); 193 } 190 } 194 191 195 if (particle==heliumDef) 192 if (particle==heliumDef) 196 { 193 { 197 SetLowEnergyLimit(lowEnergyLimit[helium]); 194 SetLowEnergyLimit(lowEnergyLimit[helium]); 198 SetHighEnergyLimit(highEnergyLimit[helium] 195 SetHighEnergyLimit(highEnergyLimit[helium]); 199 } 196 } 200 197 201 // 198 // 202 199 203 nLevels = waterExcitation.NumberOfLevels(); 200 nLevels = waterExcitation.NumberOfLevels(); 204 201 205 // 202 // 206 if( verboseLevel>0 ) 203 if( verboseLevel>0 ) 207 { 204 { 208 G4cout << "Miller & Green excitation model 205 G4cout << "Miller & Green excitation model is initialized " << G4endl 209 << "Energy range: " 206 << "Energy range: " 210 << LowEnergyLimit() / eV << " eV - 207 << LowEnergyLimit() / eV << " eV - " 211 << HighEnergyLimit() / keV << " keV 208 << HighEnergyLimit() / keV << " keV for " 212 << particle->GetParticleName() 209 << particle->GetParticleName() 213 << G4endl; 210 << G4endl; 214 } 211 } 215 212 216 // Initialize water density pointer 213 // Initialize water density pointer 217 fpMolWaterDensity = G4DNAMolecularMaterial:: 214 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER")); 218 215 219 if (isInitialised) { return; } 216 if (isInitialised) { return; } 220 fParticleChangeForGamma = GetParticleChangeF 217 fParticleChangeForGamma = GetParticleChangeForGamma(); 221 isInitialised = true; 218 isInitialised = true; 222 219 223 } 220 } 224 221 225 //....oooOO0OOooo........oooOO0OOooo........oo 222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 226 223 227 G4double G4DNAMillerGreenExcitationModel::Cros 224 G4double G4DNAMillerGreenExcitationModel::CrossSectionPerVolume(const G4Material* material, 228 225 const G4ParticleDefinition* particleDefinition, 229 226 G4double k, 230 227 G4double, 231 228 G4double) 232 { 229 { 233 if (verboseLevel > 3) 230 if (verboseLevel > 3) 234 G4cout << "Calling CrossSectionPerVolume() 231 G4cout << "Calling CrossSectionPerVolume() of G4DNAMillerGreenExcitationModel" << G4endl; 235 232 236 // Calculate total cross section for model 233 // Calculate total cross section for model 237 234 >> 235 G4DNAGenericIonsManager *instance; >> 236 instance = G4DNAGenericIonsManager::Instance(); >> 237 238 if ( 238 if ( 239 particleDefinition != protonDef << 239 particleDefinition != G4Proton::ProtonDefinition() 240 && 240 && 241 particleDefinition != hydrogenDef << 241 particleDefinition != instance->GetIon("hydrogen") 242 && 242 && 243 particleDefinition != alphaPlusPlusDef << 243 particleDefinition != instance->GetIon("alpha++") 244 && 244 && 245 particleDefinition != alphaPlusDef << 245 particleDefinition != instance->GetIon("alpha+") 246 && 246 && 247 particleDefinition != heliumDef << 247 particleDefinition != instance->GetIon("helium") 248 ) 248 ) 249 249 250 return 0; 250 return 0; 251 251 252 G4double lowLim = 0; 252 G4double lowLim = 0; 253 G4double highLim = 0; 253 G4double highLim = 0; 254 G4double crossSection = 0.; 254 G4double crossSection = 0.; 255 255 256 G4double waterDensity = (*fpMolWaterDensity) 256 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()]; 257 257 258 const G4String& particleName = particleDefin 258 const G4String& particleName = particleDefinition->GetParticleName(); 259 259 260 std::map< G4String,G4double,std::less<G4Stri 260 std::map< G4String,G4double,std::less<G4String> >::iterator pos1; 261 pos1 = lowEnergyLimit.find(particleName); 261 pos1 = lowEnergyLimit.find(particleName); 262 262 263 if (pos1 != lowEnergyLimit.end()) 263 if (pos1 != lowEnergyLimit.end()) 264 { 264 { 265 lowLim = pos1->second; 265 lowLim = pos1->second; 266 } 266 } 267 267 268 std::map< G4String,G4double,std::less<G4Stri 268 std::map< G4String,G4double,std::less<G4String> >::iterator pos2; 269 pos2 = highEnergyLimit.find(particleName); 269 pos2 = highEnergyLimit.find(particleName); 270 270 271 if (pos2 != highEnergyLimit.end()) 271 if (pos2 != highEnergyLimit.end()) 272 { 272 { 273 highLim = pos2->second; 273 highLim = pos2->second; 274 } 274 } 275 275 276 if (k >= lowLim && k <= highLim) 276 if (k >= lowLim && k <= highLim) 277 { 277 { 278 crossSection = Sum(k,particleDefinition); 278 crossSection = Sum(k,particleDefinition); 279 279 280 // add ONE or TWO electron-water excitatio 280 // add ONE or TWO electron-water excitation for alpha+ and helium 281 /* 281 /* 282 if ( particleDefinition == alphaPlusDef << 282 if ( particleDefinition == instance->GetIon("alpha+") 283 || 283 || 284 particleDefinition == heliumDef << 284 particleDefinition == instance->GetIon("helium") 285 ) 285 ) 286 { 286 { 287 287 288 G4DNAEmfietzoglouExcitationModel * excit 288 G4DNAEmfietzoglouExcitationModel * excitationXS = new G4DNAEmfietzoglouExcitationModel(); 289 excitationXS->Initialise(G4Electron: 289 excitationXS->Initialise(G4Electron::ElectronDefinition()); 290 290 291 G4double sigmaExcitation=0; 291 G4double sigmaExcitation=0; 292 G4double tmp =0.; 292 G4double tmp =0.; 293 293 294 if (k*0.511/3728 > 8.23*eV && k*0.511/37 294 if (k*0.511/3728 > 8.23*eV && k*0.511/3728 < 10*MeV ) sigmaExcitation = 295 excitationXS->CrossSectionPerVolume(ma 295 excitationXS->CrossSectionPerVolume(material,G4Electron::ElectronDefinition(),k*0.511/3728,tmp,tmp) 296 /material->GetAtomicNumDensityVector() 296 /material->GetAtomicNumDensityVector()[1]; 297 297 298 if ( particleDefinition == alphaPlusDef << 298 if ( particleDefinition == instance->GetIon("alpha+") ) 299 crossSection = crossSection + sigmaEx 299 crossSection = crossSection + sigmaExcitation ; 300 300 301 if ( particleDefinition == heliumDef ) << 301 if ( particleDefinition == instance->GetIon("helium") ) 302 crossSection = crossSection + 2*sigmaE 302 crossSection = crossSection + 2*sigmaExcitation ; 303 303 304 delete excitationXS; 304 delete excitationXS; 305 305 306 // Alternative excitation model 306 // Alternative excitation model 307 307 308 G4DNABornExcitationModel * excitatio 308 G4DNABornExcitationModel * excitationXS = new G4DNABornExcitationModel(); 309 excitationXS->Initialise(G4Electron: 309 excitationXS->Initialise(G4Electron::ElectronDefinition()); 310 310 311 G4double sigmaExcitation=0; 311 G4double sigmaExcitation=0; 312 G4double tmp=0; 312 G4double tmp=0; 313 313 314 if (k*0.511/3728 > 9*eV && k*0.511/3728 314 if (k*0.511/3728 > 9*eV && k*0.511/3728 < 1*MeV ) sigmaExcitation = 315 excitationXS->CrossSectionPerVolume(ma 315 excitationXS->CrossSectionPerVolume(material,G4Electron::ElectronDefinition(),k*0.511/3728,tmp,tmp) 316 /material->GetAtomicNumDensityVector() 316 /material->GetAtomicNumDensityVector()[1]; 317 317 318 if ( particleDefinition == alphaPlusDef << 318 if ( particleDefinition == instance->GetIon("alpha+") ) 319 crossSection = crossSection + sigmaEx 319 crossSection = crossSection + sigmaExcitation ; 320 320 321 if ( particleDefinition == heliumDef ) << 321 if ( particleDefinition == instance->GetIon("helium") ) 322 crossSection = crossSection + 2*sigmaE 322 crossSection = crossSection + 2*sigmaExcitation ; 323 323 324 delete excitationXS; 324 delete excitationXS; 325 325 326 } 326 } 327 */ 327 */ 328 328 329 } 329 } 330 330 331 if (verboseLevel > 2) 331 if (verboseLevel > 2) 332 { 332 { 333 G4cout << "_______________________________ 333 G4cout << "__________________________________" << G4endl; 334 G4cout << "G4DNAMillerGreenExcitationModel 334 G4cout << "G4DNAMillerGreenExcitationModel - XS INFO START" << G4endl; 335 G4cout << "Kinetic energy(eV)=" << k/eV << 335 G4cout << "Kinetic energy(eV)=" << k/eV << " particle : " << particleDefinition->GetParticleName() << G4endl; 336 G4cout << "Cross section per water molecul 336 G4cout << "Cross section per water molecule (cm^2)=" << crossSection/cm/cm << G4endl; 337 G4cout << "Cross section per water molecul 337 G4cout << "Cross section per water molecule (cm^-1)=" << crossSection*waterDensity/(1./cm) << G4endl; 338 // G4cout << " - Cross section per water m 338 // G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl; 339 G4cout << "G4DNAMillerGreenExcitationModel 339 G4cout << "G4DNAMillerGreenExcitationModel - XS INFO END" << G4endl; 340 } 340 } 341 341 342 return crossSection*waterDensity; 342 return crossSection*waterDensity; 343 } 343 } 344 344 345 //....oooOO0OOooo........oooOO0OOooo........oo 345 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 346 346 347 void G4DNAMillerGreenExcitationModel::SampleSe 347 void G4DNAMillerGreenExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/, 348 348 const G4MaterialCutsCouple* /*couple*/, 349 349 const G4DynamicParticle* aDynamicParticle, 350 350 G4double, 351 351 G4double) 352 { 352 { 353 353 354 if (verboseLevel > 3) 354 if (verboseLevel > 3) 355 G4cout << "Calling SampleSecondaries() of 355 G4cout << "Calling SampleSecondaries() of G4DNAMillerGreenExcitationModel" << G4endl; 356 356 357 G4double particleEnergy0 = aDynamicParticle- 357 G4double particleEnergy0 = aDynamicParticle->GetKineticEnergy(); 358 358 359 G4int level = RandomSelect(particleEnergy0,a 359 G4int level = RandomSelect(particleEnergy0,aDynamicParticle->GetDefinition()); 360 360 361 // Dingfelder's excitation levels 361 // Dingfelder's excitation levels 362 const G4double excitation[]={ 8.17*eV, 10.13 362 const G4double excitation[]={ 8.17*eV, 10.13*eV, 11.31*eV, 12.91*eV, 14.50*eV}; 363 G4double excitationEnergy = excitation[level 363 G4double excitationEnergy = excitation[level]; 364 364 365 G4double newEnergy = 0.; 365 G4double newEnergy = 0.; 366 366 367 if (!statCode) newEnergy = particleEnergy0 - 367 if (!statCode) newEnergy = particleEnergy0 - excitationEnergy; 368 368 369 else newEnergy = particleEnergy0; 369 else newEnergy = particleEnergy0; 370 370 371 if (newEnergy>0) 371 if (newEnergy>0) 372 { 372 { 373 fParticleChangeForGamma->ProposeMomentumDi 373 fParticleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection()); 374 fParticleChangeForGamma->SetProposedKineti 374 fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy); 375 fParticleChangeForGamma->ProposeLocalEnerg 375 fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy); 376 376 377 const G4Track * theIncomingTrack = fPartic 377 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack(); 378 G4DNAChemistryManager::Instance()->CreateW 378 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eExcitedMolecule, 379 level, theIncomingTrack); 379 level, theIncomingTrack); 380 380 381 } 381 } 382 382 383 } 383 } 384 384 385 //....oooOO0OOooo........oooOO0OOooo........oo 385 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 386 386 387 G4double G4DNAMillerGreenExcitationModel::GetP 387 G4double G4DNAMillerGreenExcitationModel::GetPartialCrossSection(const G4Material*, 388 G4in 388 G4int level, 389 cons 389 const G4ParticleDefinition* particleDefinition, 390 G4do 390 G4double kineticEnergy) 391 { 391 { 392 return PartialCrossSection(kineticEnergy, le 392 return PartialCrossSection(kineticEnergy, level, particleDefinition); 393 } 393 } 394 394 395 //....oooOO0OOooo........oooOO0OOooo........oo 395 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 396 396 397 G4double G4DNAMillerGreenExcitationModel::Part 397 G4double G4DNAMillerGreenExcitationModel::PartialCrossSection(G4double k, G4int excitationLevel, 398 398 const G4ParticleDefinition* particleDefinition) 399 { 399 { 400 // ( ( z * aj 400 // ( ( z * aj ) ^ omegaj ) * ( t - ej ) ^ nu 401 // sigma(t) = zEff^2 * sigma0 * ------------ 401 // sigma(t) = zEff^2 * sigma0 * -------------------------------------------- 402 // jj ^ ( omeg 402 // jj ^ ( omegaj + nu ) + t ^ ( omegaj + nu ) 403 // 403 // 404 // where t is the kinetic energy corrected b 404 // where t is the kinetic energy corrected by Helium mass over proton mass for Helium ions 405 // 405 // 406 // zEff is: 406 // zEff is: 407 // 1 for protons 407 // 1 for protons 408 // 2 for alpha++ 408 // 2 for alpha++ 409 // and 2 - c1 S_1s - c2 S_2s - c3 S_2p for 409 // and 2 - c1 S_1s - c2 S_2s - c3 S_2p for alpha+ and He 410 // 410 // 411 // Dingfelder et al., RPC 59, 255-275, 2000 411 // Dingfelder et al., RPC 59, 255-275, 2000 from Miller and Green (1973) 412 // Formula (34) and Table 2 412 // Formula (34) and Table 2 413 413 414 const G4double sigma0(1.E+8 * barn); 414 const G4double sigma0(1.E+8 * barn); 415 const G4double nu(1.); 415 const G4double nu(1.); 416 const G4double aj[]={876.*eV, 2084.* eV, 137 416 const G4double aj[]={876.*eV, 2084.* eV, 1373.*eV, 692.*eV, 900.*eV}; 417 const G4double jj[]={19820.*eV, 23490.*eV, 2 417 const G4double jj[]={19820.*eV, 23490.*eV, 27770.*eV, 30830.*eV, 33080.*eV}; 418 const G4double omegaj[]={0.85, 0.88, 0.88, 0 418 const G4double omegaj[]={0.85, 0.88, 0.88, 0.78, 0.78}; 419 419 420 // Dingfelder's excitation levels 420 // Dingfelder's excitation levels 421 const G4double Eliq[5]={ 8.17*eV, 10.13*eV, 421 const G4double Eliq[5]={ 8.17*eV, 10.13*eV, 11.31*eV, 12.91*eV, 14.50*eV}; 422 422 423 G4int particleTypeIndex = 0; 423 G4int particleTypeIndex = 0; 424 << 424 G4DNAGenericIonsManager* instance; 425 if (particleDefinition == protonDef) particl << 425 instance = G4DNAGenericIonsManager::Instance(); 426 if (particleDefinition == hydrogenDef) parti << 426 427 if (particleDefinition == alphaPlusPlusDef) << 427 if (particleDefinition == G4Proton::ProtonDefinition()) particleTypeIndex=0; 428 if (particleDefinition == alphaPlusDef) part << 428 if (particleDefinition == instance->GetIon("hydrogen")) particleTypeIndex=0; 429 if (particleDefinition == heliumDef) particl << 429 if (particleDefinition == instance->GetIon("alpha++")) particleTypeIndex=1; >> 430 if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=2; >> 431 if (particleDefinition == instance->GetIon("helium")) particleTypeIndex=3; 430 432 431 G4double tCorrected; 433 G4double tCorrected; 432 tCorrected = k * kineticEnergyCorrection[par 434 tCorrected = k * kineticEnergyCorrection[particleTypeIndex]; 433 435 434 // SI - added protection 436 // SI - added protection 435 if (tCorrected < Eliq[excitationLevel]) retu 437 if (tCorrected < Eliq[excitationLevel]) return 0; 436 // 438 // 437 439 438 G4int z = 10; 440 G4int z = 10; 439 441 440 G4double numerator; 442 G4double numerator; 441 numerator = gpow->powA(z * aj[excitationLeve << 443 numerator = std::pow(z * aj[excitationLevel], omegaj[excitationLevel]) * 442 gpow->powA(tCorrected - Eliq[exc << 444 std::pow(tCorrected - Eliq[excitationLevel], nu); 443 445 444 // H case : see S. Uehara et al. IJRB 77, 2, 446 // H case : see S. Uehara et al. IJRB 77, 2, 139-154 (2001) - section 3.3 445 447 446 if (particleDefinition == hydrogenDef) << 448 if (particleDefinition == instance->GetIon("hydrogen")) 447 numerator = gpow->powA(z * 0.75*aj[excit << 449 numerator = std::pow(z * 0.75*aj[excitationLevel], omegaj[excitationLevel]) * 448 gpow->powA(tCorrected - Eliq << 450 std::pow(tCorrected - Eliq[excitationLevel], nu); 449 451 450 452 451 G4double power; 453 G4double power; 452 power = omegaj[excitationLevel] + nu; 454 power = omegaj[excitationLevel] + nu; 453 455 454 G4double denominator; 456 G4double denominator; 455 denominator = gpow->powA(jj[excitationLevel] << 457 denominator = std::pow(jj[excitationLevel], power) + std::pow(tCorrected, power); 456 458 457 G4double zEff = particleDefinition->GetPDGCh 459 G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber(); 458 460 459 zEff -= ( sCoefficient[0][particleTypeIndex] 461 zEff -= ( sCoefficient[0][particleTypeIndex] * S_1s(k, Eliq[excitationLevel], slaterEffectiveCharge[0][particleTypeIndex], 1.) + 460 sCoefficient[1][particleTypeIndex] * 462 sCoefficient[1][particleTypeIndex] * S_2s(k, Eliq[excitationLevel], slaterEffectiveCharge[1][particleTypeIndex], 2.) + 461 sCoefficient[2][particleTypeIndex] * 463 sCoefficient[2][particleTypeIndex] * S_2p(k, Eliq[excitationLevel], slaterEffectiveCharge[2][particleTypeIndex], 2.) ); 462 464 463 if (particleDefinition == hydrogenDef) zEff << 465 if (particleDefinition == instance->GetIon("hydrogen")) zEff = 1.; 464 466 465 G4double cross = sigma0 * zEff * zEff * nume 467 G4double cross = sigma0 * zEff * zEff * numerator / denominator; 466 468 467 469 468 return cross; 470 return cross; 469 } 471 } 470 472 471 //....oooOO0OOooo........oooOO0OOooo........oo 473 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 472 474 473 G4int G4DNAMillerGreenExcitationModel::RandomS 475 G4int G4DNAMillerGreenExcitationModel::RandomSelect(G4double k,const G4ParticleDefinition* particle) 474 { 476 { 475 G4int i = nLevels; 477 G4int i = nLevels; 476 G4double value = 0.; 478 G4double value = 0.; 477 std::deque<G4double> values; 479 std::deque<G4double> values; 478 480 479 if ( particle == alphaPlusPlusDef || << 481 G4DNAGenericIonsManager *instance; 480 particle == protonDef|| << 482 instance = G4DNAGenericIonsManager::Instance(); 481 particle == hydrogenDef || << 483 482 particle == alphaPlusDef || << 484 if ( particle == instance->GetIon("alpha++") || 483 particle == heliumDef << 485 particle == G4Proton::ProtonDefinition()|| >> 486 particle == instance->GetIon("hydrogen") || >> 487 particle == instance->GetIon("alpha+") || >> 488 particle == instance->GetIon("helium") 484 ) 489 ) 485 { 490 { 486 while (i > 0) 491 while (i > 0) 487 { 492 { 488 i--; 493 i--; 489 G4double partial = PartialCrossSection(k 494 G4double partial = PartialCrossSection(k,i,particle); 490 values.push_front(partial); 495 values.push_front(partial); 491 value += partial; 496 value += partial; 492 } 497 } 493 498 494 value *= G4UniformRand(); 499 value *= G4UniformRand(); 495 500 496 i = nLevels; 501 i = nLevels; 497 502 498 while (i > 0) 503 while (i > 0) 499 { 504 { 500 i--; 505 i--; 501 if (values[i] > value) return i; 506 if (values[i] > value) return i; 502 value -= values[i]; 507 value -= values[i]; 503 } 508 } 504 } 509 } 505 510 506 /* 511 /* 507 // add ONE or TWO electron-water excitation 512 // add ONE or TWO electron-water excitation for alpha+ and helium 508 513 509 if ( particle == alphaPlusDef << 514 if ( particle == instance->GetIon("alpha+") 510 || 515 || 511 particle == heliumDef << 516 particle == instance->GetIon("helium") 512 ) 517 ) 513 { 518 { 514 while (i>0) 519 while (i>0) 515 { 520 { 516 i--; 521 i--; 517 522 518 G4DNAEmfietzoglouExcitationModel * e 523 G4DNAEmfietzoglouExcitationModel * excitationXS = new G4DNAEmfietzoglouExcitationModel(); 519 excitationXS->Initialise(G4Electron: 524 excitationXS->Initialise(G4Electron::ElectronDefinition()); 520 525 521 G4double sigmaExcitation=0; 526 G4double sigmaExcitation=0; 522 527 523 if (k*0.511/3728 > 8.23*eV && k*0.511/37 528 if (k*0.511/3728 > 8.23*eV && k*0.511/3728 < 10*MeV ) sigmaExcitation = excitationXS->PartialCrossSection(k*0.511/3728,i); 524 529 525 G4double partial = PartialCrossSection(k 530 G4double partial = PartialCrossSection(k,i,particle); 526 531 527 if (particle == alphaPlusDef) partial = << 532 if (particle == instance->GetIon("alpha+")) partial = PartialCrossSection(k,i,particle) + sigmaExcitation; 528 if (particle == heliumDef) partial = Par << 533 if (particle == instance->GetIon("helium")) partial = PartialCrossSection(k,i,particle) + 2*sigmaExcitation; 529 534 530 values.push_front(partial); 535 values.push_front(partial); 531 value += partial; 536 value += partial; 532 delete excitationXS; 537 delete excitationXS; 533 } 538 } 534 539 535 value*=G4UniformRand(); 540 value*=G4UniformRand(); 536 541 537 i=5; 542 i=5; 538 while (i>0) 543 while (i>0) 539 { 544 { 540 i--; 545 i--; 541 546 542 if (values[i]>value) return i; 547 if (values[i]>value) return i; 543 548 544 value-=values[i]; 549 value-=values[i]; 545 } 550 } 546 } 551 } 547 */ 552 */ 548 553 549 return 0; 554 return 0; 550 } 555 } 551 556 552 //....oooOO0OOooo........oooOO0OOooo........oo 557 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 553 558 554 G4double G4DNAMillerGreenExcitationModel::Sum( 559 G4double G4DNAMillerGreenExcitationModel::Sum(G4double k, const G4ParticleDefinition* particle) 555 { 560 { 556 G4double totalCrossSection = 0.; 561 G4double totalCrossSection = 0.; 557 562 558 for (G4int i=0; i<nLevels; i++) 563 for (G4int i=0; i<nLevels; i++) 559 { 564 { 560 totalCrossSection += PartialCrossSection(k 565 totalCrossSection += PartialCrossSection(k,i,particle); 561 } 566 } 562 return totalCrossSection; 567 return totalCrossSection; 563 } 568 } 564 569 565 //....oooOO0OOooo........oooOO0OOooo........oo 570 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 566 571 567 G4double G4DNAMillerGreenExcitationModel::S_1s 572 G4double G4DNAMillerGreenExcitationModel::S_1s(G4double t, 568 573 G4double energyTransferred, 569 574 G4double _slaterEffectiveCharge, 570 575 G4double shellNumber) 571 { 576 { 572 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2) 577 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2) 573 // Dingfelder, in Chattanooga 2005 proceedin 578 // Dingfelder, in Chattanooga 2005 proceedings, formula (7) 574 579 575 G4double r = R(t, energyTransferred, _slater 580 G4double r = R(t, energyTransferred, _slaterEffectiveCharge, shellNumber); 576 G4double value = 1. - G4Exp(-2 * r) * ( ( 2. 581 G4double value = 1. - G4Exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. ); 577 582 578 return value; 583 return value; 579 } 584 } 580 585 581 586 582 //....oooOO0OOooo........oooOO0OOooo........oo 587 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 583 588 584 G4double G4DNAMillerGreenExcitationModel::S_2s 589 G4double G4DNAMillerGreenExcitationModel::S_2s(G4double t, 585 590 G4double energyTransferred, 586 591 G4double _slaterEffectiveCharge, 587 592 G4double shellNumber) 588 { 593 { 589 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4) 594 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4) 590 // Dingfelder, in Chattanooga 2005 proceedin 595 // Dingfelder, in Chattanooga 2005 proceedings, formula (8) 591 596 592 G4double r = R(t, energyTransferred, _slater 597 G4double r = R(t, energyTransferred, _slaterEffectiveCharge, shellNumber); 593 G4double value = 1. - G4Exp(-2 * r) * (((2. 598 G4double value = 1. - G4Exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.); 594 599 595 return value; 600 return value; 596 601 597 } 602 } 598 603 599 //....oooOO0OOooo........oooOO0OOooo........oo 604 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 600 605 601 G4double G4DNAMillerGreenExcitationModel::S_2p 606 G4double G4DNAMillerGreenExcitationModel::S_2p(G4double t, 602 607 G4double energyTransferred, 603 608 G4double _slaterEffectiveCharge, 604 609 G4double shellNumber) 605 { 610 { 606 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^ 611 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4) 607 // Dingfelder, in Chattanooga 2005 proceedin 612 // Dingfelder, in Chattanooga 2005 proceedings, formula (9) 608 613 609 G4double r = R(t, energyTransferred, _slater 614 G4double r = R(t, energyTransferred, _slaterEffectiveCharge, shellNumber); 610 G4double value = 1. - G4Exp(-2 * r) * (((( 615 G4double value = 1. - G4Exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.); 611 616 612 return value; 617 return value; 613 } 618 } 614 619 615 //....oooOO0OOooo........oooOO0OOooo........oo 620 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 616 621 617 G4double G4DNAMillerGreenExcitationModel::R(G4 622 G4double G4DNAMillerGreenExcitationModel::R(G4double t, 618 G4 623 G4double energyTransferred, 619 G4 624 G4double _slaterEffectiveCharge, 620 G4 625 G4double shellNumber) 621 { 626 { 622 // tElectron = m_electron / m_alpha * t 627 // tElectron = m_electron / m_alpha * t 623 // Dingfelder, in Chattanooga 2005 proceedin 628 // Dingfelder, in Chattanooga 2005 proceedings, p 4 624 629 625 G4double tElectron = 0.511/3728. * t; 630 G4double tElectron = 0.511/3728. * t; 626 631 627 // The following is provided by M. Dingfelde 632 // The following is provided by M. Dingfelder 628 G4double H = 2.*13.60569172 * eV; 633 G4double H = 2.*13.60569172 * eV; 629 G4double value = std::sqrt ( 2. * tElectron 634 G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) * (_slaterEffectiveCharge/shellNumber); 630 635 631 return value; 636 return value; 632 } 637 } 633 638 634 639