Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4PenelopeRayleighModel Class Reference

#include <G4PenelopeRayleighModel.hh>

+ Inheritance diagram for G4PenelopeRayleighModel:

Public Member Functions

 G4PenelopeRayleighModel (const G4ParticleDefinition *p=0, const G4String &processName="PenRayleigh")
 
virtual ~G4PenelopeRayleighModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)
 
virtual void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel)
 
virtual G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
 
void SetVerbosityLevel (G4int lev)
 
G4int GetVerbosityLevel ()
 
void DumpFormFactorTable (const G4Material *)
 
- Public Member Functions inherited from G4VEmModel
 G4VEmModel (const G4String &nam)
 
virtual ~G4VEmModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)=0
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin=0.0, G4double tmax=DBL_MAX)=0
 
virtual void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel)
 
virtual void InitialiseForMaterial (const G4ParticleDefinition *, const G4Material *)
 
virtual void InitialiseForElement (const G4ParticleDefinition *, G4int Z)
 
virtual G4double ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
virtual G4double CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double GetPartialCrossSection (const G4Material *, G4int level, const G4ParticleDefinition *, G4double kineticEnergy)
 
virtual G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double ComputeCrossSectionPerShell (const G4ParticleDefinition *, G4int Z, G4int shellIdx, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double ChargeSquareRatio (const G4Track &)
 
virtual G4double GetChargeSquareRatio (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual G4double GetParticleCharge (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void StartTracking (G4Track *)
 
virtual void CorrectionsAlongStep (const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length)
 
virtual G4double Value (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
 
virtual G4double MinPrimaryEnergy (const G4Material *, const G4ParticleDefinition *, G4double cut=0.0)
 
virtual G4double MinEnergyCut (const G4ParticleDefinition *, const G4MaterialCutsCouple *)
 
virtual void SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void DefineForRegion (const G4Region *)
 
virtual void ModelDescription (std::ostream &outFile) const
 
void InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
 
std::vector< G4EmElementSelector * > * GetElementSelectors ()
 
void SetElementSelectors (std::vector< G4EmElementSelector * > *)
 
virtual G4double ComputeDEDX (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
G4double CrossSection (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeMeanFreePath (const G4ParticleDefinition *, G4double kineticEnergy, const G4Material *, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, const G4Element *, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4ElementSelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4ElementSelectTargetAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double logKineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4ElementSelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4int SelectRandomAtomNumber (const G4Material *)
 
G4int SelectIsotopeNumber (const G4Element *)
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=nullptr)
 
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
 
G4ElementDataGetElementData ()
 
G4PhysicsTableGetCrossSectionTable ()
 
G4VEmFluctuationModelGetModelOfFluctuations ()
 
G4VEmAngularDistributionGetAngularDistribution ()
 
G4VEmModelGetTripletModel ()
 
void SetTripletModel (G4VEmModel *)
 
void SetAngularDistribution (G4VEmAngularDistribution *)
 
G4double HighEnergyLimit () const
 
G4double LowEnergyLimit () const
 
G4double HighEnergyActivationLimit () const
 
G4double LowEnergyActivationLimit () const
 
G4double PolarAngleLimit () const
 
G4double SecondaryThreshold () const
 
G4bool LPMFlag () const
 
G4bool DeexcitationFlag () const
 
G4bool ForceBuildTableFlag () const
 
G4bool UseAngularGeneratorFlag () const
 
void SetAngularGeneratorFlag (G4bool)
 
void SetHighEnergyLimit (G4double)
 
void SetLowEnergyLimit (G4double)
 
void SetActivationHighEnergyLimit (G4double)
 
void SetActivationLowEnergyLimit (G4double)
 
G4bool IsActive (G4double kinEnergy) const
 
void SetPolarAngleLimit (G4double)
 
void SetSecondaryThreshold (G4double)
 
void SetLPMFlag (G4bool val)
 
void SetDeexcitationFlag (G4bool val)
 
void SetForceBuildTable (G4bool val)
 
void SetFluctuationFlag (G4bool val)
 
void SetMasterThread (G4bool val)
 
G4bool IsMaster () const
 
void SetUseBaseMaterials (G4bool val)
 
G4bool UseBaseMaterials () const
 
G4double MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
 
const G4StringGetName () const
 
void SetCurrentCouple (const G4MaterialCutsCouple *)
 
const G4ElementGetCurrentElement () const
 
const G4IsotopeGetCurrentIsotope () const
 
G4bool IsLocked () const
 
void SetLocked (G4bool)
 
G4VEmModeloperator= (const G4VEmModel &right)=delete
 
 G4VEmModel (const G4VEmModel &)=delete
 

Protected Attributes

G4ParticleChangeForGammafParticleChange
 
const G4ParticleDefinitionfParticle
 
- Protected Attributes inherited from G4VEmModel
G4ElementDatafElementData
 
G4VParticleChangepParticleChange
 
G4PhysicsTablexSectionTable
 
const G4MaterialpBaseMaterial
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 
size_t idxTable
 
G4bool lossFlucFlag
 
G4double inveplus
 
G4double pFactor
 

Additional Inherited Members

- Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLossGetParticleChangeForLoss ()
 
G4ParticleChangeForGammaGetParticleChangeForGamma ()
 
virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
 
const G4MaterialCutsCoupleCurrentCouple () const
 
void SetCurrentElement (const G4Element *)
 

Detailed Description

Definition at line 57 of file G4PenelopeRayleighModel.hh.

Constructor & Destructor Documentation

◆ G4PenelopeRayleighModel()

G4PenelopeRayleighModel::G4PenelopeRayleighModel ( const G4ParticleDefinition p = 0,
const G4String processName = "PenRayleigh" 
)

Definition at line 53 of file G4PenelopeRayleighModel.cc.

55 :G4VEmModel(nam),fParticleChange(0),fParticle(0),isInitialised(false),
56 logAtomicCrossSection(0),atomicFormFactor(0),logFormFactorTable(0),
57 pMaxTable(0),samplingTable(0),fLocalTable(false)
58{
59 fIntrinsicLowEnergyLimit = 100.0*eV;
60 fIntrinsicHighEnergyLimit = 100.0*GeV;
61 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
62 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
63
64 if (part)
65 SetParticle(part);
66
67 //
68 verboseLevel= 0;
69 // Verbosity scale:
70 // 0 = nothing
71 // 1 = warning for energy non-conservation
72 // 2 = details of energy budget
73 // 3 = calculation of cross sections, file openings, sampling of atoms
74 // 4 = entering in methods
75
76 //build the energy grid. It is the same for all materials
77 G4double logenergy = G4Log(fIntrinsicLowEnergyLimit/2.);
78 G4double logmaxenergy = G4Log(1.5*fIntrinsicHighEnergyLimit);
79 //finer grid below 160 keV
80 G4double logtransitionenergy = G4Log(160*keV);
81 G4double logfactor1 = G4Log(10.)/250.;
82 G4double logfactor2 = logfactor1*10;
83 logEnergyGridPMax.push_back(logenergy);
84 do{
85 if (logenergy < logtransitionenergy)
86 logenergy += logfactor1;
87 else
88 logenergy += logfactor2;
89 logEnergyGridPMax.push_back(logenergy);
90 }while (logenergy < logmaxenergy);
91}
G4double G4Log(G4double x)
Definition: G4Log.hh:226
double G4double
Definition: G4Types.hh:83
const G4ParticleDefinition * fParticle
G4ParticleChangeForGamma * fParticleChange
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:757

◆ ~G4PenelopeRayleighModel()

G4PenelopeRayleighModel::~G4PenelopeRayleighModel ( )
virtual

Definition at line 95 of file G4PenelopeRayleighModel.cc.

96{
97 if (IsMaster() || fLocalTable)
98 {
99 if (logAtomicCrossSection)
100 {
101 for (auto& item : (*logAtomicCrossSection))
102 if (item.second) delete item.second;
103 delete logAtomicCrossSection;
104 logAtomicCrossSection = nullptr;
105 }
106 if (atomicFormFactor)
107 {
108 for (auto& item : (*atomicFormFactor))
109 if (item.second) delete item.second;
110 delete atomicFormFactor;
111 atomicFormFactor = nullptr;
112 }
113 ClearTables();
114 }
115}
G4bool IsMaster() const
Definition: G4VEmModel.hh:736

Member Function Documentation

◆ ComputeCrossSectionPerAtom()

G4double G4PenelopeRayleighModel::ComputeCrossSectionPerAtom ( const G4ParticleDefinition ,
G4double  kinEnergy,
G4double  Z,
G4double  A = 0,
G4double  cut = 0,
G4double  emax = DBL_MAX 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 273 of file G4PenelopeRayleighModel.cc.

279{
280 // Cross section of Rayleigh scattering in Penelope v2008 is calculated by the EPDL97
281 // tabulation, Cuellen et al. (1997), with non-relativistic form factors from Hubbel
282 // et al. J. Phys. Chem. Ref. Data 4 (1975) 471; Erratum ibid. 6 (1977) 615.
283
284 if (verboseLevel > 3)
285 G4cout << "Calling CrossSectionPerAtom() of G4PenelopeRayleighModel" << G4endl;
286
287 G4int iZ = (G4int) Z;
288
289 //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
290 //not invoked
291 if (!logAtomicCrossSection)
292 {
293 //create a **thread-local** version of the table. Used only for G4EmCalculator and
294 //Unit Tests
295 fLocalTable = true;
296 logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>;
297 }
298 //now it should be ok
299 if (!logAtomicCrossSection->count(iZ))
300 {
301 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
302 //not filled up. This can happen in a UnitTest or via G4EmCalculator
303 if (verboseLevel > 0)
304 {
305 //Issue a G4Exception (warning) only in verbose mode
307 ed << "Unable to retrieve the cross section table for Z=" << iZ << G4endl;
308 ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl;
309 G4Exception("G4PenelopeRayleighModel::ComputeCrossSectionPerAtom()",
310 "em2040",JustWarning,ed);
311 }
312 //protect file reading via autolock
313 G4AutoLock lock(&PenelopeRayleighModelMutex);
314 ReadDataFile(iZ);
315 lock.unlock();
316 }
317
318 G4double cross = 0;
319
320 G4PhysicsFreeVector* atom = logAtomicCrossSection->find(iZ)->second;
321 if (!atom)
322 {
324 ed << "Unable to find Z=" << iZ << " in the atomic cross section table" << G4endl;
325 G4Exception("G4PenelopeRayleighModel::ComputeCrossSectionPerAtom()",
326 "em2041",FatalException,ed);
327 return 0;
328 }
329 G4double logene = G4Log(energy);
330 G4double logXS = atom->Value(logene);
331 cross = G4Exp(logXS);
332
333 if (verboseLevel > 2)
334 {
335 G4cout << "Rayleigh cross section at " << energy/keV << " keV for Z=" << Z <<
336 " = " << cross/barn << " barn" << G4endl;
337 }
338
339 return cross;
340}
@ JustWarning
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:40
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:179
int G4int
Definition: G4Types.hh:85
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
G4double Value(G4double theEnergy, std::size_t &lastidx) const
G4double energy(const ThreeVector &p, const G4double m)

◆ DumpFormFactorTable()

void G4PenelopeRayleighModel::DumpFormFactorTable ( const G4Material mat)

Definition at line 1332 of file G4PenelopeRayleighModel.cc.

1333{
1334 G4cout << "*****************************************************************" << G4endl;
1335 G4cout << "G4PenelopeRayleighModel: Form Factor Table for " << mat->GetName() << G4endl;
1336 //try to use the same format as Penelope-Fortran, namely Q (/m_e*c) and F
1337 G4cout << "Q/(m_e*c) F(Q) " << G4endl;
1338 G4cout << "*****************************************************************" << G4endl;
1339 if (!logFormFactorTable->count(mat))
1340 BuildFormFactorTable(mat);
1341
1342 G4PhysicsFreeVector* theVec = logFormFactorTable->find(mat)->second;
1343 for (size_t i=0;i<theVec->GetVectorLength();i++)
1344 {
1345 G4double logQ2 = theVec->GetLowEdgeEnergy(i);
1346 G4double Q = G4Exp(0.5*logQ2);
1347 G4double logF2 = (*theVec)[i];
1348 G4double F = G4Exp(0.5*logF2);
1349 G4cout << Q << " " << F << G4endl;
1350 }
1351 //DONE
1352 return;
1353}
const G4String & GetName() const
Definition: G4Material.hh:175
G4double GetLowEdgeEnergy(std::size_t binNumber) const
std::size_t GetVectorLength() const

◆ GetVerbosityLevel()

G4int G4PenelopeRayleighModel::GetVerbosityLevel ( )
inline

Definition at line 84 of file G4PenelopeRayleighModel.hh.

84{return verboseLevel;};

◆ Initialise()

void G4PenelopeRayleighModel::Initialise ( const G4ParticleDefinition part,
const G4DataVector  
)
virtual

Implements G4VEmModel.

Definition at line 156 of file G4PenelopeRayleighModel.cc.

158{
159 if (verboseLevel > 3)
160 G4cout << "Calling G4PenelopeRayleighModel::Initialise()" << G4endl;
161
162 SetParticle(part);
163
164 //Only the master model creates/fills/destroys the tables
165 if (IsMaster() && part == fParticle)
166 {
167 //clear tables depending on materials, not the atomic ones
168 ClearTables();
169
170 if (verboseLevel > 3)
171 G4cout << "Calling G4PenelopeRayleighModel::Initialise() [master]" << G4endl;
172
173 //create new tables
174 //
175 // logAtomicCrossSection and atomicFormFactor are created only once,
176 // since they are never cleared
177 if (!logAtomicCrossSection)
178 logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>;
179 if (!atomicFormFactor)
180 atomicFormFactor = new std::map<G4int,G4PhysicsFreeVector*>;
181
182 if (!logFormFactorTable)
183 logFormFactorTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
184 if (!pMaxTable)
185 pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
186 if (!samplingTable)
187 samplingTable = new std::map<const G4Material*,G4PenelopeSamplingData*>;
188
189
190 G4ProductionCutsTable* theCoupleTable =
192
193 for (size_t i=0;i<theCoupleTable->GetTableSize();i++)
194 {
195 const G4Material* material =
196 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
197 const G4ElementVector* theElementVector = material->GetElementVector();
198
199 for (size_t j=0;j<material->GetNumberOfElements();j++)
200 {
201 G4int iZ = theElementVector->at(j)->GetZasInt();
202 //read data files only in the master
203 if (!logAtomicCrossSection->count(iZ))
204 ReadDataFile(iZ);
205 }
206
207 //1) If the table has not been built for the material, do it!
208 if (!logFormFactorTable->count(material))
209 BuildFormFactorTable(material);
210
211 //2) retrieve or build the sampling table
212 if (!(samplingTable->count(material)))
213 InitializeSamplingAlgorithm(material);
214
215 //3) retrieve or build the pMax data
216 if (!pMaxTable->count(material))
217 GetPMaxTable(material);
218
219 }
220
221 if (verboseLevel > 1) {
222 G4cout << "Penelope Rayleigh model v2008 is initialized " << G4endl
223 << "Energy range: "
224 << LowEnergyLimit() / keV << " keV - "
225 << HighEnergyLimit() / GeV << " GeV"
226 << G4endl;
227 }
228 }
229
230 if(isInitialised) return;
232 isInitialised = true;
233}
std::vector< G4Element * > G4ElementVector
const G4Material * GetMaterial() const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:188
size_t GetNumberOfElements() const
Definition: G4Material.hh:184
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
std::size_t GetTableSize() const
static G4ProductionCutsTable * GetProductionCutsTable()
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:133
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:652
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:645

◆ InitialiseLocal()

void G4PenelopeRayleighModel::InitialiseLocal ( const G4ParticleDefinition part,
G4VEmModel masterModel 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 237 of file G4PenelopeRayleighModel.cc.

239{
240 if (verboseLevel > 3)
241 G4cout << "Calling G4PenelopeRayleighModel::InitialiseLocal()" << G4endl;
242
243 //
244 //Check that particle matches: one might have multiple master models (e.g.
245 //for e+ and e-).
246 //
247 if (part == fParticle)
248 {
249 //Get the const table pointers from the master to the workers
250 const G4PenelopeRayleighModel* theModel =
251 static_cast<G4PenelopeRayleighModel*> (masterModel);
252
253 //Copy pointers to the data tables
254 logAtomicCrossSection = theModel->logAtomicCrossSection;
255 atomicFormFactor = theModel->atomicFormFactor;
256 logFormFactorTable = theModel->logFormFactorTable;
257 pMaxTable = theModel->pMaxTable;
258 samplingTable = theModel->samplingTable;
259
260 //copy the G4DataVector with the grid
261 logQSquareGrid = theModel->logQSquareGrid;
262
263 //Same verbosity for all workers, as the master
264 verboseLevel = theModel->verboseLevel;
265 }
266
267 return;
268}

◆ SampleSecondaries()

void G4PenelopeRayleighModel::SampleSecondaries ( std::vector< G4DynamicParticle * > *  ,
const G4MaterialCutsCouple couple,
const G4DynamicParticle aDynamicGamma,
G4double  tmin,
G4double  maxEnergy 
)
virtual

Implements G4VEmModel.

Definition at line 415 of file G4PenelopeRayleighModel.cc.

420{
421 // Sampling of the Rayleigh final state (namely, scattering angle of the photon)
422 // from the Penelope2008 model. The scattering angle is sampled from the atomic
423 // cross section dOmega/d(cosTheta) from Born ("Atomic Phyisics", 1969), disregarding
424 // anomalous scattering effects. The Form Factor F(Q) function which appears in the
425 // analytical cross section is retrieved via the method GetFSquared(); atomic data
426 // are tabulated for F(Q). Form factor for compounds is calculated according to
427 // the additivity rule. The sampling from the F(Q) is made via a Rational Inverse
428 // Transform with Aliasing (RITA) algorithm; RITA parameters are calculated once
429 // for each material and managed by G4PenelopeSamplingData objects.
430 // The sampling algorithm (rejection method) has efficiency 67% at low energy, and
431 // increases with energy. For E=100 keV the efficiency is 100% and 86% for
432 // hydrogen and uranium, respectively.
433
434 if (verboseLevel > 3)
435 G4cout << "Calling SamplingSecondaries() of G4PenelopeRayleighModel" << G4endl;
436
437 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
438
439 if (photonEnergy0 <= fIntrinsicLowEnergyLimit)
440 {
444 return ;
445 }
446
447 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
448
449 const G4Material* theMat = couple->GetMaterial();
450
451
452 //1) Verify if tables are ready
453 //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
454 //not invoked
455 if (!pMaxTable || !samplingTable || !logAtomicCrossSection || !atomicFormFactor ||
456 !logFormFactorTable)
457 {
458 //create a **thread-local** version of the table. Used only for G4EmCalculator and
459 //Unit Tests
460 fLocalTable = true;
461 if (!logAtomicCrossSection)
462 logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>;
463 if (!atomicFormFactor)
464 atomicFormFactor = new std::map<G4int,G4PhysicsFreeVector*>;
465 if (!logFormFactorTable)
466 logFormFactorTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
467 if (!pMaxTable)
468 pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
469 if (!samplingTable)
470 samplingTable = new std::map<const G4Material*,G4PenelopeSamplingData*>;
471 }
472
473 if (!samplingTable->count(theMat))
474 {
475 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
476 //not filled up. This can happen in a UnitTest
477 if (verboseLevel > 0)
478 {
479 //Issue a G4Exception (warning) only in verbose mode
481 ed << "Unable to find the samplingTable data for " <<
482 theMat->GetName() << G4endl;
483 ed << "This can happen only in Unit Tests" << G4endl;
484 G4Exception("G4PenelopeRayleighModel::SampleSecondaries()",
485 "em2019",JustWarning,ed);
486 }
487 const G4ElementVector* theElementVector = theMat->GetElementVector();
488 //protect file reading via autolock
489 G4AutoLock lock(&PenelopeRayleighModelMutex);
490 for (size_t j=0;j<theMat->GetNumberOfElements();j++)
491 {
492 G4int iZ = theElementVector->at(j)->GetZasInt();
493 if (!logAtomicCrossSection->count(iZ))
494 {
495 lock.lock();
496 ReadDataFile(iZ);
497 lock.unlock();
498 }
499 }
500 lock.lock();
501 //1) If the table has not been built for the material, do it!
502 if (!logFormFactorTable->count(theMat))
503 BuildFormFactorTable(theMat);
504
505 //2) retrieve or build the sampling table
506 if (!(samplingTable->count(theMat)))
507 InitializeSamplingAlgorithm(theMat);
508
509 //3) retrieve or build the pMax data
510 if (!pMaxTable->count(theMat))
511 GetPMaxTable(theMat);
512 lock.unlock();
513 }
514
515 //Ok, restart the job
516
517 G4PenelopeSamplingData* theDataTable = samplingTable->find(theMat)->second;
518 G4PhysicsFreeVector* thePMax = pMaxTable->find(theMat)->second;
519
520 G4double cosTheta = 1.0;
521
522 //OK, ready to go!
523 G4double qmax = 2.0*photonEnergy0/electron_mass_c2; //this is non-dimensional now
524
525 if (qmax < 1e-10) //very low momentum transfer
526 {
527 G4bool loopAgain=false;
528 do
529 {
530 loopAgain = false;
531 cosTheta = 1.0-2.0*G4UniformRand();
532 G4double G = 0.5*(1+cosTheta*cosTheta);
533 if (G4UniformRand()>G)
534 loopAgain = true;
535 }while(loopAgain);
536 }
537 else //larger momentum transfer
538 {
539 size_t nData = theDataTable->GetNumberOfStoredPoints();
540 G4double LastQ2inTheTable = theDataTable->GetX(nData-1);
541 G4double q2max = std::min(qmax*qmax,LastQ2inTheTable);
542
543 G4bool loopAgain = false;
544 G4double MaxPValue = thePMax->Value(photonEnergy0);
545 G4double xx=0;
546
547 //Sampling by rejection method. The rejection function is
548 //G = 0.5*(1+cos^2(theta))
549 //
550 do{
551 loopAgain = false;
552 G4double RandomMax = G4UniformRand()*MaxPValue;
553 xx = theDataTable->SampleValue(RandomMax);
554 //xx is a random value of q^2 in (0,q2max),sampled according to
555 //F(Q^2) via the RITA algorithm
556 if (xx > q2max)
557 loopAgain = true;
558 cosTheta = 1.0-2.0*xx/q2max;
559 G4double G = 0.5*(1+cosTheta*cosTheta);
560 if (G4UniformRand()>G)
561 loopAgain = true;
562 }while(loopAgain);
563 }
564
565 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
566
567 // Scattered photon angles. ( Z - axis along the parent photon)
568 G4double phi = twopi * G4UniformRand() ;
569 G4double dirX = sinTheta*std::cos(phi);
570 G4double dirY = sinTheta*std::sin(phi);
571 G4double dirZ = cosTheta;
572
573 // Update G4VParticleChange for the scattered photon
574 G4ThreeVector photonDirection1(dirX, dirY, dirZ);
575
576 photonDirection1.rotateUz(photonDirection0);
577 fParticleChange->ProposeMomentumDirection(photonDirection1) ;
579
580 return;
581}
@ fStopAndKill
bool G4bool
Definition: G4Types.hh:86
#define G4UniformRand()
Definition: Randomize.hh:52
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4double GetX(size_t index)
G4double SampleValue(G4double rndm)
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)

◆ SetVerbosityLevel()

void G4PenelopeRayleighModel::SetVerbosityLevel ( G4int  lev)
inline

Definition at line 83 of file G4PenelopeRayleighModel.hh.

83{verboseLevel = lev;};

Member Data Documentation

◆ fParticle

const G4ParticleDefinition* G4PenelopeRayleighModel::fParticle
protected

Definition at line 91 of file G4PenelopeRayleighModel.hh.

Referenced by Initialise(), and InitialiseLocal().

◆ fParticleChange

G4ParticleChangeForGamma* G4PenelopeRayleighModel::fParticleChange
protected

Definition at line 90 of file G4PenelopeRayleighModel.hh.

Referenced by Initialise(), and SampleSecondaries().


The documentation for this class was generated from the following files: