G4LivermoreIonisationModel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // Author: Luciano Pandola
00029 //         on base of G4LowEnergyIonisation developed by A.Forti and V.Ivanchenko
00030 //
00031 // History:
00032 // --------
00033 // 12 Jan 2009   L Pandola    Migration from process to model 
00034 // 03 Mar 2009   L Pandola    Bug fix (release memory in the destructor)
00035 // 15 Apr 2009   V Ivanchenko Cleanup initialisation and generation of secondaries:
00036 //                  - apply internal high-energy limit only in constructor 
00037 //                  - do not apply low-energy limit (default is 0)
00038 //                  - simplify sampling of deexcitation by using cut in energy
00039 //                  - set activation of Auger "false"
00040 //                  - remove initialisation of element selectors
00041 // 19 May 2009   L Pandola    Explicitely set to zero pointers deleted in 
00042 //                            Initialise(), since they might be checked later on
00043 // 23 Oct 2009   L Pandola
00044 //                  - atomic deexcitation managed via G4VEmModel::DeexcitationFlag() is
00045 //                    set as "true" (default would be false)
00046 // 12 Oct 2010   L Pandola
00047 //                  - add debugging information about energy in 
00048 //                    SampleDeexcitationAlongStep()
00049 //                  - generate fluorescence SampleDeexcitationAlongStep() only above 
00050 //                    the cuts.
00051 // 01 Jun 2011   V Ivanchenko general cleanup - all old deexcitation code removed
00052 //
00053 
00054 #include "G4LivermoreIonisationModel.hh"
00055 #include "G4PhysicalConstants.hh"
00056 #include "G4SystemOfUnits.hh"
00057 #include "G4ParticleDefinition.hh"
00058 #include "G4MaterialCutsCouple.hh"
00059 #include "G4ProductionCutsTable.hh"
00060 #include "G4DynamicParticle.hh"
00061 #include "G4Element.hh"
00062 #include "G4ParticleChangeForLoss.hh"
00063 #include "G4Electron.hh"
00064 #include "G4CrossSectionHandler.hh"
00065 #include "G4VEMDataSet.hh"
00066 #include "G4eIonisationCrossSectionHandler.hh"
00067 #include "G4eIonisationSpectrum.hh"
00068 #include "G4VEnergySpectrum.hh"
00069 #include "G4SemiLogInterpolation.hh"
00070 #include "G4AtomicTransitionManager.hh"
00071 #include "G4AtomicShell.hh"
00072 
00073 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00074 
00075 
00076 G4LivermoreIonisationModel::G4LivermoreIonisationModel(const G4ParticleDefinition*,
00077                                                        const G4String& nam) : 
00078   G4VEmModel(nam), fParticleChange(0), 
00079   isInitialised(false),
00080   crossSectionHandler(0), energySpectrum(0)
00081 {
00082   fIntrinsicLowEnergyLimit = 10.0*eV;
00083   fIntrinsicHighEnergyLimit = 100.0*GeV;
00084 
00085   verboseLevel = 0;
00086 
00087   transitionManager = G4AtomicTransitionManager::Instance();
00088 }
00089 
00090 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00091 
00092 G4LivermoreIonisationModel::~G4LivermoreIonisationModel()
00093 {
00094   delete energySpectrum;
00095   delete crossSectionHandler;
00096 }
00097 
00098 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00099 
00100 void G4LivermoreIonisationModel::Initialise(const G4ParticleDefinition* particle,
00101                                             const G4DataVector& cuts)
00102 {
00103   //Check that the Livermore Ionisation is NOT attached to e+
00104   if (particle != G4Electron::Electron())
00105     {
00106       G4Exception("G4LivermoreIonisationModel::Initialise",
00107                     "em0002",FatalException,"Livermore Ionisation Model is applicable only to electrons");
00108     }
00109 
00110   //Read energy spectrum
00111   if (energySpectrum) 
00112     {
00113       delete energySpectrum;
00114       energySpectrum = 0;
00115     }
00116   energySpectrum = new G4eIonisationSpectrum();
00117   if (verboseLevel > 3)
00118     G4cout << "G4VEnergySpectrum is initialized" << G4endl;
00119 
00120   //Initialize cross section handler
00121   if (crossSectionHandler) 
00122     {
00123       delete crossSectionHandler;
00124       crossSectionHandler = 0;
00125     }
00126 
00127   const size_t nbins = 20;
00128   G4double emin = LowEnergyLimit();
00129   G4double emax = HighEnergyLimit();
00130   G4int ndec = G4int(std::log10(emax/emin) + 0.5);
00131   if(ndec <= 0) { ndec = 1; }
00132 
00133   G4VDataSetAlgorithm* interpolation = new G4SemiLogInterpolation();
00134   crossSectionHandler = new G4eIonisationCrossSectionHandler(energySpectrum,interpolation,
00135                                                              emin,emax,nbins*ndec);
00136   crossSectionHandler->Clear();
00137   crossSectionHandler->LoadShellData("ioni/ion-ss-cs-");
00138   //This is used to retrieve cross section values later on
00139   G4VEMDataSet* emdata = 
00140     crossSectionHandler->BuildMeanFreePathForMaterials(&cuts);
00141   //The method BuildMeanFreePathForMaterials() is required here only to force 
00142   //the building of an internal table: the output pointer can be deleted
00143   delete emdata;  
00144 
00145   if (verboseLevel > 0)
00146     {
00147       G4cout << "Livermore Ionisation model is initialized " << G4endl
00148              << "Energy range: "
00149              << LowEnergyLimit() / keV << " keV - "
00150              << HighEnergyLimit() / GeV << " GeV"
00151              << G4endl;
00152     }
00153 
00154   if (verboseLevel > 3)
00155     {
00156       G4cout << "Cross section data: " << G4endl; 
00157       crossSectionHandler->PrintData();
00158       G4cout << "Parameters: " << G4endl;
00159       energySpectrum->PrintData();
00160     }
00161 
00162   if(isInitialised) { return; }
00163   fParticleChange = GetParticleChangeForLoss();
00164   isInitialised = true; 
00165 }
00166 
00167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00168 
00169 G4double 
00170 G4LivermoreIonisationModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
00171                                                        G4double energy,
00172                                                        G4double Z, G4double,
00173                                                        G4double cutEnergy, 
00174                                                        G4double)
00175 {
00176   G4int iZ = (G4int) Z;
00177   if (!crossSectionHandler)
00178     {
00179       G4Exception("G4LivermoreIonisationModel::ComputeCrossSectionPerAtom",
00180                     "em1007",FatalException,"The cross section handler is not correctly initialized");
00181       return 0;
00182     }
00183   
00184   //The cut is already included in the crossSectionHandler
00185   G4double cs = 
00186     crossSectionHandler->GetCrossSectionAboveThresholdForElement(energy,cutEnergy,iZ);
00187 
00188   if (verboseLevel > 1)
00189     {
00190       G4cout << "G4LivermoreIonisationModel " << G4endl;
00191       G4cout << "Cross section for delta emission > " << cutEnergy/keV << " keV at " <<
00192         energy/keV << " keV and Z = " << iZ << " --> " << cs/barn << " barn" << G4endl;
00193     }
00194   return cs;
00195 }
00196 
00197 
00198 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00199 
00200 G4double G4LivermoreIonisationModel::ComputeDEDXPerVolume(const G4Material* material,
00201                                                           const G4ParticleDefinition* ,
00202                                                           G4double kineticEnergy,
00203                                                           G4double cutEnergy)
00204 {
00205   G4double sPower = 0.0;
00206 
00207   const G4ElementVector* theElementVector = material->GetElementVector();
00208   size_t NumberOfElements = material->GetNumberOfElements() ;
00209   const G4double* theAtomicNumDensityVector =
00210                     material->GetAtomicNumDensityVector();
00211 
00212   // loop for elements in the material
00213   for (size_t iel=0; iel<NumberOfElements; iel++ ) 
00214     {
00215       G4int iZ = (G4int)((*theElementVector)[iel]->GetZ());
00216       G4int nShells = transitionManager->NumberOfShells(iZ);
00217       for (G4int n=0; n<nShells; n++) 
00218         {
00219           G4double e = energySpectrum->AverageEnergy(iZ, 0.0,cutEnergy,
00220                                                      kineticEnergy, n);
00221           G4double cs= crossSectionHandler->FindValue(iZ,kineticEnergy, n);
00222           sPower   += e * cs * theAtomicNumDensityVector[iel];
00223         }
00224       G4double esp = energySpectrum->Excitation(iZ,kineticEnergy);
00225       sPower   += esp * theAtomicNumDensityVector[iel];
00226     }
00227 
00228   if (verboseLevel > 2)
00229     {
00230       G4cout << "G4LivermoreIonisationModel " << G4endl;
00231       G4cout << "Stopping power < " << cutEnergy/keV << " keV at " << 
00232         kineticEnergy/keV << " keV = " << sPower/(keV/mm) << " keV/mm" << G4endl;
00233     }
00234   
00235   return sPower;
00236 }
00237 
00238 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00239 
00240 void G4LivermoreIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00241                                                    const G4MaterialCutsCouple* couple,
00242                                                    const G4DynamicParticle* aDynamicParticle,
00243                                                    G4double cutE,
00244                                                    G4double maxE)
00245 {
00246   
00247   G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
00248 
00249   if (kineticEnergy <= fIntrinsicLowEnergyLimit)
00250     {
00251       fParticleChange->SetProposedKineticEnergy(0.);
00252       fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy);
00253       return;
00254     }
00255 
00256    // Select atom and shell
00257   G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy);
00258   G4int shellIndex = crossSectionHandler->SelectRandomShell(Z, kineticEnergy);
00259   const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
00260   G4double bindingEnergy = shell->BindingEnergy();
00261 
00262   // Sample delta energy using energy interval for delta-electrons 
00263   G4double energyMax = 
00264     std::min(maxE,energySpectrum->MaxEnergyOfSecondaries(kineticEnergy));
00265   G4double energyDelta = energySpectrum->SampleEnergy(Z, cutE, energyMax,
00266                                                       kineticEnergy, shellIndex);
00267 
00268   if (energyDelta == 0.) //nothing happens
00269     { return; }
00270 
00271   // Transform to shell potential
00272   G4double deltaKinE = energyDelta + 2.0*bindingEnergy;
00273   G4double primaryKinE = kineticEnergy + 2.0*bindingEnergy;
00274 
00275   // sampling of scattering angle neglecting atomic motion
00276   G4double deltaMom = std::sqrt(deltaKinE*(deltaKinE + 2.0*electron_mass_c2));
00277   G4double primaryMom = std::sqrt(primaryKinE*(primaryKinE + 2.0*electron_mass_c2));
00278 
00279   G4double cost = deltaKinE * (primaryKinE + 2.0*electron_mass_c2)
00280                             / (deltaMom * primaryMom);
00281   if (cost > 1.) { cost = 1.; }
00282   G4double sint = std::sqrt((1. - cost)*(1. + cost));
00283   G4double phi  = twopi * G4UniformRand();
00284   G4double dirx = sint * std::cos(phi);
00285   G4double diry = sint * std::sin(phi);
00286   G4double dirz = cost;
00287   
00288   // Rotate to incident electron direction
00289   G4ThreeVector primaryDirection = aDynamicParticle->GetMomentumDirection();
00290   G4ThreeVector deltaDir(dirx,diry,dirz);
00291   deltaDir.rotateUz(primaryDirection);
00292   //Updated components
00293   dirx = deltaDir.x();
00294   diry = deltaDir.y();
00295   dirz = deltaDir.z();
00296 
00297   // Take into account atomic motion del is relative momentum of the motion
00298   // kinetic energy of the motion == bindingEnergy in V.Ivanchenko model
00299   cost = 2.0*G4UniformRand() - 1.0;
00300   sint = std::sqrt(1. - cost*cost);
00301   phi  = twopi * G4UniformRand();
00302   G4double del = std::sqrt(bindingEnergy *(bindingEnergy + 2.0*electron_mass_c2))
00303                / deltaMom;
00304   dirx += del* sint * std::cos(phi);
00305   diry += del* sint * std::sin(phi);
00306   dirz += del* cost;
00307 
00308   // Find out new primary electron direction
00309   G4double finalPx = primaryMom*primaryDirection.x() - deltaMom*dirx;
00310   G4double finalPy = primaryMom*primaryDirection.y() - deltaMom*diry;
00311   G4double finalPz = primaryMom*primaryDirection.z() - deltaMom*dirz;
00312 
00313   //Ok, ready to create the delta ray
00314   G4DynamicParticle* theDeltaRay = new G4DynamicParticle();
00315   theDeltaRay->SetKineticEnergy(energyDelta);
00316   G4double norm = 1.0/std::sqrt(dirx*dirx + diry*diry + dirz*dirz);
00317   dirx *= norm;
00318   diry *= norm;
00319   dirz *= norm;
00320   theDeltaRay->SetMomentumDirection(dirx, diry, dirz);
00321   theDeltaRay->SetDefinition(G4Electron::Electron());
00322   fvect->push_back(theDeltaRay);
00323 
00324   //This is the amount of energy available for fluorescence
00325   G4double theEnergyDeposit = bindingEnergy;
00326 
00327   // fill ParticleChange
00328   // changed energy and momentum of the actual particle
00329   G4double finalKinEnergy = kineticEnergy - energyDelta - theEnergyDeposit;
00330   if(finalKinEnergy < 0.0) 
00331     {
00332       theEnergyDeposit += finalKinEnergy;
00333       finalKinEnergy    = 0.0;
00334     } 
00335   else 
00336     {
00337       G4double normLocal = 1.0/std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
00338       finalPx *= normLocal;
00339       finalPy *= normLocal;
00340       finalPz *= normLocal;
00341       fParticleChange->ProposeMomentumDirection(finalPx, finalPy, finalPz);
00342     }
00343   fParticleChange->SetProposedKineticEnergy(finalKinEnergy);
00344 
00345   if (theEnergyDeposit < 0)
00346     {
00347       G4cout <<  "G4LivermoreIonisationModel: Negative energy deposit: "
00348              << theEnergyDeposit/eV << " eV" << G4endl;
00349       theEnergyDeposit = 0.0;
00350     }
00351 
00352   //Assign local energy deposit
00353   fParticleChange->ProposeLocalEnergyDeposit(theEnergyDeposit);
00354 
00355   if (verboseLevel > 1)
00356     {
00357       G4cout << "-----------------------------------------------------------" << G4endl;
00358       G4cout << "Energy balance from G4LivermoreIonisation" << G4endl;
00359       G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
00360       G4cout << "-----------------------------------------------------------" << G4endl;
00361       G4cout << "Outgoing primary energy: " << finalKinEnergy/keV << " keV" << G4endl;
00362       G4cout << "Delta ray " << energyDelta/keV << " keV" << G4endl;
00363       G4cout << "Fluorescence: " << (bindingEnergy-theEnergyDeposit)/keV << " keV" << G4endl;
00364       G4cout << "Local energy deposit " << theEnergyDeposit/keV << " keV" << G4endl;
00365       G4cout << "Total final state: " << (finalKinEnergy+energyDelta+bindingEnergy)
00366                                           << " keV" << G4endl;
00367       G4cout << "-----------------------------------------------------------" << G4endl;
00368     }
00369   return;
00370 }
00371 
00372 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00373 

Generated on Mon May 27 17:48:48 2013 for Geant4 by  doxygen 1.4.7