Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4LivermorePolarizedComptonModel Class Reference

#include <G4LivermorePolarizedComptonModel.hh>

Inheritance diagram for G4LivermorePolarizedComptonModel:
Collaboration diagram for G4LivermorePolarizedComptonModel:

Public Member Functions

 G4LivermorePolarizedComptonModel (const G4ParticleDefinition *p=0, const G4String &nam="LivermorePolarizedCompton")
 
virtual ~G4LivermorePolarizedComptonModel ()
 
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 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 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)
 

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 *)
 
- Protected Attributes inherited from G4VEmModel
G4ElementDatafElementData
 
G4VParticleChangepParticleChange
 
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 53 of file G4LivermorePolarizedComptonModel.hh.

Constructor & Destructor Documentation

G4LivermorePolarizedComptonModel::G4LivermorePolarizedComptonModel ( const G4ParticleDefinition p = 0,
const G4String nam = "LivermorePolarizedCompton" 
)

Definition at line 70 of file G4LivermorePolarizedComptonModel.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 Polarized Compton is constructed " << G4endl;
83 
84  //Mark this model as "applicable" for atomic deexcitation
85  SetDeexcitationFlag(true);
86 
87  fParticleChange = 0;
88  fAtomDeexcitation = 0;
89 }
G4VEmModel(const G4String &nam)
Definition: G4VEmModel.cc:68
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:780

Here is the call graph for this function:

G4LivermorePolarizedComptonModel::~G4LivermorePolarizedComptonModel ( )
virtual

Definition at line 93 of file G4LivermorePolarizedComptonModel.cc.

94 {
95  if(IsMaster()) {
96  delete shellData;
97  shellData = 0;
98  delete profileData;
99  profileData = 0;
100  delete scatterFunctionData;
101  scatterFunctionData = 0;
102  for(G4int i=0; i<maxZ; ++i) {
103  if(data[i]) {
104  delete data[i];
105  data[i] = 0;
106  }
107  }
108  }
109 }
G4bool IsMaster() const
Definition: G4VEmModel.hh:717
int G4int
Definition: G4Types.hh:78
const XML_Char const XML_Char * data
Definition: expat.h:268

Here is the call graph for this function:

Member Function Documentation

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

Reimplemented from G4VEmModel.

Definition at line 253 of file G4LivermorePolarizedComptonModel.cc.

258 {
259  if (verboseLevel > 3)
260  G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePolarizedComptonModel" << G4endl;
261 
262  G4double cs = 0.0;
263 
264  if (GammaEnergy < LowEnergyLimit())
265  return 0.0;
266 
267  G4int intZ = G4lrint(Z);
268  if(intZ < 1 || intZ > maxZ) { return cs; }
269 
270  G4LPhysicsFreeVector* pv = data[intZ];
271 
272  // if element was not initialised
273  // do initialisation safely for MT mode
274  if(!pv)
275  {
276  InitialiseForElement(0, intZ);
277  pv = data[intZ];
278  if(!pv) { return cs; }
279  }
280 
281  G4int n = pv->GetVectorLength() - 1;
282  G4double e1 = pv->Energy(0);
283  G4double e2 = pv->Energy(n);
284 
285  if(GammaEnergy <= e1) { cs = GammaEnergy/(e1*e1)*pv->Value(e1); }
286  else if(GammaEnergy <= e2) { cs = pv->Value(GammaEnergy)/GammaEnergy; }
287  else if(GammaEnergy > e2) { cs = pv->Value(e2)/GammaEnergy; }
288 
289  return cs;
290 
291 }
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:640
size_t GetVectorLength() const
int G4int
Definition: G4Types.hh:78
const XML_Char const XML_Char * data
Definition: expat.h:268
G4GLOB_DLL std::ostream G4cout
G4double Energy(size_t index) const
G4double Value(G4double theEnergy, size_t &lastidx) const
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
int G4lrint(double ad)
Definition: templates.hh:163
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

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

Implements G4VEmModel.

Definition at line 113 of file G4LivermorePolarizedComptonModel.cc.

