Geant4 11.1.1
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=nullptr, const G4String &processName="PenRayleigh")
 
virtual ~G4PenelopeRayleighModel ()
 
void Initialise (const G4ParticleDefinition *, const G4DataVector &) override
 
void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel) override
 
G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX) override
 
void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
 
void SetVerbosityLevel (G4int lev)
 
G4int GetVerbosityLevel ()
 
void DumpFormFactorTable (const G4Material *)
 
G4PenelopeRayleighModeloperator= (const G4PenelopeRayleighModel &right)=delete
 
 G4PenelopeRayleighModel (const G4PenelopeRayleighModel &)=delete
 
- 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 *, const G4double &length, G4double &eloss)
 
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 FillNumberOfSecondaries (G4int &numberOfTriplets, G4int &numberOfRecoil)
 
virtual void ModelDescription (std::ostream &outFile) const
 
void InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
 
std::vector< G4EmElementSelector * > * GetElementSelectors ()
 
void SetElementSelectors (std::vector< G4EmElementSelector * > *)
 
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)
 
const G4ElementGetCurrentElement (const G4Material *mat=nullptr) const
 
G4int SelectRandomAtomNumber (const G4Material *) const
 
const G4IsotopeGetCurrentIsotope (const G4Element *elm=nullptr) const
 
G4int SelectIsotopeNumber (const G4Element *) const
 
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 *)
 
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 = nullptr
 
G4VParticleChangepParticleChange = nullptr
 
G4PhysicsTablexSectionTable = nullptr
 
const G4MaterialpBaseMaterial = nullptr
 
const std::vector< G4double > * theDensityFactor = nullptr
 
const std::vector< G4int > * theDensityIdx = nullptr
 
G4double inveplus
 
G4double pFactor = 1.0
 
size_t currentCoupleIndex = 0
 
size_t basedCoupleIndex = 0
 
G4bool lossFlucFlag = true
 

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 58 of file G4PenelopeRayleighModel.hh.

Constructor & Destructor Documentation

◆ G4PenelopeRayleighModel() [1/2]

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

Definition at line 59 of file G4PenelopeRayleighModel.cc.

61 :G4VEmModel(nam),fParticleChange(nullptr),fParticle(nullptr),
62 fLogFormFactorTable(nullptr),fPMaxTable(nullptr),fSamplingTable(nullptr),
63 fIsInitialised(false),fLocalTable(false)
64{
65 fIntrinsicLowEnergyLimit = 100.0*eV;
66 fIntrinsicHighEnergyLimit = 100.0*GeV;
67 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
68 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
69
70 if (part)
71 SetParticle(part);
72
73 //
74 fVerboseLevel= 0;
75 // Verbosity scale:
76 // 0 = nothing
77 // 1 = warning for energy non-conservation
78 // 2 = details of energy budget
79 // 3 = calculation of cross sections, file openings, sampling of atoms
80 // 4 = entering in methods
81
82 //build the energy grid. It is the same for all materials
83 G4double logenergy = G4Log(fIntrinsicLowEnergyLimit/2.);
84 G4double logmaxenergy = G4Log(1.5*fIntrinsicHighEnergyLimit);
85 //finer grid below 160 keV
86 G4double logtransitionenergy = G4Log(160*keV);
87 G4double logfactor1 = G4Log(10.)/250.;
88 G4double logfactor2 = logfactor1*10;
89 fLogEnergyGridPMax.push_back(logenergy);
90 do{
91 if (logenergy < logtransitionenergy)
92 logenergy += logfactor1;
93 else
94 logenergy += logfactor2;
95 fLogEnergyGridPMax.push_back(logenergy);
96 }while (logenergy < logmaxenergy);
97}
G4double G4Log(G4double x)
Definition: G4Log.hh:227
double G4double
Definition: G4Types.hh:83
const G4ParticleDefinition * fParticle
G4ParticleChangeForGamma * fParticleChange
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:746

◆ ~G4PenelopeRayleighModel()

G4PenelopeRayleighModel::~G4PenelopeRayleighModel ( )
virtual

