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$ >> 27 // GEANT4 tag $Name: $ 26 // 28 // 27 // ------------------------------------------- 29 // -------------------------------------------------------------- 28 // 30 // 29 // File name: G4LowEnergyIonisation 31 // File name: G4LowEnergyIonisation 30 // 32 // 31 // Author: Alessandra Forti, Vladimir I 33 // Author: Alessandra Forti, Vladimir Ivanchenko 32 // 34 // 33 // Creation date: March 1999 35 // Creation date: March 1999 34 // 36 // 35 // Modifications: 37 // Modifications: 36 // - 11.04.2000 VL 38 // - 11.04.2000 VL 37 // Changing use of float and G4float casts t 39 // Changing use of float and G4float casts to G4double casts 38 // because of problems with optimisation (bu 40 // because of problems with optimisation (bug ?) 39 // 10.04.2000 VL 41 // 10.04.2000 VL 40 // - Correcting Fluorescence transition probab 42 // - Correcting Fluorescence transition probabilities in order to take into account 41 // non-radiative transitions. No Auger elect 43 // non-radiative transitions. No Auger electron simulated yet: energy is locally deposited. 42 // 10.04.2000 VL 44 // 10.04.2000 VL 43 // - Correction of incident electron final mom 45 // - Correction of incident electron final momentum direction 44 // 07.04.2000 VL+LU 46 // 07.04.2000 VL+LU 45 // - First implementation of continuous energy 47 // - First implementation of continuous energy loss 46 // 22.03.2000 VL 48 // 22.03.2000 VL 47 // - 1 bug corrected in SelectRandomAtom metho 49 // - 1 bug corrected in SelectRandomAtom method (units) 48 // 17.02.2000 Veronique Lefebure 50 // 17.02.2000 Veronique Lefebure 49 // - 5 bugs corrected: 51 // - 5 bugs corrected: 50 // *in Fluorescence, 2 bugs affecting 52 // *in Fluorescence, 2 bugs affecting 51 // . localEnergyDeposition and 53 // . localEnergyDeposition and 52 // . number of emitted photons that was then 54 // . number of emitted photons that was then always 1 less 53 // *in EnergySampling method: 55 // *in EnergySampling method: 54 // . expon = Parms[13]+1; (instead of uncorr 56 // . expon = Parms[13]+1; (instead of uncorrect -1) 55 // . rejection /= Parms[6];(instead of uncor 57 // . rejection /= Parms[6];(instead of uncorrect Parms[7]) 56 // . Parms[6] is apparently corrupted in the 58 // . Parms[6] is apparently corrupted in the data file (often = 0) 57 // -->Compute normalisation into local var 59 // -->Compute normalisation into local variable rejectionMax 58 // and use rejectionMax in stead of Parms 60 // and use rejectionMax in stead of Parms[6] 59 // 61 // 60 // Added Livermore data table construction met 62 // Added Livermore data table construction methods A. Forti 61 // Modified BuildMeanFreePath to read new data 63 // Modified BuildMeanFreePath to read new data tables A. Forti 62 // Added EnergySampling method A. Forti 64 // Added EnergySampling method A. Forti 63 // Modified PostStepDoIt to insert sampling wi 65 // Modified PostStepDoIt to insert sampling with EEDL data A. Forti 64 // Added SelectRandomAtom A. Forti 66 // Added SelectRandomAtom A. Forti 65 // Added map of the elements A. Forti 67 // Added map of the elements A. Forti 66 // 20.09.00 V.Ivanchenko update fluctuations 68 // 20.09.00 V.Ivanchenko update fluctuations 67 // 24.04.01 V.Ivanchenko remove RogueWave 69 // 24.04.01 V.Ivanchenko remove RogueWave 68 // 22.05.01 V.Ivanchenko update calculation of 70 // 22.05.01 V.Ivanchenko update calculation of delta-ray kinematic + 69 // clean up the code 71 // clean up the code 70 // 02.08.01 V.Ivanchenko fix energy conservati 72 // 02.08.01 V.Ivanchenko fix energy conservation for small steps 71 // 18.08.01 V.Ivanchenko fix energy conservati 73 // 18.08.01 V.Ivanchenko fix energy conservation for pathalogical delta-energy 72 // 01.10.01 E. Guardincerri Replaced fluoresce 74 // 01.10.01 E. Guardincerri Replaced fluorescence generation in PostStepDoIt 73 // according to desig 75 // according to design iteration 74 // 04.10.01 MGP Minor clean-up in 76 // 04.10.01 MGP Minor clean-up in the fluo section, removal of 75 // compilation warnin 77 // compilation warnings and extra protection to 76 // prevent from acces 78 // prevent from accessing a null pointer 77 // 29.09.01 V.Ivanchenko revision based on 79 // 29.09.01 V.Ivanchenko revision based on design iteration 78 // 10.10.01 MGP Revision to improv 80 // 10.10.01 MGP Revision to improve code quality and 79 // consistency with d 81 // consistency with design 80 // 18.10.01 V.Ivanchenko Add fluorescence A 82 // 18.10.01 V.Ivanchenko Add fluorescence AlongStepDoIt 81 // 18.10.01 MGP Revision to improv 83 // 18.10.01 MGP Revision to improve code quality and 82 // consistency with d 84 // consistency with design 83 // 19.10.01 V.Ivanchenko update according t 85 // 19.10.01 V.Ivanchenko update according to new design, V.Ivanchenko 84 // 26.10.01 V.Ivanchenko clean up deexcitat 86 // 26.10.01 V.Ivanchenko clean up deexcitation 85 // 28.10.01 V.Ivanchenko update printout 87 // 28.10.01 V.Ivanchenko update printout 86 // 29.11.01 V.Ivanchenko New parametrisatio 88 // 29.11.01 V.Ivanchenko New parametrisation introduced 87 // 25.03.02 V.Ivanchneko Fix in fluorescenc 89 // 25.03.02 V.Ivanchneko Fix in fluorescence 88 // 28.03.02 V.Ivanchenko Add flag of fluore 90 // 28.03.02 V.Ivanchenko Add flag of fluorescence 89 // 28.05.02 V.Ivanchenko Remove flag fStopA 91 // 28.05.02 V.Ivanchenko Remove flag fStopAndKill 90 // 31.05.02 V.Ivanchenko Add path of Fluo + 92 // 31.05.02 V.Ivanchenko Add path of Fluo + Auger cuts to 91 // AtomicDeexcitation 93 // AtomicDeexcitation 92 // 03.06.02 MGP Restore fStopAndKi 94 // 03.06.02 MGP Restore fStopAndKill 93 // 19.06.02 VI Additional printou 95 // 19.06.02 VI Additional printout 94 // 30.07.02 VI Fix in restricted 96 // 30.07.02 VI Fix in restricted energy loss 95 // 20.09.02 VI Remove ActivateFlu 97 // 20.09.02 VI Remove ActivateFlurescence from SetCut... 96 // 21.01.03 VI Cut per region 98 // 21.01.03 VI Cut per region 97 // 12.02.03 VI Change signature f 99 // 12.02.03 VI Change signature for Deexcitation 98 // 12.04.03 V.Ivanchenko Cut per region for 100 // 12.04.03 V.Ivanchenko Cut per region for fluo AlongStep 99 // 31.08.04 V.Ivanchenko Add density correc 101 // 31.08.04 V.Ivanchenko Add density correction 100 // 102 // 101 // ------------------------------------------- 103 // -------------------------------------------------------------- 102 104 103 #include "G4LowEnergyIonisation.hh" 105 #include "G4LowEnergyIonisation.hh" 104 #include "G4PhysicalConstants.hh" 106 #include "G4PhysicalConstants.hh" 105 #include "G4SystemOfUnits.hh" 107 #include "G4SystemOfUnits.hh" 106 #include "G4RDeIonisationSpectrum.hh" 108 #include "G4RDeIonisationSpectrum.hh" 107 #include "G4RDeIonisationCrossSectionHandler.h 109 #include "G4RDeIonisationCrossSectionHandler.hh" 108 #include "G4RDAtomicTransitionManager.hh" 110 #include "G4RDAtomicTransitionManager.hh" 109 #include "G4RDAtomicShell.hh" 111 #include "G4RDAtomicShell.hh" 110 #include "G4RDVDataSetAlgorithm.hh" 112 #include "G4RDVDataSetAlgorithm.hh" 111 #include "G4RDSemiLogInterpolation.hh" 113 #include "G4RDSemiLogInterpolation.hh" 112 #include "G4RDLogLogInterpolation.hh" 114 #include "G4RDLogLogInterpolation.hh" 113 #include "G4RDEMDataSet.hh" 115 #include "G4RDEMDataSet.hh" 114 #include "G4RDVEMDataSet.hh" 116 #include "G4RDVEMDataSet.hh" 115 #include "G4RDCompositeEMDataSet.hh" 117 #include "G4RDCompositeEMDataSet.hh" 116 #include "G4EnergyLossTables.hh" 118 #include "G4EnergyLossTables.hh" 117 #include "G4RDShellVacancy.hh" 119 #include "G4RDShellVacancy.hh" 118 #include "G4UnitsTable.hh" 120 #include "G4UnitsTable.hh" 119 #include "G4Electron.hh" 121 #include "G4Electron.hh" 120 #include "G4Gamma.hh" 122 #include "G4Gamma.hh" 121 #include "G4ProductionCutsTable.hh" 123 #include "G4ProductionCutsTable.hh" 122 124 123 G4LowEnergyIonisation::G4LowEnergyIonisation(c 125 G4LowEnergyIonisation::G4LowEnergyIonisation(const G4String& nam) 124 : G4eLowEnergyLoss(nam), 126 : G4eLowEnergyLoss(nam), 125 crossSectionHandler(0), 127 crossSectionHandler(0), 126 theMeanFreePath(0), 128 theMeanFreePath(0), 127 energySpectrum(0), 129 energySpectrum(0), 128 shellVacancy(0) 130 shellVacancy(0) 129 { 131 { 130 cutForPhotons = 250.0*eV; 132 cutForPhotons = 250.0*eV; 131 cutForElectrons = 250.0*eV; 133 cutForElectrons = 250.0*eV; 132 verboseLevel = 0; 134 verboseLevel = 0; 133 } 135 } 134 136 135 137 136 G4LowEnergyIonisation::~G4LowEnergyIonisation( 138 G4LowEnergyIonisation::~G4LowEnergyIonisation() 137 { 139 { 138 delete crossSectionHandler; 140 delete crossSectionHandler; 139 delete energySpectrum; 141 delete energySpectrum; 140 delete theMeanFreePath; 142 delete theMeanFreePath; 141 delete shellVacancy; 143 delete shellVacancy; 142 } 144 } 143 145 144 146 145 void G4LowEnergyIonisation::BuildPhysicsTable( 147 void G4LowEnergyIonisation::BuildPhysicsTable(const G4ParticleDefinition& aParticleType) 146 { 148 { 147 if(verboseLevel > 0) { 149 if(verboseLevel > 0) { 148 G4cout << "G4LowEnergyIonisation::BuildPhy 150 G4cout << "G4LowEnergyIonisation::BuildPhysicsTable start" 149 << G4endl; 151 << G4endl; 150 } 152 } 151 153 152 cutForDelta.clear(); 154 cutForDelta.clear(); 153 155 154 // Create and fill IonisationParameters once 156 // Create and fill IonisationParameters once 155 if( energySpectrum != 0 ) delete energySpect 157 if( energySpectrum != 0 ) delete energySpectrum; 156 energySpectrum = new G4RDeIonisationSpectrum 158 energySpectrum = new G4RDeIonisationSpectrum(); 157 159 158 if(verboseLevel > 0) { 160 if(verboseLevel > 0) { 159 G4cout << "G4RDVEnergySpectrum is initiali 161 G4cout << "G4RDVEnergySpectrum is initialized" 160 << G4endl; 162 << G4endl; 161 } 163 } 162 164 163 // Create and fill G4RDCrossSectionHandler o 165 // Create and fill G4RDCrossSectionHandler once 164 166 165 if ( crossSectionHandler != 0 ) delete cross 167 if ( crossSectionHandler != 0 ) delete crossSectionHandler; 166 G4RDVDataSetAlgorithm* interpolation = new G 168 G4RDVDataSetAlgorithm* interpolation = new G4RDSemiLogInterpolation(); 167 G4double lowKineticEnergy = GetLowerBoundEl 169 G4double lowKineticEnergy = GetLowerBoundEloss(); 168 G4double highKineticEnergy = GetUpperBoundEl 170 G4double highKineticEnergy = GetUpperBoundEloss(); 169 G4int totBin = GetNbinEloss(); 171 G4int totBin = GetNbinEloss(); 170 crossSectionHandler = new G4RDeIonisationCro 172 crossSectionHandler = new G4RDeIonisationCrossSectionHandler(energySpectrum, 171 interpolation, 173 interpolation, 172 lowKineticEnergy, 174 lowKineticEnergy, 173 highKineticEnergy, 175 highKineticEnergy, 174 totBin); 176 totBin); 175 crossSectionHandler->LoadShellData("ioni/ion 177 crossSectionHandler->LoadShellData("ioni/ion-ss-cs-"); 176 178 177 if (verboseLevel > 0) { 179 if (verboseLevel > 0) { 178 G4cout << GetProcessName() 180 G4cout << GetProcessName() 179 << " is created; Cross section data 181 << " is created; Cross section data: " 180 << G4endl; 182 << G4endl; 181 crossSectionHandler->PrintData(); 183 crossSectionHandler->PrintData(); 182 G4cout << "Parameters: " 184 G4cout << "Parameters: " 183 << G4endl; 185 << G4endl; 184 energySpectrum->PrintData(); 186 energySpectrum->PrintData(); 185 } 187 } 186 188 187 // Build loss table for IonisationIV 189 // Build loss table for IonisationIV 188 190 189 BuildLossTable(aParticleType); 191 BuildLossTable(aParticleType); 190 192 191 if(verboseLevel > 0) { 193 if(verboseLevel > 0) { 192 G4cout << "The loss table is built" 194 G4cout << "The loss table is built" 193 << G4endl; 195 << G4endl; 194 } 196 } 195 197 196 if (&aParticleType==G4Electron::Electron()) 198 if (&aParticleType==G4Electron::Electron()) { 197 199 198 RecorderOfElectronProcess[CounterOfElectro 200 RecorderOfElectronProcess[CounterOfElectronProcess] = (*this).theLossTable; 199 CounterOfElectronProcess++; 201 CounterOfElectronProcess++; 200 PrintInfoDefinition(); 202 PrintInfoDefinition(); 201 203 202 } else { 204 } else { 203 205 204 RecorderOfPositronProcess[CounterOfPositro 206 RecorderOfPositronProcess[CounterOfPositronProcess] = (*this).theLossTable; 205 CounterOfPositronProcess++; 207 CounterOfPositronProcess++; 206 } 208 } 207 209 208 // Build mean free path data using cut value 210 // Build mean free path data using cut values 209 211 210 if( theMeanFreePath ) delete theMeanFreePath 212 if( theMeanFreePath ) delete theMeanFreePath; 211 theMeanFreePath = crossSectionHandler-> 213 theMeanFreePath = crossSectionHandler-> 212 BuildMeanFreePathForMateri 214 BuildMeanFreePathForMaterials(&cutForDelta); 213 215 214 if(verboseLevel > 0) { 216 if(verboseLevel > 0) { 215 G4cout << "The MeanFreePath table is built 217 G4cout << "The MeanFreePath table is built" 216 << G4endl; 218 << G4endl; 217 if(verboseLevel > 1) theMeanFreePath->Prin 219 if(verboseLevel > 1) theMeanFreePath->PrintData(); 218 } 220 } 219 221 220 // Build common DEDX table for all ionisatio 222 // Build common DEDX table for all ionisation processes 221 223 222 BuildDEDXTable(aParticleType); 224 BuildDEDXTable(aParticleType); 223 225 224 if (verboseLevel > 0) { 226 if (verboseLevel > 0) { 225 G4cout << "G4LowEnergyIonisation::BuildPhy 227 G4cout << "G4LowEnergyIonisation::BuildPhysicsTable end" 226 << G4endl; 228 << G4endl; 227 } 229 } 228 } 230 } 229 231 230 232 231 void G4LowEnergyIonisation::BuildLossTable(con 233 void G4LowEnergyIonisation::BuildLossTable(const G4ParticleDefinition& ) 232 { 234 { 233 // Build table for energy loss due to soft b 235 // Build table for energy loss due to soft brems 234 // the tables are built for *MATERIALS* binn 236 // the tables are built for *MATERIALS* binning is taken from LowEnergyLoss 235 237 236 G4double lowKineticEnergy = GetLowerBoundEl 238 G4double lowKineticEnergy = GetLowerBoundEloss(); 237 G4double highKineticEnergy = GetUpperBoundEl 239 G4double highKineticEnergy = GetUpperBoundEloss(); 238 size_t totBin = GetNbinEloss(); 240 size_t totBin = GetNbinEloss(); 239 241 240 // create table 242 // create table 241 243 242 if (theLossTable) { 244 if (theLossTable) { 243 theLossTable->clearAndDestroy(); 245 theLossTable->clearAndDestroy(); 244 delete theLossTable; 246 delete theLossTable; 245 } 247 } 246 const G4ProductionCutsTable* theCoupleTable= 248 const G4ProductionCutsTable* theCoupleTable= 247 G4ProductionCutsTable::GetProductionCu 249 G4ProductionCutsTable::GetProductionCutsTable(); 248 size_t numOfCouples = theCoupleTable->GetTab 250 size_t numOfCouples = theCoupleTable->GetTableSize(); 249 theLossTable = new G4PhysicsTable(numOfCoupl 251 theLossTable = new G4PhysicsTable(numOfCouples); 250 252 251 if (shellVacancy != 0) delete shellVacancy; 253 if (shellVacancy != 0) delete shellVacancy; 252 shellVacancy = new G4RDShellVacancy(); 254 shellVacancy = new G4RDShellVacancy(); 253 G4DataVector* ksi = 0; 255 G4DataVector* ksi = 0; 254 G4DataVector* energy = 0; 256 G4DataVector* energy = 0; 255 size_t binForFluo = totBin/10; 257 size_t binForFluo = totBin/10; 256 258 257 G4PhysicsLogVector* bVector = new G4PhysicsL 259 G4PhysicsLogVector* bVector = new G4PhysicsLogVector(lowKineticEnergy, 258 highKineticEner 260 highKineticEnergy, 259 binForFluo); 261 binForFluo); 260 const G4RDAtomicTransitionManager* transitio 262 const G4RDAtomicTransitionManager* transitionManager = G4RDAtomicTransitionManager::Instance(); 261 263 262 // Clean up the vector of cuts 264 // Clean up the vector of cuts 263 265 264 cutForDelta.clear(); 266 cutForDelta.clear(); 265 267 266 // Loop for materials 268 // Loop for materials 267 269 268 for (size_t m=0; m<numOfCouples; m++) { 270 for (size_t m=0; m<numOfCouples; m++) { 269 271 270 // create physics vector and fill it 272 // create physics vector and fill it 271 G4PhysicsLogVector* aVector = new G4Physic 273 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(lowKineticEnergy, 272 highKineticEnergy, 274 highKineticEnergy, 273 totBin); 275 totBin); 274 276 275 // get material parameters needed for the 277 // get material parameters needed for the energy loss calculation 276 const G4MaterialCutsCouple* couple = theCo 278 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(m); 277 const G4Material* material= couple->GetMat 279 const G4Material* material= couple->GetMaterial(); 278 280 279 // the cut cannot be below lowest limit 281 // the cut cannot be below lowest limit 280 G4double tCut = (*(theCoupleTable->GetEner 282 G4double tCut = (*(theCoupleTable->GetEnergyCutsVector(1)))[m]; 281 if(tCut > highKineticEnergy) tCut = highKi 283 if(tCut > highKineticEnergy) tCut = highKineticEnergy; 282 cutForDelta.push_back(tCut); 284 cutForDelta.push_back(tCut); 283 const G4ElementVector* theElementVector = 285 const G4ElementVector* theElementVector = material->GetElementVector(); 284 size_t NumberOfElements = material->GetNum 286 size_t NumberOfElements = material->GetNumberOfElements() ; 285 const G4double* theAtomicNumDensityVector 287 const G4double* theAtomicNumDensityVector = 286 material->GetAtomicNumDens 288 material->GetAtomicNumDensityVector(); 287 if(verboseLevel > 0) { 289 if(verboseLevel > 0) { 288 G4cout << "Energy loss for material # " 290 G4cout << "Energy loss for material # " << m 289 << " tCut(keV)= " << tCut/keV 291 << " tCut(keV)= " << tCut/keV 290 << G4endl; 292 << G4endl; 291 } 293 } 292 294 293 // now comes the loop for the kinetic ener 295 // now comes the loop for the kinetic energy values 294 for (size_t i = 0; i<totBin; i++) { 296 for (size_t i = 0; i<totBin; i++) { 295 297 296 G4double lowEdgeEnergy = aVector->GetLow 298 G4double lowEdgeEnergy = aVector->GetLowEdgeEnergy(i); 297 G4double ionloss = 0.; 299 G4double ionloss = 0.; 298 300 299 // loop for elements in the material 301 // loop for elements in the material 300 for (size_t iel=0; iel<NumberOfElements; 302 for (size_t iel=0; iel<NumberOfElements; iel++ ) { 301 303 302 G4int Z = (G4int)((*theElementVector)[ 304 G4int Z = (G4int)((*theElementVector)[iel]->GetZ()); 303 305 304 G4int nShells = transitionManager->NumberOfS 306 G4int nShells = transitionManager->NumberOfShells(Z); 305 307 306 for (G4int n=0; n<nShells; n++) { 308 for (G4int n=0; n<nShells; n++) { 307 309 308 G4double e = energySpectrum->Average 310 G4double e = energySpectrum->AverageEnergy(Z, 0.0, tCut, 309 311 lowEdgeEnergy, n); 310 G4double cs= crossSectionHandler->Fi 312 G4double cs= crossSectionHandler->FindValue(Z, lowEdgeEnergy, n); 311 ionloss += e * cs * theAtomicNumDe 313 ionloss += e * cs * theAtomicNumDensityVector[iel]; 312 314 313 if(verboseLevel > 1 || (Z == 14 && l 315 if(verboseLevel > 1 || (Z == 14 && lowEdgeEnergy>1. && lowEdgeEnergy<0.)) { 314 G4cout << "Z= " << Z 316 G4cout << "Z= " << Z 315 << " shell= " << n 317 << " shell= " << n 316 << " E(keV)= " << lowEdgeEn 318 << " E(keV)= " << lowEdgeEnergy/keV 317 << " Eav(keV)= " << e/keV 319 << " Eav(keV)= " << e/keV 318 << " cs= " << cs 320 << " cs= " << cs 319 << " loss= " << ionloss 321 << " loss= " << ionloss 320 << " rho= " << theAtomicNumDensit 322 << " rho= " << theAtomicNumDensityVector[iel] 321 << G4endl; 323 << G4endl; 322 } 324 } 323 } 325 } 324 G4double esp = energySpectrum->Excitat 326 G4double esp = energySpectrum->Excitation(Z, lowEdgeEnergy); 325 ionloss += esp * theAtomicNumDensity 327 ionloss += esp * theAtomicNumDensityVector[iel]; 326 328 327 } 329 } 328 if(verboseLevel > 1 || (m == 0 && lowEdg 330 if(verboseLevel > 1 || (m == 0 && lowEdgeEnergy>=1. && lowEdgeEnergy<=0.)) { 329 G4cout << "Sum: " 331 G4cout << "Sum: " 330 << " E(keV)= " << lowEdgeEn 332 << " E(keV)= " << lowEdgeEnergy/keV 331 << " loss(MeV/mm)= " << ionloss*m 333 << " loss(MeV/mm)= " << ionloss*mm/MeV 332 << G4endl; 334 << G4endl; 333 } 335 } 334 aVector->PutValue(i,ionloss); 336 aVector->PutValue(i,ionloss); 335 } 337 } 336 theLossTable->insert(aVector); 338 theLossTable->insert(aVector); 337 339 338 // fill data for fluorescence 340 // fill data for fluorescence 339 341 340 G4RDVDataSetAlgorithm* interp = new G4RDLo 342 G4RDVDataSetAlgorithm* interp = new G4RDLogLogInterpolation(); 341 G4RDVEMDataSet* xsis = new G4RDCompositeEM 343 G4RDVEMDataSet* xsis = new G4RDCompositeEMDataSet(interp, 1., 1.); 342 for (size_t iel=0; iel<NumberOfElements; i 344 for (size_t iel=0; iel<NumberOfElements; iel++ ) { 343 345 344 G4int Z = (G4int)((*theElementVector)[ie 346 G4int Z = (G4int)((*theElementVector)[iel]->GetZ()); 345 energy = new G4DataVector(); 347 energy = new G4DataVector(); 346 ksi = new G4DataVector(); 348 ksi = new G4DataVector(); 347 349 348 for (size_t j = 0; j<binForFluo; j++) { 350 for (size_t j = 0; j<binForFluo; j++) { 349 351 350 G4double lowEdgeEnergy = bVector->GetL 352 G4double lowEdgeEnergy = bVector->GetLowEdgeEnergy(j); 351 G4double cross = 0.; 353 G4double cross = 0.; 352 G4double eAverage= 0.; 354 G4double eAverage= 0.; 353 G4int nShells = transitionManager->NumberOfS 355 G4int nShells = transitionManager->NumberOfShells(Z); 354 356 355 for (G4int n=0; n<nShells; n++) { 357 for (G4int n=0; n<nShells; n++) { 356 358 357 G4double e = energySpectrum->Average 359 G4double e = energySpectrum->AverageEnergy(Z, 0.0, tCut, 358 360 lowEdgeEnergy, n); 359 G4double pro = energySpectrum->Proba 361 G4double pro = energySpectrum->Probability(Z, 0.0, tCut, 360 362 lowEdgeEnergy, n); 361 G4double cs= crossSectionHandler->Fi 363 G4double cs= crossSectionHandler->FindValue(Z, lowEdgeEnergy, n); 362 eAverage += e * cs * theAtomicNumD 364 eAverage += e * cs * theAtomicNumDensityVector[iel]; 363 cross += cs * pro * theAtomicNu 365 cross += cs * pro * theAtomicNumDensityVector[iel]; 364 if(verboseLevel > 1) { 366 if(verboseLevel > 1) { 365 G4cout << "Z= " << Z 367 G4cout << "Z= " << Z 366 << " shell= " << n 368 << " shell= " << n 367 << " E(keV)= " << lowEdgeEn 369 << " E(keV)= " << lowEdgeEnergy/keV 368 << " Eav(keV)= " << e/keV 370 << " Eav(keV)= " << e/keV 369 << " pro= " << pro 371 << " pro= " << pro 370 << " cs= " << cs 372 << " cs= " << cs 371 << G4endl; 373 << G4endl; 372 } 374 } 373 } 375 } 374 376 375 G4double coeff = 0.0; 377 G4double coeff = 0.0; 376 if(eAverage > 0.) { 378 if(eAverage > 0.) { 377 coeff = cross/eAverage; 379 coeff = cross/eAverage; 378 eAverage /= cross; 380 eAverage /= cross; 379 } 381 } 380 382 381 if(verboseLevel > 1) { 383 if(verboseLevel > 1) { 382 G4cout << "Ksi Coefficient for Z= 384 G4cout << "Ksi Coefficient for Z= " << Z 383 << " E(keV)= " << lowEdgeEn 385 << " E(keV)= " << lowEdgeEnergy/keV 384 << " Eav(keV)= " << eAverag 386 << " Eav(keV)= " << eAverage/keV 385 << " coeff= " << coeff 387 << " coeff= " << coeff 386 << G4endl; 388 << G4endl; 387 } 389 } 388 390 389 energy->push_back(lowEdgeEnergy); 391 energy->push_back(lowEdgeEnergy); 390 ksi->push_back(coeff); 392 ksi->push_back(coeff); 391 } 393 } 392 interp = new G4RDLogLogInterpolation(); 394 interp = new G4RDLogLogInterpolation(); 393 G4RDVEMDataSet* set = new G4RDEMDataSet( 395 G4RDVEMDataSet* set = new G4RDEMDataSet(Z,energy,ksi,interp,1.,1.); 394 xsis->AddComponent(set); 396 xsis->AddComponent(set); 395 } 397 } 396 if(verboseLevel) xsis->PrintData(); 398 if(verboseLevel) xsis->PrintData(); 397 shellVacancy->AddXsiTable(xsis); 399 shellVacancy->AddXsiTable(xsis); 398 } 400 } 399 delete bVector; 401 delete bVector; 400 } 402 } 401 403 402 404 403 G4VParticleChange* G4LowEnergyIonisation::Post 405 G4VParticleChange* G4LowEnergyIonisation::PostStepDoIt(const G4Track& track, 404 const G4Step& step) 406 const G4Step& step) 405 { 407 { 406 // Delta electron production mechanism on ba 408 // Delta electron production mechanism on base of the model 407 // J. Stepanek " A program to determine the 409 // J. Stepanek " A program to determine the radiation spectra due 408 // to a single atomic subshell ionisation by 410 // to a single atomic subshell ionisation by a particle or due to 409 // deexcitation or decay of radionuclides", 411 // deexcitation or decay of radionuclides", 410 // Comp. Phys. Comm. 1206 pp 1-19 (1997) 412 // Comp. Phys. Comm. 1206 pp 1-19 (1997) 411 413 412 aParticleChange.Initialize(track); 414 aParticleChange.Initialize(track); 413 415 414 const G4MaterialCutsCouple* couple = track.G 416 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple(); 415 G4double kineticEnergy = track.GetKineticEne 417 G4double kineticEnergy = track.GetKineticEnergy(); 416 418 417 // Select atom and shell 419 // Select atom and shell 418 420 419 G4int Z = crossSectionHandler->SelectRandomA 421 G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy); 420 G4int shell = crossSectionHandler->SelectRan 422 G4int shell = crossSectionHandler->SelectRandomShell(Z, kineticEnergy); 421 const G4RDAtomicShell* atomicShell = 423 const G4RDAtomicShell* atomicShell = 422 (G4RDAtomicTransitionManager:: 424 (G4RDAtomicTransitionManager::Instance())->Shell(Z, shell); 423 G4double bindingEnergy = atomicShell->Bindin 425 G4double bindingEnergy = atomicShell->BindingEnergy(); 424 G4int shellId = atomicShell->ShellId(); 426 G4int shellId = atomicShell->ShellId(); 425 427 426 // Sample delta energy 428 // Sample delta energy 427 429 428 G4int index = couple->GetIndex(); 430 G4int index = couple->GetIndex(); 429 G4double tCut = cutForDelta[index]; 431 G4double tCut = cutForDelta[index]; 430 G4double tmax = energySpectrum->MaxEnergyO 432 G4double tmax = energySpectrum->MaxEnergyOfSecondaries(kineticEnergy); 431 G4double tDelta = energySpectrum->SampleEner 433 G4double tDelta = energySpectrum->SampleEnergy(Z, tCut, tmax, 432 434 kineticEnergy, shell); 433 435 434 if(tDelta == 0.0) 436 if(tDelta == 0.0) 435 return G4VContinuousDiscreteProcess::PostS 437 return G4VContinuousDiscreteProcess::PostStepDoIt(track, step); 436 438 437 // Transform to shell potential 439 // Transform to shell potential 438 G4double deltaKinE = tDelta + 2.0*bindingEne 440 G4double deltaKinE = tDelta + 2.0*bindingEnergy; 439 G4double primaryKinE = kineticEnergy + 2.0*b 441 G4double primaryKinE = kineticEnergy + 2.0*bindingEnergy; 440 442 441 // sampling of scattering angle neglecting a 443 // sampling of scattering angle neglecting atomic motion 442 G4double deltaMom = std::sqrt(deltaKinE*(del 444 G4double deltaMom = std::sqrt(deltaKinE*(deltaKinE + 2.0*electron_mass_c2)); 443 G4double primaryMom = std::sqrt(primaryKinE* 445 G4double primaryMom = std::sqrt(primaryKinE*(primaryKinE + 2.0*electron_mass_c2)); 444 446 445 G4double cost = deltaKinE * (primaryKinE + 2 447 G4double cost = deltaKinE * (primaryKinE + 2.0*electron_mass_c2) 446 / (deltaMom * prim 448 / (deltaMom * primaryMom); 447 449 448 if (cost > 1.) cost = 1.; 450 if (cost > 1.) cost = 1.; 449 G4double sint = std::sqrt(1. - cost*cost); 451 G4double sint = std::sqrt(1. - cost*cost); 450 G4double phi = twopi * G4UniformRand(); 452 G4double phi = twopi * G4UniformRand(); 451 G4double dirx = sint * std::cos(phi); 453 G4double dirx = sint * std::cos(phi); 452 G4double diry = sint * std::sin(phi); 454 G4double diry = sint * std::sin(phi); 453 G4double dirz = cost; 455 G4double dirz = cost; 454 456 455 // Rotate to incident electron direction 457 // Rotate to incident electron direction 456 G4ThreeVector primaryDirection = track.GetMo 458 G4ThreeVector primaryDirection = track.GetMomentumDirection(); 457 G4ThreeVector deltaDir(dirx,diry,dirz); 459 G4ThreeVector deltaDir(dirx,diry,dirz); 458 deltaDir.rotateUz(primaryDirection); 460 deltaDir.rotateUz(primaryDirection); 459 dirx = deltaDir.x(); 461 dirx = deltaDir.x(); 460 diry = deltaDir.y(); 462 diry = deltaDir.y(); 461 dirz = deltaDir.z(); 463 dirz = deltaDir.z(); 462 464 463 465 464 // Take into account atomic motion del is re 466 // Take into account atomic motion del is relative momentum of the motion 465 // kinetic energy of the motion == bindingEn 467 // kinetic energy of the motion == bindingEnergy in V.Ivanchenko model 466 468 467 cost = 2.0*G4UniformRand() - 1.0; 469 cost = 2.0*G4UniformRand() - 1.0; 468 sint = std::sqrt(1. - cost*cost); 470 sint = std::sqrt(1. - cost*cost); 469 phi = twopi * G4UniformRand(); 471 phi = twopi * G4UniformRand(); 470 G4double del = std::sqrt(bindingEnergy *(bin 472 G4double del = std::sqrt(bindingEnergy *(bindingEnergy + 2.0*electron_mass_c2)) 471 / deltaMom; 473 / deltaMom; 472 dirx += del* sint * std::cos(phi); 474 dirx += del* sint * std::cos(phi); 473 diry += del* sint * std::sin(phi); 475 diry += del* sint * std::sin(phi); 474 dirz += del* cost; 476 dirz += del* cost; 475 477 476 // Find out new primary electron direction 478 // Find out new primary electron direction 477 G4double finalPx = primaryMom*primaryDirecti 479 G4double finalPx = primaryMom*primaryDirection.x() - deltaMom*dirx; 478 G4double finalPy = primaryMom*primaryDirecti 480 G4double finalPy = primaryMom*primaryDirection.y() - deltaMom*diry; 479 G4double finalPz = primaryMom*primaryDirecti 481 G4double finalPz = primaryMom*primaryDirection.z() - deltaMom*dirz; 480 482 481 // create G4DynamicParticle object for delta 483 // create G4DynamicParticle object for delta ray 482 G4DynamicParticle* theDeltaRay = new G4Dynam 484 G4DynamicParticle* theDeltaRay = new G4DynamicParticle(); 483 theDeltaRay->SetKineticEnergy(tDelta); 485 theDeltaRay->SetKineticEnergy(tDelta); 484 G4double norm = 1.0/std::sqrt(dirx*dirx + di 486 G4double norm = 1.0/std::sqrt(dirx*dirx + diry*diry + dirz*dirz); 485 dirx *= norm; 487 dirx *= norm; 486 diry *= norm; 488 diry *= norm; 487 dirz *= norm; 489 dirz *= norm; 488 theDeltaRay->SetMomentumDirection(dirx, diry 490 theDeltaRay->SetMomentumDirection(dirx, diry, dirz); 489 theDeltaRay->SetDefinition(G4Electron::Elect 491 theDeltaRay->SetDefinition(G4Electron::Electron()); 490 492 491 G4double theEnergyDeposit = bindingEnergy; 493 G4double theEnergyDeposit = bindingEnergy; 492 494 493 // fill ParticleChange 495 // fill ParticleChange 494 // changed energy and momentum of the actual 496 // changed energy and momentum of the actual particle 495 497 496 G4double finalKinEnergy = kineticEnergy - tD 498 G4double finalKinEnergy = kineticEnergy - tDelta - theEnergyDeposit; 497 if(finalKinEnergy < 0.0) { 499 if(finalKinEnergy < 0.0) { 498 theEnergyDeposit += finalKinEnergy; 500 theEnergyDeposit += finalKinEnergy; 499 finalKinEnergy = 0.0; 501 finalKinEnergy = 0.0; 500 aParticleChange.ProposeTrackStatus(fStopAn 502 aParticleChange.ProposeTrackStatus(fStopAndKill); 501 503 502 } else { 504 } else { 503 505 504 G4double norm = 1.0/std::sqrt(finalPx*fina 506 G4double norm = 1.0/std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz); 505 finalPx *= norm; 507 finalPx *= norm; 506 finalPy *= norm; 508 finalPy *= norm; 507 finalPz *= norm; 509 finalPz *= norm; 508 aParticleChange.ProposeMomentumDirection(f 510 aParticleChange.ProposeMomentumDirection(finalPx, finalPy, finalPz); 509 } 511 } 510 512 511 aParticleChange.ProposeEnergy(finalKinEnergy 513 aParticleChange.ProposeEnergy(finalKinEnergy); 512 514 513 // Generation of Fluorescence and Auger 515 // Generation of Fluorescence and Auger 514 size_t nSecondaries = 0; 516 size_t nSecondaries = 0; 515 size_t totalNumber = 1; 517 size_t totalNumber = 1; 516 std::vector<G4DynamicParticle*>* secondaryVe 518 std::vector<G4DynamicParticle*>* secondaryVector = 0; 517 G4DynamicParticle* aSecondary = 0; 519 G4DynamicParticle* aSecondary = 0; 518 G4ParticleDefinition* type = 0; 520 G4ParticleDefinition* type = 0; 519 521 520 // Fluorescence data start from element 6 522 // Fluorescence data start from element 6 521 523 522 if (Fluorescence() && Z > 5 && (bindingEnerg 524 if (Fluorescence() && Z > 5 && (bindingEnergy >= cutForPhotons 523 || bindingEnergy >= cutForElectro 525 || bindingEnergy >= cutForElectrons)) { 524 526 525 secondaryVector = deexcitationManager.Gene 527 secondaryVector = deexcitationManager.GenerateParticles(Z, shellId); 526 528 527 if (secondaryVector != 0) { 529 if (secondaryVector != 0) { 528 530 529 nSecondaries = secondaryVector->size(); 531 nSecondaries = secondaryVector->size(); 530 for (size_t i = 0; i<nSecondaries; i++) 532 for (size_t i = 0; i<nSecondaries; i++) { 531 533 532 aSecondary = (*secondaryVector)[i]; 534 aSecondary = (*secondaryVector)[i]; 533 if (aSecondary) { 535 if (aSecondary) { 534 536 535 G4double e = aSecondary->GetKineticE 537 G4double e = aSecondary->GetKineticEnergy(); 536 type = aSecondary->GetDefinition(); 538 type = aSecondary->GetDefinition(); 537 if (e < theEnergyDeposit && 539 if (e < theEnergyDeposit && 538 ((type == G4Gamma::Gamma() && 540 ((type == G4Gamma::Gamma() && e > cutForPhotons ) || 539 (type == G4Electron::Electron 541 (type == G4Electron::Electron() && e > cutForElectrons ))) { 540 542 541 theEnergyDeposit -= e; 543 theEnergyDeposit -= e; 542 totalNumber++; 544 totalNumber++; 543 545 544 } else { 546 } else { 545 547 546 delete aSecondary; 548 delete aSecondary; 547 (*secondaryVector)[i] = 0; 549 (*secondaryVector)[i] = 0; 548 } 550 } 549 } 551 } 550 } 552 } 551 } 553 } 552 } 554 } 553 555 554 // Save delta-electrons 556 // Save delta-electrons 555 557 556 aParticleChange.SetNumberOfSecondaries(total 558 aParticleChange.SetNumberOfSecondaries(totalNumber); 557 aParticleChange.AddSecondary(theDeltaRay); 559 aParticleChange.AddSecondary(theDeltaRay); 558 560 559 // Save Fluorescence and Auger 561 // Save Fluorescence and Auger 560 562 561 if (secondaryVector) { 563 if (secondaryVector) { 562 564 563 for (size_t l = 0; l < nSecondaries; l++) 565 for (size_t l = 0; l < nSecondaries; l++) { 564 566 565 aSecondary = (*secondaryVector)[l]; 567 aSecondary = (*secondaryVector)[l]; 566 568 567 if(aSecondary) { 569 if(aSecondary) { 568 570 569 aParticleChange.AddSecondary(aSecondar 571 aParticleChange.AddSecondary(aSecondary); 570 } 572 } 571 } 573 } 572 delete secondaryVector; 574 delete secondaryVector; 573 } 575 } 574 576 575 if(theEnergyDeposit < 0.) { 577 if(theEnergyDeposit < 0.) { 576 G4cout << "G4LowEnergyIonisation: Negative 578 G4cout << "G4LowEnergyIonisation: Negative energy deposit: " 577 << theEnergyDeposit/eV << " eV" << 579 << theEnergyDeposit/eV << " eV" << G4endl; 578 theEnergyDeposit = 0.0; 580 theEnergyDeposit = 0.0; 579 } 581 } 580 aParticleChange.ProposeLocalEnergyDeposit(th 582 aParticleChange.ProposeLocalEnergyDeposit(theEnergyDeposit); 581 583 582 return G4VContinuousDiscreteProcess::PostSte 584 return G4VContinuousDiscreteProcess::PostStepDoIt(track, step); 583 } 585 } 584 586 585 587 586 void G4LowEnergyIonisation::PrintInfoDefinitio 588 void G4LowEnergyIonisation::PrintInfoDefinition() 587 { 589 { 588 G4String comments = "Total cross sections fr 590 G4String comments = "Total cross sections from EEDL database."; 589 comments += "\n Gamma energy sampled fr 591 comments += "\n Gamma energy sampled from a parametrised formula."; 590 comments += "\n Implementation of the c 592 comments += "\n Implementation of the continuous dE/dx part."; 591 comments += "\n At present it can be us 593 comments += "\n At present it can be used for electrons "; 592 comments += "in the energy range [250eV,100G 594 comments += "in the energy range [250eV,100GeV]."; 593 comments += "\n The process must work w 595 comments += "\n The process must work with G4LowEnergyBremsstrahlung."; 594 596 595 G4cout << G4endl << GetProcessName() << ": 597 G4cout << G4endl << GetProcessName() << ": " << comments << G4endl; 596 } 598 } 597 599 598 G4bool G4LowEnergyIonisation::IsApplicable(con 600 G4bool G4LowEnergyIonisation::IsApplicable(const G4ParticleDefinition& particle) 599 { 601 { 600 return ( (&particle == G4Electron::Electron 602 return ( (&particle == G4Electron::Electron()) ); 601 } 603 } 602 604 603 std::vector<G4DynamicParticle*>* 605 std::vector<G4DynamicParticle*>* 604 G4LowEnergyIonisation::DeexciteAtom(const G4Ma 606 G4LowEnergyIonisation::DeexciteAtom(const G4MaterialCutsCouple* couple, 605 G4double incidentEnerg 607 G4double incidentEnergy, 606 G4double eLoss) 608 G4double eLoss) 607 { 609 { 608 // create vector of secondary particles 610 // create vector of secondary particles 609 const G4Material* material = couple->GetMate 611 const G4Material* material = couple->GetMaterial(); 610 612 611 std::vector<G4DynamicParticle*>* partVector 613 std::vector<G4DynamicParticle*>* partVector = 612 new std::vect 614 new std::vector<G4DynamicParticle*>; 613 615 614 if(eLoss > cutForPhotons && eLoss > cutForEl 616 if(eLoss > cutForPhotons && eLoss > cutForElectrons) { 615 617 616 const G4RDAtomicTransitionManager* transit 618 const G4RDAtomicTransitionManager* transitionManager = 617 G4RDAtomicTrans 619 G4RDAtomicTransitionManager::Instance(); 618 620 619 size_t nElements = material->GetNumberOfEl 621 size_t nElements = material->GetNumberOfElements(); 620 const G4ElementVector* theElementVector = 622 const G4ElementVector* theElementVector = material->GetElementVector(); 621 623 622 std::vector<G4DynamicParticle*>* secVector 624 std::vector<G4DynamicParticle*>* secVector = 0; 623 G4DynamicParticle* aSecondary = 0; 625 G4DynamicParticle* aSecondary = 0; 624 G4ParticleDefinition* type = 0; 626 G4ParticleDefinition* type = 0; 625 G4double e; 627 G4double e; 626 G4ThreeVector position; 628 G4ThreeVector position; 627 G4int shell, shellId; 629 G4int shell, shellId; 628 630 629 // sample secondaries 631 // sample secondaries 630 632 631 G4double eTot = 0.0; 633 G4double eTot = 0.0; 632 std::vector<G4int> n = 634 std::vector<G4int> n = 633 shellVacancy->GenerateNumberOfIonis 635 shellVacancy->GenerateNumberOfIonisations(couple, 634 636 incidentEnergy,eLoss); 635 for (size_t i=0; i<nElements; i++) { 637 for (size_t i=0; i<nElements; i++) { 636 638 637 G4int Z = (G4int)((*theElementVector)[i] 639 G4int Z = (G4int)((*theElementVector)[i]->GetZ()); 638 size_t nVacancies = n[i]; 640 size_t nVacancies = n[i]; 639 641 640 G4double maxE = transitionManager->Shell 642 G4double maxE = transitionManager->Shell(Z, 0)->BindingEnergy(); 641 643 642 if (nVacancies && Z > 5 && (maxE>cutForP 644 if (nVacancies && Z > 5 && (maxE>cutForPhotons || maxE>cutForElectrons)) { 643 645 644 for (size_t j=0; j<nVacancies; j++) { 646 for (size_t j=0; j<nVacancies; j++) { 645 647 646 shell = crossSectionHandler->SelectRandomS 648 shell = crossSectionHandler->SelectRandomShell(Z, incidentEnergy); 647 shellId = transitionManager->Shell(Z 649 shellId = transitionManager->Shell(Z, shell)->ShellId(); 648 G4double maxEShell = 650 G4double maxEShell = 649 transitionManager->Shell( 651 transitionManager->Shell(Z, shell)->BindingEnergy(); 650 652 651 if (maxEShell>cutForPhotons || maxES 653 if (maxEShell>cutForPhotons || maxEShell>cutForElectrons ) { 652 654 653 secVector = deexcitationManager.Generate 655 secVector = deexcitationManager.GenerateParticles(Z, shellId); 654 656 655 if (secVector != 0) { 657 if (secVector != 0) { 656 658 657 for (size_t l = 0; l<secVector->size() 659 for (size_t l = 0; l<secVector->size(); l++) { 658 660 659 aSecondary = (*secVector)[l]; 661 aSecondary = (*secVector)[l]; 660 if (aSecondary != 0) { 662 if (aSecondary != 0) { 661 663 662 e = aSecondary->GetKineticEnergy() 664 e = aSecondary->GetKineticEnergy(); 663 type = aSecondary->GetDefinition() 665 type = aSecondary->GetDefinition(); 664 if ( eTot + e <= eLoss && 666 if ( eTot + e <= eLoss && 665 ((type == G4Gamma::Gamma() && e 667 ((type == G4Gamma::Gamma() && e>cutForPhotons ) || 666 (type == G4Electron::Electron() 668 (type == G4Electron::Electron() && e>cutForElectrons))) { 667 669 668 eTot += e; 670 eTot += e; 669 partVector->push_bac 671 partVector->push_back(aSecondary); 670 672 671 } else { 673 } else { 672 674 673 delete aSecondary; 675 delete aSecondary; 674 676 675 } 677 } 676 } 678 } 677 } 679 } 678 delete secVector; 680 delete secVector; 679 } 681 } 680 } 682 } 681 } 683 } 682 } 684 } 683 } 685 } 684 } 686 } 685 return partVector; 687 return partVector; 686 } 688 } 687 689 688 G4double G4LowEnergyIonisation::GetMeanFreePat 690 G4double G4LowEnergyIonisation::GetMeanFreePath(const G4Track& track, 689 G4double , // previousStepSize 691 G4double , // previousStepSize 690 G4ForceCondition* cond) 692 G4ForceCondition* cond) 691 { 693 { 692 *cond = NotForced; 694 *cond = NotForced; 693 G4int index = (track.GetMaterialCutsCouple( 695 G4int index = (track.GetMaterialCutsCouple())->GetIndex(); 694 const G4RDVEMDataSet* data = theMeanFreePat 696 const G4RDVEMDataSet* data = theMeanFreePath->GetComponent(index); 695 G4double meanFreePath = data->FindValue(tra 697 G4double meanFreePath = data->FindValue(track.GetKineticEnergy()); 696 return meanFreePath; 698 return meanFreePath; 697 } 699 } 698 700 699 void G4LowEnergyIonisation::SetCutForLowEnSecP 701 void G4LowEnergyIonisation::SetCutForLowEnSecPhotons(G4double cut) 700 { 702 { 701 cutForPhotons = cut; 703 cutForPhotons = cut; 702 deexcitationManager.SetCutForSecondaryPhoton 704 deexcitationManager.SetCutForSecondaryPhotons(cut); 703 } 705 } 704 706 705 void G4LowEnergyIonisation::SetCutForLowEnSecE 707 void G4LowEnergyIonisation::SetCutForLowEnSecElectrons(G4double cut) 706 { 708 { 707 cutForElectrons = cut; 709 cutForElectrons = cut; 708 deexcitationManager.SetCutForAugerElectrons( 710 deexcitationManager.SetCutForAugerElectrons(cut); 709 } 711 } 710 712 711 void G4LowEnergyIonisation::ActivateAuger(G4bo 713 void G4LowEnergyIonisation::ActivateAuger(G4bool val) 712 { 714 { 713 deexcitationManager.ActivateAugerElectronPro 715 deexcitationManager.ActivateAugerElectronProduction(val); 714 } 716 } 715 717 716 718