Geant4  10.02.p03
G4LivermorePolarizedGammaConversionModel Class Reference

#include <G4LivermorePolarizedGammaConversionModel.hh>

Inheritance diagram for G4LivermorePolarizedGammaConversionModel:
Collaboration diagram for G4LivermorePolarizedGammaConversionModel:

Public Member Functions

 G4LivermorePolarizedGammaConversionModel (const G4ParticleDefinition *p=0, const G4String &nam="LivermorePolarizedGammaConversion")
 
virtual ~G4LivermorePolarizedGammaConversionModel ()
 
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)
 
virtual G4double MinPrimaryEnergy (const G4Material *, const G4ParticleDefinition *, G4double)
 
- 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 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 ScreenFunction1 (G4double screenVariable)
 
G4double ScreenFunction2 (G4double screenVariable)
 
G4ThreeVector GetRandomPolarization (G4ThreeVector &direction0)
 
G4ThreeVector GetPerpendicularPolarization (const G4ThreeVector &direction0, const G4ThreeVector &polarization0) const
 
G4ThreeVector SetPerpendicularVector (G4ThreeVector &a)
 
void SystemOfRefChange (G4ThreeVector &direction0, G4ThreeVector &direction1, G4ThreeVector &polarization0)
 
G4double SetPhi (G4double)
 
G4double SetPsi (G4double, G4double)
 
void SetTheta (G4double *, G4double *, G4double)
 
G4double Poli (G4double, G4double, G4double, G4double)
 
G4double Fln (G4double, G4double, G4double)
 
G4double Encu (G4double *, G4double *, G4double)
 
G4double Flor (G4double *, G4double)
 
G4double Glor (G4double *, G4double)
 
G4double Fdlor (G4double *, G4double)
 
G4double Fintlor (G4double *, G4double)
 
G4double Finvlor (G4double *, G4double, G4double)
 
G4double Ftan (G4double *, G4double)
 
G4double Fdtan (G4double *, G4double)
 
G4double Finttan (G4double *, G4double)
 
G4double Finvtan (G4double *, G4double, G4double)
 
G4LivermorePolarizedGammaConversionModeloperator= (const G4LivermorePolarizedGammaConversionModel &right)
 
 G4LivermorePolarizedGammaConversionModel (const G4LivermorePolarizedGammaConversionModel &)
 

Private Attributes

G4ParticleChangeForGamma * fParticleChange
 
G4double lowEnergyLimit
 
G4bool isInitialised
 
G4int verboseLevel
 
G4double smallEnergy
 
G4double Psi
 
G4double Phi
 

Static Private Attributes

static G4int maxZ = 99
 
static G4LPhysicsFreeVectordata [100] = {0}
 

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 43 of file G4LivermorePolarizedGammaConversionModel.hh.

Constructor & Destructor Documentation

◆ G4LivermorePolarizedGammaConversionModel() [1/2]

G4LivermorePolarizedGammaConversionModel::G4LivermorePolarizedGammaConversionModel ( const G4ParticleDefinition p = 0,
const G4String nam = "LivermorePolarizedGammaConversion" 
)

Definition at line 51 of file G4LivermorePolarizedGammaConversionModel.cc.

53  :G4VEmModel(nam), isInitialised(false),smallEnergy(2.*MeV)
54 {
55  fParticleChange = nullptr;
57 
58  Phi=0.;
59  Psi=0.;
60 
61  verboseLevel= 0;
62  // Verbosity scale:
63  // 0 = nothing
64  // 1 = calculation of cross sections, file openings, samping of atoms
65  // 2 = entering in methods
66 
67  if(verboseLevel > 0) {
68  G4cout << "Livermore Polarized GammaConversion is constructed "
69  << G4endl;
70  }
71 
72 }
static const double MeV
Definition: G4SIunits.hh:211
G4VEmModel(const G4String &nam)
Definition: G4VEmModel.cc:69
G4GLOB_DLL std::ostream G4cout
float electron_mass_c2
Definition: hepunit.py:274
#define G4endl
Definition: G4ios.hh:61

◆ ~G4LivermorePolarizedGammaConversionModel()

G4LivermorePolarizedGammaConversionModel::~G4LivermorePolarizedGammaConversionModel ( )
virtual

Definition at line 76 of file G4LivermorePolarizedGammaConversionModel.cc.

77 {
78  if(IsMaster()) {
79  for(G4int i=0; i<maxZ; ++i) {
80  if(data[i]) {
81  delete data[i];
82  data[i] = 0;
83  }
84  }
85  }
86 }
G4bool IsMaster() const
Definition: G4VEmModel.hh:718
int G4int
Definition: G4Types.hh:78
Here is the call graph for this function:

◆ G4LivermorePolarizedGammaConversionModel() [2/2]

G4LivermorePolarizedGammaConversionModel::G4LivermorePolarizedGammaConversionModel ( const G4LivermorePolarizedGammaConversionModel )
private

Member Function Documentation

◆ ComputeCrossSectionPerAtom()

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

Reimplemented from G4VEmModel.

Definition at line 221 of file G4LivermorePolarizedGammaConversionModel.cc.

