Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-31 09:22:06 

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // 
0027 // --------------------------------------------------------------
0028 //
0029 // File name:     G4LowEnergyIonisation
0030 //
0031 // Author:        Alessandra Forti, Vladimir Ivanchenko
0032 // 
0033 // Creation date: March 1999
0034 //
0035 // Modifications:
0036 // - 11.04.2000 VL
0037 //   Changing use of float and G4float casts to G4double casts 
0038 //   because of problems with optimisation (bug ?)
0039 //   10.04.2000 VL
0040 // - Correcting Fluorescence transition probabilities in order to take into account 
0041 //   non-radiative transitions. No Auger electron simulated yet: energy is locally deposited.
0042 //   10.04.2000 VL
0043 // - Correction of incident electron final momentum direction
0044 //   07.04.2000 VL+LU
0045 // - First implementation of continuous energy loss
0046 //   22.03.2000 VL
0047 // - 1 bug corrected in SelectRandomAtom method (units)
0048 //   17.02.2000 Veronique Lefebure
0049 // - 5 bugs corrected: 
0050 //   *in Fluorescence, 2 bugs affecting 
0051 //   . localEnergyDeposition and
0052 //   . number of emitted photons that was then always 1 less
0053 //   *in EnergySampling method: 
0054 //   . expon = Parms[13]+1; (instead of uncorrect -1)
0055 //   . rejection /= Parms[6];(instead of uncorrect Parms[7])
0056 //   . Parms[6] is apparently corrupted in the data file (often = 0)  
0057 //     -->Compute normalisation into local variable rejectionMax
0058 //     and use rejectionMax  in stead of Parms[6]
0059 //
0060 // Added Livermore data table construction methods A. Forti
0061 // Modified BuildMeanFreePath to read new data tables A. Forti
0062 // Added EnergySampling method A. Forti
0063 // Modified PostStepDoIt to insert sampling with EEDL data A. Forti
0064 // Added SelectRandomAtom A. Forti
0065 // Added map of the elements A. Forti
0066 // 20.09.00 V.Ivanchenko update fluctuations 
0067 // 24.04.01 V.Ivanchenko remove RogueWave 
0068 // 22.05.01 V.Ivanchenko update calculation of delta-ray kinematic + 
0069 //                       clean up the code 
0070 // 02.08.01 V.Ivanchenko fix energy conservation for small steps 
0071 // 18.08.01 V.Ivanchenko fix energy conservation for pathalogical delta-energy
0072 // 01.10.01 E. Guardincerri Replaced fluorescence generation in PostStepDoIt
0073 //                          according to design iteration
0074 // 04.10.01 MGP             Minor clean-up in the fluo section, removal of
0075 //                          compilation warnings and extra protection to
0076 //                          prevent from accessing a null pointer        
0077 // 29.09.01 V.Ivanchenko    revision based on design iteration
0078 // 10.10.01 MGP             Revision to improve code quality and 
0079 //                          consistency with design
0080 // 18.10.01 V.Ivanchenko    Add fluorescence AlongStepDoIt
0081 // 18.10.01 MGP             Revision to improve code quality and
0082 //                          consistency with design
0083 // 19.10.01 V.Ivanchenko    update according to new design, V.Ivanchenko
0084 // 26.10.01 V.Ivanchenko    clean up deexcitation
0085 // 28.10.01 V.Ivanchenko    update printout
0086 // 29.11.01 V.Ivanchenko    New parametrisation introduced
0087 // 25.03.02 V.Ivanchneko    Fix in fluorescence
0088 // 28.03.02 V.Ivanchenko    Add flag of fluorescence
0089 // 28.05.02 V.Ivanchenko    Remove flag fStopAndKill
0090 // 31.05.02 V.Ivanchenko    Add path of Fluo + Auger cuts to
0091 //                          AtomicDeexcitation
0092 // 03.06.02 MGP             Restore fStopAndKill
0093 // 19.06.02 VI              Additional printout
0094 // 30.07.02 VI              Fix in restricted energy loss
0095 // 20.09.02 VI              Remove ActivateFlurescence from SetCut...
0096 // 21.01.03 VI              Cut per region
0097 // 12.02.03 VI              Change signature for Deexcitation
0098 // 12.04.03 V.Ivanchenko    Cut per region for fluo AlongStep
0099 // 31.08.04 V.Ivanchenko    Add density correction
0100 //
0101 // --------------------------------------------------------------
0102 
0103 #include "G4LowEnergyIonisation.hh"
0104 #include "G4PhysicalConstants.hh"
0105 #include "G4SystemOfUnits.hh"
0106 #include "G4RDeIonisationSpectrum.hh"
0107 #include "G4RDeIonisationCrossSectionHandler.hh"
0108 #include "G4RDAtomicTransitionManager.hh"
0109 #include "G4RDAtomicShell.hh"
0110 #include "G4RDVDataSetAlgorithm.hh"
0111 #include "G4RDSemiLogInterpolation.hh"
0112 #include "G4RDLogLogInterpolation.hh"
0113 #include "G4RDEMDataSet.hh"
0114 #include "G4RDVEMDataSet.hh"
0115 #include "G4RDCompositeEMDataSet.hh"
0116 #include "G4EnergyLossTables.hh"
0117 #include "G4RDShellVacancy.hh"
0118 #include "G4UnitsTable.hh"
0119 #include "G4Electron.hh"
0120 #include "G4Gamma.hh"
0121 #include "G4ProductionCutsTable.hh"
0122 
0123 G4LowEnergyIonisation::G4LowEnergyIonisation(const G4String& nam)
0124   : G4eLowEnergyLoss(nam), 
0125   crossSectionHandler(0),
0126   theMeanFreePath(0),
0127   energySpectrum(0),
0128   shellVacancy(0)
0129 {
0130   cutForPhotons = 250.0*eV;
0131   cutForElectrons = 250.0*eV;
0132   verboseLevel = 0;
0133 }
0134 
0135 
0136 G4LowEnergyIonisation::~G4LowEnergyIonisation()
0137 {
0138   delete crossSectionHandler;
0139   delete energySpectrum;
0140   delete theMeanFreePath;
0141   delete shellVacancy;
0142 }
0143 
0144 
0145 void G4LowEnergyIonisation::BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
0146 {
0147   if(verboseLevel > 0) {
0148     G4cout << "G4LowEnergyIonisation::BuildPhysicsTable start"
0149            << G4endl;
0150       }
0151 
0152   cutForDelta.clear();
0153 
0154   // Create and fill IonisationParameters once
0155   if( energySpectrum != 0 ) delete energySpectrum;
0156   energySpectrum = new G4RDeIonisationSpectrum();
0157 
0158   if(verboseLevel > 0) {
0159     G4cout << "G4RDVEnergySpectrum is initialized"
0160            << G4endl;
0161       }
0162 
0163   // Create and fill G4RDCrossSectionHandler once
0164 
0165   if ( crossSectionHandler != 0 ) delete crossSectionHandler;
0166   G4RDVDataSetAlgorithm* interpolation = new G4RDSemiLogInterpolation();
0167   G4double lowKineticEnergy  = GetLowerBoundEloss();
0168   G4double highKineticEnergy = GetUpperBoundEloss();
0169   G4int    totBin = GetNbinEloss();
0170   crossSectionHandler = new G4RDeIonisationCrossSectionHandler(energySpectrum,
0171                                  interpolation,
0172                                  lowKineticEnergy,
0173                                  highKineticEnergy,
0174                                  totBin);
0175   crossSectionHandler->LoadShellData("ioni/ion-ss-cs-");
0176 
0177   if (verboseLevel > 0) {
0178     G4cout << GetProcessName()
0179            << " is created; Cross section data: "
0180            << G4endl;
0181     crossSectionHandler->PrintData();
0182     G4cout << "Parameters: "
0183            << G4endl;
0184     energySpectrum->PrintData();
0185   }
0186 
0187   // Build loss table for IonisationIV
0188 
0189   BuildLossTable(aParticleType);
0190 
0191   if(verboseLevel > 0) {
0192     G4cout << "The loss table is built"
0193            << G4endl;
0194       }
0195 
0196   if (&aParticleType==G4Electron::Electron()) {
0197 
0198     RecorderOfElectronProcess[CounterOfElectronProcess] = (*this).theLossTable;
0199     CounterOfElectronProcess++;
0200     PrintInfoDefinition();  
0201 
0202   } else {
0203 
0204     RecorderOfPositronProcess[CounterOfPositronProcess] = (*this).theLossTable;
0205     CounterOfPositronProcess++;
0206   }
0207 
0208   // Build mean free path data using cut values
0209 
0210   if( theMeanFreePath ) delete theMeanFreePath;
0211   theMeanFreePath = crossSectionHandler->
0212                     BuildMeanFreePathForMaterials(&cutForDelta);
0213 
0214   if(verboseLevel > 0) {
0215     G4cout << "The MeanFreePath table is built"
0216            << G4endl;
0217     if(verboseLevel > 1) theMeanFreePath->PrintData();
0218   }
0219 
0220   // Build common DEDX table for all ionisation processes
0221  
0222   BuildDEDXTable(aParticleType);
0223 
0224   if (verboseLevel > 0) {
0225     G4cout << "G4LowEnergyIonisation::BuildPhysicsTable end"
0226            << G4endl;
0227   }
0228 }
0229 
0230 
0231 void G4LowEnergyIonisation::BuildLossTable(const G4ParticleDefinition& )
0232 {
0233   // Build table for energy loss due to soft brems
0234   // the tables are built for *MATERIALS* binning is taken from LowEnergyLoss
0235 
0236   G4double lowKineticEnergy  = GetLowerBoundEloss();
0237   G4double highKineticEnergy = GetUpperBoundEloss();
0238   size_t   totBin = GetNbinEloss();
0239  
0240   //  create table
0241 
0242   if (theLossTable) { 
0243       theLossTable->clearAndDestroy();
0244       delete theLossTable;
0245   }
0246   const G4ProductionCutsTable* theCoupleTable=
0247         G4ProductionCutsTable::GetProductionCutsTable();
0248   size_t numOfCouples = theCoupleTable->GetTableSize();
0249   theLossTable = new G4PhysicsTable(numOfCouples);
0250 
0251   if (shellVacancy != 0) delete shellVacancy;
0252   shellVacancy = new G4RDShellVacancy();
0253   G4DataVector* ksi = 0;
0254   G4DataVector* energy = 0;
0255   size_t binForFluo = totBin/10;
0256 
0257   G4PhysicsLogVector* bVector = new G4PhysicsLogVector(lowKineticEnergy,
0258                                        highKineticEnergy,
0259                                binForFluo);
0260   const G4RDAtomicTransitionManager* transitionManager = G4RDAtomicTransitionManager::Instance();
0261   
0262   // Clean up the vector of cuts
0263 
0264   cutForDelta.clear();
0265 
0266   // Loop for materials
0267 
0268   for (size_t m=0; m<numOfCouples; m++) {
0269 
0270     // create physics vector and fill it
0271     G4PhysicsLogVector* aVector = new G4PhysicsLogVector(lowKineticEnergy,
0272                                  highKineticEnergy,
0273                              totBin);
0274 
0275     // get material parameters needed for the energy loss calculation
0276     const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(m);
0277     const G4Material* material= couple->GetMaterial();
0278 
0279     // the cut cannot be below lowest limit
0280     G4double tCut = (*(theCoupleTable->GetEnergyCutsVector(1)))[m];
0281     if(tCut > highKineticEnergy) tCut = highKineticEnergy;
0282     cutForDelta.push_back(tCut);
0283     const G4ElementVector* theElementVector = material->GetElementVector();
0284     size_t NumberOfElements = material->GetNumberOfElements() ;
0285     const G4double* theAtomicNumDensityVector =
0286                     material->GetAtomicNumDensityVector();
0287     if(verboseLevel > 0) {
0288       G4cout << "Energy loss for material # " << m
0289              << " tCut(keV)= " << tCut/keV
0290              << G4endl;
0291       }
0292 
0293     // now comes the loop for the kinetic energy values
0294     for (size_t i = 0; i<totBin; i++) {
0295 
0296       G4double lowEdgeEnergy = aVector->GetLowEdgeEnergy(i);
0297       G4double ionloss = 0.;
0298 
0299       // loop for elements in the material
0300       for (size_t iel=0; iel<NumberOfElements; iel++ ) {
0301 
0302         G4int Z = (G4int)((*theElementVector)[iel]->GetZ());
0303 
0304     G4int nShells = transitionManager->NumberOfShells(Z);
0305 
0306         for (G4int n=0; n<nShells; n++) {
0307 
0308           G4double e = energySpectrum->AverageEnergy(Z, 0.0, tCut,
0309                                                              lowEdgeEnergy, n);
0310           G4double cs= crossSectionHandler->FindValue(Z, lowEdgeEnergy, n);
0311           ionloss   += e * cs * theAtomicNumDensityVector[iel];
0312 
0313           if(verboseLevel > 1 || (Z == 14 && lowEdgeEnergy>1. && lowEdgeEnergy<0.)) {
0314             G4cout << "Z= " << Z
0315                    << " shell= " << n
0316                    << " E(keV)= " << lowEdgeEnergy/keV
0317                    << " Eav(keV)= " << e/keV
0318                    << " cs= " << cs
0319                << " loss= " << ionloss
0320                << " rho= " << theAtomicNumDensityVector[iel]
0321                    << G4endl;
0322           }
0323         }
0324         G4double esp = energySpectrum->Excitation(Z, lowEdgeEnergy);
0325         ionloss   += esp * theAtomicNumDensityVector[iel];
0326 
0327       }
0328       if(verboseLevel > 1 || (m == 0 && lowEdgeEnergy>=1. && lowEdgeEnergy<=0.)) {
0329             G4cout << "Sum: "
0330                    << " E(keV)= " << lowEdgeEnergy/keV
0331                << " loss(MeV/mm)= " << ionloss*mm/MeV
0332                    << G4endl;
0333       }
0334       aVector->PutValue(i,ionloss);
0335     }
0336     theLossTable->insert(aVector);
0337 
0338     // fill data for fluorescence
0339 
0340     G4RDVDataSetAlgorithm* interp = new G4RDLogLogInterpolation();
0341     G4RDVEMDataSet* xsis = new G4RDCompositeEMDataSet(interp, 1., 1.);
0342     for (size_t iel=0; iel<NumberOfElements; iel++ ) {
0343 
0344       G4int Z = (G4int)((*theElementVector)[iel]->GetZ());
0345       energy = new G4DataVector();
0346       ksi    = new G4DataVector();
0347 
0348       for (size_t j = 0; j<binForFluo; j++) {
0349 
0350         G4double lowEdgeEnergy = bVector->GetLowEdgeEnergy(j);
0351         G4double cross   = 0.;
0352         G4double eAverage= 0.;
0353     G4int nShells = transitionManager->NumberOfShells(Z);
0354 
0355         for (G4int n=0; n<nShells; n++) {
0356 
0357           G4double e = energySpectrum->AverageEnergy(Z, 0.0, tCut,
0358                                                              lowEdgeEnergy, n);
0359           G4double pro = energySpectrum->Probability(Z, 0.0, tCut,
0360                                                              lowEdgeEnergy, n);
0361           G4double cs= crossSectionHandler->FindValue(Z, lowEdgeEnergy, n);
0362           eAverage   += e * cs * theAtomicNumDensityVector[iel];
0363           cross      += cs * pro * theAtomicNumDensityVector[iel];
0364           if(verboseLevel > 1) {
0365             G4cout << "Z= " << Z
0366                    << " shell= " << n
0367                    << " E(keV)= " << lowEdgeEnergy/keV
0368                    << " Eav(keV)= " << e/keV
0369                    << " pro= " << pro
0370                    << " cs= " << cs
0371                    << G4endl;
0372           }
0373     }
0374 
0375         G4double coeff = 0.0;
0376         if(eAverage > 0.) {
0377           coeff = cross/eAverage;
0378           eAverage /= cross;
0379     }
0380 
0381         if(verboseLevel > 1) {
0382             G4cout << "Ksi Coefficient for Z= " << Z
0383                    << " E(keV)= " << lowEdgeEnergy/keV
0384                    << " Eav(keV)= " << eAverage/keV
0385                    << " coeff= " << coeff
0386                    << G4endl;
0387         }
0388 
0389         energy->push_back(lowEdgeEnergy);
0390         ksi->push_back(coeff);
0391       }
0392       interp = new G4RDLogLogInterpolation();
0393       G4RDVEMDataSet* set = new G4RDEMDataSet(Z,energy,ksi,interp,1.,1.);
0394       xsis->AddComponent(set);
0395     }
0396     if(verboseLevel) xsis->PrintData();
0397     shellVacancy->AddXsiTable(xsis);
0398   }
0399   delete bVector;
0400 }
0401 
0402 
0403 G4VParticleChange* G4LowEnergyIonisation::PostStepDoIt(const G4Track& track,
0404                                    const G4Step&  step)
0405 {
0406   // Delta electron production mechanism on base of the model
0407   // J. Stepanek " A program to determine the radiation spectra due
0408   // to a single atomic subshell ionisation by a particle or due to
0409   // deexcitation or decay of radionuclides",
0410   // Comp. Phys. Comm. 1206 pp 1-19 (1997)
0411 
0412   aParticleChange.Initialize(track);
0413 
0414   const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
0415   G4double kineticEnergy = track.GetKineticEnergy();
0416 
0417   // Select atom and shell
0418 
0419   G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy);
0420   G4int shell = crossSectionHandler->SelectRandomShell(Z, kineticEnergy);
0421   const G4RDAtomicShell* atomicShell =
0422                 (G4RDAtomicTransitionManager::Instance())->Shell(Z, shell);
0423   G4double bindingEnergy = atomicShell->BindingEnergy();
0424   G4int shellId = atomicShell->ShellId();
0425 
0426   // Sample delta energy
0427 
0428   G4int    index  = couple->GetIndex();
0429   G4double tCut   = cutForDelta[index];
0430   G4double tmax   = energySpectrum->MaxEnergyOfSecondaries(kineticEnergy);
0431   G4double tDelta = energySpectrum->SampleEnergy(Z, tCut, tmax,
0432                                                  kineticEnergy, shell);
0433 
0434   if(tDelta == 0.0)
0435     return G4VContinuousDiscreteProcess::PostStepDoIt(track, step);
0436 
0437   // Transform to shell potential
0438   G4double deltaKinE = tDelta + 2.0*bindingEnergy;
0439   G4double primaryKinE = kineticEnergy + 2.0*bindingEnergy;
0440 
0441   // sampling of scattering angle neglecting atomic motion
0442   G4double deltaMom = std::sqrt(deltaKinE*(deltaKinE + 2.0*electron_mass_c2));
0443   G4double primaryMom = std::sqrt(primaryKinE*(primaryKinE + 2.0*electron_mass_c2));
0444 
0445   G4double cost = deltaKinE * (primaryKinE + 2.0*electron_mass_c2)
0446                             / (deltaMom * primaryMom);
0447 
0448   if (cost > 1.) cost = 1.;
0449   G4double sint = std::sqrt(1. - cost*cost);
0450   G4double phi  = twopi * G4UniformRand();
0451   G4double dirx = sint * std::cos(phi);
0452   G4double diry = sint * std::sin(phi);
0453   G4double dirz = cost;
0454 
0455   // Rotate to incident electron direction
0456   G4ThreeVector primaryDirection = track.GetMomentumDirection();
0457   G4ThreeVector deltaDir(dirx,diry,dirz);
0458   deltaDir.rotateUz(primaryDirection);
0459   dirx = deltaDir.x();
0460   diry = deltaDir.y();
0461   dirz = deltaDir.z();
0462 
0463 
0464   // Take into account atomic motion del is relative momentum of the motion
0465   // kinetic energy of the motion == bindingEnergy in V.Ivanchenko model
0466 
0467   cost = 2.0*G4UniformRand() - 1.0;
0468   sint = std::sqrt(1. - cost*cost);
0469   phi  = twopi * G4UniformRand();
0470   G4double del = std::sqrt(bindingEnergy *(bindingEnergy + 2.0*electron_mass_c2))
0471                / deltaMom;
0472   dirx += del* sint * std::cos(phi);
0473   diry += del* sint * std::sin(phi);
0474   dirz += del* cost;
0475 
0476   // Find out new primary electron direction
0477   G4double finalPx = primaryMom*primaryDirection.x() - deltaMom*dirx;
0478   G4double finalPy = primaryMom*primaryDirection.y() - deltaMom*diry;
0479   G4double finalPz = primaryMom*primaryDirection.z() - deltaMom*dirz;
0480 
0481   // create G4DynamicParticle object for delta ray
0482   G4DynamicParticle* theDeltaRay = new G4DynamicParticle();
0483   theDeltaRay->SetKineticEnergy(tDelta);
0484   G4double norm = 1.0/std::sqrt(dirx*dirx + diry*diry + dirz*dirz);
0485   dirx *= norm;
0486   diry *= norm;
0487   dirz *= norm;
0488   theDeltaRay->SetMomentumDirection(dirx, diry, dirz);
0489   theDeltaRay->SetDefinition(G4Electron::Electron());
0490 
0491   G4double theEnergyDeposit = bindingEnergy;
0492 
0493   // fill ParticleChange
0494   // changed energy and momentum of the actual particle
0495 
0496   G4double finalKinEnergy = kineticEnergy - tDelta - theEnergyDeposit;
0497   if(finalKinEnergy < 0.0) {
0498     theEnergyDeposit += finalKinEnergy;
0499     finalKinEnergy    = 0.0;
0500     aParticleChange.ProposeTrackStatus(fStopAndKill);
0501 
0502   } else {
0503 
0504     G4double norm = 1.0/std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
0505     finalPx *= norm;
0506     finalPy *= norm;
0507     finalPz *= norm;
0508     aParticleChange.ProposeMomentumDirection(finalPx, finalPy, finalPz);
0509   }
0510 
0511   aParticleChange.ProposeEnergy(finalKinEnergy);
0512 
0513   // Generation of Fluorescence and Auger
0514   size_t nSecondaries = 0;
0515   size_t totalNumber  = 1;
0516   std::vector<G4DynamicParticle*>* secondaryVector = 0;
0517   G4DynamicParticle* aSecondary = 0;
0518   G4ParticleDefinition* type = 0;
0519 
0520   // Fluorescence data start from element 6
0521 
0522   if (Fluorescence() && Z > 5 && (bindingEnergy >= cutForPhotons
0523             ||  bindingEnergy >= cutForElectrons)) {
0524 
0525     secondaryVector = deexcitationManager.GenerateParticles(Z, shellId);
0526 
0527     if (secondaryVector != 0) {
0528 
0529       nSecondaries = secondaryVector->size();
0530       for (size_t i = 0; i<nSecondaries; i++) {
0531 
0532         aSecondary = (*secondaryVector)[i];
0533         if (aSecondary) {
0534 
0535           G4double e = aSecondary->GetKineticEnergy();
0536           type = aSecondary->GetDefinition();
0537           if (e < theEnergyDeposit &&
0538                 ((type == G4Gamma::Gamma() && e > cutForPhotons ) ||
0539                  (type == G4Electron::Electron() && e > cutForElectrons ))) {
0540 
0541              theEnergyDeposit -= e;
0542              totalNumber++;
0543 
0544       } else {
0545 
0546              delete aSecondary;
0547              (*secondaryVector)[i] = 0;
0548       }
0549     }
0550       }
0551     }
0552   }
0553 
0554   // Save delta-electrons
0555 
0556   aParticleChange.SetNumberOfSecondaries(totalNumber);
0557   aParticleChange.AddSecondary(theDeltaRay);
0558 
0559   // Save Fluorescence and Auger
0560 
0561   if (secondaryVector) {
0562 
0563     for (size_t l = 0; l < nSecondaries; l++) {
0564 
0565       aSecondary = (*secondaryVector)[l];
0566 
0567       if(aSecondary) {
0568 
0569         aParticleChange.AddSecondary(aSecondary);
0570       }
0571     }
0572     delete secondaryVector;
0573   }
0574 
0575   if(theEnergyDeposit < 0.) {
0576     G4cout << "G4LowEnergyIonisation: Negative energy deposit: "
0577            << theEnergyDeposit/eV << " eV" << G4endl;
0578     theEnergyDeposit = 0.0;
0579   }
0580   aParticleChange.ProposeLocalEnergyDeposit(theEnergyDeposit);
0581 
0582   return G4VContinuousDiscreteProcess::PostStepDoIt(track, step);
0583 }
0584 
0585 
0586 void G4LowEnergyIonisation::PrintInfoDefinition()
0587 {
0588   G4String comments = "Total cross sections from EEDL database.";
0589   comments += "\n      Gamma energy sampled from a parametrised formula.";
0590   comments += "\n      Implementation of the continuous dE/dx part.";
0591   comments += "\n      At present it can be used for electrons ";
0592   comments += "in the energy range [250eV,100GeV].";
0593   comments += "\n      The process must work with G4LowEnergyBremsstrahlung.";
0594 
0595   G4cout << G4endl << GetProcessName() << ":  " << comments << G4endl;
0596 }
0597 
0598 G4bool G4LowEnergyIonisation::IsApplicable(const G4ParticleDefinition& particle)
0599 {
0600    return ( (&particle == G4Electron::Electron()) );
0601 }
0602 
0603 std::vector<G4DynamicParticle*>*
0604 G4LowEnergyIonisation::DeexciteAtom(const G4MaterialCutsCouple* couple,
0605                               G4double incidentEnergy,
0606                               G4double eLoss)
0607 {
0608   // create vector of secondary particles
0609   const G4Material* material = couple->GetMaterial();
0610 
0611   std::vector<G4DynamicParticle*>* partVector =
0612                                  new std::vector<G4DynamicParticle*>;
0613 
0614   if(eLoss > cutForPhotons && eLoss > cutForElectrons) {
0615 
0616     const G4RDAtomicTransitionManager* transitionManager =
0617                                G4RDAtomicTransitionManager::Instance();
0618 
0619     size_t nElements = material->GetNumberOfElements();
0620     const G4ElementVector* theElementVector = material->GetElementVector();
0621 
0622     std::vector<G4DynamicParticle*>* secVector = 0;
0623     G4DynamicParticle* aSecondary = 0;
0624     G4ParticleDefinition* type = 0;
0625     G4double e;
0626     G4ThreeVector position;
0627     G4int shell, shellId;
0628 
0629     // sample secondaries
0630 
0631     G4double eTot = 0.0;
0632     std::vector<G4int> n =
0633            shellVacancy->GenerateNumberOfIonisations(couple,
0634                                                      incidentEnergy,eLoss);
0635     for (size_t i=0; i<nElements; i++) {
0636 
0637       G4int Z = (G4int)((*theElementVector)[i]->GetZ());
0638       size_t nVacancies = n[i];
0639 
0640       G4double maxE = transitionManager->Shell(Z, 0)->BindingEnergy();
0641 
0642       if (nVacancies && Z > 5 && (maxE>cutForPhotons || maxE>cutForElectrons)) {
0643 
0644     for (size_t j=0; j<nVacancies; j++) {
0645 
0646       shell = crossSectionHandler->SelectRandomShell(Z, incidentEnergy);
0647           shellId = transitionManager->Shell(Z, shell)->ShellId();
0648       G4double maxEShell =
0649                      transitionManager->Shell(Z, shell)->BindingEnergy();
0650 
0651           if (maxEShell>cutForPhotons || maxEShell>cutForElectrons ) {
0652 
0653         secVector = deexcitationManager.GenerateParticles(Z, shellId);
0654 
0655         if (secVector != 0) {
0656 
0657           for (size_t l = 0; l<secVector->size(); l++) {
0658 
0659             aSecondary = (*secVector)[l];
0660             if (aSecondary != 0) {
0661 
0662               e = aSecondary->GetKineticEnergy();
0663               type = aSecondary->GetDefinition();
0664               if ( eTot + e <= eLoss &&
0665                  ((type == G4Gamma::Gamma() && e>cutForPhotons ) ||
0666                  (type == G4Electron::Electron() && e>cutForElectrons))) {
0667 
0668               eTot += e;
0669                           partVector->push_back(aSecondary);
0670 
0671           } else {
0672 
0673                            delete aSecondary;
0674 
0675               }
0676             }
0677           }
0678               delete secVector;
0679         }
0680       }
0681     }
0682       }
0683     }
0684   }
0685   return partVector;
0686 }
0687 
0688 G4double G4LowEnergyIonisation::GetMeanFreePath(const G4Track& track,
0689                         G4double , // previousStepSize
0690                         G4ForceCondition* cond)
0691 {
0692    *cond = NotForced;
0693    G4int index = (track.GetMaterialCutsCouple())->GetIndex();
0694    const G4RDVEMDataSet* data = theMeanFreePath->GetComponent(index);
0695    G4double meanFreePath = data->FindValue(track.GetKineticEnergy());
0696    return meanFreePath;
0697 }
0698 
0699 void G4LowEnergyIonisation::SetCutForLowEnSecPhotons(G4double cut)
0700 {
0701   cutForPhotons = cut;
0702   deexcitationManager.SetCutForSecondaryPhotons(cut);
0703 }
0704 
0705 void G4LowEnergyIonisation::SetCutForLowEnSecElectrons(G4double cut)
0706 {
0707   cutForElectrons = cut;
0708   deexcitationManager.SetCutForAugerElectrons(cut);
0709 }
0710 
0711 void G4LowEnergyIonisation::ActivateAuger(G4bool val)
0712 {
0713   deexcitationManager.ActivateAugerElectronProduction(val);
0714 }
0715
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!