66 :
G4VEmModel(nam),fParticleChange(0),isInitialised(false),
67 scatterFunctionData(0),
68 crossSectionHandler(0),fAtomDeexcitation(0)
79 G4cout <<
"Livermore Modified Compton model is constructed " <<
G4endl;
100 G4cout <<
"Calling G4LivermoreComptonModifiedModel::Initialise()" <<
G4endl;
112 G4String crossSectionFile =
"comp/ce-cs-";
116 G4String scatterFile =
"comp/ce-sf-";
128 G4cout <<
"Loaded cross section files for Livermore Modified Compton model" <<
G4endl;
139 G4cout <<
"Livermore modified Compton model is initialized " <<
G4endl
156 G4cout <<
"Calling ComputeCrossSectionPerAtom() of G4LivermoreComptonModifiedModel" <<
G4endl;
190 G4cout <<
"G4LivermoreComptonModifiedModel::SampleSecondaries() E(MeV)= "
208 G4double epsilon0Local = 1. / (1. + 2. * e0m);
209 G4double epsilon0Sq = epsilon0Local * epsilon0Local;
210 G4double alpha1 = -std::log(epsilon0Local);
211 G4double alpha2 = 0.5 * (1. - epsilon0Sq);
232 epsilonSq = epsilon0Sq + (1. - epsilon0Sq) *
G4UniformRand();
233 epsilon = std::sqrt(epsilonSq);
236 oneCosT = (1. -
epsilon) / ( epsilon * e0m);
237 sinT2 = oneCosT * (2. - oneCosT);
238 G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/
cm);
240 gReject = (1. - epsilon * sinT2 / (1. + epsilonSq)) * scatteringFunction;
245 G4double sinTheta = std::sqrt (sinT2);
247 G4double dirx = sinTheta * std::cos(phi);
248 G4double diry = sinTheta * std::sin(phi);
257 G4int maxDopplerIterations = 1000;
259 G4double photonEoriginal = epsilon * photonEnergy0;
266 G4double momentum_au_to_nat = 1.992851740*std::pow(10.,-24.);
267 G4double e_mass_kg = 9.10938188 * std::pow(10.,-31.);
292 }
while(Alpha >= (
pi/2.0));
294 ePAU = pSample / std::cos(Alpha);
298 G4double ePSI = ePAU * momentum_au_to_nat;
299 G4double u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / ((e_mass_kg*e_mass_kg)*(vel_c*vel_c)+(ePSI*ePSI)))/vel_c;
303 systemE = eEIncident+photonEnergy0;
307 G4double pDoppler2 = pDoppler * pDoppler;
309 G4double var3 = var2*var2 - pDoppler2;
310 G4double var4 = var2 - pDoppler2 * cosTheta;
311 G4double var = var4*var4 - var3 + pDoppler2 * var3;
317 if (
G4UniformRand() < 0.5) { photonE = (var4 - varSqrt) * scale; }
318 else { photonE = (var4 + varSqrt) * scale; }
324 }
while ( iteration <= maxDopplerIterations &&
325 (photonE < 0. || photonE > eMax ) );
336 if(eKineticEnergy < 0.0) {
337 G4cout <<
"Error, kinetic energy of electron less than zero" <<
G4endl;
344 G4double E_num = photonEnergy0 - photonE*cosTheta;
345 G4double E_dom = sqrt(photonEnergy0*photonEnergy0 + photonE*photonE -2*photonEnergy0*photonE*cosTheta);
347 G4double sinThetaE = -sqrt((1. - cosThetaE) * (1. + cosThetaE));
349 eDirX = sinThetaE * std::cos(phi);
350 eDirY = sinThetaE * std::sin(phi);
354 eDirection.
rotateUz(photonDirection0);
356 eDirection,eKineticEnergy) ;
357 fvect->push_back(dp);
363 if (iteration >= maxDopplerIterations)
365 photonE = photonEoriginal;
372 photonDirection1.
rotateUz(photonDirection0);
377 if (photonEnergy1 > 0.)
381 if (iteration < maxDopplerIterations)
384 eDirection.
rotateUz(photonDirection0);
386 eDirection,eKineticEnergy) ;
387 fvect->push_back(dp);
402 size_t nbefore = fvect->size();
406 size_t nafter = fvect->size();
407 if(nafter > nbefore) {
408 for (
size_t i=nbefore; i<nafter; ++i) {
409 bindingE -= ((*fvect)[i])->GetKineticEnergy();
414 if(bindingE < 0.0) { bindingE = 0.0; }
G4double FindValue(G4int Z, G4double e) const
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
void LoadData(const G4String &dataFile)
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4VAtomDeexcitation * AtomDeexcitation()
static constexpr double MeV
const G4ThreeVector & GetMomentumDirection() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Hep3Vector & rotateUz(const Hep3Vector &)
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4VAtomDeexcitation * fAtomDeexcitation
G4LivermoreComptonModifiedModel(const G4ParticleDefinition *p=0, const G4String &nam="LivermoreModifiedCompton")
static constexpr double h_Planck
G4double LowEnergyLimit() const
const G4String & GetName() const
G4double BindingEnergy(G4int Z, G4int shellIndex) const
void SetDeexcitationFlag(G4bool val)
G4ParticleDefinition * GetDefinition() const
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
static constexpr double electron_mass_c2
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
static constexpr double twopi
G4int SelectRandomShell(G4int Z) const
void LoadData(const G4String &fileName)
static constexpr double eV
G4VCrossSectionHandler * crossSectionHandler
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
static G4Electron * Electron()
double epsilon(double density, double temperature)
virtual G4bool LoadData(const G4String &fileName)=0
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
static constexpr double c_light
virtual ~G4LivermoreComptonModifiedModel()
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
G4DopplerProfile profileData
G4double RandomSelectMomentum(G4int Z, G4int shellIndex) const
G4double GetKineticEnergy() const
static constexpr double cm
G4GLOB_DLL std::ostream G4cout
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
static G4LossTableManager * Instance()
const G4Material * GetMaterial() const
static constexpr double pi
G4ParticleChangeForGamma * GetParticleChangeForGamma()
G4VEMDataSet * scatterFunctionData
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static constexpr double fine_structure_const
G4double HighEnergyLimit() const
void ProposeTrackStatus(G4TrackStatus status)
static constexpr double GeV
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
G4ParticleChangeForGamma * fParticleChange