Definition at line 101 of file G4PenelopeRayleighModel.cc.

102{
103 if (IsMaster() || fLocalTable)
104 {
105
106 for(G4int i=0; i<=fMaxZ; ++i)
107 {
108 if(fLogAtomicCrossSection[i])
109 {
110 delete fLogAtomicCrossSection[i];
111 fLogAtomicCrossSection[i] = nullptr;
112 }
113 if(fAtomicFormFactor[i])
114 {
115 delete fAtomicFormFactor[i];
116 fAtomicFormFactor[i] = nullptr;
117 }
118 }
119 ClearTables();
120 }
121}
int G4int
Definition: G4Types.hh:85
G4bool IsMaster() const
Definition: G4VEmModel.hh:725

◆ G4PenelopeRayleighModel() [2/2]

G4PenelopeRayleighModel::G4PenelopeRayleighModel ( const G4PenelopeRayleighModel )
delete

Member Function Documentation

◆ ComputeCrossSectionPerAtom()

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

Reimplemented from G4VEmModel.

Definition at line 261 of file G4PenelopeRayleighModel.cc.

267{
268 // Cross section of Rayleigh scattering in Penelope v2008 is calculated by the EPDL97
269 // tabulation, Cuellen et al. (1997), with non-relativistic form factors from Hubbel
270 // et al. J. Phys. Chem. Ref. Data 4 (1975) 471; Erratum ibid. 6 (1977) 615.
271
272 if (fVerboseLevel > 3)
273 G4cout << "Calling CrossSectionPerAtom() of G4PenelopeRayleighModel" << G4endl;
274
275 G4int iZ = G4int(Z);
276
277 if (!fLogAtomicCrossSection[iZ])
278 {
279 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
280 //not filled up. This can happen in a UnitTest or via G4EmCalculator
281 if (fVerboseLevel > 0)
282 {
283 //Issue a G4Exception (warning) only in verbose mode
285 ed << "Unable to retrieve the cross section table for Z=" << iZ << G4endl;
286 ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl;
287 G4Exception("G4PenelopeRayleighModel::ComputeCrossSectionPerAtom()",
288 "em2040",JustWarning,ed);
289 }
290 //protect file reading via autolock
291 G4AutoLock lock(&PenelopeRayleighModelMutex);
292 ReadDataFile(iZ);
293 lock.unlock();
294 }
295
296 G4double cross = 0;
297 G4PhysicsFreeVector* atom = fLogAtomicCrossSection[iZ];
298 if (!atom)
299 {
301 ed << "Unable to find Z=" << iZ << " in the atomic cross section table" << G4endl;
302 G4Exception("G4PenelopeRayleighModel::ComputeCrossSectionPerAtom()",
303 "em2041",FatalException,ed);
304 return 0;
305 }
306 G4double logene = G4Log(energy);
307 G4double logXS = atom->Value(logene);
308 cross = G4Exp(logXS);
309
310 if (fVerboseLevel > 2)
311 {
312 G4cout << "Rayleigh cross section at " << energy/keV << " keV for Z=" << Z <<
313 " = " << cross/barn << " barn" << G4endl;
314 }
315 return cross;
316}
@ JustWarning
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:59
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:40
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
const G4int Z[17]
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
G4double Value(const G4double energy, std::size_t &lastidx) const
G4double energy(const ThreeVector &p, const G4double m)

◆ DumpFormFactorTable()

void G4PenelopeRayleighModel::DumpFormFactorTable ( const G4Material mat)

Definition at line 1264 of file G4PenelopeRayleighModel.cc.