115 {
116  if (verboseLevel > 1)
117  G4cout << "Calling G4LivermorePolarizedComptonModel::Initialise()" << G4endl;
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 
130  G4int numOfCouples = theCoupleTable->GetTableSize();
131 
132  for(G4int i=0; i<numOfCouples; ++i) {
133  const G4Material* material =
134  theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
135  const G4ElementVector* theElementVector = material->GetElementVector();
136  G4int nelm = material->GetNumberOfElements();
137 
138  for (G4int j=0; j<nelm; ++j) {
139  G4int Z = G4lrint((*theElementVector)[j]->GetZ());
140  if(Z < 1) { Z = 1; }
141  else if(Z > maxZ){ Z = maxZ; }
142 
143  if( (!data[Z]) ) { ReadData(Z, path); }
144  }
145  }
146 
147  // For Doppler broadening
148  if(!shellData) {
149  shellData = new G4ShellData();
150  shellData->SetOccupancyData();
151  G4String file = "/doppler/shell-doppler";
152  shellData->LoadData(file);
153  }
154  if(!profileData) { profileData = new G4DopplerProfile(); }
155 
156  // Scattering Function
157 
158  if(!scatterFunctionData)
159  {
160 
161  G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation;
162  G4String scatterFile = "comp/ce-sf-";
163  scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.);
164  scatterFunctionData->LoadData(scatterFile);
165  }
166 
167  InitialiseElementSelectors(particle, cuts);
168  }
169 
170  if (verboseLevel > 2) {
171  G4cout << "Loaded cross section files" << G4endl;
172  }
173 
174  if( verboseLevel>1 ) {
175  G4cout << "G4LivermoreComptonModel is initialized " << G4endl
176  << "Energy range: "
177  << LowEnergyLimit() / eV << " eV - "
178  << HighEnergyLimit() / GeV << " GeV"
179  << G4endl;
180  }
181  //
182  if(isInitialised) { return; }
183 
184  fParticleChange = GetParticleChangeForGamma();
185  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
186  isInitialised = true;
187 }
void SetOccupancyData()
Definition: G4ShellData.hh:70
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:640
static G4LossTableManager * Instance()
std::vector< G4Element * > G4ElementVector
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:146
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:633
G4bool IsMaster() const
Definition: G4VEmModel.hh:717
virtual G4bool LoadData(const G4String &fileName)
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
void LoadData(const G4String &fileName)
Definition: G4ShellData.cc:234
const XML_Char const XML_Char * data
Definition: expat.h:268
G4GLOB_DLL std::ostream G4cout
static constexpr double eV
Definition: G4SIunits.hh:215
static G4ProductionCutsTable * GetProductionCutsTable()
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
int G4lrint(double ad)
Definition: templates.hh:163
static constexpr double GeV
Definition: G4SIunits.hh:217
#define G4endl
Definition: G4ios.hh:61
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
G4VAtomDeexcitation * AtomDeexcitation()
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:132
const G4Material * GetMaterial() const

Here is the call graph for this function:

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

Reimplemented from G4VEmModel.

Definition at line 875 of file G4LivermorePolarizedComptonModel.cc.

877 {
878  G4AutoLock l(&LivermorePolarizedComptonModelMutex);
879  // G4cout << "G4LivermoreComptonModel::InitialiseForElement Z= "
880  // << Z << G4endl;
881  if(!data[Z]) { ReadData(Z); }
882  l.unlock();
883 }
const XML_Char const XML_Char * data
Definition: expat.h:268

Here is the call graph for this function:

Here is the caller graph for this function:

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

Reimplemented from G4VEmModel.

Definition at line 190 of file G4LivermorePolarizedComptonModel.cc.

192 {
194 }
std::vector< G4EmElementSelector * > * GetElementSelectors()
Definition: G4VEmModel.hh:801
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
Definition: G4VEmModel.hh:809

Here is the call graph for this function:

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

Implements G4VEmModel.

Definition at line 295 of file G4LivermorePolarizedComptonModel.cc.

