50 :
G4VEmModel(
"G4JAEAPolarizedElasticScatteringModel"),isInitialised(false)
66 G4cout <<
"G4JAEAPolarizedElasticScatteringModel is constructed " <<
G4endl;
96 G4cout <<
"Calling Initialise() of G4JAEAPolarizedElasticScatteringModel." <<
G4endl
109 char* path = std::getenv(
"G4LEDATA");
114 for(
G4int i=0; i<numOfCouples; ++i)
122 for (
G4int j=0; j<nelm; ++j)
152 G4cout <<
"Calling ReadData() of G4JAEAPolarizedElasticScatteringModel"
158 const char* datadir = path;
161 datadir = std::getenv(
"G4LEDATA");
164 G4Exception(
"G4JAEAPolarizedElasticScatteringModel::ReadData()",
"em0006",
166 "Environment variable G4LEDATA not defined");
171 std::ostringstream ostCS;
172 ostCS << datadir <<
"/JAEAESData/amp_Z_" <<
Z ;
173 std::ifstream ES_Data_Buffer(ostCS.str().c_str(),ios::binary);
174 if( !ES_Data_Buffer.is_open() )
177 ed <<
"G4JAEAPolarizedElasticScattering Model data file <" << ostCS.str().c_str()
178 <<
"> is not opened!" <<
G4endl;
180 ed,
"G4LEDATA version should be G4EMLOW7.11 or later. Polarized Elastic Scattering Data are not loaded");
186 G4cout <<
"File " << ostCS.str()
187 <<
" is opened by G4JAEAPolarizedElasticScatteringModel" <<
G4endl;
196 while (ES_Data_Buffer.read(
reinterpret_cast<char*
>(&buffer_var),
sizeof(
float)))
203 for (
G4int i=0;i<300;++i)
224 G4cout <<
"G4JAEAPolarizedElasticScatteringModel::ComputeCrossSectionPerAtom()"
234 if(intZ < 1 || intZ >
maxZ) {
return xs; }
243 if(!pv) {
return xs; }
251 }
else if(e >= pv->
Energy(0)) {
257 G4cout <<
"****** DEBUG: tcs value for Z=" <<
Z <<
" at energy (MeV)="
259 G4cout <<
" cs (Geant4 internal unit)=" << xs <<
G4endl;
260 G4cout <<
" -> first E*E*cs value in CS data file (iu) =" << (*pv)[0]
262 G4cout <<
" -> last E*E*cs value in CS data file (iu) =" << (*pv)[
n]
264 G4cout <<
"*********************************************************"
274 std::vector<G4DynamicParticle*>*,
281 G4cout <<
"Calling SampleSecondaries() of G4JAEAPolarizedElasticScatteringModel."
300 G4int energyindex=round(100*photonEnergy0)-1;
302 for (
G4int i=0;i<=180;++i)
320 G4double am1=0,am2=0,am3=0,am4=0,aparaSquare=0,aperpSquare=0,
321 apara_aper_Asterisk=0,img_apara_aper_Asterisk=0;
326 aparaSquare=am1*am1+am2*am2;
327 aperpSquare=am3*am3+am4*am4;
328 apara_aper_Asterisk=2*a1*a3+2*a2*a4;
329 img_apara_aper_Asterisk=2*a1*a4-2*a2*a3;
341 G4double Xi1=0, Xi2=0, Xi3=0, Xi1_Prime=0,Xi2_Prime=0,Xi3_Prime=0;
345 Xi1=gammaPolarization0.
x();
346 Xi2=gammaPolarization0.
y();
347 Xi3=gammaPolarization0.
z();
350 if ((gammaPolarization0.
mag())>1.05 || (Xi1*Xi1>1.05) || (Xi2*Xi2>1.05) || (Xi3*Xi3>1.05))
352 G4Exception(
"G4JAEAPolarizedElasticScatteringModel::SampleSecondaries()",
"em1006",
354 "WARNING: G4JAEAPolarizedElasticScatteringModel is only compatible with a unit polarization vector.");
358 if (Xi1==0 && Xi2==0 && Xi3==0)
365 Direction_Unpolarized.
setX(sin(theta)*cos(Phi_Unpolarized));
366 Direction_Unpolarized.
setY(sin(theta)*sin(Phi_Unpolarized));
367 Direction_Unpolarized.
setZ(cos(theta));
369 Xi1_Prime=(aparaSquare-aperpSquare)/(aparaSquare+aperpSquare);
370 Polarization_Unpolarized.
setX(Xi1_Prime);
371 Polarization_Unpolarized.
setY(0.);
372 Polarization_Unpolarized.
setZ(0.);
380 if(InitialAzimuth<0) InitialAzimuth=InitialAzimuth+
CLHEP::twopi;
385 aparaSquare+aperpSquare-Xi1*(aparaSquare-aperpSquare),InitialAzimuth);
387 Xi1_Prime=((aparaSquare-aperpSquare)+Xi1*(aparaSquare+aperpSquare)*cos(2*Phi_Linear1))/
388 ((aparaSquare+aperpSquare)+Xi1*(aparaSquare-aperpSquare)*cos(2*Phi_Linear1));
389 Xi2_Prime=(-Xi1*apara_aper_Asterisk*sin(2*Phi_Linear1))/
390 ((aparaSquare+aperpSquare)+Xi1*(aparaSquare-aperpSquare)*cos(2*Phi_Linear1));
391 Xi3_Prime=(-Xi1*img_apara_aper_Asterisk*sin(2*Phi_Linear1))/
392 ((aparaSquare+aperpSquare)+Xi1*(aparaSquare-aperpSquare)*cos(2*Phi_Linear1));
394 Direction_Linear1.
setX(sin(theta)*cos(Phi_Linear1));
395 Direction_Linear1.
setY(sin(theta)*sin(Phi_Linear1));
396 Direction_Linear1.
setZ(cos(theta));
397 Polarization_Linear1.
setX(Xi1_Prime);
398 Polarization_Linear1.
setY(Xi2_Prime);
399 Polarization_Linear1.
setZ(Xi3_Prime);
407 if(abs(Xi1)>0.0) dsigmaL1=0.25*((aparaSquare+aperpSquare)*
408 (1+Xi1*Xi1_Prime*cos(2*Phi_Linear1))+
409 (aparaSquare-aperpSquare)*(Xi1*cos(2*Phi_Linear1)+Xi1_Prime)
410 -Xi1*Xi2_Prime*apara_aper_Asterisk*sin(2*Phi_Linear1)-
411 Xi1*Xi3_Prime*img_apara_aper_Asterisk*sin(2*Phi_Linear1));
417 InitialAzimuth=InitialAzimuth-
CLHEP::pi/4.;
418 if(InitialAzimuth<0) InitialAzimuth=InitialAzimuth+
CLHEP::twopi;
421 ,aparaSquare+aperpSquare-Xi1*(aparaSquare-aperpSquare),InitialAzimuth);
423 Xi1_Prime=((aparaSquare-aperpSquare)+Xi2*(aparaSquare+aperpSquare)*sin(2*Phi_Linear2))/
424 ((aparaSquare+aperpSquare)+Xi2*(aparaSquare-aperpSquare)*sin(2*Phi_Linear2));
425 Xi2_Prime=(Xi2*apara_aper_Asterisk*cos(2*Phi_Linear2))/
426 ((aparaSquare+aperpSquare)+Xi2*(aparaSquare-aperpSquare)*sin(2*Phi_Linear2));
427 Xi3_Prime=(Xi2*img_apara_aper_Asterisk*cos(2*Phi_Linear2))/
428 ((aparaSquare+aperpSquare)+Xi2*(aparaSquare-aperpSquare)*sin(2*Phi_Linear2));
430 Direction_Linear2.
setX(sin(theta)*cos(Phi_Linear2));
431 Direction_Linear2.
setY(sin(theta)*sin(Phi_Linear2));
432 Direction_Linear2.
setZ(cos(theta));
433 Polarization_Linear2.
setX(Xi1_Prime);
434 Polarization_Linear2.
setY(Xi2_Prime);
435 Polarization_Linear2.
setZ(Xi3_Prime);
444 dsigmaL2=0.25*((aparaSquare+aperpSquare)*(1+Xi2*Xi1_Prime*sin(2*Phi_Linear2))+
445 (aparaSquare-aperpSquare)*(Xi2*sin(2*Phi_Linear2)+Xi1_Prime)
446 +Xi2*Xi2_Prime*apara_aper_Asterisk*cos(2*Phi_Linear2)-
447 Xi2*Xi3_Prime*img_apara_aper_Asterisk*cos(2*Phi_Linear2));
453 Xi1_Prime=(aparaSquare-aperpSquare)/(aparaSquare+aperpSquare);
454 Xi2_Prime=(-Xi3*img_apara_aper_Asterisk)/(aparaSquare+aperpSquare);
455 Xi3_Prime=(Xi3*apara_aper_Asterisk)/(aparaSquare+aperpSquare);
457 Polarization_Circular.
setX(Xi1_Prime);
458 Polarization_Circular.
setY(Xi2_Prime);
459 Polarization_Circular.
setZ(Xi3_Prime);
468 dsigmaC=0.25*(aparaSquare+aperpSquare+Xi1_Prime*(aparaSquare-aperpSquare)-
469 Xi3*Xi2_Prime*img_apara_aper_Asterisk
470 +Xi3*Xi3_Prime*apara_aper_Asterisk);
472 if (abs(Xi3)==0.0 && abs(Xi1_Prime)==0.0)
474 Direction_Circular.
setX(sin(theta)*cos(Phi_Circular));
475 Direction_Circular.
setY(sin(theta)*sin(Phi_Circular));
476 Direction_Circular.
setZ(cos(theta));
481 for (
G4int i=0;i<=180;++i)
488 Xi1_Prime*(c1*c1+c2*c2-c3*c3-c4*c4)-
489 Xi3*Xi2_Prime*(2*c1*c4-2*c2*c3)
490 +Xi3*Xi3_Prime*(2*c1*c4-2*c2*c3));
494 Direction_Circular.
setX(sin(Theta_Circular)*cos(Phi_Circular));
495 Direction_Circular.
setY(sin(Theta_Circular)*sin(Phi_Circular));
496 Direction_Circular.
setZ(cos(Theta_Circular));
500 G4double totalSigma= dsigmaL1+dsigmaL2+dsigmaC;
506 if (abs(probc - dsigmaC/totalSigma)>=0.0001)
508 G4Exception(
"G4JAEAPolarizedElasticScatteringModel::SampleSecondaries()",
"em1007",
510 "WARNING: Polarization mixing might be incorrect.");
521 finaldirection.
setX(Direction_Linear1.
x());
522 finaldirection.
setY(Direction_Linear1.
y());
523 finaldirection.
setZ(Direction_Linear1.
z());
524 outcomingPhotonPolarization.
setX(Polarization_Linear1.
x());
525 outcomingPhotonPolarization.
setY(Polarization_Linear1.
y());
526 outcomingPhotonPolarization.
setZ(Polarization_Linear1.
z());
528 else if ((polmix>prob1) && (polmix<=prob1+prob2))
530 finaldirection.
setX(Direction_Linear2.
x());
531 finaldirection.
setY(Direction_Linear2.
y());
532 finaldirection.
setZ(Direction_Linear2.
z());
533 outcomingPhotonPolarization.
setX(Polarization_Linear2.
x());
534 outcomingPhotonPolarization.
setY(Polarization_Linear2.
y());
535 outcomingPhotonPolarization.
setZ(Polarization_Linear2.
z());
537 else if (polmix>prob1+prob2)
539 finaldirection.
setX(Direction_Circular.
x());
540 finaldirection.
setY(Direction_Circular.
y());
541 finaldirection.
setZ(Direction_Circular.
z());
542 outcomingPhotonPolarization.
setX(Polarization_Circular.
x());
543 outcomingPhotonPolarization.
setY(Polarization_Circular.
y());
544 outcomingPhotonPolarization.
setZ(Polarization_Circular.
z());
563 G4double Probability=Sigma_perp/(Sigma_para+Sigma_perp);
569 phiProbability = cos(phi+initial_Pol_Plane)*cos(phi+initial_Pol_Plane);
579 phiProbability = sin(phi+initial_Pol_Plane)*sin(phi+initial_Pol_Plane);
std::vector< const G4Element * > G4ElementVector
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
static constexpr double keV
static constexpr double eV
static constexpr double GeV
#define G4MUTEX_INITIALIZER
G4GLOB_DLL std::ostream G4cout
Hep3Vector & rotateUz(const Hep3Vector &)
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
const G4ThreeVector & GetPolarization() const
G4JAEAPolarizedElasticScatteringModel()
void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel) override
G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX) override
G4bool fCircularPolarizationSensitvity
void ReadData(size_t Z, const char *path=0)
G4bool fLinearPolarizationSensitvity2
static G4DataVector * Polarized_ES_Data[maxZ+1]
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
void InitialiseForElement(const G4ParticleDefinition *, G4int Z) override
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
virtual ~G4JAEAPolarizedElasticScatteringModel()
G4bool fLinearPolarizationSensitvity1
G4ParticleChangeForGamma * fParticleChange
static G4PhysicsFreeVector * dataCS[maxZ+1]
G4double cdistribution[181]
G4double GeneratePolarizedPhi(G4double Sigma_para, G4double Sigma_perp, G4double initial_Pol_Plane)
G4double distribution[181]
const G4Material * GetMaterial() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposePolarization(const G4ThreeVector &dir)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4double Energy(const std::size_t index) const
G4double Value(const G4double energy, std::size_t &lastidx) const
std::size_t GetVectorLength() const
void FillSecondDerivatives(const G4SplineType=G4SplineType::Base, const G4double dir1=0.0, const G4double dir2=0.0)
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
std::size_t GetTableSize() const
static G4ProductionCutsTable * GetProductionCutsTable()
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
G4double LowEnergyLimit() const
std::vector< G4EmElementSelector * > * GetElementSelectors()
G4double HighEnergyLimit() const
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static constexpr double twopi
static constexpr double pi
G4Mutex G4JAEAPolarizedElasticScatteringModelMutex