1265{
1266 G4cout << "*****************************************************************" << G4endl;
1267 G4cout << "G4PenelopeRayleighModel: Form Factor Table for " << mat->GetName() << G4endl;
1268 //try to use the same format as Penelope-Fortran, namely Q (/m_e*c) and F
1269 G4cout << "Q/(m_e*c) F(Q) " << G4endl;
1270 G4cout << "*****************************************************************" << G4endl;
1271 if (!fLogFormFactorTable->count(mat))
1272 BuildFormFactorTable(mat);
1273
1274 G4PhysicsFreeVector* theVec = fLogFormFactorTable->find(mat)->second;
1275 for (std::size_t i=0;i<theVec->GetVectorLength();++i)
1276 {
1277 G4double logQ2 = theVec->GetLowEdgeEnergy(i);
1278 G4double Q = G4Exp(0.5*logQ2);
1279 G4double logF2 = (*theVec)[i];
1280 G4double F = G4Exp(0.5*logF2);
1281 G4cout << Q << " " << F << G4endl;
1282 }
1283 //DONE
1284 return;
1285}
const G4String & GetName() const
Definition: G4Material.hh:172
G4double GetLowEdgeEnergy(const std::size_t index) const
std::size_t GetVectorLength() const

◆ GetVerbosityLevel()

G4int G4PenelopeRayleighModel::GetVerbosityLevel ( )
inline

Definition at line 85 of file G4PenelopeRayleighModel.hh.

85{return fVerboseLevel;};

◆ Initialise()

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

Implements G4VEmModel.

Definition at line 152 of file G4PenelopeRayleighModel.cc.

154{
155 if (fVerboseLevel > 3)
156 G4cout << "Calling G4PenelopeRayleighModel::Initialise()" << G4endl;
157
158 SetParticle(part);
159
160 //Only the master model creates/fills/destroys the tables
161 if (IsMaster() && part == fParticle)
162 {
163 //clear tables depending on materials, not the atomic ones
164 ClearTables();
165
166 if (fVerboseLevel > 3)
167 G4cout << "Calling G4PenelopeRayleighModel::Initialise() [master]" << G4endl;
168
169 //create new tables
170 if (!fLogFormFactorTable)
171 fLogFormFactorTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
172 if (!fPMaxTable)
173 fPMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
174 if (!fSamplingTable)
175 fSamplingTable = new std::map<const G4Material*,G4PenelopeSamplingData*>;
176
177 G4ProductionCutsTable* theCoupleTable =
179
180 for (G4int i=0;i<(G4int)theCoupleTable->GetTableSize();++i)
181 {
182 const G4Material* material =
183 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
184 const G4ElementVector* theElementVector = material->GetElementVector();
185
186 for (std::size_t j=0;j<material->GetNumberOfElements();++j)
187 {
188 G4int iZ = theElementVector->at(j)->GetZasInt();
189 //read data files only in the master
190 if (!fLogAtomicCrossSection[iZ])
191 ReadDataFile(iZ);
192 }
193
194 //1) If the table has not been built for the material, do it!
195 if (!fLogFormFactorTable->count(material))
196 BuildFormFactorTable(material);
197
198 //2) retrieve or build the sampling table
199 if (!(fSamplingTable->count(material)))
200 InitializeSamplingAlgorithm(material);
201
202 //3) retrieve or build the pMax data
203 if (!fPMaxTable->count(material))
204 GetPMaxTable(material);
205 }
206
207 if (fVerboseLevel > 1) {
208 G4cout << "Penelope Rayleigh model v2008 is initialized " << G4endl
209 << "Energy range: "
210 << LowEnergyLimit() / keV << " keV - "
211 << HighEnergyLimit() / GeV << " GeV"
212 << G4endl;
213 }
214 }
215
216 if(fIsInitialised) return;
218 fIsInitialised = true;
219}
std::vector< const G4Element * > G4ElementVector
const G4Material * GetMaterial() const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:185
size_t GetNumberOfElements() const
Definition: G4Material.hh:181
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
std::size_t GetTableSize() const
static G4ProductionCutsTable * GetProductionCutsTable()
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:124
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:634

◆ InitialiseLocal()

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

Reimplemented from G4VEmModel.

Definition at line 223 of file G4PenelopeRayleighModel.cc.

