Geant4  10.02.p03
G4LivermoreComptonModel Class Reference

#include <G4LivermoreComptonModel.hh>

Inheritance diagram for G4LivermoreComptonModel:
Collaboration diagram for G4LivermoreComptonModel:

Public Member Functions

 G4LivermoreComptonModel (const G4ParticleDefinition *p=0, const G4String &nam="LivermoreCompton")
 
virtual ~G4LivermoreComptonModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)
 
virtual void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel)
 
virtual void InitialiseForElement (const G4ParticleDefinition *, G4int Z)
 
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)
 
- Public Member Functions inherited from G4VEmModel
 G4VEmModel (const G4String &nam)
 
virtual ~G4VEmModel ()
 
virtual void InitialiseForMaterial (const G4ParticleDefinition *, const G4Material *)
 
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, const G4ParticleDefinition *, G4double)
 
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 *> *)
 
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=0)
 
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 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)
 

Private Member Functions

void ReadData (size_t Z, const char *path=0)
 
G4double ComputeScatteringFunction (G4double x, G4int Z)
 
G4LivermoreComptonModeloperator= (const G4LivermoreComptonModel &right)
 
 G4LivermoreComptonModel (const G4LivermoreComptonModel &)
 

Private Attributes

G4bool isInitialised
 
G4int verboseLevel
 
G4ParticleChangeForGamma * fParticleChange
 
G4VAtomDeexcitationfAtomDeexcitation
 

Static Private Attributes

static G4ShellDatashellData = 0
 
static G4DopplerProfileprofileData = 0
 
static G4int maxZ = 99
 
static G4LPhysicsFreeVectordata [100] = {0}
 
static const G4double ScatFuncFitParam [101][9]
 

Additional Inherited Members

- Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLoss * GetParticleChangeForLoss ()
 
G4ParticleChangeForGamma * GetParticleChangeForGamma ()
 
virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
 
const G4MaterialCutsCoupleCurrentCouple () const
 
void SetCurrentElement (const G4Element *)
 
- Protected Attributes inherited from G4VEmModel
G4ElementDatafElementData
 
G4VParticleChange * pParticleChange
 
G4PhysicsTablexSectionTable
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 
size_t idxTable
 
- Static Protected Attributes inherited from G4VEmModel
static const G4double inveplus = 1.0/CLHEP::eplus
 

Detailed Description

Definition at line 48 of file G4LivermoreComptonModel.hh.

Constructor & Destructor Documentation

◆ G4LivermoreComptonModel() [1/2]

G4LivermoreComptonModel::G4LivermoreComptonModel ( const G4ParticleDefinition p = 0,
const G4String nam = "LivermoreCompton" 
)

Definition at line 69 of file G4LivermoreComptonModel.cc.

71  : G4VEmModel(nam),isInitialised(false)
72 {
73  verboseLevel=1 ;
74  // Verbosity scale:
75  // 0 = nothing
76  // 1 = warning for energy non-conservation
77  // 2 = details of energy budget
78  // 3 = calculation of cross sections, file openings, sampling of atoms
79  // 4 = entering in methods
80 
81  if( verboseLevel>1 ) {
82  G4cout << "Livermore Compton model is constructed " << G4endl;
83  }
84 
85  //Mark this model as "applicable" for atomic deexcitation
86  SetDeexcitationFlag(true);
87 
88  fParticleChange = 0;
90 }
G4VEmModel(const G4String &nam)
Definition: G4VEmModel.cc:69
G4GLOB_DLL std::ostream G4cout
G4VAtomDeexcitation * fAtomDeexcitation
G4ParticleChangeForGamma * fParticleChange
#define G4endl
Definition: G4ios.hh:61
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:781
Here is the call graph for this function:

◆ ~G4LivermoreComptonModel()

G4LivermoreComptonModel::~G4LivermoreComptonModel ( )
virtual

Definition at line 94 of file G4LivermoreComptonModel.cc.

