Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // 27 // Author: Luciano Pandola 28 // on base of G4LowEnergyIonisation developed by A.Forti and V.Ivanchenko 29 // 30 // History: 31 // -------- 32 // 12 Jan 2009 L Pandola Migration from process to model 33 // 03 Mar 2009 L Pandola Bug fix (release memory in the destructor) 34 // 15 Apr 2009 V Ivanchenko Cleanup initialisation and generation of secondaries: 35 // - apply internal high-energy limit only in constructor 36 // - do not apply low-energy limit (default is 0) 37 // - simplify sampling of deexcitation by using cut in energy 38 // - set activation of Auger "false" 39 // - remove initialisation of element selectors 40 // 19 May 2009 L Pandola Explicitely set to zero pointers deleted in 41 // Initialise(), since they might be checked later on 42 // 23 Oct 2009 L Pandola 43 // - atomic deexcitation managed via G4VEmModel::DeexcitationFlag() is 44 // set as "true" (default would be false) 45 // 12 Oct 2010 L Pandola 46 // - add debugging information about energy in 47 // SampleDeexcitationAlongStep() 48 // - generate fluorescence SampleDeexcitationAlongStep() only above 49 // the cuts. 50 // 01 Jun 2011 V Ivanchenko general cleanup - all old deexcitation code removed 51 // 52 53 #include "G4LivermoreIonisationModel.hh" 54 #include "G4PhysicalConstants.hh" 55 #include "G4SystemOfUnits.hh" 56 #include "G4ParticleDefinition.hh" 57 #include "G4MaterialCutsCouple.hh" 58 #include "G4ProductionCutsTable.hh" 59 #include "G4DynamicParticle.hh" 60 #include "G4Element.hh" 61 #include "G4ParticleChangeForLoss.hh" 62 #include "G4Electron.hh" 63 #include "G4CrossSectionHandler.hh" 64 #include "G4VEMDataSet.hh" 65 #include "G4eIonisationCrossSectionHandler.hh" 66 #include "G4eIonisationSpectrum.hh" 67 #include "G4VEnergySpectrum.hh" 68 #include "G4SemiLogInterpolation.hh" 69 #include "G4AtomicTransitionManager.hh" 70 #include "G4AtomicShell.hh" 71 #include "G4DeltaAngle.hh" 72 73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 74 75 76 G4LivermoreIonisationModel::G4LivermoreIonisationModel(const G4ParticleDefinition*, 77 const G4String& nam) : 78 G4VEmModel(nam), fParticleChange(nullptr), 79 crossSectionHandler(nullptr), energySpectrum(nullptr), 80 isInitialised(false) 81 { 82 fIntrinsicLowEnergyLimit = 12.*eV; 83 fIntrinsicHighEnergyLimit = 100.0*GeV; 84 85 verboseLevel = 0; 86 SetAngularDistribution(new G4DeltaAngle()); 87 88 transitionManager = G4AtomicTransitionManager::Instance(); 89 } 90 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 92 93 G4LivermoreIonisationModel::~G4LivermoreIonisationModel() 94 { 95 delete energySpectrum; 96 delete crossSectionHandler; 97 } 98 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 100 101 void G4LivermoreIonisationModel::Initialise(const G4ParticleDefinition* particle, 102 const G4DataVector& cuts) 103 { 104 //Check that the Livermore Ionisation is NOT attached to e+ 105 if (particle != G4Electron::Electron()) 106 { 107 G4Exception("G4LivermoreIonisationModel::Initialise", 108 "em0002",FatalException, 109 "Livermore Ionisation Model is applicable only to electrons"); 110 } 111 transitionManager->Initialise(); 112 113 //Read energy spectrum 114 if (energySpectrum) 115 { 116 delete energySpectrum; 117 energySpectrum = nullptr; 118 } 119 energySpectrum = new G4eIonisationSpectrum(); 120 if (verboseLevel > 3) 121 G4cout << "G4VEnergySpectrum is initialized" << G4endl; 122 123 //Initialize cross section handler 124 if (crossSectionHandler) 125 { 126 delete crossSectionHandler; 127 crossSectionHandler = nullptr; 128 } 129 130 const size_t nbins = 20; 131 G4double emin = LowEnergyLimit(); 132 G4double emax = HighEnergyLimit(); 133 G4int ndec = G4int(std::log10(emax/emin) + 0.5); 134 if(ndec <= 0) { ndec = 1; } 135 136 G4VDataSetAlgorithm* interpolation = new G4SemiLogInterpolation(); 137 crossSectionHandler = 138 new G4eIonisationCrossSectionHandler(energySpectrum,interpolation, 139 emin,emax,nbins*ndec); 140 crossSectionHandler->Clear(); 141 crossSectionHandler->LoadShellData("ioni/ion-ss-cs-"); 142 //This is used to retrieve cross section values later on 143 G4VEMDataSet* emdata = 144 crossSectionHandler->BuildMeanFreePathForMaterials(&cuts); 145 //The method BuildMeanFreePathForMaterials() is required here only to force 146 //the building of an internal table: the output pointer can be deleted 147 delete emdata; 148 149 if (verboseLevel > 0) 150 { 151 G4cout << "Livermore Ionisation model is initialized " << G4endl 152 << "Energy range: " 153 << LowEnergyLimit() / keV << " keV - " 154 << HighEnergyLimit() / GeV << " GeV" 155 << G4endl; 156 } 157 158 if (verboseLevel > 3) 159 { 160 G4cout << "Cross section data: " << G4endl; 161 crossSectionHandler->PrintData(); 162 G4cout << "Parameters: " << G4endl; 163 energySpectrum->PrintData(); 164 } 165 166 if(isInitialised) { return; } 167 fParticleChange = GetParticleChangeForLoss(); 168 isInitialised = true; 169 } 170 171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 172 173 G4double 174 G4LivermoreIonisationModel::ComputeCrossSectionPerAtom( 175 const G4ParticleDefinition*, 176 G4double energy, 177 G4double Z, G4double, 178 G4double cutEnergy, 179 G4double) 180 { 181 G4int iZ = G4int(Z); 182 if (!crossSectionHandler) 183 { 184 G4Exception("G4LivermoreIonisationModel::ComputeCrossSectionPerAtom", 185 "em1007",FatalException, 186 "The cross section handler is not correctly initialized"); 187 return 0; 188 } 189 190 //The cut is already included in the crossSectionHandler 191 G4double cs = 192 crossSectionHandler->GetCrossSectionAboveThresholdForElement(energy, 193 cutEnergy, 194 iZ); 195 196 if (verboseLevel > 1) 197 { 198 G4cout << "G4LivermoreIonisationModel " << G4endl; 199 G4cout << "Cross section for delta emission > " 200 << cutEnergy/keV << " keV at " 201 << energy/keV << " keV and Z = " << iZ << " --> " 202 << cs/barn << " barn" << G4endl; 203 } 204 return cs; 205 } 206 207 208 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 209 210 G4double 211 G4LivermoreIonisationModel::ComputeDEDXPerVolume(const G4Material* material, 212 const G4ParticleDefinition*, 213 G4double kineticEnergy, 214 G4double cutEnergy) 215 { 216 G4double sPower = 0.0; 217 218 const G4ElementVector* theElementVector = material->GetElementVector(); 219 size_t NumberOfElements = material->GetNumberOfElements() ; 220 const G4double* theAtomicNumDensityVector = 221 material->GetAtomicNumDensityVector(); 222 223 // loop for elements in the material 224 for (size_t iel=0; iel<NumberOfElements; iel++ ) 225 { 226 G4int iZ = (G4int)((*theElementVector)[iel]->GetZ()); 227 G4int nShells = transitionManager->NumberOfShells(iZ); 228 for (G4int n=0; n<nShells; n++) 229 { 230 G4double e = energySpectrum->AverageEnergy(iZ, 0.0,cutEnergy, 231 kineticEnergy, n); 232 G4double cs= crossSectionHandler->FindValue(iZ,kineticEnergy, n); 233 sPower += e * cs * theAtomicNumDensityVector[iel]; 234 } 235 G4double esp = energySpectrum->Excitation(iZ,kineticEnergy); 236 sPower += esp * theAtomicNumDensityVector[iel]; 237 } 238 239 if (verboseLevel > 2) 240 { 241 G4cout << "G4LivermoreIonisationModel " << G4endl; 242 G4cout << "Stopping power < " << cutEnergy/keV 243 << " keV at " << kineticEnergy/keV << " keV = " 244 << sPower/(keV/mm) << " keV/mm" << G4endl; 245 } 246 247 return sPower; 248 } 249 250 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 251 252 void G4LivermoreIonisationModel::SampleSecondaries( 253 std::vector<G4DynamicParticle*>* fvect, 254 const G4MaterialCutsCouple* couple, 255 const G4DynamicParticle* aDynamicParticle, 256 G4double cutE, 257 G4double maxE) 258 { 259 260 G4double kineticEnergy = aDynamicParticle->GetKineticEnergy(); 261 262 if (kineticEnergy <= fIntrinsicLowEnergyLimit) 263 { 264 fParticleChange->SetProposedKineticEnergy(0.); 265 fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy); 266 return; 267 } 268 269 // Select atom and shell 270 G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy); 271 G4int shellIndex = crossSectionHandler->SelectRandomShell(Z, kineticEnergy); 272 const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex); 273 G4double bindingEnergy = shell->BindingEnergy(); 274 275 // Sample delta energy using energy interval for delta-electrons 276 G4double energyMax = 277 std::min(maxE,energySpectrum->MaxEnergyOfSecondaries(kineticEnergy)); 278 G4double energyDelta = energySpectrum->SampleEnergy(Z, cutE, energyMax, 279 kineticEnergy, shellIndex); 280 281 if (energyDelta == 0.) //nothing happens 282 { return; } 283 284 const G4ParticleDefinition* electron = G4Electron::Electron(); 285 G4DynamicParticle* delta = new G4DynamicParticle(electron, 286 GetAngularDistribution()->SampleDirectionForShell(aDynamicParticle, energyDelta, 287 Z, shellIndex, 288 couple->GetMaterial()), 289 energyDelta); 290 291 fvect->push_back(delta); 292 293 // Change kinematics of primary particle 294 G4ThreeVector direction = aDynamicParticle->GetMomentumDirection(); 295 G4double totalMomentum = std::sqrt(kineticEnergy*(kineticEnergy + 2*electron_mass_c2)); 296 297 G4ThreeVector finalP = totalMomentum*direction - delta->GetMomentum(); 298 finalP = finalP.unit(); 299 300 //This is the amount of energy available for fluorescence 301 G4double theEnergyDeposit = bindingEnergy; 302 303 // fill ParticleChange 304 // changed energy and momentum of the actual particle 305 G4double finalKinEnergy = kineticEnergy - energyDelta - theEnergyDeposit; 306 if(finalKinEnergy < 0.0) 307 { 308 theEnergyDeposit += finalKinEnergy; 309 finalKinEnergy = 0.0; 310 } 311 else 312 { 313 fParticleChange->ProposeMomentumDirection(finalP); 314 } 315 fParticleChange->SetProposedKineticEnergy(finalKinEnergy); 316 317 if (theEnergyDeposit < 0) 318 { 319 G4cout << "G4LivermoreIonisationModel: Negative energy deposit: " 320 << theEnergyDeposit/eV << " eV" << G4endl; 321 theEnergyDeposit = 0.0; 322 } 323 324 //Assign local energy deposit 325 fParticleChange->ProposeLocalEnergyDeposit(theEnergyDeposit); 326 327 if (verboseLevel > 1) 328 { 329 G4cout << "-----------------------------------------------------------" << G4endl; 330 G4cout << "Energy balance from G4LivermoreIonisation" << G4endl; 331 G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl; 332 G4cout << "-----------------------------------------------------------" << G4endl; 333 G4cout << "Outgoing primary energy: " << finalKinEnergy/keV << " keV" << G4endl; 334 G4cout << "Delta ray " << energyDelta/keV << " keV" << G4endl; 335 G4cout << "Fluorescence: " << (bindingEnergy-theEnergyDeposit)/keV << " keV" << G4endl; 336 G4cout << "Local energy deposit " << theEnergyDeposit/keV << " keV" << G4endl; 337 G4cout << "Total final state: " << (finalKinEnergy+energyDelta+bindingEnergy) 338 << " keV" << G4endl; 339 G4cout << "-----------------------------------------------------------" << G4endl; 340 } 341 return; 342 } 343 344 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 345 346