300 {
301  // The scattered gamma energy is sampled according to Klein - Nishina formula.
302  // The random number techniques of Butcher & Messel are used (Nuc Phys 20(1960),15).
303  // GEANT4 internal units
304  //
305  // Note : Effects due to binding of atomic electrons are negliged.
306 
307  if (verboseLevel > 3)
308  G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedComptonModel" << G4endl;
309 
310  G4double gammaEnergy0 = aDynamicGamma->GetKineticEnergy();
311 
312  // do nothing below the threshold
313  // should never get here because the XS is zero below the limit
314  if (gammaEnergy0 < LowEnergyLimit())
315  return ;
316 
317 
318  G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
319 
320  // Protection: a polarisation parallel to the
321  // direction causes problems;
322  // in that case find a random polarization
323 
324  G4ThreeVector gammaDirection0 = aDynamicGamma->GetMomentumDirection();
325 
326  // Make sure that the polarization vector is perpendicular to the
327  // gamma direction. If not
328 
329  if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
330  { // only for testing now
331  gammaPolarization0 = GetRandomPolarization(gammaDirection0);
332  }
333  else
334  {
335  if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
336  {
337  gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
338  }
339  }
340 
341  // End of Protection
342 
343  G4double E0_m = gammaEnergy0 / electron_mass_c2 ;
344 
345  // Select randomly one element in the current material
346  //G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0);
347  const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
348  const G4Element* elm = SelectRandomAtom(couple,particle,gammaEnergy0);
349  G4int Z = (G4int)elm->GetZ();
350 
351  // Sample the energy and the polarization of the scattered photon
352 
353  G4double epsilon, epsilonSq, onecost, sinThetaSqr, greject ;
354 
355  G4double epsilon0Local = 1./(1. + 2*E0_m);
356  G4double epsilon0Sq = epsilon0Local*epsilon0Local;
357  G4double alpha1 = - std::log(epsilon0Local);
358  G4double alpha2 = 0.5*(1.- epsilon0Sq);
359 
360  G4double wlGamma = h_Planck*c_light/gammaEnergy0;
361  G4double gammaEnergy1;
362  G4ThreeVector gammaDirection1;
363 
364  do {
365  if ( alpha1/(alpha1+alpha2) > G4UniformRand() )
366  {
367  epsilon = G4Exp(-alpha1*G4UniformRand());
368  epsilonSq = epsilon*epsilon;
369  }
370  else
371  {
372  epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4UniformRand();
373  epsilon = std::sqrt(epsilonSq);
374  }
375 
376  onecost = (1.- epsilon)/(epsilon*E0_m);
377  sinThetaSqr = onecost*(2.-onecost);
378 
379  // Protection
380  if (sinThetaSqr > 1.)
381  {
382  G4cout
383  << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
384  << "sin(theta)**2 = "
385  << sinThetaSqr
386  << "; set to 1"
387  << G4endl;
388  sinThetaSqr = 1.;
389  }
390  if (sinThetaSqr < 0.)
391  {
392  G4cout
393  << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
394  << "sin(theta)**2 = "
395  << sinThetaSqr
396  << "; set to 0"
397  << G4endl;
398  sinThetaSqr = 0.;
399  }
400  // End protection
401 
402  G4double x = std::sqrt(onecost/2.) / (wlGamma/cm);;
403  G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
404  greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction;
405 
406  } while(greject < G4UniformRand()*Z);
407 
408 
409  // ****************************************************
410  // Phi determination
411  // ****************************************************
412 
413  G4double phi = SetPhi(epsilon,sinThetaSqr);
414 
415  //
416  // scattered gamma angles. ( Z - axis along the parent gamma)
417  //
418 
419  G4double cosTheta = 1. - onecost;
420 
421  // Protection
422 
423  if (cosTheta > 1.)
424  {
425  G4cout
426  << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
427  << "cosTheta = "
428  << cosTheta
429  << "; set to 1"
430  << G4endl;
431  cosTheta = 1.;
432  }
433  if (cosTheta < -1.)
434  {
435  G4cout
436  << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
437  << "cosTheta = "
438  << cosTheta
439  << "; set to -1"
440  << G4endl;
441  cosTheta = -1.;
442  }
443  // End protection
444 
445 
446  G4double sinTheta = std::sqrt (sinThetaSqr);
447 
448  // Protection
449  if (sinTheta > 1.)
450  {
451  G4cout
452  << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
453  << "sinTheta = "
454  << sinTheta
455  << "; set to 1"
456  << G4endl;
457  sinTheta = 1.;
458  }
459  if (sinTheta < -1.)
460  {
461  G4cout
462  << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
463  << "sinTheta = "
464  << sinTheta
465  << "; set to -1"
466  << G4endl;
467  sinTheta = -1.;
468  }
469  // End protection
470 
471 
472  G4double dirx = sinTheta*std::cos(phi);
473  G4double diry = sinTheta*std::sin(phi);
474  G4double dirz = cosTheta ;
475 
476 
477  // oneCosT , eom
478 
479  // Doppler broadening - Method based on:
480  // Y. Namito, S. Ban and H. Hirayama,
481  // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code"
482  // NIM A 349, pp. 489-494, 1994
483 
484  // Maximum number of sampling iterations
485 
486  static G4int maxDopplerIterations = 1000;
487  G4double bindingE = 0.;
488  G4double photonEoriginal = epsilon * gammaEnergy0;
489  G4double photonE = -1.;
490  G4int iteration = 0;
491  G4double eMax = gammaEnergy0;
492 
493  G4int shellIdx = 0;
494 
495  if (verboseLevel > 3) {
496  G4cout << "Started loop to sample broading" << G4endl;
497  }
498 
499  do
500  {
501  iteration++;
502  // Select shell based on shell occupancy
503  shellIdx = shellData->SelectRandomShell(Z);
504  bindingE = shellData->BindingEnergy(Z,shellIdx);
505 
506  if (verboseLevel > 3) {
507  G4cout << "Shell ID= " << shellIdx
508  << " Ebind(keV)= " << bindingE/keV << G4endl;
509  }
510  eMax = gammaEnergy0 - bindingE;
511 
512  // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
513  G4double pSample = profileData->RandomSelectMomentum(Z,shellIdx);
514 
515  if (verboseLevel > 3) {
516  G4cout << "pSample= " << pSample << G4endl;
517  }
518  // Rescale from atomic units
519  G4double pDoppler = pSample * fine_structure_const;
520  G4double pDoppler2 = pDoppler * pDoppler;
521  G4double var2 = 1. + onecost * E0_m;
522  G4double var3 = var2*var2 - pDoppler2;
523  G4double var4 = var2 - pDoppler2 * cosTheta;
524  G4double var = var4*var4 - var3 + pDoppler2 * var3;
525  if (var > 0.)
526  {
527  G4double varSqrt = std::sqrt(var);
528  G4double scale = gammaEnergy0 / var3;
529  // Random select either root
530  if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;
531  else photonE = (var4 + varSqrt) * scale;
532  }
533  else
534  {
535  photonE = -1.;
536  }
537  } while ( iteration <= maxDopplerIterations &&
538  (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) );
539  //while (iteration <= maxDopplerIterations && photonE > eMax); ???
540 
541 
542  // End of recalculation of photon energy with Doppler broadening
543  // Revert to original if maximum number of iterations threshold has been reached
544  if (iteration >= maxDopplerIterations)
545  {
546  photonE = photonEoriginal;
547  bindingE = 0.;
548  }
549 
550  gammaEnergy1 = photonE;
551 
552  //
553  // update G4VParticleChange for the scattered photon
554  //
555 
556 
557 
558  // gammaEnergy1 = epsilon*gammaEnergy0;
559 
560 
561  // New polarization
562 
563  G4ThreeVector gammaPolarization1 = SetNewPolarization(epsilon,
564  sinThetaSqr,
565  phi,
566  cosTheta);
567 
568  // Set new direction
569  G4ThreeVector tmpDirection1( dirx,diry,dirz );
570  gammaDirection1 = tmpDirection1;
571 
572  // Change reference frame.
573 
574  SystemOfRefChange(gammaDirection0,gammaDirection1,
575  gammaPolarization0,gammaPolarization1);
576 
577  if (gammaEnergy1 > 0.)
578  {
579  fParticleChange->SetProposedKineticEnergy( gammaEnergy1 ) ;
580  fParticleChange->ProposeMomentumDirection( gammaDirection1 );
581  fParticleChange->ProposePolarization( gammaPolarization1 );
582  }
583  else
584  {
585  gammaEnergy1 = 0.;
586  fParticleChange->SetProposedKineticEnergy(0.) ;
587  fParticleChange->ProposeTrackStatus(fStopAndKill);
588  }
589 
590  //
591  // kinematic of the scattered electron
592  //
593 
594  G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;
595 
596  // SI -protection against negative final energy: no e- is created
597  // like in G4LivermoreComptonModel.cc
598  if(ElecKineEnergy < 0.0) {
599  fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0 - gammaEnergy1);
600  return;
601  }
602 
603  // SI - Removed range test
604 
605  G4double ElecMomentum = std::sqrt(ElecKineEnergy*(ElecKineEnergy+2.*electron_mass_c2));
606 
607  G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 -
608  gammaEnergy1 * gammaDirection1) * (1./ElecMomentum));
609 
610  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),ElecDirection.unit(),ElecKineEnergy) ;
611  fvect->push_back(dp);
612 
613  // sample deexcitation
614  //
615 
616  if (verboseLevel > 3) {
617  G4cout << "Started atomic de-excitation " << fAtomDeexcitation << G4endl;
618  }
619 
620  if(fAtomDeexcitation && iteration < maxDopplerIterations) {
621  G4int index = couple->GetIndex();
622  if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
623  size_t nbefore = fvect->size();
625  const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
626  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
627  size_t nafter = fvect->size();
628  if(nafter > nbefore) {
629  for (size_t i=nbefore; i<nafter; ++i) {
630  //Check if there is enough residual energy
631  if (bindingE >= ((*fvect)[i])->GetKineticEnergy())
632  {
633  //Ok, this is a valid secondary: keep it
634  bindingE -= ((*fvect)[i])->GetKineticEnergy();
635  }
636  else
637  {
638  //Invalid secondary: not enough energy to create it!
639  //Keep its energy in the local deposit
640  delete (*fvect)[i];
641  (*fvect)[i]=0;
642  }
643  }
644  }
645  }
646  }
647  //This should never happen
648  if(bindingE < 0.0)
649  G4Exception("G4LivermoreComptonModel::SampleSecondaries()",
650  "em2050",FatalException,"Negative local energy deposit");
651 
652  fParticleChange->ProposeLocalEnergyDeposit(bindingE);
653 
654 }
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:640
static constexpr double h_Planck
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
G4double GetKineticEnergy() const
G4double GetZ() const
Definition: G4Element.hh:131
G4ParticleDefinition * GetDefinition() const
bool isOrthogonal(const Hep3Vector &v, double epsilon=tolerance) const
Definition: SpaceVector.cc:237
int G4int
Definition: G4Types.hh:78
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4int SelectRandomShell(G4int Z) const
Definition: G4ShellData.cc:363
double howOrthogonal(const Hep3Vector &v) const
Definition: SpaceVector.cc:219
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static constexpr double electron_mass_c2
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetMomentumDirection() const
G4double BindingEnergy(G4int Z, G4int shellIndex) const
Definition: G4ShellData.cc:166
static constexpr double cm
Definition: G4SIunits.hh:119
void ProposePolarization(const G4ThreeVector &dir)
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
virtual G4double FindValue(G4double x, G4int componentId=0) const
static constexpr double c_light
const G4ThreeVector & GetPolarization() const
static G4Electron * Electron()
Definition: G4Electron.cc:94
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
static constexpr double fine_structure_const
double mag() const
static constexpr double keV
Definition: G4SIunits.hh:216
double epsilon(double density, double temperature)
G4double RandomSelectMomentum(G4int Z, G4int shellIndex) const
G4AtomicShellEnumerator
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:541

Here is the call graph for this function:


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