226 {
227  if (verboseLevel > 1) {
228  G4cout << "G4LivermorePolarizedGammaConversionModel::ComputeCrossSectionPerAtom()"
229  << G4endl;
230  }
231  if (GammaEnergy < lowEnergyLimit) { return 0.0; }
232 
233  G4double xs = 0.0;
234 
235  G4int intZ=G4int(Z);
236 
237  if(intZ < 1 || intZ > maxZ) { return xs; }
238 
239  G4LPhysicsFreeVector* pv = data[intZ];
240 
241  // if element was not initialised
242  // do initialisation safely for MT mode
243  if(!pv)
244  {
245  InitialiseForElement(0, intZ);
246  pv = data[intZ];
247  if(!pv) { return xs; }
248  }
249  // x-section is taken from the table
250  xs = pv->Value(GammaEnergy);
251 
252  if(verboseLevel > 0)
253  {
254  G4int n = pv->GetVectorLength() - 1;
255  G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)="
256  << GammaEnergy/MeV << G4endl;
257  G4cout << " cs (Geant4 internal unit)=" << xs << G4endl;
258  G4cout << " -> first cs value in EADL data file (iu) =" << (*pv)[0] << G4endl;
259  G4cout << " -> last cs value in EADL data file (iu) =" << (*pv)[n] << G4endl;
260  G4cout << "*********************************************************" << G4endl;
261  }
262 
263  return xs;
264 }
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
static const double MeV
Definition: G4SIunits.hh:211
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
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:

◆ Encu()

G4double G4LivermorePolarizedGammaConversionModel::Encu ( G4double p_p1,
G4double p_p2,
G4double  x0 
)
private

Definition at line 856 of file G4LivermorePolarizedGammaConversionModel.cc.

857 {
858  G4int i=0;
859  G4double fx = 1.;
860  G4double x = x0;
861  const G4double xmax = 3.0;
862 
863  for(i=0; i<100; ++i)
864  {
865  fx = (Flor(p_p1,x)*Glor(p_p1,x) - Ftan(p_p2, x))/
866  (Fdlor(p_p1,x) - Fdtan(p_p2,x));
867  x -= fx;
868  if(x > xmax) { return xmax; }
869  // x -= (Flor(p_p1, x)*Glor(p_p1,x) - Ftan(p_p2, x))/
870  // (Fdlor(p_p1,x) - Fdtan(p_p2,x));
871  // fx = Flor(p_p1,x)*Glor(p_p1,x) - Ftan(p_p2, x);
872  // G4cout << std::fabs(fx) << " " << i << " " << x << "dentro ENCU " << G4endl;
873  if(std::fabs(fx) <= x*1.0e-6) { break; }
874  }
875 
876  if(x < 0.0) { x = 0.0; }
877  return x;
878 }
int G4int
Definition: G4Types.hh:78
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Fdlor()

G4double G4LivermorePolarizedGammaConversionModel::Fdlor ( G4double p_p1,
G4double  x 
)
private

Definition at line 907 of file G4LivermorePolarizedGammaConversionModel.cc.

908 {
909  G4double value =0.;
910  //G4double y0 = p_p1[0];
911  G4double A = p_p1[1];
912  G4double w = p_p1[2];
913  G4double xc = p_p1[3];
914 
915  value = (-16.*A*w*(x-xc))/
916  (pi*(w*w+4.*(x-xc)*(x-xc))*(w*w+4.*(x-xc)*(x-xc)));
917  return value;
918 }
double A(double temperature)
static const double pi
Definition: G4SIunits.hh:74
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Fdtan()

G4double G4LivermorePolarizedGammaConversionModel::Fdtan ( G4double p_p1,
G4double  x 
)
private

Definition at line 961 of file G4LivermorePolarizedGammaConversionModel.cc.

962 {
963  G4double value =0.;
964  G4double a = p_p1[0];
965  G4double b = p_p1[1];
966 
967  value = -1.*a / ((x-b)*(x-b));
968  return value;
969 }
double G4double
Definition: G4Types.hh:76
Here is the caller graph for this function:

◆ Fintlor()

G4double G4LivermorePolarizedGammaConversionModel::Fintlor ( G4double p_p1,
G4double  x 
)
private

Definition at line 921 of file G4LivermorePolarizedGammaConversionModel.cc.

922 {
923  G4double value =0.;
924  G4double y0 = p_p1[0];
925  G4double A = p_p1[1];
926  G4double w = p_p1[2];
927  G4double xc = p_p1[3];
928 
929  value = y0*x + A*atan( 2*(x-xc)/w) / pi;
930  return value;
931 }
double A(double temperature)
static const double pi
Definition: G4SIunits.hh:74
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Finttan()

G4double G4LivermorePolarizedGammaConversionModel::Finttan ( G4double p_p1,
G4double  x 
)
private

Definition at line 972 of file G4LivermorePolarizedGammaConversionModel.cc.

973 {
974  G4double value =0.;
975  G4double a = p_p1[0];
976  G4double b = p_p1[1];
977 
978 
979  value = a*log(b-x);
980  return value;
981 }
double G4double
Definition: G4Types.hh:76
Here is the caller graph for this function:

◆ Finvlor()

G4double G4LivermorePolarizedGammaConversionModel::Finvlor ( G4double p_p1,
G4double  x,
G4double  r 
)
private

Definition at line 934 of file G4LivermorePolarizedGammaConversionModel.cc.

935 {
936  G4double value = 0.;
937  G4double nor = 0.;
938  //G4double y0 = p_p1[0];
939  // G4double A = p_p1[1];
940  G4double w = p_p1[2];
941  G4double xc = p_p1[3];
942 
943  nor = atan(2.*(pi-xc)/w)/(2.*pi*w) - atan(2.*(x-xc)/w)/(2.*pi*w);
944  value = xc - (w/2.)*tan(-2.*r*nor*pi*w+atan(2*(xc-x)/w));
945 
946  return value;
947 }
static const double pi
Definition: G4SIunits.hh:74
double G4double
Definition: G4Types.hh:76
Here is the caller graph for this function:

◆ Finvtan()

G4double G4LivermorePolarizedGammaConversionModel::Finvtan ( G4double p_p1,
G4double  cnor,
G4double  r 
)
private

Definition at line 984 of file G4LivermorePolarizedGammaConversionModel.cc.

985 {
986  G4double value =0.;
987  G4double a = p_p1[0];
988  G4double b = p_p1[1];
989 
990  value = b*(1-G4Exp(r*cnor/a));
991 
992  return value;
993 }
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Fln()