95 {
96  if(IsMaster()) {
97  delete shellData;
98  shellData = 0;
99  delete profileData;
100  profileData = 0;
101  for(G4int i=0; i<maxZ; ++i) {
102  if(data[i]) {
103  delete data[i];
104  data[i] = 0;
105  }
106  }
107  }
108 }
G4bool IsMaster() const
Definition: G4VEmModel.hh:718
int G4int
Definition: G4Types.hh:78
static G4DopplerProfile * profileData
static G4ShellData * shellData
static G4LPhysicsFreeVector * data[100]
Here is the call graph for this function:

◆ G4LivermoreComptonModel() [2/2]

G4LivermoreComptonModel::G4LivermoreComptonModel ( const G4LivermoreComptonModel )
private

Member Function Documentation

◆ ComputeCrossSectionPerAtom()

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

Reimplemented from G4VEmModel.

Definition at line 240 of file G4LivermoreComptonModel.cc.

244 {
245  if (verboseLevel > 3) {
246  G4cout << "G4LivermoreComptonModel::ComputeCrossSectionPerAtom()"
247  << G4endl;
248  }
249  G4double cs = 0.0;
250 
251  if (GammaEnergy < LowEnergyLimit()) { return 0.0; }
252 
253  G4int intZ = G4lrint(Z);
254  if(intZ < 1 || intZ > maxZ) { return cs; }
255 
256  G4LPhysicsFreeVector* pv = data[intZ];
257 
258  // if element was not initialised
259  // do initialisation safely for MT mode
260  if(!pv)
261  {
262  InitialiseForElement(0, intZ);
263  pv = data[intZ];
264  if(!pv) { return cs; }
265  }
266 
267  G4int n = pv->GetVectorLength() - 1;
268  G4double e1 = pv->Energy(0);
269  G4double e2 = pv->Energy(n);
270 
271  if(GammaEnergy <= e1) { cs = GammaEnergy/(e1*e1)*pv->Value(e1); }
272  else if(GammaEnergy <= e2) { cs = pv->Value(GammaEnergy)/GammaEnergy; }
273  else if(GammaEnergy > e2) { cs = pv->Value(e2)/GammaEnergy; }
274 
275  return cs;
276 }
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
static const G4double e2
int G4int
Definition: G4Types.hh:78
Char_t n[5]
G4GLOB_DLL std::ostream G4cout
Float_t Z
size_t GetVectorLength() const
G4double Value(G4double theEnergy, size_t &lastidx) const
static const G4double e1
int G4lrint(double ad)
Definition: templates.hh:163
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
#define G4endl
Definition: G4ios.hh:61
G4double Energy(size_t index) const
double G4double
Definition: G4Types.hh:76
static G4LPhysicsFreeVector * data[100]
Here is the call graph for this function:

◆ ComputeScatteringFunction()

G4double G4LivermoreComptonModel::ComputeScatteringFunction ( G4double  x,
G4int  Z 
)
private

Definition at line 539 of file G4LivermoreComptonModel.cc.

540 {
541  G4double value = Z;
542  if (x <= ScatFuncFitParam[Z][2]) {
543 
544  G4double lgq = G4Log(x)/ln10;
545 
546  if (lgq < ScatFuncFitParam[Z][1]) {
547  value = ScatFuncFitParam[Z][3] + lgq*ScatFuncFitParam[Z][4];
548  } else {
549  value = ScatFuncFitParam[Z][5] + lgq*ScatFuncFitParam[Z][6] +
550  lgq*lgq*ScatFuncFitParam[Z][7] + lgq*lgq*lgq*ScatFuncFitParam[Z][8];
551  }
552  value = G4Exp(value*ln10);
553  }
554  //G4cout << " value= " << value << G4endl;
555  return value;
556 }
static const G4double ln10
Float_t Z
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
double G4double
Definition: G4Types.hh:76
static const G4double ScatFuncFitParam[101][9]
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Initialise()

void G4LivermoreComptonModel::Initialise ( const G4ParticleDefinition particle,
const G4DataVector cuts 
)
virtual

