81 G4cout <<
"Livermore Polarized Compton is constructed " <<
G4endl;
116 G4cout <<
"Calling G4LivermorePolarizedComptonModel::Initialise()" <<
G4endl;
124 char* path = std::getenv(
"G4LEDATA");
131 for(
G4int i=0; i<numOfCouples; ++i) {
137 for (
G4int j=0; j<nelm; ++j) {
161 G4String scatterFile =
"comp/ce-sf-";
174 G4cout <<
"G4LivermoreComptonModel is initialized " << G4endl
201 G4cout <<
"G4LivermorePolarizedComptonModel::ReadData()"
204 if(
data[Z]) {
return; }
205 const char* datadir = path;
208 datadir = std::getenv(
"G4LEDATA");
211 G4Exception(
"G4LivermorePolarizedComptonModel::ReadData()",
213 "Environment variable G4LEDATA not defined");
221 data[
Z]->SetSpline(
false);
223 std::ostringstream ost;
224 ost << datadir <<
"/livermore/comp/ce-cs-" << Z <<
".dat";
225 std::ifstream
fin(ost.str().c_str());
230 ed <<
"G4LivermorePolarizedComptonModel data file <" << ost.str().c_str()
231 <<
"> is not opened!" <<
G4endl;
234 ed,
"G4LEDATA version should be G4EMLOW6.34 or later");
238 G4cout <<
"File " << ost.str()
239 <<
" is opened by G4LivermorePolarizedComptonModel" <<
G4endl;
259 G4cout <<
"Calling ComputeCrossSectionPerAtom() of G4LivermorePolarizedComptonModel" <<
G4endl;
267 if(intZ < 1 || intZ >
maxZ) {
return cs; }
277 if(!pv) {
return cs; }
284 if(GammaEnergy <= e1) { cs = GammaEnergy/(e1*
e1)*pv->
Value(e1); }
285 else if(GammaEnergy <= e2) { cs = pv->
Value(GammaEnergy)/GammaEnergy; }
286 else if(GammaEnergy > e2) { cs = pv->
Value(e2)/GammaEnergy; }
307 G4cout <<
"Calling SampleSecondaries() of G4LivermorePolarizedComptonModel" <<
G4endl;
328 if(!(gammaPolarization0.
isOrthogonal(gammaDirection0, 1
e-6))||(gammaPolarization0.
mag()==0))
354 G4double epsilon0Local = 1./(1. + 2*E0_m);
355 G4double epsilon0Sq = epsilon0Local*epsilon0Local;
356 G4double alpha1 = - std::log(epsilon0Local);
357 G4double alpha2 = 0.5*(1.- epsilon0Sq);
372 epsilon = std::sqrt(epsilonSq);
375 onecost = (1.-
epsilon)/(epsilon*E0_m);
376 sinThetaSqr = onecost*(2.-onecost);
379 if (sinThetaSqr > 1.)
382 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
383 <<
"sin(theta)**2 = "
389 if (sinThetaSqr < 0.)
392 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
393 <<
"sin(theta)**2 = "
401 G4double x = std::sqrt(onecost/2.) / (wlGamma/
cm);;
403 greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction;
425 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
435 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
445 G4double sinTheta = std::sqrt (sinThetaSqr);
451 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
461 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
471 G4double dirx = sinTheta*std::cos(phi);
472 G4double diry = sinTheta*std::sin(phi);
485 static G4int maxDopplerIterations = 1000;
487 G4double photonEoriginal = epsilon * gammaEnergy0;
506 G4cout <<
"Shell ID= " << shellIdx
507 <<
" Ebind(keV)= " << bindingE/
keV <<
G4endl;
509 eMax = gammaEnergy0 - bindingE;
519 G4double pDoppler2 = pDoppler * pDoppler;
520 G4double var2 = 1. + onecost * E0_m;
521 G4double var3 = var2*var2 - pDoppler2;
522 G4double var4 = var2 - pDoppler2 * cosTheta;
523 G4double var = var4*var4 - var3 + pDoppler2 * var3;
529 if (
G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;
530 else photonE = (var4 + varSqrt) * scale;
536 }
while ( iteration <= maxDopplerIterations &&
537 (photonE < 0. || photonE > eMax || photonE < eMax*
G4UniformRand()) );
543 if (iteration >= maxDopplerIterations)
545 photonE = photonEoriginal;
549 gammaEnergy1 = photonE;
569 gammaDirection1 = tmpDirection1;
574 gammaPolarization0,gammaPolarization1);
576 if (gammaEnergy1 > 0.)
593 G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;
597 if(ElecKineEnergy < 0.0) {
606 G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 -
607 gammaEnergy1 * gammaDirection1) * (1./ElecMomentum));
610 fvect->push_back(dp);
622 size_t nbefore = fvect->size();
626 size_t nafter = fvect->size();
627 if(nafter > nbefore) {
628 for (
size_t i=nbefore; i<nafter; ++i) {
630 if (bindingE >= ((*fvect)[i])->GetKineticEnergy())
633 bindingE -= ((*fvect)[i])->GetKineticEnergy();
648 G4Exception(
"G4LivermoreComptonModel::SampleSecondaries()",
674 b = energyRate + 1/energyRate;
676 phiProbability = 1 - (a/
b)*(std::cos(phi)*std::cos(phi));
681 while ( rand2 > phiProbability );
718 c.
setX(std::cos(angle)*(a0.
x())+std::sin(angle)*b0.
x());
719 c.
setY(std::cos(angle)*(a0.
y())+std::sin(angle)*b0.
y());
720 c.
setZ(std::cos(angle)*(a0.
z())+std::sin(angle)*b0.
z());
744 return gammaPolarization - gammaPolarization.
dot(gammaDirection)/gammaDirection.
dot(gammaDirection) * gammaDirection;
758 G4double sinTheta = std::sqrt(sinSqrTh);
762 G4double normalisation = std::sqrt(1. - cosSqrPhi*sinSqrTh);
801 if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsilon+1.0/epsilon)-4.0*sinSqrTh*cosSqrPhi))
816 G4double sinBeta = std::sqrt(1-cosBeta*cosBeta);
820 G4double xParallel = normalisation*cosBeta;
821 G4double yParallel = -(sinSqrTh*cosPhi*sinPhi)*cosBeta/normalisation;
822 G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation;
824 G4double yPerpendicular = (costheta)*sinBeta/normalisation;
825 G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation;
827 G4double xTotal = (xParallel + xPerpendicular);
828 G4double yTotal = (yParallel + yPerpendicular);
829 G4double zTotal = (zParallel + zPerpendicular);
831 gammaPolarization1.
setX(xTotal);
832 gammaPolarization1.
setY(yTotal);
833 gammaPolarization1.
setZ(zTotal);
835 return gammaPolarization1;
857 direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
862 polarization1 = (polarization_x*Axis_X0 + polarization_y*Axis_Y0 + polarization_z*Axis_Z0).unit();
877 G4AutoLock l(&LivermorePolarizedComptonModelMutex);