G4double G4LivermorePolarizedGammaConversionModel::Fln ( G4double  a,
G4double  b,
G4double  x 
)
private

Definition at line 840 of file G4LivermorePolarizedGammaConversionModel.cc.

841 {
842  G4double value=0.;
843  if(x>0.)
844  {
845  value =(a*log(x)-b);
846  }
847  else
848  {
849  //G4cout << "ERROR in Fln! " << G4endl;
850  }
851  return value;
852 }
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Flor()

G4double G4LivermorePolarizedGammaConversionModel::Flor ( G4double p_p1,
G4double  x 
)
private

Definition at line 881 of file G4LivermorePolarizedGammaConversionModel.cc.

882 {
883  G4double value =0.;
884  // G4double y0 = p_p1[0];
885  // G4double A = p_p1[1];
886  G4double w = p_p1[2];
887  G4double xc = p_p1[3];
888 
889  value = 1./(pi*(w*w + 4.*(x-xc)*(x-xc)));
890  return value;
891 }
static const double pi
Definition: G4SIunits.hh:74
double G4double
Definition: G4Types.hh:76
Here is the caller graph for this function:

◆ Ftan()

G4double G4LivermorePolarizedGammaConversionModel::Ftan ( G4double p_p1,
G4double  x 
)
private

Definition at line 950 of file G4LivermorePolarizedGammaConversionModel.cc.

951 {
952  G4double value =0.;
953  G4double a = p_p1[0];
954  G4double b = p_p1[1];
955 
956  value = a /(x-b);
957  return value;
958 }
double G4double
Definition: G4Types.hh:76
Here is the caller graph for this function:

◆ GetPerpendicularPolarization()

G4ThreeVector G4LivermorePolarizedGammaConversionModel::GetPerpendicularPolarization ( const G4ThreeVector direction0,
const G4ThreeVector polarization0 
) const
private

Definition at line 1043 of file G4LivermorePolarizedGammaConversionModel.cc.

1044 {
1045 
1046  //
1047  // The polarization of a photon is always perpendicular to its momentum direction.
1048  // Therefore this function removes those vector component of gammaPolarization, which
1049  // points in direction of gammaDirection
1050  //
1051  // Mathematically we search the projection of the vector a on the plane E, where n is the
1052  // plains normal vector.
1053  // The basic equation can be found in each geometry book (e.g. Bronstein):
1054  // p = a - (a o n)/(n o n)*n
1055 
1056  return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
1057 }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetRandomPolarization()

G4ThreeVector G4LivermorePolarizedGammaConversionModel::GetRandomPolarization ( G4ThreeVector direction0)
private

Definition at line 1017 of file G4LivermorePolarizedGammaConversionModel.cc.

1018 {
1019  G4ThreeVector d0 = direction0.unit();
1020  G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
1021  G4ThreeVector a0 = a1.unit(); // unit vector
1022 
1023  G4double rand1 = G4UniformRand();
1024 
1025  G4double angle = twopi*rand1; // random polar angle
1026  G4ThreeVector b0 = d0.cross(a0); // cross product
1027 
1028  G4ThreeVector c;
1029 
1030  c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
1031  c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
1032  c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
1033 
1034  G4ThreeVector c0 = c.unit();
1035 
1036  return c0;
1037 
1038 }
const G4double a0
static const G4double a1
static G4double angle[DIM]
void setY(double)
void setZ(double)
void setX(double)
#define G4UniformRand()
Definition: Randomize.hh:97
Hep3Vector cross(const Hep3Vector &) const
Hep3Vector unit() const
static const double twopi
Definition: G4SIunits.hh:75
double x() const
double y() const
static const G4double c0
double z() const
double G4double
Definition: G4Types.hh:76
static const G4double b0
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Glor()

G4double G4LivermorePolarizedGammaConversionModel::Glor ( G4double p_p1,
G4double  x 
)
private

Definition at line 894 of file G4LivermorePolarizedGammaConversionModel.cc.

895 {
896  G4double value =0.;
897  G4double y0 = p_p1[0];
898  G4double A = p_p1[1];
899  G4double w = p_p1[2];
900  G4double xc = p_p1[3];
901 
902  value = (y0 *pi*(w*w + 4.*(x-xc)*(x-xc)) + 2.*A*w);
903  return value;
904 }
double A(double temperature)
static const double pi
Definition: G4SIunits.hh:74
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Initialise()

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

Implements G4VEmModel.

Definition at line 90 of file G4LivermorePolarizedGammaConversionModel.cc.

92 {
93  if (verboseLevel > 1)
94  {
95  G4cout << "Calling1 G4LivermorePolarizedGammaConversionModel::Initialise()"
96  << G4endl
97  << "Energy range: "
98  << LowEnergyLimit() / MeV << " MeV - "
99  << HighEnergyLimit() / GeV << " GeV"
100  << G4endl;
101  }
102 
103 
104  if(IsMaster())
105  {
106 
107  // Initialise element selector
108 
109  InitialiseElementSelectors(particle, cuts);
110 
111  // Access to elements
112 
113  char* path = getenv("G4LEDATA");
114 
115  G4ProductionCutsTable* theCoupleTable =
117 
118  G4int numOfCouples = theCoupleTable->GetTableSize();
119 
120  for(G4int i=0; i<numOfCouples; ++i)
121  {
122  const G4Material* material =
123  theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
124  const G4ElementVector* theElementVector = material->GetElementVector();
125  G4int nelm = material->GetNumberOfElements();
126 
127  for (G4int j=0; j<nelm; ++j)
128  {
129  G4int Z = (G4int)(*theElementVector)[j]->GetZ();
130  if(Z < 1) { Z = 1; }
131  else if(Z > maxZ) { Z = maxZ; }
132  if(!data[Z]) { ReadData(Z, path); }
133  }
134  }
135  }
136  if(isInitialised) { return; }
138  isInitialised = true;
139 
140 }
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
static const double MeV
Definition: G4SIunits.hh:211
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
int G4int
Definition: G4Types.hh:78
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
static G4ProductionCutsTable * GetProductionCutsTable()
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
#define G4endl
Definition: G4ios.hh:61
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:134
Here is the call graph for this function:

◆ InitialiseForElement()

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

Reimplemented from G4VEmModel.

Definition at line 1087 of file G4LivermorePolarizedGammaConversionModel.cc.

1090 {
1091  G4AutoLock l(&LivermorePolarizedGammaConversionModelMutex);
1092  // G4cout << "G4LivermorePolarizedGammaConversionModel::InitialiseForElement Z= "
1093  // << Z << G4endl;
1094  if(!data[Z]) { ReadData(Z); }
1095  l.unlock();
1096 }
Float_t Z
Here is the call graph for this function:
Here is the caller graph for this function:

◆ InitialiseLocal()

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

Reimplemented from G4VEmModel.

Definition at line 145 of file G4LivermorePolarizedGammaConversionModel.cc.

147 {
149 }
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:

◆ MinPrimaryEnergy()

G4double G4LivermorePolarizedGammaConversionModel::MinPrimaryEnergy ( const G4Material ,
const G4ParticleDefinition ,
G4double   
)
virtual

Reimplemented from G4VEmModel.

Definition at line 153 of file G4LivermorePolarizedGammaConversionModel.cc.

◆ operator=()

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

◆ Poli()

G4double G4LivermorePolarizedGammaConversionModel::Poli ( G4double  a,
G4double  b,
G4double  c,
G4double  x 
)
private

Definition at line 826 of file G4LivermorePolarizedGammaConversionModel.cc.

827 {
828  G4double value=0.;
829  if(x>0.)
830  {
831  value =(a + b/x + c/(x*x*x));
832  }
833  else
834  {
835  //G4cout << "ERROR in Poli! " << G4endl;
836  }
837  return value;
838 }
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ReadData()

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

Definition at line 161 of file G4LivermorePolarizedGammaConversionModel.cc.

162 {
163  if (verboseLevel > 1)
164  {
165  G4cout << "Calling ReadData() of G4LivermorePolarizedGammaConversionModel"
166  << G4endl;
167  }
168 
169  if(data[Z]) { return; }
170 
171  const char* datadir = path;
172 
173  if(!datadir)
174  {
175  datadir = getenv("G4LEDATA");
176  if(!datadir)
177  {
178  G4Exception("G4LivermorePolarizedGammaConversionModel::ReadData()",
179  "em0006",FatalException,
180  "Environment variable G4LEDATA not defined");
181  return;
182  }
183  }
184 
185  //
186 
187  data[Z] = new G4LPhysicsFreeVector();
188 
189  //
190 
191  std::ostringstream ost;
192  ost << datadir << "/livermore/pair/pp-cs-" << Z <<".dat";
193  std::ifstream fin(ost.str().c_str());
194 
195  if( !fin.is_open())
196  {
198  ed << "G4LivermorePolarizedGammaConversionModel data file <" << ost.str().c_str()
199  << "> is not opened!" << G4endl;
200  G4Exception("G4LivermorePolarizedGammaConversionModel::ReadData()",
201  "em0003",FatalException,
202  ed,"G4LEDATA version should be G4EMLOW6.27 or later.");
203  return;
204  }
205  else
206  {
207 
208  if(verboseLevel > 3) { G4cout << "File " << ost.str()
209  << " is opened by G4LivermorePolarizedGammaConversionModel" << G4endl;}
210 
211  data[Z]->Retrieve(fin, true);
212  }
213 
214  // Activation of spline interpolation
215  data[Z] ->SetSpline(true);
216 
217 }
TString fin
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
void SetSpline(G4bool)
G4GLOB_DLL std::ostream G4cout
Float_t Z
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
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SampleSecondaries()

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

Implements G4VEmModel.

Definition at line 269 of file G4LivermorePolarizedGammaConversionModel.cc.