Implements G4VEmModel.

Definition at line 112 of file G4LivermoreComptonModel.cc.

114 {
115  if (verboseLevel > 1) {
116  G4cout << "Calling G4LivermoreComptonModel::Initialise()" << G4endl;
117  }
118 
119  // Initialise element selector
120 
121  if(IsMaster()) {
122 
123  // Access to elements
124 
125  char* path = getenv("G4LEDATA");
126 
127  G4ProductionCutsTable* theCoupleTable =
129  G4int numOfCouples = theCoupleTable->GetTableSize();
130 
131  for(G4int i=0; i<numOfCouples; ++i) {
132  const G4Material* material =
133  theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
134  const G4ElementVector* theElementVector = material->GetElementVector();
135  G4int nelm = material->GetNumberOfElements();
136 
137  for (G4int j=0; j<nelm; ++j) {
138  G4int Z = G4lrint((*theElementVector)[j]->GetZ());
139  if(Z < 1) { Z = 1; }
140  else if(Z > maxZ){ Z = maxZ; }
141 
142  if( (!data[Z]) ) { ReadData(Z, path); }
143  }
144  }
145 
146  // For Doppler broadening
147  if(!shellData) {
148  shellData = new G4ShellData();
150  G4String file = "/doppler/shell-doppler";
151  shellData->LoadData(file);
152  }
153  if(!profileData) { profileData = new G4DopplerProfile(); }
154 
155  InitialiseElementSelectors(particle, cuts);
156  }
157 
158  if (verboseLevel > 2) {
159  G4cout << "Loaded cross section files" << G4endl;
160  }
161 
162  if( verboseLevel>1 ) {
163  G4cout << "G4LivermoreComptonModel is initialized " << G4endl
164  << "Energy range: "
165  << LowEnergyLimit() / eV << " eV - "
166  << HighEnergyLimit() / GeV << " GeV"
167  << G4endl;
168  }
169  //
170  if(isInitialised) { return; }
171 
174  isInitialised = true;
175 }
void SetOccupancyData()
Definition: G4ShellData.hh:70
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
void ReadData(size_t Z, const char *path=0)
static G4LossTableManager * Instance()
std::vector< G4Element * > G4ElementVector
const G4Material * GetMaterial() const
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:148
G4bool IsMaster() const
Definition: G4VEmModel.hh:718
TFile * file
int G4int
Definition: G4Types.hh:78
static G4DopplerProfile * profileData
static G4ShellData * shellData
void LoadData(const G4String &fileName)
Definition: G4ShellData.cc:234
string material
Definition: eplot.py:19
G4GLOB_DLL std::ostream G4cout
Float_t Z
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:634
static const double GeV
Definition: G4SIunits.hh:214
G4VAtomDeexcitation * fAtomDeexcitation
static G4ProductionCutsTable * GetProductionCutsTable()
static const double eV
Definition: G4SIunits.hh:212
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
int G4lrint(double ad)
Definition: templates.hh:163
G4ParticleChangeForGamma * fParticleChange
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
#define G4endl
Definition: G4ios.hh:61
G4VAtomDeexcitation * AtomDeexcitation()
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:134
static G4LPhysicsFreeVector * data[100]
Here is the call graph for this function:

◆ InitialiseForElement()

