Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // 27 // Author: Luciano Pandola 28 // 29 // History: 30 // -------- 31 // 27 Jul 2010 L Pandola First complete implementation 32 // 18 Jan 2011 L.Pandola Stricter check on production of sub-treshold delta-rays. 33 // Should never happen now 34 // 01 Feb 2011 L Pandola Suppress fake energy-violation warning when Auger is active. 35 // Make sure that fluorescence/Auger is generated only if 36 // above threshold 37 // 25 May 2011 L Pandola Renamed (make v2008 as default Penelope) 38 // 26 Jan 2012 L Pandola Migration of AtomicDeexcitation to the new interface 39 // 09 Mar 2012 L Pandola Moved the management and calculation of 40 // cross sections to a separate class. Use a different method to 41 // get normalized shell cross sections 42 // 07 Oct 2013 L. Pandola Migration to MT 43 // 23 Jun 2015 L. Pandola Keep track of the PIXE flag, to avoid double-production of 44 // atomic de-excitation (bug #1761) 45 // 29 Aug 2018 L. Pandola Fix bug causing energy non-conservation 46 // 47 48 #include "G4PenelopeIonisationModel.hh" 49 #include "G4PhysicalConstants.hh" 50 #include "G4SystemOfUnits.hh" 51 #include "G4ParticleDefinition.hh" 52 #include "G4MaterialCutsCouple.hh" 53 #include "G4ProductionCutsTable.hh" 54 #include "G4DynamicParticle.hh" 55 #include "G4AtomicTransitionManager.hh" 56 #include "G4AtomicShell.hh" 57 #include "G4Gamma.hh" 58 #include "G4Electron.hh" 59 #include "G4Positron.hh" 60 #include "G4PenelopeOscillatorManager.hh" 61 #include "G4PenelopeOscillator.hh" 62 #include "G4PenelopeCrossSection.hh" 63 #include "G4PhysicsFreeVector.hh" 64 #include "G4PhysicsLogVector.hh" 65 #include "G4LossTableManager.hh" 66 #include "G4PenelopeIonisationXSHandler.hh" 67 #include "G4EmParameters.hh" 68 #include "G4AutoLock.hh" 69 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 71 namespace { G4Mutex PenelopeIonisationModelMutex = G4MUTEX_INITIALIZER; } 72 73 G4PenelopeIonisationModel::G4PenelopeIonisationModel(const G4ParticleDefinition* part, 74 const G4String& nam) 75 :G4VEmModel(nam),fParticleChange(nullptr),fParticle(nullptr), 76 fCrossSectionHandler(nullptr), 77 fAtomDeexcitation(nullptr), fKineticEnergy1(0.*eV), 78 fCosThetaPrimary(1.0),fEnergySecondary(0.*eV), 79 fCosThetaSecondary(0.0),fTargetOscillator(-1), 80 fIsInitialised(false),fPIXEflag(false),fLocalTable(false) 81 { 82 fIntrinsicLowEnergyLimit = 100.0*eV; 83 fIntrinsicHighEnergyLimit = 100.0*GeV; 84 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit); 85 SetHighEnergyLimit(fIntrinsicHighEnergyLimit); 86 fNBins = 200; 87 88 if (part) 89 SetParticle(part); 90 91 // 92 fOscManager = G4PenelopeOscillatorManager::GetOscillatorManager(); 93 // 94 fVerboseLevel= 0; 95 // Verbosity scale: 96 // 0 = nothing 97 // 1 = warning for energy non-conservation 98 // 2 = details of energy budget 99 // 3 = calculation of cross sections, file openings, sampling of atoms 100 // 4 = entering in methods 101 102 // Atomic deexcitation model activated by default 103 SetDeexcitationFlag(true); 104 } 105 106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 107 108 G4PenelopeIonisationModel::~G4PenelopeIonisationModel() 109 { 110 if (IsMaster() || fLocalTable) 111 { 112 if (fCrossSectionHandler) 113 delete fCrossSectionHandler; 114 } 115 } 116 117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 118 119 void G4PenelopeIonisationModel::Initialise(const G4ParticleDefinition* particle, 120 const G4DataVector& theCuts) 121 { 122 if (fVerboseLevel > 3) 123 G4cout << "Calling G4PenelopeIonisationModel::Initialise()" << G4endl; 124 125 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation(); 126 //Issue warning if the AtomicDeexcitation has not been declared 127 if (!fAtomDeexcitation) 128 { 129 G4cout << G4endl; 130 G4cout << "WARNING from G4PenelopeIonisationModel " << G4endl; 131 G4cout << "Atomic de-excitation module is not instantiated, so there will not be "; 132 G4cout << "any fluorescence/Auger emission." << G4endl; 133 G4cout << "Please make sure this is intended" << G4endl; 134 } 135 136 if (fAtomDeexcitation) 137 fPIXEflag = fAtomDeexcitation->IsPIXEActive(); 138 139 //If the PIXE flag is active, the PIXE interface will take care of the 140 //atomic de-excitation. The model does not need to do that. 141 //Issue warnings here 142 if (fPIXEflag && IsMaster() && particle==G4Electron::Electron()) 143 { 144 G4String theModel = G4EmParameters::Instance()->PIXEElectronCrossSectionModel(); 145 G4cout << "======================================================================" << G4endl; 146 G4cout << "The G4PenelopeIonisationModel is being used with the PIXE flag ON." << G4endl; 147 G4cout << "Atomic de-excitation will be produced statistically by the PIXE " << G4endl; 148 G4cout << "interface by using the shell cross section --> " << theModel << G4endl; 149 G4cout << "The built-in model procedure for atomic de-excitation is disabled. " << G4endl; 150 G4cout << "*Please be sure this is intended*, or disable PIXE by" << G4endl; 151 G4cout << "/process/em/pixe false" << G4endl; 152 G4cout << "======================================================================" << G4endl; 153 } 154 155 SetParticle(particle); 156 157 //Only the master model creates/manages the tables. All workers get the 158 //pointer to the table, and use it as readonly 159 if (IsMaster() && particle == fParticle) 160 { 161 //Set the number of bins for the tables. 20 points per decade 162 fNBins = (std::size_t) (20*std::log10(HighEnergyLimit()/LowEnergyLimit())); 163 fNBins = std::max(fNBins,(std::size_t)100); 164 165 //Clear and re-build the tables 166 if (fCrossSectionHandler) 167 { 168 delete fCrossSectionHandler; 169 fCrossSectionHandler = 0; 170 } 171 fCrossSectionHandler = new G4PenelopeIonisationXSHandler(fNBins); 172 fCrossSectionHandler->SetVerboseLevel(fVerboseLevel); 173 174 //Build tables for all materials 175 G4ProductionCutsTable* theCoupleTable = 176 G4ProductionCutsTable::GetProductionCutsTable(); 177 for (G4int i=0;i<(G4int)theCoupleTable->GetTableSize();++i) 178 { 179 const G4Material* theMat = 180 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial(); 181 //Forces the building of the cross section tables 182 fCrossSectionHandler->BuildXSTable(theMat,theCuts.at(i),particle, 183 IsMaster()); 184 } 185 186 if (fVerboseLevel > 2) { 187 G4cout << "Penelope Ionisation model v2008 is initialized " << G4endl 188 << "Energy range: " 189 << LowEnergyLimit() / keV << " keV - " 190 << HighEnergyLimit() / GeV << " GeV. Using " 191 << fNBins << " bins." 192 << G4endl; 193 } 194 } 195 196 if(fIsInitialised) 197 return; 198 fParticleChange = GetParticleChangeForLoss(); 199 fIsInitialised = true; 200 201 return; 202 } 203 204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 205 206 void G4PenelopeIonisationModel::InitialiseLocal(const G4ParticleDefinition* part, 207 G4VEmModel *masterModel) 208 { 209 if (fVerboseLevel > 3) 210 G4cout << "Calling G4PenelopeIonisationModel::InitialiseLocal()" << G4endl; 211 // 212 //Check that particle matches: one might have multiple master models (e.g. 213 //for e+ and e-). 214 // 215 if (part == fParticle) 216 { 217 //Get the const table pointers from the master to the workers 218 const G4PenelopeIonisationModel* theModel = 219 static_cast<G4PenelopeIonisationModel*> (masterModel); 220 221 //Copy pointers to the data tables 222 fCrossSectionHandler = theModel->fCrossSectionHandler; 223 224 //copy data 225 fNBins = theModel->fNBins; 226 227 //Same verbosity for all workers, as the master 228 fVerboseLevel = theModel->fVerboseLevel; 229 } 230 return; 231 } 232 233 234 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 235 236 G4double G4PenelopeIonisationModel::CrossSectionPerVolume(const G4Material* material, 237 const G4ParticleDefinition* 238 theParticle, 239 G4double energy, 240 G4double cutEnergy, 241 G4double) 242 { 243 // Penelope model v2008 to calculate the cross section for inelastic collisions above the 244 // threshold. It makes use of the Generalised Oscillator Strength (GOS) model from 245 // D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567 246 // 247 // The total cross section is calculated analytically by taking 248 // into account the atomic oscillators coming into the play for a given threshold. 249 // 250 // For incident e- the maximum energy allowed for the delta-rays is energy/2. 251 // because particles are undistinghishable. 252 // 253 // The contribution is splitted in three parts: distant longitudinal collisions, 254 // distant transverse collisions and close collisions. Each term is described by 255 // its own analytical function. 256 // Fermi density correction is calculated analytically according to 257 // U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1 258 // 259 if (fVerboseLevel > 3) 260 G4cout << "Calling CrossSectionPerVolume() of G4PenelopeIonisationModel" << G4endl; 261 262 SetupForMaterial(theParticle, material, energy); 263 264 G4double totalCross = 0.0; 265 G4double crossPerMolecule = 0.; 266 267 //Either Initialize() was not called, or we are in a slave and InitializeLocal() was 268 //not invoked 269 if (!fCrossSectionHandler) 270 { 271 //create a **thread-local** version of the table. Used only for G4EmCalculator and 272 //Unit Tests 273 fLocalTable = true; 274 fCrossSectionHandler = new G4PenelopeIonisationXSHandler(fNBins); 275 } 276 277 const G4PenelopeCrossSection* theXS = 278 fCrossSectionHandler->GetCrossSectionTableForCouple(theParticle, 279 material, 280 cutEnergy); 281 if (!theXS) 282 { 283 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was 284 //not filled up. This can happen in a UnitTest or via G4EmCalculator 285 if (fVerboseLevel > 0) 286 { 287 //Issue a G4Exception (warning) only in verbose mode 288 G4ExceptionDescription ed; 289 ed << "Unable to retrieve the cross section table for " << 290 theParticle->GetParticleName() << 291 " in " << material->GetName() << ", cut = " << cutEnergy/keV << " keV " << G4endl; 292 ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl; 293 G4Exception("G4PenelopeIonisationModel::CrossSectionPerVolume()", 294 "em2038",JustWarning,ed); 295 } 296 //protect file reading via autolock 297 G4AutoLock lock(&PenelopeIonisationModelMutex); 298 fCrossSectionHandler->BuildXSTable(material,cutEnergy,theParticle); 299 lock.unlock(); 300 //now it should be ok 301 theXS = 302 fCrossSectionHandler->GetCrossSectionTableForCouple(theParticle, 303 material, 304 cutEnergy); 305 } 306 307 if (theXS) 308 crossPerMolecule = theXS->GetHardCrossSection(energy); 309 310 G4double atomDensity = material->GetTotNbOfAtomsPerVolume(); 311 G4double atPerMol = fOscManager->GetAtomsPerMolecule(material); 312 313 if (fVerboseLevel > 3) 314 G4cout << "Material " << material->GetName() << " has " << atPerMol << 315 "atoms per molecule" << G4endl; 316 317 G4double moleculeDensity = 0.; 318 if (atPerMol) 319 moleculeDensity = atomDensity/atPerMol; 320 G4double crossPerVolume = crossPerMolecule*moleculeDensity; 321 322 if (fVerboseLevel > 2) 323 { 324 G4cout << "G4PenelopeIonisationModel " << G4endl; 325 G4cout << "Mean free path for delta emission > " << cutEnergy/keV << " keV at " << 326 energy/keV << " keV = " << (1./crossPerVolume)/mm << " mm" << G4endl; 327 if (theXS) 328 totalCross = (theXS->GetTotalCrossSection(energy))*moleculeDensity; 329 G4cout << "Total free path for ionisation (no threshold) at " << 330 energy/keV << " keV = " << (1./totalCross)/mm << " mm" << G4endl; 331 } 332 return crossPerVolume; 333 } 334 335 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 336 337 //This is a dummy method. Never inkoved by the tracking, it just issues 338 //a warning if one tries to get Cross Sections per Atom via the 339 //G4EmCalculator. 340 G4double G4PenelopeIonisationModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*, 341 G4double, 342 G4double, 343 G4double, 344 G4double, 345 G4double) 346 { 347 G4cout << "*** G4PenelopeIonisationModel -- WARNING ***" << G4endl; 348 G4cout << "Penelope Ionisation model v2008 does not calculate cross section _per atom_ " << G4endl; 349 G4cout << "so the result is always zero. For physics values, please invoke " << G4endl; 350 G4cout << "GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" << G4endl; 351 return 0; 352 } 353 354 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 355 356 G4double G4PenelopeIonisationModel::ComputeDEDXPerVolume(const G4Material* material, 357 const G4ParticleDefinition* theParticle, 358 G4double kineticEnergy, 359 G4double cutEnergy) 360 { 361 // Penelope model v2008 to calculate the stopping power for soft inelastic collisions 362 // below the threshold. It makes use of the Generalised Oscillator Strength (GOS) 363 // model from 364 // D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567 365 // 366 // The stopping power is calculated analytically using the dsigma/dW cross 367 // section from the GOS models, which includes separate contributions from 368 // distant longitudinal collisions, distant transverse collisions and 369 // close collisions. Only the atomic oscillators that come in the play 370 // (according to the threshold) are considered for the calculation. 371 // Differential cross sections have a different form for e+ and e-. 372 // 373 // Fermi density correction is calculated analytically according to 374 // U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1 375 376 if (fVerboseLevel > 3) 377 G4cout << "Calling ComputeDEDX() of G4PenelopeIonisationModel" << G4endl; 378 379 //Either Initialize() was not called, or we are in a slave and InitializeLocal() was 380 //not invoked 381 if (!fCrossSectionHandler) 382 { 383 //create a **thread-local** version of the table. Used only for G4EmCalculator and 384 //Unit Tests 385 fLocalTable = true; 386 fCrossSectionHandler = new G4PenelopeIonisationXSHandler(fNBins); 387 } 388 389 const G4PenelopeCrossSection* theXS = 390 fCrossSectionHandler->GetCrossSectionTableForCouple(theParticle,material, 391 cutEnergy); 392 if (!theXS) 393 { 394 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was 395 //not filled up. This can happen in a UnitTest or via G4EmCalculator 396 if (fVerboseLevel > 0) 397 { 398 //Issue a G4Exception (warning) only in verbose mode 399 G4ExceptionDescription ed; 400 ed << "Unable to retrieve the cross section table for " << 401 theParticle->GetParticleName() << 402 " in " << material->GetName() << ", cut = " << cutEnergy/keV << " keV " << G4endl; 403 ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl; 404 G4Exception("G4PenelopeIonisationModel::ComputeDEDXPerVolume()", 405 "em2038",JustWarning,ed); 406 } 407 //protect file reading via autolock 408 G4AutoLock lock(&PenelopeIonisationModelMutex); 409 fCrossSectionHandler->BuildXSTable(material,cutEnergy,theParticle); 410 lock.unlock(); 411 //now it should be ok 412 theXS = 413 fCrossSectionHandler->GetCrossSectionTableForCouple(theParticle, 414 material, 415 cutEnergy); 416 } 417 418 G4double sPowerPerMolecule = 0.0; 419 if (theXS) 420 sPowerPerMolecule = theXS->GetSoftStoppingPower(kineticEnergy); 421 422 G4double atomDensity = material->GetTotNbOfAtomsPerVolume(); 423 G4double atPerMol = fOscManager->GetAtomsPerMolecule(material); 424 425 G4double moleculeDensity = 0.; 426 if (atPerMol) 427 moleculeDensity = atomDensity/atPerMol; 428 G4double sPowerPerVolume = sPowerPerMolecule*moleculeDensity; 429 430 if (fVerboseLevel > 2) 431 { 432 G4cout << "G4PenelopeIonisationModel " << G4endl; 433 G4cout << "Stopping power < " << cutEnergy/keV << " keV at " << 434 kineticEnergy/keV << " keV = " << 435 sPowerPerVolume/(keV/mm) << " keV/mm" << G4endl; 436 } 437 return sPowerPerVolume; 438 } 439 440 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 441 442 G4double G4PenelopeIonisationModel::MinEnergyCut(const G4ParticleDefinition*, 443 const G4MaterialCutsCouple*) 444 { 445 return fIntrinsicLowEnergyLimit; 446 } 447 448 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 449 450 void G4PenelopeIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect, 451 const G4MaterialCutsCouple* couple, 452 const G4DynamicParticle* aDynamicParticle, 453 G4double cutE, G4double) 454 { 455 // Penelope model v2008 to sample the final state following an hard inelastic interaction. 456 // It makes use of the Generalised Oscillator Strength (GOS) model from 457 // D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567 458 // 459 // The GOS model is used to calculate the individual cross sections for all 460 // the atomic oscillators coming in the play, taking into account the three 461 // contributions (distant longitudinal collisions, distant transverse collisions and 462 // close collisions). Then the target shell and the interaction channel are 463 // sampled. Final state of the delta-ray (energy, angle) are generated according 464 // to the analytical distributions (dSigma/dW) for the selected interaction 465 // channels. 466 // For e-, the maximum energy for the delta-ray is initialEnergy/2. (because 467 // particles are indistinghusbable), while it is the full initialEnergy for 468 // e+. 469 // The efficiency of the random sampling algorithm (e.g. for close collisions) 470 // decreases when initial and cutoff energy increase (e.g. 87% for 10-keV primary 471 // and 1 keV threshold, 99% for 10-MeV primary and 10-keV threshold). 472 // Differential cross sections have a different form for e+ and e-. 473 // 474 // WARNING: The model provides an _average_ description of inelastic collisions. 475 // Anyway, the energy spectrum associated to distant excitations of a given 476 // atomic shell is approximated as a single resonance. The simulated energy spectra 477 // show _unphysical_ narrow peaks at energies that are multiple of the shell 478 // resonance energies. The spurious speaks are automatically smoothed out after 479 // multiple inelastic collisions. 480 // 481 // The model determines also the original shell from which the delta-ray is expelled, 482 // in order to produce fluorescence de-excitation (from G4DeexcitationManager) 483 // 484 // Fermi density correction is calculated analytically according to 485 // U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1 486 487 if (fVerboseLevel > 3) 488 G4cout << "Calling SamplingSecondaries() of G4PenelopeIonisationModel" << G4endl; 489 490 G4double kineticEnergy0 = aDynamicParticle->GetKineticEnergy(); 491 const G4ParticleDefinition* theParticle = aDynamicParticle->GetDefinition(); 492 493 if (kineticEnergy0 <= fIntrinsicLowEnergyLimit) 494 { 495 fParticleChange->SetProposedKineticEnergy(0.); 496 fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy0); 497 return ; 498 } 499 500 const G4Material* material = couple->GetMaterial(); 501 const G4PenelopeOscillatorTable* theTable = fOscManager->GetOscillatorTableIonisation(material); 502 503 G4ParticleMomentum particleDirection0 = aDynamicParticle->GetMomentumDirection(); 504 505 //Initialise final-state variables. The proper values will be set by the methods 506 // SampleFinalStateElectron() and SampleFinalStatePositron() 507 fKineticEnergy1=kineticEnergy0; 508 fCosThetaPrimary=1.0; 509 fEnergySecondary=0.0; 510 fCosThetaSecondary=1.0; 511 fTargetOscillator = -1; 512 513 if (theParticle == G4Electron::Electron()) 514 SampleFinalStateElectron(material,cutE,kineticEnergy0); 515 else if (theParticle == G4Positron::Positron()) 516 SampleFinalStatePositron(material,cutE,kineticEnergy0); 517 else 518 { 519 G4ExceptionDescription ed; 520 ed << "Invalid particle " << theParticle->GetParticleName() << G4endl; 521 G4Exception("G4PenelopeIonisationModel::SamplingSecondaries()", 522 "em0001",FatalException,ed); 523 524 } 525 if (fEnergySecondary == 0) return; 526 527 if (fVerboseLevel > 3) 528 { 529 G4cout << "G4PenelopeIonisationModel::SamplingSecondaries() for " << 530 theParticle->GetParticleName() << G4endl; 531 G4cout << "Final eKin = " << fKineticEnergy1 << " keV" << G4endl; 532 G4cout << "Final cosTheta = " << fCosThetaPrimary << G4endl; 533 G4cout << "Delta-ray eKin = " << fEnergySecondary << " keV" << G4endl; 534 G4cout << "Delta-ray cosTheta = " << fCosThetaSecondary << G4endl; 535 G4cout << "Oscillator: " << fTargetOscillator << G4endl; 536 } 537 538 //Update the primary particle 539 G4double sint = std::sqrt(1. - fCosThetaPrimary*fCosThetaPrimary); 540 G4double phiPrimary = twopi * G4UniformRand(); 541 G4double dirx = sint * std::cos(phiPrimary); 542 G4double diry = sint * std::sin(phiPrimary); 543 G4double dirz = fCosThetaPrimary; 544 545 G4ThreeVector electronDirection1(dirx,diry,dirz); 546 electronDirection1.rotateUz(particleDirection0); 547 548 if (fKineticEnergy1 > 0) 549 { 550 fParticleChange->ProposeMomentumDirection(electronDirection1); 551 fParticleChange->SetProposedKineticEnergy(fKineticEnergy1); 552 } 553 else 554 fParticleChange->SetProposedKineticEnergy(0.); 555 556 //Generate the delta ray 557 G4double ionEnergyInPenelopeDatabase = 558 (*theTable)[fTargetOscillator]->GetIonisationEnergy(); 559 560 //Now, try to handle fluorescence 561 //Notice: merged levels are indicated with Z=0 and flag=30 562 G4int shFlag = (*theTable)[fTargetOscillator]->GetShellFlag(); 563 G4int Z = (G4int) (*theTable)[fTargetOscillator]->GetParentZ(); 564 565 //initialize here, then check photons created by Atomic-Deexcitation, and the final state e- 566 const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance(); 567 G4double bindingEnergy = 0.*eV; 568 569 const G4AtomicShell* shell = nullptr; 570 //Real level 571 if (Z > 0 && shFlag<30) 572 { 573 shell = transitionManager->Shell(Z,shFlag-1); 574 bindingEnergy = shell->BindingEnergy(); 575 //shellId = shell->ShellId(); 576 } 577 578 //correct the fEnergySecondary to account for the fact that the Penelope 579 //database of ionisation energies is in general (slightly) different 580 //from the fluorescence database used in Geant4. 581 fEnergySecondary += ionEnergyInPenelopeDatabase-bindingEnergy; 582 583 G4double localEnergyDeposit = bindingEnergy; 584 //testing purposes only 585 G4double energyInFluorescence = 0; 586 G4double energyInAuger = 0; 587 588 if (fEnergySecondary < 0) 589 { 590 //It means that there was some problem/mismatch between the two databases. 591 //In this case, the available energy is ok to excite the level according 592 //to the Penelope database, but not according to the Geant4 database 593 //Full residual energy is deposited locally 594 localEnergyDeposit += fEnergySecondary; 595 fEnergySecondary = 0.0; 596 } 597 598 //Notice: shell might be nullptr (invalid!) if shFlag=30. Must be protected 599 //Disable the built-in de-excitation of the PIXE flag is active. In this 600 //case, the PIXE interface takes care (statistically) of producing the 601 //de-excitation. 602 //Now, take care of fluorescence, if required 603 if (fAtomDeexcitation && !fPIXEflag && shell) 604 { 605 G4int index = couple->GetIndex(); 606 if (fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) 607 { 608 std::size_t nBefore = fvect->size(); 609 fAtomDeexcitation->GenerateParticles(fvect,shell,Z,index); 610 std::size_t nAfter = fvect->size(); 611 612 if (nAfter>nBefore) //actual production of fluorescence 613 { 614 for (std::size_t j=nBefore;j<nAfter;++j) //loop on products 615 { 616 G4double itsEnergy = ((*fvect)[j])->GetKineticEnergy(); 617 if (itsEnergy < localEnergyDeposit) // valid secondary, generate it 618 { 619 localEnergyDeposit -= itsEnergy; 620 if (((*fvect)[j])->GetParticleDefinition() == G4Gamma::Definition()) 621 energyInFluorescence += itsEnergy; 622 else if (((*fvect)[j])->GetParticleDefinition() == G4Electron::Definition()) 623 energyInAuger += itsEnergy; 624 } 625 else //invalid secondary: takes more than the available energy: delete it 626 { 627 delete (*fvect)[j]; 628 (*fvect)[j] = nullptr; 629 } 630 } 631 } 632 } 633 } 634 635 // Generate the delta ray --> to be done only if above cut 636 if (fEnergySecondary > cutE) 637 { 638 G4DynamicParticle* electron = nullptr; 639 G4double sinThetaE = std::sqrt(1.-fCosThetaSecondary*fCosThetaSecondary); 640 G4double phiEl = phiPrimary+pi; //pi with respect to the primary electron/positron 641 G4double xEl = sinThetaE * std::cos(phiEl); 642 G4double yEl = sinThetaE * std::sin(phiEl); 643 G4double zEl = fCosThetaSecondary; 644 G4ThreeVector eDirection(xEl,yEl,zEl); //electron direction 645 eDirection.rotateUz(particleDirection0); 646 electron = new G4DynamicParticle (G4Electron::Electron(), 647 eDirection,fEnergySecondary) ; 648 fvect->push_back(electron); 649 } 650 else 651 { 652 localEnergyDeposit += fEnergySecondary; 653 fEnergySecondary = 0; 654 } 655 656 if (localEnergyDeposit < 0) //Should not be: issue a G4Exception (warning) 657 { 658 G4Exception("G4PenelopeIonisationModel::SampleSecondaries()", 659 "em2099",JustWarning,"WARNING: Negative local energy deposit"); 660 localEnergyDeposit=0.; 661 } 662 fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit); 663 664 if (fVerboseLevel > 1) 665 { 666 G4cout << "-----------------------------------------------------------" << G4endl; 667 G4cout << "Energy balance from G4PenelopeIonisation" << G4endl; 668 G4cout << "Incoming primary energy: " << kineticEnergy0/keV << " keV" << G4endl; 669 G4cout << "-----------------------------------------------------------" << G4endl; 670 G4cout << "Outgoing primary energy: " << fKineticEnergy1/keV << " keV" << G4endl; 671 G4cout << "Delta ray " << fEnergySecondary/keV << " keV" << G4endl; 672 if (energyInFluorescence) 673 G4cout << "Fluorescence x-rays: " << energyInFluorescence/keV << " keV" << G4endl; 674 if (energyInAuger) 675 G4cout << "Auger electrons: " << energyInAuger/keV << " keV" << G4endl; 676 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl; 677 G4cout << "Total final state: " << (fEnergySecondary+energyInFluorescence+fKineticEnergy1+ 678 localEnergyDeposit+energyInAuger)/keV << 679 " keV" << G4endl; 680 G4cout << "-----------------------------------------------------------" << G4endl; 681 } 682 683 if (fVerboseLevel > 0) 684 { 685 G4double energyDiff = std::fabs(fEnergySecondary+energyInFluorescence+fKineticEnergy1+ 686 localEnergyDeposit+energyInAuger-kineticEnergy0); 687 if (energyDiff > 0.05*keV) 688 G4cout << "Warning from G4PenelopeIonisation: problem with energy conservation: " << 689 (fEnergySecondary+energyInFluorescence+fKineticEnergy1+localEnergyDeposit+energyInAuger)/keV << 690 " keV (final) vs. " << 691 kineticEnergy0/keV << " keV (initial)" << G4endl; 692 } 693 694 } 695 696 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 697 void G4PenelopeIonisationModel::SampleFinalStateElectron(const G4Material* mat, 698 G4double cutEnergy, 699 G4double kineticEnergy) 700 { 701 // This method sets the final ionisation parameters 702 // fKineticEnergy1, fCosThetaPrimary (= updates of the primary e-) 703 // fEnergySecondary, fCosThetaSecondary (= info of delta-ray) 704 // fTargetOscillator (= ionised oscillator) 705 // 706 // The method implements SUBROUTINE EINa of Penelope 707 // 708 709 G4PenelopeOscillatorTable* theTable = fOscManager->GetOscillatorTableIonisation(mat); 710 std::size_t numberOfOscillators = theTable->size(); 711 const G4PenelopeCrossSection* theXS = 712 fCrossSectionHandler->GetCrossSectionTableForCouple(G4Electron::Electron(),mat, 713 cutEnergy); 714 G4double delta = fCrossSectionHandler->GetDensityCorrection(mat,kineticEnergy); 715 716 // Selection of the active oscillator 717 G4double TST = G4UniformRand(); 718 fTargetOscillator = G4int(numberOfOscillators-1); //initialization, last oscillator 719 G4double XSsum = 0.; 720 721 for (std::size_t i=0;i<numberOfOscillators-1;++i) 722 { 723 XSsum += theXS->GetNormalizedShellCrossSection(i,kineticEnergy); 724 725 if (XSsum > TST) 726 { 727 fTargetOscillator = (G4int) i; 728 break; 729 } 730 } 731 732 if (fVerboseLevel > 3) 733 { 734 G4cout << "SampleFinalStateElectron: sampled oscillator #" << 735 fTargetOscillator << "." << G4endl; 736 G4cout << "Ionisation energy: " << 737 (*theTable)[fTargetOscillator]->GetIonisationEnergy()/eV << 738 " eV " << G4endl; 739 G4cout << "Resonance energy: : " << 740 (*theTable)[fTargetOscillator]->GetResonanceEnergy()/eV << " eV " 741 << G4endl; 742 } 743 //Constants 744 G4double rb = kineticEnergy + 2.0*electron_mass_c2; 745 G4double gam = 1.0+kineticEnergy/electron_mass_c2; 746 G4double gam2 = gam*gam; 747 G4double beta2 = (gam2-1.0)/gam2; 748 G4double amol = ((gam-1.0)/gam)*((gam-1.0)/gam); 749 750 //Partial cross section of the active oscillator 751 G4double resEne = (*theTable)[fTargetOscillator]->GetResonanceEnergy(); 752 G4double invResEne = 1.0/resEne; 753 G4double ionEne = (*theTable)[fTargetOscillator]->GetIonisationEnergy(); 754 G4double cutoffEne = (*theTable)[fTargetOscillator]->GetCutoffRecoilResonantEnergy(); 755 G4double XHDL = 0.; 756 G4double XHDT = 0.; 757 G4double QM = 0.; 758 G4double cps = 0.; 759 G4double cp = 0.; 760 761 //Distant excitations 762 if (resEne > cutEnergy && resEne < kineticEnergy) 763 { 764 cps = kineticEnergy*rb; 765 cp = std::sqrt(cps); 766 G4double XHDT0 = std::max(G4Log(gam2)-beta2-delta,0.); 767 if (resEne > 1.0e-6*kineticEnergy) 768 { 769 G4double cpp = std::sqrt((kineticEnergy-resEne)*(kineticEnergy-resEne+2.0*electron_mass_c2)); 770 QM = std::sqrt((cp-cpp)*(cp-cpp)+electron_mass_c2*electron_mass_c2)-electron_mass_c2; 771 } 772 else 773 { 774 QM = resEne*resEne/(beta2*2.0*electron_mass_c2); 775 QM *= (1.0-QM*0.5/electron_mass_c2); 776 } 777 if (QM < cutoffEne) 778 { 779 XHDL = G4Log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2))) 780 *invResEne; 781 XHDT = XHDT0*invResEne; 782 } 783 else 784 { 785 QM = cutoffEne; 786 XHDL = 0.; 787 XHDT = 0.; 788 } 789 } 790 else 791 { 792 QM = cutoffEne; 793 cps = 0.; 794 cp = 0.; 795 XHDL = 0.; 796 XHDT = 0.; 797 } 798 799 //Close collisions 800 G4double EE = kineticEnergy + ionEne; 801 G4double wmaxc = 0.5*EE; 802 G4double wcl = std::max(cutEnergy,cutoffEne); 803 G4double rcl = wcl/EE; 804 G4double XHC = 0.; 805 if (wcl < wmaxc) 806 { 807 G4double rl1 = 1.0-rcl; 808 G4double rrl1 = 1.0/rl1; 809 XHC = (amol*(0.5-rcl)+1.0/rcl-rrl1+ 810 (1.0-amol)*G4Log(rcl*rrl1))/EE; 811 } 812 813 //Total cross section per molecule for the active shell, in cm2 814 G4double XHTOT = XHC + XHDL + XHDT; 815 816 //very small cross section, do nothing 817 if (XHTOT < 1.e-14*barn) 818 { 819 fKineticEnergy1=kineticEnergy; 820 fCosThetaPrimary=1.0; 821 fEnergySecondary=0.0; 822 fCosThetaSecondary=1.0; 823 fTargetOscillator = G4int(numberOfOscillators-1); 824 return; 825 } 826 827 //decide which kind of interaction we'll have 828 TST = XHTOT*G4UniformRand(); 829 830 // Hard close collision 831 G4double TS1 = XHC; 832 833 if (TST < TS1) 834 { 835 G4double A = 5.0*amol; 836 G4double ARCL = A*0.5*rcl; 837 G4double rk=0.; 838 G4bool loopAgain = false; 839 do 840 { 841 loopAgain = false; 842 G4double fb = (1.0+ARCL)*G4UniformRand(); 843 if (fb < 1) 844 rk = rcl/(1.0-fb*(1.0-(rcl+rcl))); 845 else 846 rk = rcl + (fb-1.0)*(0.5-rcl)/ARCL; 847 G4double rk2 = rk*rk; 848 G4double rkf = rk/(1.0-rk); 849 G4double phi = 1.0+rkf*rkf-rkf+amol*(rk2+rkf); 850 if (G4UniformRand()*(1.0+A*rk2) > phi) 851 loopAgain = true; 852 }while(loopAgain); 853 //energy and scattering angle (primary electron) 854 G4double deltaE = rk*EE; 855 856 fKineticEnergy1 = kineticEnergy - deltaE; 857 fCosThetaPrimary = std::sqrt(fKineticEnergy1*rb/(kineticEnergy*(rb-deltaE))); 858 //energy and scattering angle of the delta ray 859 fEnergySecondary = deltaE - ionEne; //subtract ionisation energy 860 fCosThetaSecondary= std::sqrt(deltaE*rb/(kineticEnergy*(deltaE+2.0*electron_mass_c2))); 861 if (fVerboseLevel > 3) 862 G4cout << "SampleFinalStateElectron: sampled close collision " << G4endl; 863 return; 864 } 865 866 //Hard distant longitudinal collisions 867 TS1 += XHDL; 868 G4double deltaE = resEne; 869 fKineticEnergy1 = kineticEnergy - deltaE; 870 871 if (TST < TS1) 872 { 873 G4double QS = QM/(1.0+QM*0.5/electron_mass_c2); 874 G4double Q = QS/(std::pow((QS/cutoffEne)*(1.0+cutoffEne*0.5/electron_mass_c2),G4UniformRand()) 875 - (QS*0.5/electron_mass_c2)); 876 G4double QTREV = Q*(Q+2.0*electron_mass_c2); 877 G4double cpps = fKineticEnergy1*(fKineticEnergy1+2.0*electron_mass_c2); 878 fCosThetaPrimary = (cpps+cps-QTREV)/(2.0*cp*std::sqrt(cpps)); 879 if (fCosThetaPrimary > 1.) 880 fCosThetaPrimary = 1.0; 881 //energy and emission angle of the delta ray 882 fEnergySecondary = deltaE - ionEne; 883 fCosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV); 884 if (fCosThetaSecondary > 1.0) 885 fCosThetaSecondary = 1.0; 886 if (fVerboseLevel > 3) 887 G4cout << "SampleFinalStateElectron: sampled distant longitudinal collision " << G4endl; 888 return; 889 } 890 891 //Hard distant transverse collisions 892 fCosThetaPrimary = 1.0; 893 //energy and emission angle of the delta ray 894 fEnergySecondary = deltaE - ionEne; 895 fCosThetaSecondary = 0.5; 896 if (fVerboseLevel > 3) 897 G4cout << "SampleFinalStateElectron: sampled distant transverse collision " << G4endl; 898 899 return; 900 } 901 902 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 903 904 void G4PenelopeIonisationModel::SampleFinalStatePositron(const G4Material* mat, 905 G4double cutEnergy, 906 G4double kineticEnergy) 907 { 908 // This method sets the final ionisation parameters 909 // fKineticEnergy1, fCosThetaPrimary (= updates of the primary e-) 910 // fEnergySecondary, fCosThetaSecondary (= info of delta-ray) 911 // fTargetOscillator (= ionised oscillator) 912 // 913 // The method implements SUBROUTINE PINa of Penelope 914 // 915 916 G4PenelopeOscillatorTable* theTable = fOscManager->GetOscillatorTableIonisation(mat); 917 std::size_t numberOfOscillators = theTable->size(); 918 const G4PenelopeCrossSection* theXS = 919 fCrossSectionHandler->GetCrossSectionTableForCouple(G4Positron::Positron(),mat, 920 cutEnergy); 921 G4double delta = fCrossSectionHandler->GetDensityCorrection(mat,kineticEnergy); 922 923 // Selection of the active oscillator 924 G4double TST = G4UniformRand(); 925 fTargetOscillator = G4int(numberOfOscillators-1); //initialization, last oscillator 926 G4double XSsum = 0.; 927 for (std::size_t i=0;i<numberOfOscillators-1;++i) 928 { 929 XSsum += theXS->GetNormalizedShellCrossSection(i,kineticEnergy); 930 if (XSsum > TST) 931 { 932 fTargetOscillator = (G4int) i; 933 break; 934 } 935 } 936 937 if (fVerboseLevel > 3) 938 { 939 G4cout << "SampleFinalStatePositron: sampled oscillator #" << 940 fTargetOscillator << "." << G4endl; 941 G4cout << "Ionisation energy: " << (*theTable)[fTargetOscillator]->GetIonisationEnergy()/eV 942 << " eV " << G4endl; 943 G4cout << "Resonance energy: : " << (*theTable)[fTargetOscillator]->GetResonanceEnergy()/eV 944 << " eV " << G4endl; 945 } 946 947 //Constants 948 G4double rb = kineticEnergy + 2.0*electron_mass_c2; 949 G4double gam = 1.0+kineticEnergy/electron_mass_c2; 950 G4double gam2 = gam*gam; 951 G4double beta2 = (gam2-1.0)/gam2; 952 G4double g12 = (gam+1.0)*(gam+1.0); 953 G4double amol = ((gam-1.0)/gam)*((gam-1.0)/gam); 954 //Bhabha coefficients 955 G4double bha1 = amol*(2.0*g12-1.0)/(gam2-1.0); 956 G4double bha2 = amol*(3.0+1.0/g12); 957 G4double bha3 = amol*2.0*gam*(gam-1.0)/g12; 958 G4double bha4 = amol*(gam-1.0)*(gam-1.0)/g12; 959 960 // 961 //Partial cross section of the active oscillator 962 // 963 G4double resEne = (*theTable)[fTargetOscillator]->GetResonanceEnergy(); 964 G4double invResEne = 1.0/resEne; 965 G4double ionEne = (*theTable)[fTargetOscillator]->GetIonisationEnergy(); 966 G4double cutoffEne = (*theTable)[fTargetOscillator]->GetCutoffRecoilResonantEnergy(); 967 968 G4double XHDL = 0.; 969 G4double XHDT = 0.; 970 G4double QM = 0.; 971 G4double cps = 0.; 972 G4double cp = 0.; 973 974 //Distant excitations XS (same as for electrons) 975 if (resEne > cutEnergy && resEne < kineticEnergy) 976 { 977 cps = kineticEnergy*rb; 978 cp = std::sqrt(cps); 979 G4double XHDT0 = std::max(G4Log(gam2)-beta2-delta,0.); 980 if (resEne > 1.0e-6*kineticEnergy) 981 { 982 G4double cpp = std::sqrt((kineticEnergy-resEne)*(kineticEnergy-resEne+2.0*electron_mass_c2)); 983 QM = std::sqrt((cp-cpp)*(cp-cpp)+electron_mass_c2*electron_mass_c2)-electron_mass_c2; 984 } 985 else 986 { 987 QM = resEne*resEne/(beta2*2.0*electron_mass_c2); 988 QM *= (1.0-QM*0.5/electron_mass_c2); 989 } 990 if (QM < cutoffEne) 991 { 992 XHDL = G4Log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2))) 993 *invResEne; 994 XHDT = XHDT0*invResEne; 995 } 996 else 997 { 998 QM = cutoffEne; 999 XHDL = 0.; 1000 XHDT = 0.; 1001 } 1002 } 1003 else 1004 { 1005 QM = cutoffEne; 1006 cps = 0.; 1007 cp = 0.; 1008 XHDL = 0.; 1009 XHDT = 0.; 1010 } 1011 //Close collisions (Bhabha) 1012 G4double wmaxc = kineticEnergy; 1013 G4double wcl = std::max(cutEnergy,cutoffEne); 1014 G4double rcl = wcl/kineticEnergy; 1015 G4double XHC = 0.; 1016 if (wcl < wmaxc) 1017 { 1018 G4double rl1 = 1.0-rcl; 1019 XHC = ((1.0/rcl-1.0)+bha1*G4Log(rcl)+bha2*rl1 1020 + (bha3/2.0)*(rcl*rcl-1.0) 1021 + (bha4/3.0)*(1.0-rcl*rcl*rcl))/kineticEnergy; 1022 } 1023 1024 //Total cross section per molecule for the active shell, in cm2 1025 G4double XHTOT = XHC + XHDL + XHDT; 1026 1027 //very small cross section, do nothing 1028 if (XHTOT < 1.e-14*barn) 1029 { 1030 fKineticEnergy1=kineticEnergy; 1031 fCosThetaPrimary=1.0; 1032 fEnergySecondary=0.0; 1033 fCosThetaSecondary=1.0; 1034 fTargetOscillator = G4int(numberOfOscillators-1); 1035 return; 1036 } 1037 1038 //decide which kind of interaction we'll have 1039 TST = XHTOT*G4UniformRand(); 1040 1041 // Hard close collision 1042 G4double TS1 = XHC; 1043 if (TST < TS1) 1044 { 1045 G4double rl1 = 1.0-rcl; 1046 G4double rk=0.; 1047 G4bool loopAgain = false; 1048 do 1049 { 1050 loopAgain = false; 1051 rk = rcl/(1.0-G4UniformRand()*rl1); 1052 G4double phi = 1.0-rk*(bha1-rk*(bha2-rk*(bha3-bha4*rk))); 1053 if (G4UniformRand() > phi) 1054 loopAgain = true; 1055 }while(loopAgain); 1056 //energy and scattering angle (primary electron) 1057 G4double deltaE = rk*kineticEnergy; 1058 fKineticEnergy1 = kineticEnergy - deltaE; 1059 fCosThetaPrimary = std::sqrt(fKineticEnergy1*rb/(kineticEnergy*(rb-deltaE))); 1060 //energy and scattering angle of the delta ray 1061 fEnergySecondary = deltaE - ionEne; //subtract ionisation energy 1062 fCosThetaSecondary= std::sqrt(deltaE*rb/(kineticEnergy*(deltaE+2.0*electron_mass_c2))); 1063 if (fVerboseLevel > 3) 1064 G4cout << "SampleFinalStatePositron: sampled close collision " << G4endl; 1065 return; 1066 } 1067 1068 //Hard distant longitudinal collisions 1069 TS1 += XHDL; 1070 G4double deltaE = resEne; 1071 fKineticEnergy1 = kineticEnergy - deltaE; 1072 if (TST < TS1) 1073 { 1074 G4double QS = QM/(1.0+QM*0.5/electron_mass_c2); 1075 G4double Q = QS/(std::pow((QS/cutoffEne)*(1.0+cutoffEne*0.5/electron_mass_c2),G4UniformRand()) 1076 - (QS*0.5/electron_mass_c2)); 1077 G4double QTREV = Q*(Q+2.0*electron_mass_c2); 1078 G4double cpps = fKineticEnergy1*(fKineticEnergy1+2.0*electron_mass_c2); 1079 fCosThetaPrimary = (cpps+cps-QTREV)/(2.0*cp*std::sqrt(cpps)); 1080 if (fCosThetaPrimary > 1.) 1081 fCosThetaPrimary = 1.0; 1082 //energy and emission angle of the delta ray 1083 fEnergySecondary = deltaE - ionEne; 1084 fCosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV); 1085 if (fCosThetaSecondary > 1.0) 1086 fCosThetaSecondary = 1.0; 1087 if (fVerboseLevel > 3) 1088 G4cout << "SampleFinalStatePositron: sampled distant longitudinal collision " << G4endl; 1089 return; 1090 } 1091 1092 //Hard distant transverse collisions 1093 fCosThetaPrimary = 1.0; 1094 //energy and emission angle of the delta ray 1095 fEnergySecondary = deltaE - ionEne; 1096 fCosThetaSecondary = 0.5; 1097 1098 if (fVerboseLevel > 3) 1099 G4cout << "SampleFinalStatePositron: sampled distant transverse collision " << G4endl; 1100 1101 return; 1102 } 1103 1104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo... 1105 1106 void G4PenelopeIonisationModel::SetParticle(const G4ParticleDefinition* p) 1107 { 1108 if(!fParticle) { 1109 fParticle = p; 1110 } 1111 } 1112