274 {
275 
276  // Fluorescence generated according to:
277  // J. Stepanek ,"A program to determine the radiation spectra due to a single atomic
278  // subshell ionisation by a particle or due to deexcitation or decay of radionuclides",
279  // Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
280 
281  if (verboseLevel > 3)
282  G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedGammaConversionModel" << G4endl;
283 
284  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
285  // Within energy limit?
286 
287  if(photonEnergy <= lowEnergyLimit)
288  {
289  fParticleChange->ProposeTrackStatus(fStopAndKill);
290  fParticleChange->SetProposedKineticEnergy(0.);
291  return;
292  }
293 
294 
295  G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
296  G4ThreeVector gammaDirection0 = aDynamicGamma->GetMomentumDirection();
297 
298  // Make sure that the polarization vector is perpendicular to the
299  // gamma direction. If not
300 
301  if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
302  { // only for testing now
303  gammaPolarization0 = GetRandomPolarization(gammaDirection0);
304  }
305  else
306  {
307  if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
308  {
309  gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
310  }
311  }
312 
313  // End of Protection
314 
315 
316  G4double epsilon ;
317  G4double epsilon0Local = electron_mass_c2 / photonEnergy ;
318 
319  // Do it fast if photon energy < 2. MeV
320 
321  if (photonEnergy < smallEnergy )
322  {
323  epsilon = epsilon0Local + (0.5 - epsilon0Local) * G4UniformRand();
324  }
325  else
326  {
327 
328  // Select randomly one element in the current material
329 
330  const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
331  const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy);
332 
333 
334  if (element == 0)
335  {
336  G4cout << "G4LivermorePolarizedGammaConversionModel::SampleSecondaries - element = 0" << G4endl;
337  return;
338  }
339 
340 
341  G4IonisParamElm* ionisation = element->GetIonisation();
342 
343  if (ionisation == 0)
344  {
345  G4cout << "G4LivermorePolarizedGammaConversionModel::SampleSecondaries - ionisation = 0" << G4endl;
346  return;
347  }
348 
349  // Extract Coulomb factor for this Element
350 
351  G4double fZ = 8. * (ionisation->GetlogZ3());
352  if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb());
353 
354  // Limits of the screening variable
355  G4double screenFactor = 136. * epsilon0Local / (element->GetIonisation()->GetZ3()) ;
356  G4double screenMax = G4Exp ((42.24 - fZ)/8.368) - 0.952 ;
357  G4double screenMin = std::min(4.*screenFactor,screenMax) ;
358 
359  // Limits of the energy sampling
360  G4double epsilon1 = 0.5 - 0.5 * sqrt(1. - screenMin / screenMax) ;
361  G4double epsilonMin = std::max(epsilon0Local,epsilon1);
362  G4double epsilonRange = 0.5 - epsilonMin ;
363 
364  // Sample the energy rate of the created electron (or positron)
365  G4double screen;
366  G4double gReject ;
367 
368  G4double f10 = ScreenFunction1(screenMin) - fZ;
369  G4double f20 = ScreenFunction2(screenMin) - fZ;
370  G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
371  G4double normF2 = std::max(1.5 * f20,0.);
372 
373  do {
374  if (normF1 / (normF1 + normF2) > G4UniformRand() )
375  {
376  epsilon = 0.5 - epsilonRange * pow(G4UniformRand(), 0.3333) ;
377  screen = screenFactor / (epsilon * (1. - epsilon));
378  gReject = (ScreenFunction1(screen) - fZ) / f10 ;
379  }
380  else
381  {
382  epsilon = epsilonMin + epsilonRange * G4UniformRand();
383  screen = screenFactor / (epsilon * (1 - epsilon));
384  gReject = (ScreenFunction2(screen) - fZ) / f20 ;
385 
386 
387  }
388  } while ( gReject < G4UniformRand() );
389 
390  } // End of epsilon sampling
391 
392  // Fix charges randomly
393 
394  G4double electronTotEnergy;
395  G4double positronTotEnergy;
396 
397 
398  // if (G4int(2*G4UniformRand()))
399  if (G4UniformRand() > 0.5)
400  {
401  electronTotEnergy = (1. - epsilon) * photonEnergy;
402  positronTotEnergy = epsilon * photonEnergy;
403  }
404  else
405  {
406  positronTotEnergy = (1. - epsilon) * photonEnergy;
407  electronTotEnergy = epsilon * photonEnergy;
408  }
409 
410  // Scattered electron (positron) angles. ( Z - axis along the parent photon)
411  // Universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211),
412  // derived from Tsai distribution (Rev. Mod. Phys. 49, 421 (1977)
413 
414 /*
415  G4double u;
416  const G4double a1 = 0.625;
417  G4double a2 = 3. * a1;
418 
419  if (0.25 > G4UniformRand())
420  {
421  u = - log(G4UniformRand() * G4UniformRand()) / a1 ;
422  }
423  else
424  {
425  u = - log(G4UniformRand() * G4UniformRand()) / a2 ;
426  }
427 */
428 
429  G4double Ene = electronTotEnergy/electron_mass_c2; // Normalized energy
430 
431  G4double cosTheta = 0.;
432  G4double sinTheta = 0.;
433 
434  SetTheta(&cosTheta,&sinTheta,Ene);
435 
436  // G4double theta = u * electron_mass_c2 / photonEnergy ;
437  // G4double phi = twopi * G4UniformRand() ;
438 
439  G4double phi,psi=0.;
440 
441  //corrected e+ e- angular angular distribution //preliminary!
442 
443  // if(photonEnergy>50*MeV)
444  // {
445  phi = SetPhi(photonEnergy);
446  psi = SetPsi(photonEnergy,phi);
447  // }
448  //else
449  // {
450  //psi = G4UniformRand()*2.*pi;
451  //phi = pi; // coplanar
452  // }
453 
454  Psi = psi;
455  Phi = phi;
456  //G4cout << "PHI " << phi << G4endl;
457  //G4cout << "PSI " << psi << G4endl;
458 
459  G4double phie, phip;
460  G4double choice, choice2;
461  choice = G4UniformRand();
462  choice2 = G4UniformRand();
463 
464  if (choice2 <= 0.5)
465  {
466  // do nothing
467  // phi = phi;
468  }
469  else
470  {
471  phi = -phi;
472  }
473 
474  if (choice <= 0.5)
475  {
476  phie = psi; //azimuthal angle for the electron
477  phip = phie+phi; //azimuthal angle for the positron
478  }
479  else
480  {
481  // opzione 1 phie / phip equivalenti
482 
483  phip = psi; //azimuthal angle for the positron
484  phie = phip + phi; //azimuthal angle for the electron
485  }
486 
487 
488  // Electron Kinematics
489 
490  G4double dirX = sinTheta*cos(phie);
491  G4double dirY = sinTheta*sin(phie);
492  G4double dirZ = cosTheta;
493  G4ThreeVector electronDirection(dirX,dirY,dirZ);
494 
495  // Kinematics of the created pair:
496  // the electron and positron are assumed to have a symetric angular
497  // distribution with respect to the Z axis along the parent photon
498 
499  //G4double localEnergyDeposit = 0. ;
500 
501  G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
502 
503  SystemOfRefChange(gammaDirection0,electronDirection,
504  gammaPolarization0);
505 
507  electronDirection,
508  electronKineEnergy);
509 
510  // The e+ is always created (even with kinetic energy = 0) for further annihilation
511 
512  Ene = positronTotEnergy/electron_mass_c2; // Normalized energy
513 
514  cosTheta = 0.;
515  sinTheta = 0.;
516 
517  SetTheta(&cosTheta,&sinTheta,Ene);
518 
519  // Positron Kinematics
520 
521  dirX = sinTheta*cos(phip);
522  dirY = sinTheta*sin(phip);
523  dirZ = cosTheta;
524  G4ThreeVector positronDirection(dirX,dirY,dirZ);
525 
526  G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
527  SystemOfRefChange(gammaDirection0,positronDirection,
528  gammaPolarization0);
529 
530  // Create G4DynamicParticle object for the particle2
532  positronDirection, positronKineEnergy);
533 
534 
535  fvect->push_back(particle1);
536  fvect->push_back(particle2);
537 
538  // Kill the incident photon
539  fParticleChange->SetProposedKineticEnergy(0.);
540  fParticleChange->ProposeTrackStatus(fStopAndKill);
541 
542 }
static const double MeV
Definition: G4SIunits.hh:211
G4double GetlogZ3() const
double howOrthogonal(const Hep3Vector &v) const
Definition: SpaceVector.cc:219
G4double GetZ3() const
G4double GetKineticEnergy() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double mag() const
G4IonisParamElm * GetIonisation() const
Definition: G4Element.hh:198
float electron_mass_c2
Definition: hepunit.py:274
G4double GetfCoulomb() const
Definition: G4Element.hh:190
bool isOrthogonal(const Hep3Vector &v, double epsilon=tolerance) const
Definition: SpaceVector.cc:237
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
void SystemOfRefChange(G4ThreeVector &direction0, G4ThreeVector &direction1, G4ThreeVector &polarization0)
static G4Positron * Positron()
Definition: G4Positron.cc:94
const G4ThreeVector & GetMomentumDirection() const
const G4ThreeVector & GetPolarization() const
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4ParticleDefinition * GetDefinition() const
#define G4endl
Definition: G4ios.hh:61
G4ThreeVector GetPerpendicularPolarization(const G4ThreeVector &direction0, const G4ThreeVector &polarization0) const
double G4double
Definition: G4Types.hh:76
G4ThreeVector GetRandomPolarization(G4ThreeVector &direction0)
double epsilon(double density, double temperature)
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:

