53 :
G4VEmModel(nam),fParticleChange(0),isInitialised(false),
54 meanFreePathTable(0),scatterFunctionData(0),crossSectionHandler(0)
70 G4cout <<
"Livermore Polarized Compton is constructed " <<
G4endl
93 G4cout <<
"Calling G4LivermorePolarizedComptonModel::Initialise()" <<
G4endl;
105 G4String crossSectionFile =
"comp/ce-cs-";
112 G4String scatterFile =
"comp/ce-sf-";
118 G4String file =
"/doppler/shell-doppler";
122 G4cout <<
"Loaded cross section files for Livermore Polarized Compton model" <<
G4endl;
127 G4cout <<
"Livermore Polarized Compton model is initialized " << G4endl
150 G4cout <<
"Calling ComputeCrossSectionPerAtom() of G4LivermorePolarizedComptonModel" <<
G4endl;
152 if (GammaEnergy < lowEnergyLimit || GammaEnergy >
highEnergyLimit)
return 0.0;
173 G4cout <<
"Calling SampleSecondaries() of G4LivermorePolarizedComptonModel" <<
G4endl;
187 if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
193 if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
211 G4double E0_m = gammaEnergy0 / electron_mass_c2 ;
221 G4double epsilon, epsilonSq, onecost, sinThetaSqr, greject ;
223 G4double epsilon0Local = 1./(1. + 2*E0_m);
224 G4double epsilon0Sq = epsilon0Local*epsilon0Local;
225 G4double alpha1 = - std::log(epsilon0Local);
226 G4double alpha2 = 0.5*(1.- epsilon0Sq);
228 G4double wlGamma = h_Planck*c_light/gammaEnergy0;
236 epsilonSq = epsilon*epsilon;
241 epsilon = std::sqrt(epsilonSq);
244 onecost = (1.- epsilon)/(epsilon*E0_m);
245 sinThetaSqr = onecost*(2.-onecost);
248 if (sinThetaSqr > 1.)
251 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
252 <<
"sin(theta)**2 = "
258 if (sinThetaSqr < 0.)
261 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
262 <<
"sin(theta)**2 = "
270 G4double x = std::sqrt(onecost/2.) / (wlGamma/
cm);;
272 greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction;
294 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
304 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
314 G4double sinTheta = std::sqrt (sinThetaSqr);
320 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
330 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
340 G4double dirx = sinTheta*std::cos(phi);
341 G4double diry = sinTheta*std::sin(phi);
354 G4int maxDopplerIterations = 1000;
356 G4double photonEoriginal = epsilon * gammaEnergy0;
368 eMax = gammaEnergy0 - bindingE;
373 G4double pDoppler = pSample * fine_structure_const;
374 G4double pDoppler2 = pDoppler * pDoppler;
375 G4double var2 = 1. + onecost * E0_m;
376 G4double var3 = var2*var2 - pDoppler2;
377 G4double var4 = var2 - pDoppler2 * cosTheta;
378 G4double var = var4*var4 - var3 + pDoppler2 * var3;
382 G4double scale = gammaEnergy0 / var3;
384 if (
G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;
385 else photonE = (var4 + varSqrt) * scale;
391 }
while ( iteration <= maxDopplerIterations &&
392 (photonE < 0. || photonE > eMax || photonE < eMax*
G4UniformRand()) );
396 if (iteration >= maxDopplerIterations)
398 photonE = photonEoriginal;
402 gammaEnergy1 = photonE;
420 gammaDirection1 = tmpDirection1;
425 gammaPolarization0,gammaPolarization1);
427 if (gammaEnergy1 > 0.)
444 G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;
448 if(ElecKineEnergy < 0.0) {
455 G4double ElecMomentum = std::sqrt(ElecKineEnergy*(ElecKineEnergy+2.*electron_mass_c2));
457 G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 -
458 gammaEnergy1 * gammaDirection1) * (1./ElecMomentum));
463 fvect->push_back(dp);
486 b = energyRate + 1/energyRate;
488 phiProbability = 1 - (a/b)*(std::cos(phi)*std::cos(phi));
493 while ( rand2 > phiProbability );
530 c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
531 c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
532 c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
556 return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
570 G4double sinTheta = std::sqrt(sinSqrTh);
574 G4double normalisation = std::sqrt(1. - cosSqrPhi*sinSqrTh);
613 if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsilon+1.0/epsilon)-4.0*sinSqrTh*cosSqrPhi))
628 G4double sinBeta = std::sqrt(1-cosBeta*cosBeta);
632 G4double xParallel = normalisation*cosBeta;
633 G4double yParallel = -(sinSqrTh*cosPhi*sinPhi)*cosBeta/normalisation;
634 G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation;
636 G4double yPerpendicular = (costheta)*sinBeta/normalisation;
637 G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation;
639 G4double xTotal = (xParallel + xPerpendicular);
640 G4double yTotal = (yParallel + yPerpendicular);
641 G4double zTotal = (zParallel + zPerpendicular);
643 gammaPolarization1.setX(xTotal);
644 gammaPolarization1.setY(yTotal);
645 gammaPolarization1.setZ(zTotal);
647 return gammaPolarization1;
665 G4double direction_x = direction1.getX();
666 G4double direction_y = direction1.getY();
667 G4double direction_z = direction1.getZ();
669 direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
670 G4double polarization_x = polarization1.getX();
671 G4double polarization_y = polarization1.getY();
672 G4double polarization_z = polarization1.getZ();
674 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)
void SetHighEnergyLimit(G4double)
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)