225{
226 if (fVerboseLevel > 3)
227 G4cout << "Calling G4PenelopeRayleighModel::InitialiseLocal()" << G4endl;
228 //
229 //Check that particle matches: one might have multiple master models (e.g.
230 //for e+ and e-).
231 //
232 if (part == fParticle)
233 {
234 //Get the const table pointers from the master to the workers
235 const G4PenelopeRayleighModel* theModel =
236 static_cast<G4PenelopeRayleighModel*> (masterModel);
237
238 //Copy pointers to the data tables
239 for(G4int i=0; i<=fMaxZ; ++i)
240 {
241 fLogAtomicCrossSection[i] = theModel->fLogAtomicCrossSection[i];
242 fAtomicFormFactor[i] = theModel->fAtomicFormFactor[i];
243 }
244 fLogFormFactorTable = theModel->fLogFormFactorTable;
245 fPMaxTable = theModel->fPMaxTable;
246 fSamplingTable = theModel->fSamplingTable;
247
248 //copy the G4DataVector with the grid
249 fLogQSquareGrid = theModel->fLogQSquareGrid;
250
251 //Same verbosity for all workers, as the master
252 fVerboseLevel = theModel->fVerboseLevel;
253 }
254
255 return;
256}

◆ operator=()

G4PenelopeRayleighModel & G4PenelopeRayleighModel::operator= ( const G4PenelopeRayleighModel right)
delete

◆ SampleSecondaries()

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

Implements G4VEmModel.

Definition at line 382 of file G4PenelopeRayleighModel.cc.