◆ ScreenFunction1()

G4double G4LivermorePolarizedGammaConversionModel::ScreenFunction1 ( G4double  screenVariable)
private

Definition at line 546 of file G4LivermorePolarizedGammaConversionModel.cc.

547 {
548  // Compute the value of the screening function 3*phi1 - phi2
549 
550  G4double value;
551 
552  if (screenVariable > 1.)
553  value = 42.24 - 8.368 * log(screenVariable + 0.952);
554  else
555  value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
556 
557  return value;
558 }
double G4double
Definition: G4Types.hh:76
Here is the caller graph for this function:

◆ ScreenFunction2()

G4double G4LivermorePolarizedGammaConversionModel::ScreenFunction2 ( G4double  screenVariable)
private

Definition at line 562 of file G4LivermorePolarizedGammaConversionModel.cc.

563 {
564  // Compute the value of the screening function 1.5*phi1 - 0.5*phi2
565 
566  G4double value;
567 
568  if (screenVariable > 1.)
569  value = 42.24 - 8.368 * log(screenVariable + 0.952);
570  else
571  value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
572 
573  return value;
574 }
double G4double
Definition: G4Types.hh:76
Here is the caller graph for this function:

◆ SetPerpendicularVector()

G4ThreeVector G4LivermorePolarizedGammaConversionModel::SetPerpendicularVector ( G4ThreeVector a)
private

Definition at line 1000 of file G4LivermorePolarizedGammaConversionModel.cc.

1001 {
1002  G4double dx = a.x();
1003  G4double dy = a.y();
1004  G4double dz = a.z();
1005  G4double x = dx < 0.0 ? -dx : dx;
1006  G4double y = dy < 0.0 ? -dy : dy;
1007  G4double z = dz < 0.0 ? -dz : dz;
1008  if (x < y) {
1009  return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
1010  }else{
1011  return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
1012  }
1013 }
CLHEP::Hep3Vector G4ThreeVector
Double_t y
double x() const
double y() const
double z() const
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SetPhi()

G4double G4LivermorePolarizedGammaConversionModel::SetPhi ( G4double  Energy)
private

Definition at line 592 of file G4LivermorePolarizedGammaConversionModel.cc.

