// // ******************************************************************** // * License and Disclaimer * // * * // * The Geant4 software is copyright of the Copyright Holders of * // * the Geant4 Collaboration. It is provided under the terms and * // * conditions of the Geant4 Software License, included in the file * // * LICENSE and available at http://cern.ch/geant4/license . These * // * include a list of copyright holders. * // * * // * Neither the authors of this software system, nor their employing * // * institutes,nor the agencies providing financial support for this * // * work make any representation or warranty, express or implied, * // * regarding this software system or assume any liability for its * // * use. Please see the license in the file LICENSE and URL above * // * for the full disclaimer and the limitation of liability. * // * * // * This code implementation is the result of the scientific and * // * technical work of the GEANT4 collaboration. * // * By using, copying, modifying or distributing the software (or * // * any work based on the software) you agree to acknowledge its * // * use in resulting scientific publications, and indicate your * // * acceptance of all terms of the Geant4 Software license. * // ******************************************************************** // // $Id: G4LivermoreIonisationModel.cc,v 1.1 2009/02/10 16:15:48 pandola Exp $ // GEANT4 tag $Name: geant4-09-02-ref-02 $ // // Author: Luciano Pandola // // History: // -------- // 12 Jan 2009 L Pandola Migration from process to model // #include "G4LivermoreIonisationModel.hh" #include "G4ParticleDefinition.hh" #include "G4MaterialCutsCouple.hh" #include "G4ProductionCutsTable.hh" #include "G4DynamicParticle.hh" #include "G4Element.hh" #include "G4AtomicTransitionManager.hh" #include "G4AtomicDeexcitation.hh" #include "G4AtomicShell.hh" #include "G4Gamma.hh" #include "G4Electron.hh" #include "G4CrossSectionHandler.hh" #include "G4AtomicDeexcitation.hh" #include "G4ProcessManager.hh" #include "G4VEMDataSet.hh" #include "G4EMDataSet.hh" #include "G4eIonisationCrossSectionHandler.hh" #include "G4eIonisationSpectrum.hh" #include "G4VDataSetAlgorithm.hh" #include "G4SemiLogInterpolation.hh" #include "G4ShellVacancy.hh" #include "G4VDataSetAlgorithm.hh" #include "G4LogLogInterpolation.hh" #include "G4CompositeEMDataSet.hh" //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... G4LivermoreIonisationModel::G4LivermoreIonisationModel(const G4ParticleDefinition*, const G4String& nam) :G4VEmModel(nam),isInitialised(false),crossSectionHandler(0), energySpectrum(0),shellVacancy(0) { fIntrinsicLowEnergyLimit = 10.0*eV; fIntrinsicHighEnergyLimit = 100.0*GeV; fNBinEnergyLoss = 360; SetLowEnergyLimit(fIntrinsicLowEnergyLimit); SetHighEnergyLimit(fIntrinsicHighEnergyLimit); // verboseLevel = 0; // fUseAtomicDeexcitation = true; // // Notice: the fluorescence along step is generated only if it is // set by the PROCESS (e.g. G4eIonisation) via the command // process->ActivateDeexcitation(true); // } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... G4LivermoreIonisationModel::~G4LivermoreIonisationModel() {;} //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... void G4LivermoreIonisationModel::Initialise(const G4ParticleDefinition* particle, const G4DataVector& cuts) { //Check that the Livermore Ionisation is NOT attached to e+ if (particle != G4Electron::Electron()) { G4cout << "ERROR: Livermore Ionisation Model is applicable only to electrons" << G4endl; G4cout << "It cannot be registered to " << particle->GetParticleName() << G4endl; G4Exception(); } // if (LowEnergyLimit() < fIntrinsicLowEnergyLimit) { G4cout << "G4LivermoreIonisationModel: low energy limit increased from " << LowEnergyLimit()/eV << " eV to " << fIntrinsicLowEnergyLimit/eV << " eV" << G4endl; SetLowEnergyLimit(fIntrinsicLowEnergyLimit); } if (HighEnergyLimit() > fIntrinsicHighEnergyLimit) { G4cout << "G4LivermoreIonisationModel: high energy limit decreased from " << HighEnergyLimit()/GeV << " GeV to " << fIntrinsicHighEnergyLimit/GeV << " GeV" << G4endl; SetHighEnergyLimit(fIntrinsicHighEnergyLimit); } //Read energy spectrum if (energySpectrum) delete energySpectrum; energySpectrum = new G4eIonisationSpectrum(); if (verboseLevel > 0) G4cout << "G4VEnergySpectrum is initialized" << G4endl; //Initialize cross section handler if (crossSectionHandler) delete crossSectionHandler; G4VDataSetAlgorithm* interpolation = new G4SemiLogInterpolation(); crossSectionHandler = new G4eIonisationCrossSectionHandler(energySpectrum,interpolation, LowEnergyLimit(),HighEnergyLimit(), fNBinEnergyLoss); crossSectionHandler->Clear(); crossSectionHandler->LoadShellData("ioni/ion-ss-cs-"); //This is used to retrieve cross section values later on crossSectionHandler->BuildMeanFreePathForMaterials(&cuts); //Fluorescence data transitionManager = G4AtomicTransitionManager::Instance(); if (shellVacancy) delete shellVacancy; shellVacancy = new G4ShellVacancy(); InitialiseFluorescence(); //InitialiseElementSelectors(particle,cuts); G4cout << "Livermore Ionisation model is initialized " << G4endl << "Energy range: " << LowEnergyLimit() / keV << " keV - " << HighEnergyLimit() / GeV << " GeV" << G4endl; if (verboseLevel > 1) { G4cout << "Cross section data: " << G4endl; crossSectionHandler->PrintData(); G4cout << "Parameters: " << G4endl; energySpectrum->PrintData(); } if(isInitialised) return; if(pParticleChange) fParticleChange = reinterpret_cast(pParticleChange); else { fParticleChange = new G4ParticleChangeForLoss(); pParticleChange = fParticleChange; } isInitialised = true; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... G4double G4LivermoreIonisationModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*, G4double energy, G4double Z, G4double, G4double cutEnergy, G4double) { G4int iZ = (G4int) Z; if (!crossSectionHandler) { G4cout << "G4LivermoreIonisationModel::ComputeCrossSectionPerAtom" << G4endl; G4cout << "The cross section handler is not correctly initialized" << G4endl; G4Exception(); } //The cut is already included in the crossSectionHandler G4double cs = crossSectionHandler->GetCrossSectionAboveThresholdForElement(energy,cutEnergy,iZ); if (verboseLevel > 1) { G4cout << "G4LivermoreIonisationModel " << G4endl; G4cout << "Cross section for delta emission > " << cutEnergy/keV << " keV at " << energy/keV << " keV and Z = " << iZ << " --> " << cs/barn << " barn" << G4endl; } return cs; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... G4double G4LivermoreIonisationModel::ComputeDEDXPerVolume(const G4Material* material, const G4ParticleDefinition* , G4double kineticEnergy, G4double cutEnergy) { G4double sPower = 0.0; const G4ElementVector* theElementVector = material->GetElementVector(); size_t NumberOfElements = material->GetNumberOfElements() ; const G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector(); // loop for elements in the material for (size_t iel=0; ielGetZ()); G4int nShells = transitionManager->NumberOfShells(iZ); for (G4int n=0; nAverageEnergy(iZ, 0.0,cutEnergy, kineticEnergy, n); G4double cs= crossSectionHandler->FindValue(iZ,kineticEnergy, n); sPower += e * cs * theAtomicNumDensityVector[iel]; } G4double esp = energySpectrum->Excitation(iZ,kineticEnergy); sPower += esp * theAtomicNumDensityVector[iel]; } if (verboseLevel > 2) { G4cout << "G4LivermoreIonisationModel " << G4endl; G4cout << "Stopping power < " << cutEnergy/keV << " keV at " << kineticEnergy/keV << " keV = " << sPower/(keV/mm) << " keV/mm" << G4endl; } return sPower; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... void G4LivermoreIonisationModel::SampleSecondaries(std::vector* fvect, const G4MaterialCutsCouple* couple, const G4DynamicParticle* aDynamicParticle, G4double, G4double) { G4double kineticEnergy = aDynamicParticle->GetKineticEnergy(); const G4ParticleDefinition* theParticle = aDynamicParticle->GetDefinition(); if (kineticEnergy <= LowEnergyLimit()) { fParticleChange->SetProposedKineticEnergy(0.); fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy); //Check if there are AtRest processes if (theParticle->GetProcessManager()->GetAtRestProcessVector()->size()) //In this case there is at least one AtRest process fParticleChange->ProposeTrackStatus(fStopButAlive); else fParticleChange->ProposeTrackStatus(fStopAndKill); return ; } // Select atom and shell G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy); G4int shell = crossSectionHandler->SelectRandomShell(Z, kineticEnergy); const G4AtomicShell* atomicShell = (G4AtomicTransitionManager::Instance())->Shell(Z, shell); G4double bindingEnergy = atomicShell->BindingEnergy(); G4int shellId = atomicShell->ShellId(); G4ProductionCutsTable* theCoupleTable= G4ProductionCutsTable::GetProductionCutsTable(); // Retrieve cuts for delta-rays and for gammas G4double cutForLowEnergySecondaryParticles = 250.0*eV; G4double cutE = (*theCoupleTable->GetEnergyCutsVector(1))[couple->GetIndex()]; G4double cutG = (*theCoupleTable->GetEnergyCutsVector(0))[couple->GetIndex()]; //Production cut for delta-rays (electrons) cutE = std::max(cutForLowEnergySecondaryParticles,cutE); //Production cut for gamma (fluorescence) cutG = std::max(cutForLowEnergySecondaryParticles,cutG); G4double energyCut = cutE; // Sample delta energy G4double energyMax = energySpectrum->MaxEnergyOfSecondaries(kineticEnergy); G4double energyDelta= energySpectrum->SampleEnergy(Z, energyCut,energyMax, kineticEnergy, shell); if (energyDelta == 0.) //nothing happens return; // Transform to shell potential G4double deltaKinE = energyDelta + 2.0*bindingEnergy; G4double primaryKinE = kineticEnergy + 2.0*bindingEnergy; // sampling of scattering angle neglecting atomic motion G4double deltaMom = std::sqrt(deltaKinE*(deltaKinE + 2.0*electron_mass_c2)); G4double primaryMom = std::sqrt(primaryKinE*(primaryKinE + 2.0*electron_mass_c2)); G4double cost = deltaKinE * (primaryKinE + 2.0*electron_mass_c2) / (deltaMom * primaryMom); if (cost > 1.) cost = 1.; G4double sint = std::sqrt(1. - cost*cost); G4double phi = twopi * G4UniformRand(); G4double dirx = sint * std::cos(phi); G4double diry = sint * std::sin(phi); G4double dirz = cost; // Rotate to incident electron direction G4ThreeVector primaryDirection = aDynamicParticle->GetMomentumDirection(); G4ThreeVector deltaDir(dirx,diry,dirz); deltaDir.rotateUz(primaryDirection); //Updated components dirx = deltaDir.x(); diry = deltaDir.y(); dirz = deltaDir.z(); // Take into account atomic motion del is relative momentum of the motion // kinetic energy of the motion == bindingEnergy in V.Ivanchenko model cost = 2.0*G4UniformRand() - 1.0; sint = std::sqrt(1. - cost*cost); phi = twopi * G4UniformRand(); G4double del = std::sqrt(bindingEnergy *(bindingEnergy + 2.0*electron_mass_c2)) / deltaMom; dirx += del* sint * std::cos(phi); diry += del* sint * std::sin(phi); dirz += del* cost; // Find out new primary electron direction G4double finalPx = primaryMom*primaryDirection.x() - deltaMom*dirx; G4double finalPy = primaryMom*primaryDirection.y() - deltaMom*diry; G4double finalPz = primaryMom*primaryDirection.z() - deltaMom*dirz; //Ok, ready to create the delta ray G4DynamicParticle* theDeltaRay = new G4DynamicParticle(); theDeltaRay->SetKineticEnergy(energyDelta); G4double norm = 1.0/std::sqrt(dirx*dirx + diry*diry + dirz*dirz); dirx *= norm; diry *= norm; dirz *= norm; theDeltaRay->SetMomentumDirection(dirx, diry, dirz); theDeltaRay->SetDefinition(G4Electron::Electron()); fvect->push_back(theDeltaRay); //This is the amount of energy available for fluorescence G4double theEnergyDeposit = bindingEnergy; // fill ParticleChange // changed energy and momentum of the actual particle G4double finalKinEnergy = kineticEnergy - energyDelta - theEnergyDeposit; if(finalKinEnergy < 0.0) { theEnergyDeposit += finalKinEnergy; finalKinEnergy = 0.0; if (theParticle->GetProcessManager()->GetAtRestProcessVector()->size()) //In this case there is at least one AtRest process fParticleChange->ProposeTrackStatus(fStopButAlive); else fParticleChange->ProposeTrackStatus(fStopAndKill); } else { G4double norm = 1.0/std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz); finalPx *= norm; finalPy *= norm; finalPz *= norm; fParticleChange->ProposeMomentumDirection(finalPx, finalPy, finalPz); } fParticleChange->SetProposedKineticEnergy(finalKinEnergy); // Generation of Fluorescence and Auger std::vector* secondaryVector = 0; G4DynamicParticle* aSecondary = 0; // Fluorescence data start from element 6 //Notice: in LowEnergyIonisation, fluorescence is always generated above 250 eV //not above the tracking cut. if (fUseAtomicDeexcitation && Z > 5 && bindingEnergy >= cutForLowEnergySecondaryParticles) { secondaryVector = deexcitationManager.GenerateParticles(Z, shellId); if (secondaryVector) { for (size_t i = 0; isize(); i++) { aSecondary = (*secondaryVector)[i]; //Check if it is a valid secondary if (aSecondary) { G4double e = aSecondary->GetKineticEnergy(); G4ParticleDefinition* type = aSecondary->GetDefinition(); if (e < theEnergyDeposit && ((type == G4Gamma::Gamma() && e > cutForLowEnergySecondaryParticles) || (type == G4Electron::Electron() && e > cutForLowEnergySecondaryParticles ))) { theEnergyDeposit -= e; } else { delete aSecondary; (*secondaryVector)[i] = 0; } } } } } G4double totalEnergyInFluorescence = 0.0; //testing purposes // Save Fluorescence and Auger if (secondaryVector) { for (size_t l = 0; l < secondaryVector->size(); l++) { aSecondary = (*secondaryVector)[l]; if(aSecondary) { fvect->push_back(aSecondary); totalEnergyInFluorescence += aSecondary->GetKineticEnergy(); } } delete secondaryVector; } if (theEnergyDeposit < 0) { G4cout << "G4LivermoreIonisationModel: Negative energy deposit: " << theEnergyDeposit/eV << " eV" << G4endl; theEnergyDeposit = 0.0; } //Assign local energy deposit fParticleChange->ProposeLocalEnergyDeposit(theEnergyDeposit); if (verboseLevel > 1) { G4cout << "-----------------------------------------------------------" << G4endl; G4cout << "Energy balance from G4LivermoreIonisation" << G4endl; G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl; G4cout << "-----------------------------------------------------------" << G4endl; G4cout << "Outgoing primary energy: " << finalKinEnergy/keV << " keV" << G4endl; G4cout << "Delta ray " << energyDelta/keV << " keV" << G4endl; G4cout << "Fluorescence: " << totalEnergyInFluorescence/keV << " keV" << G4endl; G4cout << "Local energy deposit " << theEnergyDeposit/keV << " keV" << G4endl; G4cout << "Total final state: " << (finalKinEnergy+energyDelta+totalEnergyInFluorescence+ theEnergyDeposit)/keV << " keV" << G4endl; G4cout << "-----------------------------------------------------------" << G4endl; } if (verboseLevel > 0) { G4double energyDiff = std::fabs(finalKinEnergy+energyDelta+totalEnergyInFluorescence+ theEnergyDeposit-kineticEnergy); if (energyDiff > 0.05*keV) G4cout << "Warning from G4LivermoreIonisation: problem with energy conservation: " << (finalKinEnergy+energyDelta+totalEnergyInFluorescence+theEnergyDeposit)/keV << " keV (final) vs. " << kineticEnergy/keV << " keV (initial)" << G4endl; } return; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... void G4LivermoreIonisationModel::SampleDeexcitationAlongStep(const G4Material* theMaterial, const G4Track& theTrack, G4double& eloss) { //No call if there is no deexcitation along step if (!fUseAtomicDeexcitation) return; //This method gets the energy loss calculated "Along the Step" and //(including fluctuations) and produces explicit fluorescence/Auger //secondaries. The eloss value is updated. G4double energyLossBefore = eloss; if (verboseLevel > 2) G4cout << "Energy loss along step before deexcitation : " << energyLossBefore/keV << " keV" << G4endl; G4double incidentEnergy = theTrack.GetDynamicParticle()->GetKineticEnergy(); const G4MaterialCutsCouple* couple = theTrack.GetMaterialCutsCouple(); //Notice: in LowEnergyIonisation, fluorescence is always generated above 250 eV //not above the tracking cut. G4double cutForLowEnergySecondaryParticles = 250.0*eV; std::vector* deexcitationProducts = new std::vector; if(eloss > cutForLowEnergySecondaryParticles) { const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance(); size_t nElements = theMaterial->GetNumberOfElements(); const G4ElementVector* theElementVector = theMaterial->GetElementVector(); std::vector* secVector = 0; G4DynamicParticle* aSecondary = 0; G4ParticleDefinition* type = 0; G4double e; G4ThreeVector position; G4int shell, shellId; // sample secondaries G4double eTot = 0.0; std::vector n = shellVacancy->GenerateNumberOfIonisations(couple, incidentEnergy,eloss); for (size_t i=0; iGetZ()); size_t nVacancies = n[i]; G4double maxE = transitionManager->Shell(Z, 0)->BindingEnergy(); if (nVacancies && Z > 5 && (maxE> cutForLowEnergySecondaryParticles)) { for (size_t j=0; jSelectRandomShell(Z, incidentEnergy); shellId = transitionManager->Shell(Z, shell)->ShellId(); G4double maxEShell = transitionManager->Shell(Z, shell)->BindingEnergy(); if (maxEShell > cutForLowEnergySecondaryParticles ) { secVector = deexcitationManager.GenerateParticles(Z, shellId); if (secVector) { for (size_t l = 0; lsize(); l++) { aSecondary = (*secVector)[l]; if (aSecondary) { e = aSecondary->GetKineticEnergy(); type = aSecondary->GetDefinition(); if ( eTot + e <= eloss && ((type == G4Gamma::Gamma() && e> cutForLowEnergySecondaryParticles ) || (type == G4Electron::Electron() && e> cutForLowEnergySecondaryParticles))) { eTot += e; deexcitationProducts->push_back(aSecondary); } else delete aSecondary; } } } delete secVector; } } } } } size_t nSecondaries = 0; if (deexcitationProducts) nSecondaries = deexcitationProducts->size(); fParticleChange->SetNumberOfSecondaries(nSecondaries); if (nSecondaries > 0) { const G4StepPoint* preStep = theTrack.GetStep()->GetPreStepPoint(); const G4StepPoint* postStep = theTrack.GetStep()->GetPostStepPoint(); G4ThreeVector r = preStep->GetPosition(); G4ThreeVector deltaR = postStep->GetPosition(); deltaR -= r; G4double t = preStep->GetGlobalTime(); G4double deltaT = postStep->GetGlobalTime(); deltaT -= t; G4double time, q; G4ThreeVector position; for (size_t i=0; iGetKineticEnergy(); eloss -= eSecondary; if (eloss > 0.) { q = G4UniformRand(); time = deltaT*q + t; position = deltaR*q; position += r; G4Track* newTrack = new G4Track(part, time, position); pParticleChange->AddSecondary(newTrack); } else { eloss += eSecondary; delete part; part = 0; } } } } delete deexcitationProducts; if (verboseLevel > 2) G4cout << "Energy loss along step after deexcitation : " << eloss/keV << " keV" << G4endl; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... void G4LivermoreIonisationModel::InitialiseFluorescence() { G4DataVector* ksi = 0; G4DataVector* energyVector = 0; size_t binForFluo = fNBinEnergyLoss/10; G4PhysicsLogVector* eVector = new G4PhysicsLogVector(LowEnergyLimit(),HighEnergyLimit(), binForFluo); const G4ProductionCutsTable* theCoupleTable= G4ProductionCutsTable::GetProductionCutsTable(); size_t numOfCouples = theCoupleTable->GetTableSize(); // Loop on couples for (size_t m=0; mGetMaterialCutsCouple(m); const G4Material* material= couple->GetMaterial(); const G4ElementVector* theElementVector = material->GetElementVector(); size_t NumberOfElements = material->GetNumberOfElements() ; const G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector(); G4VDataSetAlgorithm* interp = new G4LogLogInterpolation(); G4VEMDataSet* xsis = new G4CompositeEMDataSet(interp, 1., 1.); //loop on elements G4double energyCut = (*(theCoupleTable->GetEnergyCutsVector(1)))[m]; for (size_t iel=0; ielGetZ()); energyVector = new G4DataVector(); ksi = new G4DataVector(); //Loop on energy for (size_t j = 0; jGetLowEdgeEnergy(j); G4double cross = 0.; G4double eAverage= 0.; G4int nShells = transitionManager->NumberOfShells(iZ); for (G4int n=0; nAverageEnergy(iZ, 0.0,energyCut, energy, n); G4double pro = energySpectrum->Probability(iZ, 0.0,energyCut, energy, n); G4double cs= crossSectionHandler->FindValue(iZ, energy, n); eAverage += e * cs * theAtomicNumDensityVector[iel]; cross += cs * pro * theAtomicNumDensityVector[iel]; if(verboseLevel > 1) { G4cout << "Z= " << iZ << " shell= " << n << " E(keV)= " << energy/keV << " Eav(keV)= " << e/keV << " pro= " << pro << " cs= " << cs << G4endl; } } G4double coeff = 0.0; if(eAverage > 0.) { coeff = cross/eAverage; eAverage /= cross; } if(verboseLevel > 1) { G4cout << "Ksi Coefficient for Z= " << iZ << " E(keV)= " << energy/keV << " Eav(keV)= " << eAverage/keV << " coeff= " << coeff << G4endl; } energyVector->push_back(energy); ksi->push_back(coeff); } G4VEMDataSet* set = new G4EMDataSet(iZ,energyVector,ksi,interp,1.,1.); xsis->AddComponent(set); } if(verboseLevel) xsis->PrintData(); shellVacancy->AddXsiTable(xsis); } } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... void G4LivermoreIonisationModel::ActivateAuger(G4bool val) { deexcitationManager.ActivateAugerElectronProduction(val); }