Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PenelopeRayleighModel Class Reference

#include <G4PenelopeRayleighModel.hh>

Inheritance diagram for G4PenelopeRayleighModel:
Collaboration 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 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 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)
 
G4int SelectIsotopeNumber (const G4Element *)
 
const G4ElementSelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, 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 *)
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=nullptr)
 
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
 
G4ElementDataGetElementData ()
 
G4PhysicsTableGetCrossSectionTable ()
 
G4VEmFluctuationModelGetModelOfFluctuations ()
 
G4VEmAngularDistributionGetAngularDistribution ()
 
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)
 
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
 
G4double MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
 
const G4StringGetName () const
 
void SetCurrentCouple (const G4MaterialCutsCouple *)
 
const G4ElementGetCurrentElement () const
 
const G4IsotopeGetCurrentIsotope () const
 
G4bool IsLocked () const
 
void SetLocked (G4bool)
 

Protected Attributes

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

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 *)
 
- Static Protected Attributes inherited from G4VEmModel
static const G4double inveplus = 1.0/CLHEP::eplus
 

Detailed Description

Definition at line 58 of file G4PenelopeRayleighModel.hh.

Constructor & Destructor Documentation

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

Definition at line 54 of file G4PenelopeRayleighModel.cc.

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

Here is the call graph for this function:

G4PenelopeRayleighModel::~G4PenelopeRayleighModel ( )
virtual

Definition at line 96 of file G4PenelopeRayleighModel.cc.

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

Here is the call graph for this function:

Member Function Documentation

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 274 of file G4PenelopeRayleighModel.cc.

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

Here is the call graph for this function:

void G4PenelopeRayleighModel::DumpFormFactorTable ( const G4Material mat)

Definition at line 1333 of file G4PenelopeRayleighModel.cc.

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

Here is the call graph for this function:

G4int G4PenelopeRayleighModel::GetVerbosityLevel ( )
inline

Definition at line 85 of file G4PenelopeRayleighModel.hh.

85 {return verboseLevel;};
void G4PenelopeRayleighModel::Initialise ( const G4ParticleDefinition part,
const G4DataVector  
)
virtual

Implements G4VEmModel.

Definition at line 157 of file G4PenelopeRayleighModel.cc.

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

Here is the call graph for this function:

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

Reimplemented from G4VEmModel.

Definition at line 238 of file G4PenelopeRayleighModel.cc.

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

Implements G4VEmModel.

Definition at line 416 of file G4PenelopeRayleighModel.cc.

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

Here is the call graph for this function:

void G4PenelopeRayleighModel::SetVerbosityLevel ( G4int  lev)
inline

Definition at line 84 of file G4PenelopeRayleighModel.hh.

84 {verboseLevel = lev;};

Member Data Documentation

const G4ParticleDefinition* G4PenelopeRayleighModel::fParticle
protected

Definition at line 92 of file G4PenelopeRayleighModel.hh.

G4ParticleChangeForGamma* G4PenelopeRayleighModel::fParticleChange
protected

Definition at line 91 of file G4PenelopeRayleighModel.hh.


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