387{
388 // Sampling of the Rayleigh final state (namely, scattering angle of the photon)
389 // from the Penelope2008 model. The scattering angle is sampled from the atomic
390 // cross section dOmega/d(cosTheta) from Born ("Atomic Phyisics", 1969), disregarding
391 // anomalous scattering effects. The Form Factor F(Q) function which appears in the
392 // analytical cross section is retrieved via the method GetFSquared(); atomic data
393 // are tabulated for F(Q). Form factor for compounds is calculated according to
394 // the additivity rule. The sampling from the F(Q) is made via a Rational Inverse
395 // Transform with Aliasing (RITA) algorithm; RITA parameters are calculated once
396 // for each material and managed by G4PenelopeSamplingData objects.
397 // The sampling algorithm (rejection method) has efficiency 67% at low energy, and
398 // increases with energy. For E=100 keV the efficiency is 100% and 86% for
399 // hydrogen and uranium, respectively.
400
401 if (fVerboseLevel > 3)
402 G4cout << "Calling SamplingSecondaries() of G4PenelopeRayleighModel" << G4endl;
403
404 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
405
406 if (photonEnergy0 <= fIntrinsicLowEnergyLimit)
407 {
411 return ;
412 }
413
414 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
415
416 const G4Material* theMat = couple->GetMaterial();
417
418 //1) Verify if tables are ready
419 //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
420 //not invoked
421 if (!fPMaxTable || !fSamplingTable || !fLogFormFactorTable)
422 {
423 //create a **thread-local** version of the table. Used only for G4EmCalculator and
424 //Unit Tests
425 fLocalTable = true;
426 if (!fLogFormFactorTable)
427 fLogFormFactorTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
428 if (!fPMaxTable)
429 fPMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
430 if (!fSamplingTable)
431 fSamplingTable = new std::map<const G4Material*,G4PenelopeSamplingData*>;
432 }
433
434 if (!fSamplingTable->count(theMat))
435 {
436 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
437 //not filled up. This can happen in a UnitTest
438 if (fVerboseLevel > 0)
439 {
440 //Issue a G4Exception (warning) only in verbose mode
442 ed << "Unable to find the fSamplingTable data for " <<
443 theMat->GetName() << G4endl;
444 ed << "This can happen only in Unit Tests" << G4endl;
445 G4Exception("G4PenelopeRayleighModel::SampleSecondaries()",
446 "em2019",JustWarning,ed);
447 }
448 const G4ElementVector* theElementVector = theMat->GetElementVector();
449 //protect file reading via autolock
450 G4AutoLock lock(&PenelopeRayleighModelMutex);
451 for (std::size_t j=0;j<theMat->GetNumberOfElements();++j)
452 {
453 G4int iZ = theElementVector->at(j)->GetZasInt();
454 if (!fLogAtomicCrossSection[iZ])
455 {
456 lock.lock();
457 ReadDataFile(iZ);
458 lock.unlock();
459 }
460 }
461 lock.lock();
462 //1) If the table has not been built for the material, do it!
463 if (!fLogFormFactorTable->count(theMat))
464 BuildFormFactorTable(theMat);
465
466 //2) retrieve or build the sampling table
467 if (!(fSamplingTable->count(theMat)))
468 InitializeSamplingAlgorithm(theMat);
469
470 //3) retrieve or build the pMax data
471 if (!fPMaxTable->count(theMat))
472 GetPMaxTable(theMat);
473 lock.unlock();
474 }
475
476 //Ok, restart the job
477 G4PenelopeSamplingData* theDataTable = fSamplingTable->find(theMat)->second;
478 G4PhysicsFreeVector* thePMax = fPMaxTable->find(theMat)->second;
479
480 G4double cosTheta = 1.0;
481
482 //OK, ready to go!
483 G4double qmax = 2.0*photonEnergy0/electron_mass_c2; //this is non-dimensional now
484
485 if (qmax < 1e-10) //very low momentum transfer
486 {
487 G4bool loopAgain=false;
488 do
489 {
490 loopAgain = false;
491 cosTheta = 1.0-2.0*G4UniformRand();
492 G4double G = 0.5*(1+cosTheta*cosTheta);
493 if (G4UniformRand()>G)
494 loopAgain = true;
495 }while(loopAgain);
496 }
497 else //larger momentum transfer
498 {
499 std::size_t nData = theDataTable->GetNumberOfStoredPoints();
500 G4double LastQ2inTheTable = theDataTable->GetX(nData-1);
501 G4double q2max = std::min(qmax*qmax,LastQ2inTheTable);
502
503 G4bool loopAgain = false;
504 G4double MaxPValue = thePMax->Value(photonEnergy0);
505 G4double xx=0;
506
507 //Sampling by rejection method. The rejection function is
508 //G = 0.5*(1+cos^2(theta))
509 //
510 do{
511 loopAgain = false;
512 G4double RandomMax = G4UniformRand()*MaxPValue;
513 xx = theDataTable->SampleValue(RandomMax);
514 //xx is a random value of q^2 in (0,q2max),sampled according to
515 //F(Q^2) via the RITA algorithm
516 if (xx > q2max)
517 loopAgain = true;
518 cosTheta = 1.0-2.0*xx/q2max;
519 G4double G = 0.5*(1+cosTheta*cosTheta);
520 if (G4UniformRand()>G)
521 loopAgain = true;
522 }while(loopAgain);
523 }
524
525 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
526
527 // Scattered photon angles. ( Z - axis along the parent photon)
528 G4double phi = twopi * G4UniformRand() ;
529 G4double dirX = sinTheta*std::cos(phi);
530 G4double dirY = sinTheta*std::sin(phi);
531 G4double dirZ = cosTheta;
532
533 // Update G4VParticleChange for the scattered photon
534 G4ThreeVector photonDirection1(dirX, dirY, dirZ);
535
536 photonDirection1.rotateUz(photonDirection0);
537 fParticleChange->ProposeMomentumDirection(photonDirection1) ;
539
540 return;
541}
@ 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(const G4ThreeVector &Pfinal)
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 84 of file G4PenelopeRayleighModel.hh.

84{fVerboseLevel = lev;};

Member Data Documentation

◆ fParticle

const G4ParticleDefinition* G4PenelopeRayleighModel::fParticle
protected

Definition at line 95 of file G4PenelopeRayleighModel.hh.

Referenced by Initialise(), and InitialiseLocal().

◆ fParticleChange

G4ParticleChangeForGamma* G4PenelopeRayleighModel::fParticleChange
protected

Definition at line 94 of file G4PenelopeRayleighModel.hh.

Referenced by Initialise(), and SampleSecondaries().


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