65 :
G4VEmModel(nam),fParticleChange(0),isInitialised(false),
66 scatterFunctionData(0),
67 crossSectionHandler(0),fAtomDeexcitation(0)
81 G4cout <<
"Livermore Modified Compton model is constructed " <<
G4endl
107 G4cout <<
"Calling G4LivermoreComptonModifiedModel::Initialise()" <<
G4endl;
119 G4String crossSectionFile =
"comp/ce-cs-";
123 G4String scatterFile =
"comp/ce-sf-";
129 G4String file =
"/doppler/shell-doppler";
135 G4cout <<
"Loaded cross section files for Livermore Modified Compton model" <<
G4endl;
146 G4cout <<
"Livermore Compton model is initialized " <<
G4endl
163 G4cout <<
"Calling ComputeCrossSectionPerAtom() of G4LivermoreComptonModifiedModel" <<
G4endl;
165 if (GammaEnergy < lowEnergyLimit || GammaEnergy >
highEnergyLimit) {
return 0.0; }
196 G4cout <<
"G4LivermoreComptonModifiedModel::SampleSecondaries() E(MeV)= "
210 G4double e0m = photonEnergy0 / electron_mass_c2 ;
218 G4double epsilon0Local = 1. / (1. + 2. * e0m);
219 G4double epsilon0Sq = epsilon0Local * epsilon0Local;
220 G4double alpha1 = -std::log(epsilon0Local);
221 G4double alpha2 = 0.5 * (1. - epsilon0Sq);
223 G4double wlPhoton = h_Planck*c_light/photonEnergy0;
238 epsilonSq = epsilon * epsilon;
242 epsilonSq = epsilon0Sq + (1. - epsilon0Sq) *
G4UniformRand();
243 epsilon = std::sqrt(epsilonSq);
246 oneCosT = (1. - epsilon) / ( epsilon * e0m);
247 sinT2 = oneCosT * (2. - oneCosT);
248 G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/
cm);
250 gReject = (1. - epsilon * sinT2 / (1. + epsilonSq)) * scatteringFunction;
255 G4double sinTheta = std::sqrt (sinT2);
257 G4double dirx = sinTheta * std::cos(phi);
258 G4double diry = sinTheta * std::sin(phi);
267 G4int maxDopplerIterations = 1000;
269 G4double photonEoriginal = epsilon * photonEnergy0;
276 G4double momentum_au_to_nat = 1.992851740*std::pow(10.,-24.);
277 G4double e_mass_kg = 9.10938188 * std::pow(10.,-31.);
302 }
while(Alpha >= (
pi/2.0));
304 ePAU = pSample / std::cos(Alpha);
308 G4double ePSI = ePAU * momentum_au_to_nat;
309 G4double u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / ((e_mass_kg*e_mass_kg)*(vel_c*vel_c)+(ePSI*ePSI)))/vel_c;
310 G4double eEIncident = electron_mass_c2 / sqrt( 1 - (u_temp*u_temp));
313 systemE = eEIncident+photonEnergy0;
315 eMax = systemE - bindingE - electron_mass_c2;
316 G4double pDoppler = pSample * fine_structure_const;
317 G4double pDoppler2 = pDoppler * pDoppler;
319 G4double var3 = var2*var2 - pDoppler2;
320 G4double var4 = var2 - pDoppler2 * cosTheta;
321 G4double var = var4*var4 - var3 + pDoppler2 * var3;
325 G4double scale = photonEnergy0 / var3;
327 if (
G4UniformRand() < 0.5) { photonE = (var4 - varSqrt) * scale; }
328 else { photonE = (var4 + varSqrt) * scale; }
334 }
while ( iteration <= maxDopplerIterations &&
335 (photonE < 0. || photonE > eMax ) );
339 G4double eKineticEnergy = systemE - photonE - bindingE - electron_mass_c2;
346 if(eKineticEnergy < 0.0) {
347 G4cout <<
"Error, kinetic energy of electron less than zero" <<
G4endl;
354 G4double E_num = photonEnergy0 - photonE*cosTheta;
355 G4double E_dom = sqrt(photonEnergy0*photonEnergy0 + photonE*photonE -2*photonEnergy0*photonE*cosTheta);
357 G4double sinThetaE = -sqrt((1. - cosThetaE) * (1. + cosThetaE));
359 eDirX = sinThetaE * std::cos(phi);
360 eDirY = sinThetaE * std::sin(phi);
364 eDirection.rotateUz(photonDirection0);
366 eDirection,eKineticEnergy) ;
367 fvect->push_back(dp);
373 if (iteration >= maxDopplerIterations)
375 photonE = photonEoriginal;
382 photonDirection1.rotateUz(photonDirection0);
387 if (photonEnergy1 > 0.)
391 if (iteration < maxDopplerIterations)
394 eDirection.rotateUz(photonDirection0);
396 eDirection,eKineticEnergy) ;
397 fvect->push_back(dp);
412 size_t nbefore = fvect->size();
416 size_t nafter = fvect->size();
417 if(nafter > nbefore) {
418 for (
size_t i=nbefore; i<nafter; ++i) {
419 bindingE -= ((*fvect)[i])->GetKineticEnergy();
424 if(bindingE < 0.0) { bindingE = 0.0; }
G4double LowEnergyLimit() const
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
static G4LossTableManager * Instance()
G4VCrossSectionHandler * crossSectionHandler
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
G4double HighEnergyLimit() const
const G4String & GetName() const
G4LivermoreComptonModifiedModel(const G4ParticleDefinition *p=0, const G4String &nam="LivermoreModifiedCompton")
G4ParticleChangeForGamma * fParticleChange
G4ParticleDefinition * GetDefinition() const
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4int SelectRandomShell(G4int Z) const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void LoadData(const G4String &fileName)
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4double FindValue(G4int Z, G4double e) const
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetMomentumDirection() const
virtual ~G4LivermoreComptonModifiedModel()
G4double BindingEnergy(G4int Z, G4int shellIndex) const
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4DopplerProfile profileData
G4VAtomDeexcitation * fAtomDeexcitation
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
void LoadData(const G4String &dataFile)
virtual G4bool LoadData(const G4String &fileName)=0
static G4Electron * Electron()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
G4VAtomDeexcitation * AtomDeexcitation()
void ProposeTrackStatus(G4TrackStatus status)
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void SetDeexcitationFlag(G4bool val)
G4ThreeVector G4ParticleMomentum
G4double RandomSelectMomentum(G4int Z, G4int shellIndex) const
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
const G4Material * GetMaterial() const
G4VEMDataSet * scatterFunctionData