593 {
594 
595 
596  G4double value = 0.;
597  G4double Ene = Energy/MeV;
598 
599  G4double pl[4];
600 
601 
602  G4double pt[2];
603  G4double xi = 0;
604  G4double xe = 0.;
605  G4double n1=0.;
606  G4double n2=0.;
607 
608 
609  if (Ene>=50.)
610  {
611  const G4double ay0=5.6, by0=18.6, aa0=2.9, ba0 = 8.16E-3;
612  const G4double aw = 0.0151, bw = 10.7, cw = -410.;
613 
614  const G4double axc = 3.1455, bxc = -1.11, cxc = 310.;
615 
616  pl[0] = Fln(ay0,by0,Ene);
617  pl[1] = aa0 + ba0*(Ene);
618  pl[2] = Poli(aw,bw,cw,Ene);
619  pl[3] = Poli(axc,bxc,cxc,Ene);
620 
621  const G4double abf = 3.1216, bbf = 2.68;
622  pt[0] = -1.4;
623  pt[1] = abf + bbf/Ene;
624 
625 
626 
627  //G4cout << "PL > 50. "<< pl[0] << " " << pl[1] << " " << pl[2] << " " <<pl[3] << " " << G4endl;
628 
629  xi = 3.0;
630  xe = Encu(pl,pt,xi);
631  //G4cout << "ENCU "<< xe << G4endl;
632  n1 = Fintlor(pl,pi) - Fintlor(pl,xe);
633  n2 = Finttan(pt,xe) - Finttan(pt,0.);
634  }
635  else
636  {
637  const G4double ay0=0.144, by0=0.11;
638  const G4double aa0=2.7, ba0 = 2.74;
639  const G4double aw = 0.21, bw = 10.8, cw = -58.;
640  const G4double axc = 3.17, bxc = -0.87, cxc = -6.;
641 
642  pl[0] = Fln(ay0, by0, Ene);
643  pl[1] = Fln(aa0, ba0, Ene);
644  pl[2] = Poli(aw,bw,cw,Ene);
645  pl[3] = Poli(axc,bxc,cxc,Ene);
646 
647  //G4cout << "PL < 50."<< pl[0] << " " << pl[1] << " " << pl[2] << " " <<pl[3] << " " << G4endl;
648  //G4cout << "ENCU "<< xe << G4endl;
649  n1 = Fintlor(pl,pi) - Fintlor(pl,xe);
650 
651  }
652 
653 
654  G4double n=0.;
655  n = n1+n2;
656 
657  G4double c1 = 0.;
658  c1 = Glor(pl, xe);
659 
660 /*
661  G4double xm = 0.;
662  xm = Flor(pl,pl[3])*Glor(pl,pl[3]);
663 */
664 
665  G4double r1,r2,r3;
666  G4double xco=0.;
667 
668  if (Ene>=50.)
669  {
670  r1= G4UniformRand();
671  if( r1>=n2/n)
672  {
673  do
674  {
675  r2 = G4UniformRand();
676  value = Finvlor(pl,xe,r2);
677  xco = Glor(pl,value)/c1;
678  r3 = G4UniformRand();
679  } while(r3>=xco);
680  }
681  else
682  {
683  value = Finvtan(pt,n,r1);
684  }
685  }
686  else
687  {
688  do
689  {
690  r2 = G4UniformRand();
691  value = Finvlor(pl,xe,r2);
692  xco = Glor(pl,value)/c1;
693  r3 = G4UniformRand();
694  } while(r3>=xco);
695  }
696 
697  // G4cout << "PHI = " <<value << G4endl;
698  return value;
699 }
G4double Poli(G4double, G4double, G4double, G4double)
static const double MeV
Definition: G4SIunits.hh:211
pl
Definition: readPY.py:5
Char_t n[5]
#define G4UniformRand()
Definition: Randomize.hh:97
TMarker * pt
Definition: egs.C:25
static const double pi
Definition: G4SIunits.hh:74
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SetPsi()

G4double G4LivermorePolarizedGammaConversionModel::SetPsi ( G4double  Energy,
G4double  PhiLocal 
)
private

Definition at line 700 of file G4LivermorePolarizedGammaConversionModel.cc.

701 {
702 
703  G4double value = 0.;
704  G4double Ene = Energy/MeV;
705 
706  G4double p0l[4];
707  G4double ppml[4];
708  G4double p0t[2];
709  G4double ppmt[2];
710 
711  G4double xi = 0.;
712  G4double xe0 = 0.;
713  G4double xepm = 0.;
714 
715  if (Ene>=50.)
716  {
717  const G4double ay00 = 3.4, by00 = 9.8, aa00 = 1.34, ba00 = 5.3;
718  const G4double aw0 = 0.014, bw0 = 9.7, cw0 = -2.E4;
719  const G4double axc0 = 3.1423, bxc0 = -2.35, cxc0 = 0.;
720  const G4double ay0p = 1.53, by0p = 3.2, aap = 0.67, bap = 8.5E-3;
721  const G4double awp = 6.9E-3, bwp = 12.6, cwp = -3.8E4;
722  const G4double axcp = 2.8E-3,bxcp = -3.133;
723  const G4double abf0 = 3.1213, bbf0 = 2.61;
724  const G4double abfpm = 3.1231, bbfpm = 2.84;
725 
726  p0l[0] = Fln(ay00, by00, Ene);
727  p0l[1] = Fln(aa00, ba00, Ene);
728  p0l[2] = Poli(aw0, bw0, cw0, Ene);
729  p0l[3] = Poli(axc0, bxc0, cxc0, Ene);
730 
731  ppml[0] = Fln(ay0p, by0p, Ene);
732  ppml[1] = aap + bap*(Ene);
733  ppml[2] = Poli(awp, bwp, cwp, Ene);
734  ppml[3] = Fln(axcp,bxcp,Ene);
735 
736  p0t[0] = -0.81;
737  p0t[1] = abf0 + bbf0/Ene;
738  ppmt[0] = -0.6;
739  ppmt[1] = abfpm + bbfpm/Ene;
740 
741  //G4cout << "P0L > 50"<< p0l[0] << " " << p0l[1] << " " << p0l[2] << " " <<p0l[3] << " " << G4endl;
742  //G4cout << "PPML > 50"<< ppml[0] << " " << ppml[1] << " " << ppml[2] << " " <<ppml[3] << " " << G4endl;
743 
744  xi = 3.0;
745  xe0 = Encu(p0l, p0t, xi);
746  //G4cout << "ENCU1 "<< xe0 << G4endl;
747  xepm = Encu(ppml, ppmt, xi);
748  //G4cout << "ENCU2 "<< xepm << G4endl;
749  }
750  else
751  {
752  const G4double ay00 = 2.82, by00 = 6.35;
753  const G4double aa00 = -1.75, ba00 = 0.25;
754 
755  const G4double aw0 = 0.028, bw0 = 5., cw0 = -50.;
756  const G4double axc0 = 3.14213, bxc0 = -2.3, cxc0 = 5.7;
757  const G4double ay0p = 1.56, by0p = 3.6;
758  const G4double aap = 0.86, bap = 8.3E-3;
759  const G4double awp = 0.022, bwp = 7.4, cwp = -51.;
760  const G4double xcp = 3.1486;
761 
762  p0l[0] = Fln(ay00, by00, Ene);
763  p0l[1] = aa00+pow(Ene, ba00);
764  p0l[2] = Poli(aw0, bw0, cw0, Ene);
765  p0l[3] = Poli(axc0, bxc0, cxc0, Ene);
766  ppml[0] = Fln(ay0p, by0p, Ene);
767  ppml[1] = aap + bap*(Ene);
768  ppml[2] = Poli(awp, bwp, cwp, Ene);
769  ppml[3] = xcp;
770 
771  }
772 
773  G4double a,b=0.;
774 
775  if (Ene>=50.)
776  {
777  if (PhiLocal>xepm)
778  {
779  b = (ppml[0]+2*ppml[1]*ppml[2]*Flor(ppml,PhiLocal));
780  }
781  else
782  {
783  b = Ftan(ppmt,PhiLocal);
784  }
785  if (PhiLocal>xe0)
786  {
787  a = (p0l[0]+2*p0l[1]*p0l[2]*Flor(p0l,PhiLocal));
788  }
789  else
790  {
791  a = Ftan(p0t,PhiLocal);
792  }
793  }
794  else
795  {
796  b = (ppml[0]+2*ppml[1]*ppml[2]*Flor(ppml,PhiLocal));
797  a = (p0l[0]+2*p0l[1]*p0l[2]*Flor(p0l,PhiLocal));
798  }
799  G4double nr =0.;
800 
801  if (b>a)
802  {
803  nr = 1./b;
804  }
805  else
806  {
807  nr = 1./a;
808  }
809 
810  G4double r1,r2=0.;
811  G4double r3 =-1.;
812  do
813  {
814  r1 = G4UniformRand();
815  r2 = G4UniformRand();
816  //value = r2*pi;
817  value = 2.*r2*pi;
818  r3 = nr*(a*cos(value)*cos(value) + b*sin(value)*sin(value));
819  }while(r1>r3);
820 
821  return value;
822 }
G4double Poli(G4double, G4double, G4double, G4double)
static const double MeV
Definition: G4SIunits.hh:211
#define G4UniformRand()
Definition: Randomize.hh:97
static const double pi
Definition: G4SIunits.hh:74
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SetTheta()