void G4LivermoreComptonModel::InitialiseForElement ( const G4ParticleDefinition ,
G4int  Z 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 565 of file G4LivermoreComptonModel.cc.

567 {
568  G4AutoLock l(&LivermoreComptonModelMutex);
569  // G4cout << "G4LivermoreComptonModel::InitialiseForElement Z= "
570  // << Z << G4endl;
571  if(!data[Z]) { ReadData(Z); }
572  l.unlock();
573 }
void ReadData(size_t Z, const char *path=0)
Float_t Z
static G4LPhysicsFreeVector * data[100]
Here is the call graph for this function:
Here is the caller graph for this function:

◆ InitialiseLocal()

void G4LivermoreComptonModel::InitialiseLocal ( const G4ParticleDefinition ,
G4VEmModel masterModel 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 179 of file G4LivermoreComptonModel.cc.

181 {
183 }
void SetElementSelectors(std::vector< G4EmElementSelector *> *)
Definition: G4VEmModel.hh:810
std::vector< G4EmElementSelector * > * GetElementSelectors()
Definition: G4VEmModel.hh:802
Here is the call graph for this function:

◆ operator=()

G4LivermoreComptonModel& G4LivermoreComptonModel::operator= ( const G4LivermoreComptonModel right)
private

◆ ReadData()

void G4LivermoreComptonModel::ReadData ( size_t  Z,
const char *  path = 0 
)
private

Definition at line 187 of file G4LivermoreComptonModel.cc.

188 {
189  if (verboseLevel > 1)
190  {
191  G4cout << "G4LivermoreComptonModel::ReadData()"
192  << G4endl;
193  }
194  if(data[Z]) { return; }
195  const char* datadir = path;
196  if(!datadir)
197  {
198  datadir = getenv("G4LEDATA");
199  if(!datadir)
200  {
201  G4Exception("G4LivermoreComptonModel::ReadData()",
202  "em0006",FatalException,
203  "Environment variable G4LEDATA not defined");
204  return;
205  }
206  }
207 
208  data[Z] = new G4LPhysicsFreeVector();
209 
210  // Activation of spline interpolation
211  data[Z]->SetSpline(false);
212 
213  std::ostringstream ost;
214  ost << datadir << "/livermore/comp/ce-cs-" << Z <<".dat";
215  std::ifstream fin(ost.str().c_str());
216 
217  if( !fin.is_open())
218  {
220  ed << "G4LivermoreComptonModel data file <" << ost.str().c_str()
221  << "> is not opened!" << G4endl;
222  G4Exception("G4LivermoreComptonModel::ReadData()",
223  "em0003",FatalException,
224  ed,"G4LEDATA version should be G4EMLOW6.34 or later");
225  return;
226  } else {
227  if(verboseLevel > 3) {
228  G4cout << "File " << ost.str()
229  << " is opened by G4LivermoreComptonModel" << G4endl;
230  }
231  data[Z]->Retrieve(fin, true);
233  }
234  fin.close();
235 }
TString fin
static const double MeV
Definition: G4SIunits.hh:211
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
void SetSpline(G4bool)
G4GLOB_DLL std::ostream G4cout
Float_t Z
virtual void ScaleVector(G4double factorE, G4double factorV)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
virtual G4bool Retrieve(std::ifstream &fIn, G4bool ascii=false)
#define G4endl
Definition: G4ios.hh:61
static const double barn
Definition: G4SIunits.hh:104
static G4LPhysicsFreeVector * data[100]
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SampleSecondaries()

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

Implements G4VEmModel.

Definition at line 281 of file G4LivermoreComptonModel.cc.

286 {
287 
288  // The scattered gamma energy is sampled according to Klein - Nishina
289  // formula then accepted or rejected depending on the Scattering Function
290  // multiplied by factor from Klein - Nishina formula.
291  // Expression of the angular distribution as Klein Nishina
292  // angular and energy distribution and Scattering fuctions is taken from
293  // D. E. Cullen "A simple model of photon transport" Nucl. Instr. Meth.
294  // Phys. Res. B 101 (1995). Method of sampling with form factors is different
295  // data are interpolated while in the article they are fitted.
296  // Reference to the article is from J. Stepanek New Photon, Positron
297  // and Electron Interaction Data for GEANT in Energy Range from 1 eV to 10
298  // TeV (draft).
299  // The random number techniques of Butcher & Messel are used
300  // (Nucl Phys 20(1960),15).
301 
302  G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
303 
304  if (verboseLevel > 3) {
305  G4cout << "G4LivermoreComptonModel::SampleSecondaries() E(MeV)= "
306  << photonEnergy0/MeV << " in " << couple->GetMaterial()->GetName()
307  << G4endl;
308  }
309 
310  // do nothing below the threshold
311  // should never get here because the XS is zero below the limit
312  if (photonEnergy0 < LowEnergyLimit())
313  return ;
314 
315  G4double e0m = photonEnergy0 / electron_mass_c2 ;
316  G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
317 
318  // Select randomly one element in the current material
319  const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
320  const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
321 
322  G4int Z = G4lrint(elm->GetZ());
323 
324  G4double epsilon0Local = 1. / (1. + 2. * e0m);
325  G4double epsilon0Sq = epsilon0Local * epsilon0Local;
326  G4double alpha1 = -G4Log(epsilon0Local);
327  G4double alpha2 = 0.5 * (1. - epsilon0Sq);
328 
329  G4double wlPhoton = h_Planck*c_light/photonEnergy0;
330 
331  // Sample the energy of the scattered photon
333  G4double epsilonSq;
334  G4double oneCosT;
335  G4double sinT2;
336  G4double gReject;
337 
338  if (verboseLevel > 3) {
339  G4cout << "Started loop to sample gamma energy" << G4endl;
340  }
341 
342  do {
343  if ( alpha1/(alpha1+alpha2) > G4UniformRand())
344  {
345  epsilon = G4Exp(-alpha1 * G4UniformRand());
346  epsilonSq = epsilon * epsilon;
347  }
348  else
349  {
350  epsilonSq = epsilon0Sq + (1. - epsilon0Sq) * G4UniformRand();
351  epsilon = std::sqrt(epsilonSq);
352  }
353 
354  oneCosT = (1. - epsilon) / ( epsilon * e0m);
355  sinT2 = oneCosT * (2. - oneCosT);
356  G4double x = std::sqrt(oneCosT/2.) * cm / wlPhoton;
357  G4double scatteringFunction = ComputeScatteringFunction(x, Z);
358  gReject = (1. - epsilon * sinT2 / (1. + epsilonSq)) * scatteringFunction;
359 
360  } while(gReject < G4UniformRand()*Z);
361 
362  G4double cosTheta = 1. - oneCosT;
363  G4double sinTheta = std::sqrt (sinT2);
364  G4double phi = twopi * G4UniformRand() ;
365  G4double dirx = sinTheta * std::cos(phi);
366  G4double diry = sinTheta * std::sin(phi);
367  G4double dirz = cosTheta ;
368 
369  // Doppler broadening - Method based on:
370  // Y. Namito, S. Ban and H. Hirayama,
371  // "Implementation of the Doppler Broadening of a Compton-Scattered Photon
372  // into the EGS4 Code", NIM A 349, pp. 489-494, 1994
373 
374  // Maximum number of sampling iterations
375  static G4int maxDopplerIterations = 1000;
376  G4double bindingE = 0.;
377  G4double photonEoriginal = epsilon * photonEnergy0;
378  G4double photonE = -1.;
379  G4int iteration = 0;
380  G4double eMax = photonEnergy0;
381 
382  G4int shellIdx = 0;
383 
384  if (verboseLevel > 3) {
385  G4cout << "Started loop to sample broading" << G4endl;
386  }
387 
388  do
389  {
390  ++iteration;
391  // Select shell based on shell occupancy
392  shellIdx = shellData->SelectRandomShell(Z);
393  bindingE = shellData->BindingEnergy(Z,shellIdx);
394 
395  if (verboseLevel > 3) {
396  G4cout << "Shell ID= " << shellIdx
397  << " Ebind(keV)= " << bindingE/keV << G4endl;
398  }
399 
400  eMax = photonEnergy0 - bindingE;
401 
402  // Randomly sample bound electron momentum
403  // (memento: the data set is in Atomic Units)
404  G4double pSample = profileData->RandomSelectMomentum(Z,shellIdx);
405  if (verboseLevel > 3) {
406  G4cout << "pSample= " << pSample << G4endl;
407  }
408  // Rescale from atomic units
409  G4double pDoppler = pSample * fine_structure_const;
410  G4double pDoppler2 = pDoppler * pDoppler;
411  G4double var2 = 1. + oneCosT * e0m;
412  G4double var3 = var2*var2 - pDoppler2;
413  G4double var4 = var2 - pDoppler2 * cosTheta;
414  G4double var = var4*var4 - var3 + pDoppler2 * var3;
415  if (var > 0.)
416  {
417  G4double varSqrt = std::sqrt(var);
418  G4double scale = photonEnergy0 / var3;
419  // Random select either root
420  if (G4UniformRand() < 0.5) { photonE = (var4 - varSqrt) * scale; }
421  else { photonE = (var4 + varSqrt) * scale; }
422  }
423  else
424  {
425  photonE = -1.;
426  }
427  } while (iteration <= maxDopplerIterations && photonE > eMax);
428  // (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) );
429 
430 
431  // End of recalculation of photon energy with Doppler broadening
432  // Revert to original if maximum number of iterations threshold
433  // has been reached
434 
435  if (iteration >= maxDopplerIterations)
436  {
437  photonE = photonEoriginal;
438  bindingE = 0.;
439  }
440 
441  // Update G4VParticleChange for the scattered photon
442 
443  G4ThreeVector photonDirection1(dirx,diry,dirz);
444  photonDirection1.rotateUz(photonDirection0);
445  fParticleChange->ProposeMomentumDirection(photonDirection1) ;
446 
447  G4double photonEnergy1 = photonE;
448 
449  if (photonEnergy1 > 0.) {
450  fParticleChange->SetProposedKineticEnergy(photonEnergy1) ;
451 
452  } else {
453  // photon absorbed
454  photonEnergy1 = 0.;
455  fParticleChange->SetProposedKineticEnergy(0.) ;
456  fParticleChange->ProposeTrackStatus(fStopAndKill);
457  fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
458  return;
459  }
460 
461  // Kinematics of the scattered electron
462  G4double eKineticEnergy = photonEnergy0 - photonEnergy1 - bindingE;
463 
464  // protection against negative final energy: no e- is created
465  if(eKineticEnergy < 0.0) {
466  fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0 - photonEnergy1);
467  return;
468  }
469 
470  G4double eTotalEnergy = eKineticEnergy + electron_mass_c2;
471 
472  G4double electronE = photonEnergy0 * (1. - epsilon) + electron_mass_c2;
473  G4double electronP2 =
474  electronE*electronE - electron_mass_c2*electron_mass_c2;
475  G4double sinThetaE = -1.;
476  G4double cosThetaE = 0.;
477  if (electronP2 > 0.)
478  {
479  cosThetaE = (eTotalEnergy + photonEnergy1 )*
480  (1. - epsilon) / std::sqrt(electronP2);
481  sinThetaE = -1. * sqrt(1. - cosThetaE * cosThetaE);
482  }
483 
484  G4double eDirX = sinThetaE * std::cos(phi);
485  G4double eDirY = sinThetaE * std::sin(phi);
486  G4double eDirZ = cosThetaE;
487 
488  G4ThreeVector eDirection(eDirX,eDirY,eDirZ);
489  eDirection.rotateUz(photonDirection0);
491  eDirection,eKineticEnergy) ;
492  fvect->push_back(dp);
493 
494  // sample deexcitation
495  //
496 
497  if (verboseLevel > 3) {
498  G4cout << "Started atomic de-excitation " << fAtomDeexcitation << G4endl;
499  }
500 
501  if(fAtomDeexcitation && iteration < maxDopplerIterations) {
502  G4int index = couple->GetIndex();
504  size_t nbefore = fvect->size();
506  const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
507  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
508  size_t nafter = fvect->size();
509  if(nafter > nbefore) {
510  for (size_t i=nbefore; i<nafter; ++i) {
511  //Check if there is enough residual energy
512  if (bindingE >= ((*fvect)[i])->GetKineticEnergy())
513  {
514  //Ok, this is a valid secondary: keep it
515  bindingE -= ((*fvect)[i])->GetKineticEnergy();
516  }
517  else
518  {
519  //Invalid secondary: not enough energy to create it!
520  //Keep its energy in the local deposit
521  delete (*fvect)[i];
522  (*fvect)[i]=0;
523  }
524  }
525  }
526  }
527  }
528  //This should never happen
529  if(bindingE < 0.0)
530  G4Exception("G4LivermoreComptonModel::SampleSecondaries()",
531  "em2050",FatalException,"Negative local energy deposit");
532 
533  fParticleChange->ProposeLocalEnergyDeposit(bindingE);
534 }
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
static const double cm
Definition: G4SIunits.hh:118
static const double MeV
Definition: G4SIunits.hh:211
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
G4double ComputeScatteringFunction(G4double x, G4int Z)
const G4Material * GetMaterial() const
Int_t index
float h_Planck
Definition: hepunit.py:263
int G4int
Definition: G4Types.hh:78
static G4DopplerProfile * profileData
static G4ShellData * shellData
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
int fine_structure_const
Definition: hepunit.py:287
G4double RandomSelectMomentum(G4int Z, G4int shellIndex) const
G4double GetKineticEnergy() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
Float_t Z
Double_t scale
static const double twopi
Definition: G4SIunits.hh:75
float electron_mass_c2
Definition: hepunit.py:274
G4double BindingEnergy(G4int Z, G4int shellIndex) const
Definition: G4ShellData.cc:166
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4VAtomDeexcitation * fAtomDeexcitation
int G4lrint(double ad)
Definition: templates.hh:163
const G4ThreeVector & GetMomentumDirection() const
G4ParticleChangeForGamma * fParticleChange
void GenerateParticles(std::vector< G4DynamicParticle *> *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4ParticleDefinition * GetDefinition() const
#define G4endl
Definition: G4ios.hh:61
static const double keV
Definition: G4SIunits.hh:213
double G4double
Definition: G4Types.hh:76
const G4String & GetName() const
Definition: G4Material.hh:178
G4int SelectRandomShell(G4int Z) const
Definition: G4ShellData.cc:363
G4double GetZ() const
Definition: G4Element.hh:131
double epsilon(double density, double temperature)
float c_light
Definition: hepunit.py:257
G4AtomicShellEnumerator
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:544
Here is the call graph for this function:

Member Data Documentation

◆ data

G4LPhysicsFreeVector * G4LivermoreComptonModel::data = {0}
staticprivate

Definition at line 96 of file G4LivermoreComptonModel.hh.

◆ fAtomDeexcitation

G4VAtomDeexcitation* G4LivermoreComptonModel::fAtomDeexcitation
private

Definition at line 90 of file G4LivermoreComptonModel.hh.

◆ fParticleChange

G4ParticleChangeForGamma* G4LivermoreComptonModel::fParticleChange
private

Definition at line 89 of file G4LivermoreComptonModel.hh.

◆ isInitialised

G4bool G4LivermoreComptonModel::isInitialised
private

Definition at line 86 of file G4LivermoreComptonModel.hh.

◆ maxZ

G4int G4LivermoreComptonModel::maxZ = 99
staticprivate

Definition at line 95 of file G4LivermoreComptonModel.hh.

◆ profileData

G4DopplerProfile * G4LivermoreComptonModel::profileData = 0
staticprivate

Definition at line 93 of file G4LivermoreComptonModel.hh.

◆ ScatFuncFitParam

const G4double G4LivermoreComptonModel::ScatFuncFitParam
staticprivate

Definition at line 98 of file G4LivermoreComptonModel.hh.

◆ shellData

G4ShellData * G4LivermoreComptonModel::shellData = 0
staticprivate

Definition at line 92 of file G4LivermoreComptonModel.hh.

◆ verboseLevel

G4int G4LivermoreComptonModel::verboseLevel
private

Definition at line 87 of file G4LivermoreComptonModel.hh.


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