53 :
G4VEmModel(nam),fParticleChange(0),isInitialised(false),
54 meanFreePathTable(0),scatterFunctionData(0),crossSectionHandler(0)
65 G4cout <<
"Livermore Polarized Compton is constructed " <<
G4endl;
84 G4cout <<
"Calling G4LivermorePolarizedComptonModel::Initialise()" <<
G4endl;
96 G4String crossSectionFile =
"comp/ce-cs-";
103 G4String scatterFile =
"comp/ce-sf-";
109 G4String file =
"/doppler/shell-doppler";
113 G4cout <<
"Loaded cross section files for Livermore Polarized Compton model" <<
G4endl;
118 G4cout <<
"Livermore Polarized Compton model is initialized " << G4endl
141 G4cout <<
"Calling ComputeCrossSectionPerAtom() of G4LivermorePolarizedComptonModel" <<
G4endl;
165 G4cout <<
"Calling SampleSecondaries() of G4LivermorePolarizedComptonModel" <<
G4endl;
186 if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
192 if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
200 G4double E0_m = gammaEnergy0 / electron_mass_c2 ;
210 G4double epsilon, epsilonSq, onecost, sinThetaSqr, greject ;
212 G4double epsilon0Local = 1./(1. + 2*E0_m);
213 G4double epsilon0Sq = epsilon0Local*epsilon0Local;
214 G4double alpha1 = - std::log(epsilon0Local);
215 G4double alpha2 = 0.5*(1.- epsilon0Sq);
217 G4double wlGamma = h_Planck*c_light/gammaEnergy0;
225 epsilonSq = epsilon*epsilon;
230 epsilon = std::sqrt(epsilonSq);
233 onecost = (1.- epsilon)/(epsilon*E0_m);
234 sinThetaSqr = onecost*(2.-onecost);
237 if (sinThetaSqr > 1.)
240 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
241 <<
"sin(theta)**2 = "
247 if (sinThetaSqr < 0.)
250 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
251 <<
"sin(theta)**2 = "
259 G4double x = std::sqrt(onecost/2.) / (wlGamma/
cm);;
261 greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction;
283 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
293 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
303 G4double sinTheta = std::sqrt (sinThetaSqr);
309 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
319 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
329 G4double dirx = sinTheta*std::cos(phi);
330 G4double diry = sinTheta*std::sin(phi);
343 G4int maxDopplerIterations = 1000;
345 G4double photonEoriginal = epsilon * gammaEnergy0;
357 eMax = gammaEnergy0 - bindingE;
362 G4double pDoppler = pSample * fine_structure_const;
363 G4double pDoppler2 = pDoppler * pDoppler;
364 G4double var2 = 1. + onecost * E0_m;
365 G4double var3 = var2*var2 - pDoppler2;
366 G4double var4 = var2 - pDoppler2 * cosTheta;
367 G4double var = var4*var4 - var3 + pDoppler2 * var3;
371 G4double scale = gammaEnergy0 / var3;
373 if (
G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;
374 else photonE = (var4 + varSqrt) * scale;
380 }
while ( iteration <= maxDopplerIterations &&
381 (photonE < 0. || photonE > eMax || photonE < eMax*
G4UniformRand()) );
385 if (iteration >= maxDopplerIterations)
387 photonE = photonEoriginal;
391 gammaEnergy1 = photonE;
409 gammaDirection1 = tmpDirection1;
414 gammaPolarization0,gammaPolarization1);
416 if (gammaEnergy1 > 0.)
433 G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;
437 if(ElecKineEnergy < 0.0) {
444 G4double ElecMomentum = std::sqrt(ElecKineEnergy*(ElecKineEnergy+2.*electron_mass_c2));
446 G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 -
447 gammaEnergy1 * gammaDirection1) * (1./ElecMomentum));
452 fvect->push_back(dp);
475 b = energyRate + 1/energyRate;
477 phiProbability = 1 - (a/b)*(std::cos(phi)*std::cos(phi));
482 while ( rand2 > phiProbability );
519 c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
520 c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
521 c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
545 return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
559 G4double sinTheta = std::sqrt(sinSqrTh);
563 G4double normalisation = std::sqrt(1. - cosSqrPhi*sinSqrTh);
602 if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsilon+1.0/epsilon)-4.0*sinSqrTh*cosSqrPhi))
617 G4double sinBeta = std::sqrt(1-cosBeta*cosBeta);
621 G4double xParallel = normalisation*cosBeta;
622 G4double yParallel = -(sinSqrTh*cosPhi*sinPhi)*cosBeta/normalisation;
623 G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation;
625 G4double yPerpendicular = (costheta)*sinBeta/normalisation;
626 G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation;
628 G4double xTotal = (xParallel + xPerpendicular);
629 G4double yTotal = (yParallel + yPerpendicular);
630 G4double zTotal = (zParallel + zPerpendicular);
632 gammaPolarization1.setX(xTotal);
633 gammaPolarization1.setY(yTotal);
634 gammaPolarization1.setZ(zTotal);
636 return gammaPolarization1;
654 G4double direction_x = direction1.getX();
655 G4double direction_y = direction1.getY();
656 G4double direction_z = direction1.getZ();
658 direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
659 G4double polarization_x = polarization1.getX();
660 G4double polarization_y = polarization1.getY();
661 G4double polarization_z = polarization1.getZ();
663 polarization1 = (polarization_x*Axis_X0 + polarization_y*Axis_Y0 + polarization_z*Axis_Z0).unit();
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4double LowEnergyLimit() const
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
G4DopplerProfile profileData
virtual ~G4LivermorePolarizedComptonModel()
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
G4double HighEnergyLimit() const
G4ParticleDefinition * GetDefinition() const
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4VCrossSectionHandler * crossSectionHandler
G4int SelectRandomShell(G4int Z) const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void LoadData(const G4String &fileName)
G4double FindValue(G4int Z, G4double e) const
G4GLOB_DLL std::ostream G4cout
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
const G4ThreeVector & GetMomentumDirection() const
G4double BindingEnergy(G4int Z, G4int shellIndex) const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
void ProposePolarization(const G4ThreeVector &dir)
G4VEMDataSet * BuildMeanFreePathForMaterials(const G4DataVector *energyCuts=0)
G4ThreeVector SetNewPolarization(G4double epsilon, G4double sinSqrTheta, G4double phi, G4double cosTheta)
G4VEMDataSet * scatterFunctionData
G4ThreeVector GetPerpendicularPolarization(const G4ThreeVector &direction0, const G4ThreeVector &polarization0) const
G4ParticleChangeForGamma * fParticleChange
G4LivermorePolarizedComptonModel(const G4ParticleDefinition *p=0, const G4String &nam="LivermorePolarizedCompton")
const G4ThreeVector & GetPolarization() const
void LoadData(const G4String &dataFile)
void SystemOfRefChange(G4ThreeVector &direction0, G4ThreeVector &direction1, G4ThreeVector &polarization0, G4ThreeVector &polarization1)
virtual G4bool LoadData(const G4String &fileName)=0
G4double SetPhi(G4double, G4double)
G4ThreeVector SetPerpendicularVector(G4ThreeVector &a)
static G4Electron * Electron()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeTrackStatus(G4TrackStatus status)
G4VEMDataSet * meanFreePathTable
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()
G4ThreeVector GetRandomPolarization(G4ThreeVector &direction0)