void G4LivermorePolarizedGammaConversionModel::SetTheta ( G4double p_cosTheta,
G4double p_sinTheta,
G4double  Energy 
)
private

Definition at line 577 of file G4LivermorePolarizedGammaConversionModel.cc.

578 {
579 
580  // to avoid computational errors since Theta could be very small
581  // Energy in Normalized Units (!)
582 
583  G4double Momentum = sqrt(Energy*Energy -1);
584  G4double Rand = G4UniformRand();
585 
586  *p_cosTheta = (Energy*((2*Rand)- 1) + Momentum)/((Momentum*(2*Rand-1))+Energy);
587  *p_sinTheta = (2*sqrt(Rand*(1-Rand)))/(Momentum*(2*Rand-1)+Energy);
588 }
#define G4UniformRand()
Definition: Randomize.hh:97
double G4double
Definition: G4Types.hh:76
Here is the caller graph for this function:

◆ SystemOfRefChange()

void G4LivermorePolarizedGammaConversionModel::SystemOfRefChange ( G4ThreeVector direction0,
G4ThreeVector direction1,
G4ThreeVector polarization0 
)
private

Definition at line 1063 of file G4LivermorePolarizedGammaConversionModel.cc.

1065 {
1066  // direction0 is the original photon direction ---> z
1067  // polarization0 is the original photon polarization ---> x
1068  // need to specify y axis in the real reference frame ---> y
1069  G4ThreeVector Axis_Z0 = direction0.unit();
1070  G4ThreeVector Axis_X0 = polarization0.unit();
1071  G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
1072 
1073  G4double direction_x = direction1.getX();
1074  G4double direction_y = direction1.getY();
1075  G4double direction_z = direction1.getZ();
1076 
1077  direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
1078 
1079 }
double getY() const
Hep3Vector cross(const Hep3Vector &) const
double getX() const
Hep3Vector unit() const
double getZ() const
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

Member Data Documentation

◆ data

G4LPhysicsFreeVector * G4LivermorePolarizedGammaConversionModel::data = {0}
staticprivate

Definition at line 138 of file G4LivermorePolarizedGammaConversionModel.hh.

◆ fParticleChange

G4ParticleChangeForGamma* G4LivermorePolarizedGammaConversionModel::fParticleChange
private

Definition at line 80 of file G4LivermorePolarizedGammaConversionModel.hh.

◆ isInitialised

G4bool G4LivermorePolarizedGammaConversionModel::isInitialised
private

Definition at line 85 of file G4LivermorePolarizedGammaConversionModel.hh.

◆ lowEnergyLimit

G4double G4LivermorePolarizedGammaConversionModel::lowEnergyLimit
private

Definition at line 82 of file G4LivermorePolarizedGammaConversionModel.hh.

◆ maxZ

G4int G4LivermorePolarizedGammaConversionModel::maxZ = 99
staticprivate

Definition at line 137 of file G4LivermorePolarizedGammaConversionModel.hh.

◆ Phi

G4double G4LivermorePolarizedGammaConversionModel::Phi
private

Definition at line 131 of file G4LivermorePolarizedGammaConversionModel.hh.

◆ Psi

G4double G4LivermorePolarizedGammaConversionModel::Psi
private

Definition at line 131 of file G4LivermorePolarizedGammaConversionModel.hh.

◆ smallEnergy

G4double G4LivermorePolarizedGammaConversionModel::smallEnergy
private

Definition at line 130 of file G4LivermorePolarizedGammaConversionModel.hh.

◆ verboseLevel

G4int G4LivermorePolarizedGammaConversionModel::verboseLevel
private

Definition at line 86 of file G4LivermorePolarizedGammaConversionModel.hh.


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