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 // $Id$ 26 // 27 // 27 28 28 #include "G4DNADingfelderChargeIncreaseModel.h 29 #include "G4DNADingfelderChargeIncreaseModel.hh" 29 #include "G4PhysicalConstants.hh" 30 #include "G4PhysicalConstants.hh" 30 #include "G4SystemOfUnits.hh" 31 #include "G4SystemOfUnits.hh" 31 #include "G4DNAMolecularMaterial.hh" 32 #include "G4DNAMolecularMaterial.hh" 32 #include "G4Log.hh" << 33 #include "G4Pow.hh" << 34 #include "G4Alpha.hh" << 35 << 36 static G4Pow * gpow = G4Pow::GetInstance(); << 37 33 38 //....oooOO0OOooo........oooOO0OOooo........oo 34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 39 35 40 using namespace std; 36 using namespace std; 41 37 42 //....oooOO0OOooo........oooOO0OOooo........oo 38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 43 39 44 G4DNADingfelderChargeIncreaseModel::G4DNADingf 40 G4DNADingfelderChargeIncreaseModel::G4DNADingfelderChargeIncreaseModel(const G4ParticleDefinition*, 45 41 const G4String& nam) : 46 G4VEmModel(nam) << 42 G4VEmModel(nam), isInitialised(false) 47 { 43 { >> 44 // nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER"); >> 45 fpMolWaterDensity = 0; >> 46 >> 47 numberOfPartialCrossSections[0] = 0; >> 48 numberOfPartialCrossSections[1] = 0; >> 49 >> 50 verboseLevel = 0; >> 51 // Verbosity scale: >> 52 // 0 = nothing >> 53 // 1 = warning for energy non-conservation >> 54 // 2 = details of energy budget >> 55 // 3 = calculation of cross sections, file openings, sampling of atoms >> 56 // 4 = entering in methods >> 57 48 if (verboseLevel > 0) 58 if (verboseLevel > 0) 49 { 59 { 50 G4cout << "Dingfelder charge increase mode 60 G4cout << "Dingfelder charge increase model is constructed " << G4endl; 51 } 61 } >> 62 fParticleChangeForGamma = 0; >> 63 >> 64 // Selection of stationary mode >> 65 >> 66 statCode = false; 52 } 67 } 53 68 54 //....oooOO0OOooo........oooOO0OOooo........oo 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 55 70 >> 71 G4DNADingfelderChargeIncreaseModel::~G4DNADingfelderChargeIncreaseModel() >> 72 {} >> 73 >> 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 75 56 void G4DNADingfelderChargeIncreaseModel::Initi 76 void G4DNADingfelderChargeIncreaseModel::Initialise(const G4ParticleDefinition* particle, 57 77 const G4DataVector& /*cuts*/) 58 { 78 { 59 79 60 if (verboseLevel > 3) 80 if (verboseLevel > 3) 61 { 81 { 62 G4cout << "Calling G4DNADingfelderChargeIn 82 G4cout << "Calling G4DNADingfelderChargeIncreaseModel::Initialise()" 63 << G4endl; 83 << G4endl; 64 } 84 } 65 85 66 // Energy limits 86 // Energy limits 67 87 68 G4DNAGenericIonsManager *instance; 88 G4DNAGenericIonsManager *instance; 69 instance = G4DNAGenericIonsManager::Instance 89 instance = G4DNAGenericIonsManager::Instance(); 70 hydrogenDef = instance->GetIon("hydrogen"); << 90 G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen"); 71 alphaPlusPlusDef = G4Alpha::Alpha(); << 91 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+"); 72 alphaPlusDef = instance->GetIon("alpha+"); << 92 G4ParticleDefinition* heliumDef = instance->GetIon("helium"); 73 heliumDef = instance->GetIon("helium"); << 74 93 75 G4String hydrogen; 94 G4String hydrogen; 76 G4String alphaPlus; 95 G4String alphaPlus; 77 G4String helium; 96 G4String helium; 78 97 79 // Limits << 98 // LIMITS 80 99 81 hydrogen = hydrogenDef->GetParticleName(); 100 hydrogen = hydrogenDef->GetParticleName(); 82 lowEnergyLimit[hydrogen] = 100. * eV; 101 lowEnergyLimit[hydrogen] = 100. * eV; 83 highEnergyLimit[hydrogen] = 100. * MeV; 102 highEnergyLimit[hydrogen] = 100. * MeV; 84 103 85 alphaPlus = alphaPlusDef->GetParticleName(); 104 alphaPlus = alphaPlusDef->GetParticleName(); 86 lowEnergyLimit[alphaPlus] = 1. * keV; 105 lowEnergyLimit[alphaPlus] = 1. * keV; 87 highEnergyLimit[alphaPlus] = 400. * MeV; 106 highEnergyLimit[alphaPlus] = 400. * MeV; 88 107 89 helium = heliumDef->GetParticleName(); 108 helium = heliumDef->GetParticleName(); 90 lowEnergyLimit[helium] = 1. * keV; 109 lowEnergyLimit[helium] = 1. * keV; 91 highEnergyLimit[helium] = 400. * MeV; 110 highEnergyLimit[helium] = 400. * MeV; 92 111 93 // 112 // 94 113 95 if (particle==hydrogenDef) 114 if (particle==hydrogenDef) 96 { 115 { 97 SetLowEnergyLimit(lowEnergyLimit[hydrogen] 116 SetLowEnergyLimit(lowEnergyLimit[hydrogen]); 98 SetHighEnergyLimit(highEnergyLimit[hydroge 117 SetHighEnergyLimit(highEnergyLimit[hydrogen]); 99 } 118 } 100 119 101 if (particle==alphaPlusDef) 120 if (particle==alphaPlusDef) 102 { 121 { 103 SetLowEnergyLimit(lowEnergyLimit[alphaPlus 122 SetLowEnergyLimit(lowEnergyLimit[alphaPlus]); 104 SetHighEnergyLimit(highEnergyLimit[alphaPl 123 SetHighEnergyLimit(highEnergyLimit[alphaPlus]); 105 } 124 } 106 125 107 if (particle==heliumDef) 126 if (particle==heliumDef) 108 { 127 { 109 SetLowEnergyLimit(lowEnergyLimit[helium]); 128 SetLowEnergyLimit(lowEnergyLimit[helium]); 110 SetHighEnergyLimit(highEnergyLimit[helium] 129 SetHighEnergyLimit(highEnergyLimit[helium]); 111 } 130 } 112 131 113 // Final state 132 // Final state 114 133 115 //ALPHA+ 134 //ALPHA+ 116 135 117 f0[0][0]=1.; 136 f0[0][0]=1.; 118 a0[0][0]=2.25; 137 a0[0][0]=2.25; 119 a1[0][0]=-0.75; 138 a1[0][0]=-0.75; 120 b0[0][0]=-32.10; 139 b0[0][0]=-32.10; 121 c0[0][0]=0.600; 140 c0[0][0]=0.600; 122 d0[0][0]=2.40; 141 d0[0][0]=2.40; 123 x0[0][0]=4.60; 142 x0[0][0]=4.60; 124 143 125 x1[0][0]=-1.; 144 x1[0][0]=-1.; 126 b1[0][0]=-1.; 145 b1[0][0]=-1.; 127 146 128 numberOfPartialCrossSections[0]=1; 147 numberOfPartialCrossSections[0]=1; 129 148 130 //HELIUM 149 //HELIUM 131 150 132 f0[0][1]=1.; 151 f0[0][1]=1.; 133 a0[0][1]=2.25; 152 a0[0][1]=2.25; 134 a1[0][1]=-0.75; 153 a1[0][1]=-0.75; 135 b0[0][1]=-30.93; 154 b0[0][1]=-30.93; 136 c0[0][1]=0.590; 155 c0[0][1]=0.590; 137 d0[0][1]=2.35; 156 d0[0][1]=2.35; 138 x0[0][1]=4.29; 157 x0[0][1]=4.29; 139 158 140 f0[1][1]=1.; 159 f0[1][1]=1.; 141 a0[1][1]=2.25; 160 a0[1][1]=2.25; 142 a1[1][1]=-0.75; 161 a1[1][1]=-0.75; 143 b0[1][1]=-32.61; 162 b0[1][1]=-32.61; 144 c0[1][1]=0.435; 163 c0[1][1]=0.435; 145 d0[1][1]=2.70; 164 d0[1][1]=2.70; 146 x0[1][1]=4.45; 165 x0[1][1]=4.45; 147 166 148 x1[0][1]=-1.; 167 x1[0][1]=-1.; 149 b1[0][1]=-1.; 168 b1[0][1]=-1.; 150 169 151 x1[1][1]=-1.; 170 x1[1][1]=-1.; 152 b1[1][1]=-1.; 171 b1[1][1]=-1.; 153 172 154 numberOfPartialCrossSections[1]=2; 173 numberOfPartialCrossSections[1]=2; 155 174 156 // 175 // 157 176 158 if( verboseLevel>0 ) 177 if( verboseLevel>0 ) 159 { 178 { 160 G4cout << "Dingfelder charge increase mode 179 G4cout << "Dingfelder charge increase model is initialized " << G4endl 161 << "Energy range: " 180 << "Energy range: " 162 << LowEnergyLimit() / keV << " keV - " 181 << LowEnergyLimit() / keV << " keV - " 163 << HighEnergyLimit() / MeV << " MeV for " 182 << HighEnergyLimit() / MeV << " MeV for " 164 << particle->GetParticleName() 183 << particle->GetParticleName() 165 << G4endl; 184 << G4endl; 166 } 185 } 167 186 168 // Initialize water density pointer 187 // Initialize water density pointer 169 fpMolWaterDensity = G4DNAMolecularMaterial:: 188 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER")); 170 189 171 if (isInitialised) return; << 190 if (isInitialised) 172 << 191 { return;} 173 fParticleChangeForGamma = GetParticleChangeF 192 fParticleChangeForGamma = GetParticleChangeForGamma(); 174 isInitialised = true; 193 isInitialised = true; 175 } 194 } 176 195 177 //....oooOO0OOooo........oooOO0OOooo........oo 196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 178 197 179 G4double G4DNADingfelderChargeIncreaseModel::C 198 G4double G4DNADingfelderChargeIncreaseModel::CrossSectionPerVolume(const G4Material* material, 180 199 const G4ParticleDefinition* particleDefinition, 181 200 G4double k, 182 201 G4double, 183 202 G4double) 184 { 203 { 185 if (verboseLevel > 3) 204 if (verboseLevel > 3) 186 { 205 { 187 G4cout 206 G4cout 188 << "Calling CrossSectionPerVolume() of 207 << "Calling CrossSectionPerVolume() of G4DNADingfelderChargeIncreaseModel" 189 << G4endl; 208 << G4endl; 190 } 209 } 191 210 192 // Calculate total cross section for model 211 // Calculate total cross section for model 193 212 >> 213 G4DNAGenericIonsManager *instance; >> 214 instance = G4DNAGenericIonsManager::Instance(); >> 215 194 if ( 216 if ( 195 particleDefinition != hydrogenDef << 217 particleDefinition != instance->GetIon("hydrogen") 196 && 218 && 197 particleDefinition != alphaPlusDef << 219 particleDefinition != instance->GetIon("alpha+") 198 && 220 && 199 particleDefinition != heliumDef << 221 particleDefinition != instance->GetIon("helium") 200 ) 222 ) 201 223 202 return 0; 224 return 0; 203 225 204 G4double lowLim = 0; 226 G4double lowLim = 0; 205 G4double highLim = 0; 227 G4double highLim = 0; 206 G4double totalCrossSection = 0.; 228 G4double totalCrossSection = 0.; 207 229 208 G4double waterDensity = (*fpMolWaterDensity) 230 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()]; 209 231 210 const G4String& particleName = particleDefin << 232 if(waterDensity!= 0.0) >> 233 // if (material == nistwater || material->GetBaseMaterial() == nistwater) >> 234 { >> 235 const G4String& particleName = particleDefinition->GetParticleName(); 211 236 212 std::map< G4String,G4double,std::less<G4Stri << 237 std::map< G4String,G4double,std::less<G4String> >::iterator pos1; 213 pos1 = lowEnergyLimit.find(particleName); << 238 pos1 = lowEnergyLimit.find(particleName); 214 239 215 if (pos1 != lowEnergyLimit.end()) << 240 if (pos1 != lowEnergyLimit.end()) 216 { << 241 { 217 lowLim = pos1->second; << 242 lowLim = pos1->second; 218 } << 243 } 219 244 220 std::map< G4String,G4double,std::less<G4Stri << 245 std::map< G4String,G4double,std::less<G4String> >::iterator pos2; 221 pos2 = highEnergyLimit.find(particleName); << 246 pos2 = highEnergyLimit.find(particleName); 222 247 223 if (pos2 != highEnergyLimit.end()) << 248 if (pos2 != highEnergyLimit.end()) 224 { << 249 { 225 highLim = pos2->second; << 250 highLim = pos2->second; 226 } << 251 } 227 252 228 if (k >= lowLim && k <= highLim) << 253 if (k >= lowLim && k < highLim) 229 { << 230 //HYDROGEN << 231 if (particleDefinition == hydrogenDef) << 232 { 254 { 233 const G4double aa = 2.835; << 255 //HYDROGEN 234 const G4double bb = 0.310; << 256 if (particleDefinition == instance->GetIon("hydrogen")) 235 const G4double cc = 2.100; << 257 { 236 const G4double dd = 0.760; << 258 const G4double aa = 2.835; 237 const G4double fac = 1.0e-18; << 259 const G4double bb = 0.310; 238 const G4double rr = 13.606 * eV; << 260 const G4double cc = 2.100; 239 << 261 const G4double dd = 0.760; 240 G4double t = k / (proton_mass_c2/electro << 262 const G4double fac = 1.0e-18; 241 G4double x = t / rr; << 263 const G4double rr = 13.606 * eV; 242 G4double temp = 4.0 * pi * Bohr_radius/n << 264 243 G4double sigmal = temp * cc * (gpow->pow << 265 G4double t = k / (proton_mass_c2/electron_mass_c2); 244 G4double sigmah = temp * (aa * G4Log(1.0 << 266 G4double x = t / rr; 245 totalCrossSection = 1.0/(1.0/sigmal + 1. << 267 G4double temp = 4.0 * pi * Bohr_radius/nm * Bohr_radius/nm * fac; >> 268 G4double sigmal = temp * cc * (std::pow(x,dd)); >> 269 G4double sigmah = temp * (aa * std::log(1.0 + x) + bb) / x; >> 270 totalCrossSection = 1.0/(1.0/sigmal + 1.0/sigmah) *m*m; >> 271 } >> 272 else >> 273 { >> 274 totalCrossSection = Sum(k,particleDefinition); >> 275 } 246 } 276 } 247 else << 277 >> 278 if (verboseLevel > 2) 248 { 279 { 249 totalCrossSection = Sum(k,particleDefini << 280 G4cout << "__________________________________" << G4endl; >> 281 G4cout << "G4DNADingfelderChargeIncreaseModel - XS INFO START" << G4endl; >> 282 G4cout << "Kinetic energy(eV)=" << k/eV << " particle : " << particleName << G4endl; >> 283 G4cout << "Cross section per water molecule (cm^2)=" << totalCrossSection/cm/cm << G4endl; >> 284 G4cout << "Cross section per water molecule (cm^-1)=" << totalCrossSection*waterDensity/(1./cm) << G4endl; >> 285 // G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl; >> 286 G4cout << "G4DNADingfelderChargeIncreaseModel - XS INFO END" << G4endl; 250 } 287 } 251 } << 252 288 253 if (verboseLevel > 2) << 254 { << 255 G4cout << "_______________________________ << 256 G4cout << "G4DNADingfelderChargeIncreaseMo << 257 G4cout << "Kinetic energy(eV)=" << k/eV << << 258 G4cout << "Cross section per water molecul << 259 G4cout << "Cross section per water molecul << 260 // G4cout << " - Cross section per water << 261 // << sigma*material->GetAtomicNumDensity << 262 G4cout << "G4DNADingfelderChargeIncreaseMo << 263 } 289 } 264 290 265 return totalCrossSection*waterDensity; 291 return totalCrossSection*waterDensity; 266 // return totalCrossSection*material->GetAto << 292 // return totalCrossSection*material->GetAtomicNumDensityVector()[1]; 267 293 268 } 294 } 269 295 270 //....oooOO0OOooo........oooOO0OOooo........oo 296 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 271 297 272 void G4DNADingfelderChargeIncreaseModel::Sampl 298 void G4DNADingfelderChargeIncreaseModel::SampleSecondaries(std::vector< 273 299 G4DynamicParticle*>* fvect, 274 300 const G4MaterialCutsCouple* /*couple*/, 275 301 const G4DynamicParticle* aDynamicParticle, 276 302 G4double, 277 303 G4double) 278 { 304 { 279 if (verboseLevel > 3) 305 if (verboseLevel > 3) 280 { 306 { 281 G4cout 307 G4cout 282 << "Calling SampleSecondaries() of G4D 308 << "Calling SampleSecondaries() of G4DNADingfelderChargeIncreaseModel" 283 << G4endl; 309 << G4endl; 284 } 310 } 285 311 286 if (!statCode) fParticleChangeForGamma->Prop 312 if (!statCode) fParticleChangeForGamma->ProposeLocalEnergyDeposit(0.); 287 313 288 G4ParticleDefinition* definition = aDynamicP 314 G4ParticleDefinition* definition = aDynamicParticle->GetDefinition(); 289 315 290 G4double particleMass = definition->GetPDGMa 316 G4double particleMass = definition->GetPDGMass(); 291 317 292 G4double inK = aDynamicParticle->GetKineticE 318 G4double inK = aDynamicParticle->GetKineticEnergy(); 293 319 294 G4int finalStateIndex = RandomSelect(inK,def 320 G4int finalStateIndex = RandomSelect(inK,definition); 295 321 296 G4int n = NumberOfFinalStates(definition,fin 322 G4int n = NumberOfFinalStates(definition,finalStateIndex); 297 323 298 G4double outK = 0.; 324 G4double outK = 0.; 299 325 300 if (!statCode) outK = inK - IncomingParticle 326 if (!statCode) outK = inK - IncomingParticleBindingEnergyConstant(definition,finalStateIndex); 301 327 302 else outK = inK; 328 else outK = inK; 303 329 304 if (statCode) 330 if (statCode) 305 fParticleChangeForGamma-> << 331 fParticleChangeForGamma->ProposeLocalEnergyDeposit(IncomingParticleBindingEnergyConstant(definition,finalStateIndex)); 306 ProposeLocalEnergyDeposit(IncomingPartic << 307 332 308 fParticleChangeForGamma->ProposeTrackStatus( 333 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill); 309 334 >> 335 G4DNAGenericIonsManager* instance; >> 336 instance = G4DNAGenericIonsManager::Instance(); >> 337 310 G4double electronK; 338 G4double electronK; 311 if (definition == hydrogenDef) electronK = i << 339 if (definition == instance->GetIon("hydrogen")) electronK = inK*electron_mass_c2/proton_mass_c2; 312 else electronK = inK*electron_mass_c2/(parti 340 else electronK = inK*electron_mass_c2/(particleMass); 313 341 314 if (outK<0) 342 if (outK<0) 315 { 343 { 316 G4Exception("G4DNADingfelderChargeIncrease 344 G4Exception("G4DNADingfelderChargeIncreaseModel::SampleSecondaries","em0004", 317 FatalException,"Final kinetic energy i 345 FatalException,"Final kinetic energy is negative."); 318 } 346 } 319 347 320 auto dp = new G4DynamicParticle(OutgoingPart << 348 G4DynamicParticle* dp = new G4DynamicParticle(OutgoingParticleDefinition(definition,finalStateIndex), 321 aDynamicParticle->GetMomentumDirection() 349 aDynamicParticle->GetMomentumDirection(), 322 outK); 350 outK); 323 351 324 fvect->push_back(dp); 352 fvect->push_back(dp); 325 353 326 n = n - 1; 354 n = n - 1; 327 355 328 while (n>0) 356 while (n>0) 329 { 357 { 330 n--; 358 n--; 331 fvect->push_back(new G4DynamicParticle 359 fvect->push_back(new G4DynamicParticle 332 (G4Electron::Electron(), aDynamicParti 360 (G4Electron::Electron(), aDynamicParticle->GetMomentumDirection(), electronK) ); 333 } 361 } 334 } 362 } 335 363 336 //....oooOO0OOooo........oooOO0OOooo........oo 364 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 337 365 338 G4int G4DNADingfelderChargeIncreaseModel::Numb 366 G4int G4DNADingfelderChargeIncreaseModel::NumberOfFinalStates(G4ParticleDefinition* particleDefinition, 339 367 G4int finalStateIndex) 340 368 341 { 369 { 342 << 370 G4DNAGenericIonsManager* instance; 343 if (particleDefinition == hydrogenDef) << 371 instance = G4DNAGenericIonsManager::Instance(); >> 372 if (particleDefinition == instance->GetIon("hydrogen")) 344 return 2; 373 return 2; 345 << 374 if (particleDefinition == instance->GetIon("alpha+")) 346 if (particleDefinition == alphaPlusDef) << 347 return 2; 375 return 2; 348 376 349 if (particleDefinition == heliumDef) << 377 if (particleDefinition == instance->GetIon("helium")) 350 { 378 { 351 if (finalStateIndex == 0) 379 if (finalStateIndex == 0) 352 return 2; 380 return 2; 353 return 3; 381 return 3; 354 } 382 } 355 << 356 return 0; 383 return 0; 357 } 384 } 358 385 359 //....oooOO0OOooo........oooOO0OOooo........oo 386 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 360 387 361 G4ParticleDefinition* G4DNADingfelderChargeInc 388 G4ParticleDefinition* G4DNADingfelderChargeIncreaseModel::OutgoingParticleDefinition(G4ParticleDefinition* particleDefinition, 362 389 G4int finalStateIndex) 363 { 390 { 364 << 391 G4DNAGenericIonsManager * instance(G4DNAGenericIonsManager::Instance()); 365 if (particleDefinition == hydrogenDef) << 392 if (particleDefinition == instance->GetIon("hydrogen")) 366 return G4Proton::Proton(); 393 return G4Proton::Proton(); >> 394 if (particleDefinition == instance->GetIon("alpha+")) >> 395 return instance->GetIon("alpha++"); 367 396 368 if (particleDefinition == alphaPlusDef) << 397 if (particleDefinition == instance->GetIon("helium")) 369 return alphaPlusPlusDef; << 370 << 371 if (particleDefinition == heliumDef) << 372 { 398 { 373 if (finalStateIndex == 0) 399 if (finalStateIndex == 0) 374 return alphaPlusDef; << 400 return instance->GetIon("alpha+"); 375 return alphaPlusPlusDef; << 401 return instance->GetIon("alpha++"); 376 } 402 } 377 << 403 return 0; 378 return nullptr; << 379 } 404 } 380 405 381 //....oooOO0OOooo........oooOO0OOooo........oo 406 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 382 407 383 G4double G4DNADingfelderChargeIncreaseModel::I 408 G4double G4DNADingfelderChargeIncreaseModel::IncomingParticleBindingEnergyConstant(G4ParticleDefinition* particleDefinition, 384 409 G4int finalStateIndex) 385 { 410 { 386 << 411 G4DNAGenericIonsManager * instance(G4DNAGenericIonsManager::Instance()); 387 if (particleDefinition == hydrogenDef) << 412 if (particleDefinition == instance->GetIon("hydrogen")) 388 return 13.6 * eV; 413 return 13.6 * eV; 389 414 390 if (particleDefinition == alphaPlusDef) << 415 if (particleDefinition == instance->GetIon("alpha+")) 391 { 416 { 392 // Binding energy for He+ -> He++ + e- 417 // Binding energy for He+ -> He++ + e- 54.509 eV 393 // Binding energy for He -> He+ + e- 418 // Binding energy for He -> He+ + e- 24.587 eV 394 return 54.509 * eV; 419 return 54.509 * eV; 395 } 420 } 396 421 397 if (particleDefinition == heliumDef) << 422 if (particleDefinition == instance->GetIon("helium")) 398 { 423 { 399 // Binding energy for He+ -> He++ + e- 424 // Binding energy for He+ -> He++ + e- 54.509 eV 400 // Binding energy for He -> He+ + e- 425 // Binding energy for He -> He+ + e- 24.587 eV 401 426 402 if (finalStateIndex == 0) 427 if (finalStateIndex == 0) 403 return 24.587 * eV; 428 return 24.587 * eV; 404 return (54.509 + 24.587) * eV; 429 return (54.509 + 24.587) * eV; 405 } 430 } 406 431 407 return 0; 432 return 0; 408 } 433 } 409 434 410 //....oooOO0OOooo........oooOO0OOooo........oo 435 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 411 436 412 G4double G4DNADingfelderChargeIncreaseModel::P << 437 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 413 << 438 >> 439 G4double G4DNADingfelderChargeIncreaseModel::PartialCrossSection(G4double k, >> 440 G4int index, 414 441 const G4ParticleDefinition* particleDefinition) 415 { 442 { 416 G4int particleTypeIndex = 0; 443 G4int particleTypeIndex = 0; >> 444 G4DNAGenericIonsManager *instance; >> 445 instance = G4DNAGenericIonsManager::Instance(); 417 446 418 if (particleDefinition == alphaPlusDef) << 447 if (particleDefinition == instance->GetIon("alpha+")) 419 particleTypeIndex = 0; 448 particleTypeIndex = 0; 420 << 449 if (particleDefinition == instance->GetIon("helium")) 421 if (particleDefinition == heliumDef) << 422 particleTypeIndex = 1; 450 particleTypeIndex = 1; 423 451 424 // 452 // 425 // sigma(T) = f0 10 ^ y(log10(T/eV)) 453 // sigma(T) = f0 10 ^ y(log10(T/eV)) 426 // 454 // 427 // / a0 x + b0 x 455 // / a0 x + b0 x < x0 428 // | 456 // | 429 // y(x) = < a0 x + b0 - c0 (x - x0)^d0 x 457 // y(x) = < a0 x + b0 - c0 (x - x0)^d0 x0 <= x < x1 430 // | 458 // | 431 // \ a1 x + b1 x 459 // \ a1 x + b1 x >= x1 432 // 460 // 433 // 461 // 434 // f0, a0, a1, b0, b1, c0, d0, x0, x1 are pa 462 // f0, a0, a1, b0, b1, c0, d0, x0, x1 are parameters that change for protons and helium (0, +, ++) 435 // 463 // 436 // f0 has been added to the code in order to << 464 // f0 has been added to the code in order to manage partial (shell-dependent) cross sections (if no shell dependence is present. f0=1. Sum of f0 over the considered shells should give 1) 437 // (if no shell dependence is present. f0=1. << 438 // 465 // 439 // From Rad. Phys. and Chem. 59 (2000) 255-2 466 // From Rad. Phys. and Chem. 59 (2000) 255-275, M. Dingfelder et al. 440 // Inelastic-collision cross sections of liq 467 // Inelastic-collision cross sections of liquid water for interactions of energetic proton 441 // 468 // 442 469 443 if (x1[index][particleTypeIndex] < x0[index] 470 if (x1[index][particleTypeIndex] < x0[index][particleTypeIndex]) 444 { 471 { 445 // 472 // 446 // if x1 < x0 means that x1 and b1 will be << 473 // if x1 < x0 means that x1 and b1 will be calculated with the following formula (this piece of code is run on all alphas and not on protons) 447 // (this piece of code is run on all alpha << 448 // 474 // 449 // x1 = x0 + ((a0 - a1)/(c0 * d0)) ^ (1 / 475 // x1 = x0 + ((a0 - a1)/(c0 * d0)) ^ (1 / (d0 - 1)) 450 // 476 // 451 // b1 = (a0 - a1) * x1 + b0 - c0 * (x1 - x 477 // b1 = (a0 - a1) * x1 + b0 - c0 * (x1 - x0) ^ d0 452 // 478 // 453 479 454 x1[index][particleTypeIndex] = x0[index][p 480 x1[index][particleTypeIndex] = x0[index][particleTypeIndex] 455 + gpow->powA((a0[index][particleTypeIn << 481 + std::pow((a0[index][particleTypeIndex] - a1[index][particleTypeIndex]) 456 / (c0[index][particleTy 482 / (c0[index][particleTypeIndex] 457 * d0[index][particl 483 * d0[index][particleTypeIndex]), 458 1. / (d0[index][particleTyp 484 1. / (d0[index][particleTypeIndex] - 1.)); 459 b1[index][particleTypeIndex] = (a0[index][ 485 b1[index][particleTypeIndex] = (a0[index][particleTypeIndex] 460 - a1[index][particleTypeIndex]) * x1[i 486 - a1[index][particleTypeIndex]) * x1[index][particleTypeIndex] 461 + b0[index][particleTypeIndex] 487 + b0[index][particleTypeIndex] 462 - c0[index][particleTypeIndex] 488 - c0[index][particleTypeIndex] 463 * gpow->powA(x1[index][particleTyp << 489 * std::pow(x1[index][particleTypeIndex] 464 - x0[index][particl 490 - x0[index][particleTypeIndex], 465 d0[index][particleTypeI 491 d0[index][particleTypeIndex]); 466 } 492 } 467 493 468 G4double x(G4Log(k / eV)/gpow->logZ(10)); << 494 G4double x(std::log10(k / eV)); 469 G4double y; 495 G4double y; 470 496 471 if (x < x0[index][particleTypeIndex]) 497 if (x < x0[index][particleTypeIndex]) 472 y = a0[index][particleTypeIndex] * x + b0[ 498 y = a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex]; 473 else if (x < x1[index][particleTypeIndex]) 499 else if (x < x1[index][particleTypeIndex]) 474 y = a0[index][particleTypeIndex] * x + b0[ 500 y = a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex] 475 - c0[index][particleTypeIndex] 501 - c0[index][particleTypeIndex] 476 * gpow->powA(x - x0[index][particl << 502 * std::pow(x - x0[index][particleTypeIndex], 477 d0[index][particleTypeI 503 d0[index][particleTypeIndex]); 478 else 504 else 479 y = a1[index][particleTypeIndex] * x + b1[ 505 y = a1[index][particleTypeIndex] * x + b1[index][particleTypeIndex]; 480 506 481 return f0[index][particleTypeIndex] * gpow-> << 507 return f0[index][particleTypeIndex] * std::pow(10., y) * m * m; 482 508 483 } 509 } 484 510 485 //....oooOO0OOooo........oooOO0OOooo........oo 511 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 486 512 487 G4int G4DNADingfelderChargeIncreaseModel::Rand << 513 G4int G4DNADingfelderChargeIncreaseModel::RandomSelect(G4double k, 488 514 const G4ParticleDefinition* particleDefinition) 489 { 515 { 490 G4int particleTypeIndex = 0; 516 G4int particleTypeIndex = 0; >> 517 G4DNAGenericIonsManager *instance; >> 518 instance = G4DNAGenericIonsManager::Instance(); 491 519 492 if (particleDefinition == hydrogenDef) << 520 if (particleDefinition == instance->GetIon("hydrogen")) 493 return 0; 521 return 0; 494 << 522 if (particleDefinition == instance->GetIon("alpha+")) 495 if (particleDefinition == alphaPlusDef) << 496 particleTypeIndex = 0; 523 particleTypeIndex = 0; 497 << 524 if (particleDefinition == instance->GetIon("helium")) 498 if (particleDefinition == heliumDef) << 499 particleTypeIndex = 1; 525 particleTypeIndex = 1; 500 526 501 const G4int n = numberOfPartialCrossSections 527 const G4int n = numberOfPartialCrossSections[particleTypeIndex]; 502 auto values(new G4double[n]); << 528 G4double* values(new G4double[n]); 503 G4double value = 0; 529 G4double value = 0; 504 G4int i = n; 530 G4int i = n; 505 531 506 while (i > 0) 532 while (i > 0) 507 { 533 { 508 i--; 534 i--; 509 values[i] = PartialCrossSection(k, i, part 535 values[i] = PartialCrossSection(k, i, particleDefinition); 510 value += values[i]; 536 value += values[i]; 511 } 537 } 512 538 513 value *= G4UniformRand(); 539 value *= G4UniformRand(); 514 540 515 i = n; 541 i = n; 516 while (i > 0) 542 while (i > 0) 517 { 543 { 518 i--; 544 i--; 519 545 520 if (values[i] > value) 546 if (values[i] > value) 521 break; 547 break; 522 548 523 value -= values[i]; 549 value -= values[i]; 524 } 550 } 525 551 526 delete[] values; 552 delete[] values; 527 553 528 return i; 554 return i; 529 } 555 } 530 556 531 //....oooOO0OOooo........oooOO0OOooo........oo 557 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 532 558 533 G4double G4DNADingfelderChargeIncreaseModel::S << 559 G4double G4DNADingfelderChargeIncreaseModel::Sum(G4double k, 534 560 const G4ParticleDefinition* particleDefinition) 535 { 561 { 536 G4int particleTypeIndex = 0; 562 G4int particleTypeIndex = 0; >> 563 G4DNAGenericIonsManager *instance; >> 564 instance = G4DNAGenericIonsManager::Instance(); 537 565 538 if (particleDefinition == alphaPlusDef) << 566 if (particleDefinition == instance->GetIon("alpha+")) 539 particleTypeIndex = 0; 567 particleTypeIndex = 0; 540 << 568 if (particleDefinition == instance->GetIon("helium")) 541 if (particleDefinition == heliumDef) << 542 particleTypeIndex = 1; 569 particleTypeIndex = 1; 543 570 544 G4double totalCrossSection = 0.; 571 G4double totalCrossSection = 0.; 545 572 546 for (G4int i = 0; i < numberOfPartialCrossSe 573 for (G4int i = 0; i < numberOfPartialCrossSections[particleTypeIndex]; i++) 547 { 574 { 548 totalCrossSection += PartialCrossSection(k 575 totalCrossSection += PartialCrossSection(k, i, particleDefinition); 549 } 576 } 550 << 551 return totalCrossSection; 577 return totalCrossSection; 552 } 578 } 553 579