141 if (verboseLevel >= 1) {
142 G4cout<<
"G4PolarizedComptonModel::SampleSecondaries in "
154 theTargetPolarization =
165 if (targetIsPolarized) {
166 theTargetPolarization.
rotateUz(gamDirection0);
191 G4double alpha2 = alpha1 + 0.5*(1.- epsilon0sq);
194 theBeamPolarization.
p3()*theTargetPolarization.
p3();
207 PrintWarning(aDynamicGamma, nloop, greject, onecost, Phi,
208 "too many iterations");
215 if ( alpha1 > alpha2*rndm[0]) {
216 epsilon =
G4Exp(-alpha1*rndm[1]);
218 epsilon = std::sqrt(epsilon0sq + (1.- epsilon0sq)*rndm[1]);
221 onecost = (1.-
epsilon)/(epsilon*E0_m);
222 sint2 = onecost*(2.-onecost);
225 G4double gdist = 1./epsilon + epsilon - sint2
226 - polarization*(1./epsilon-
epsilon)*(1.-onecost);
228 greject = gdist/gdiced;
231 PrintWarning(aDynamicGamma, nloop, greject, onecost, Phi,
232 "theta majoranta wrong");
235 }
while (greject < rndm[2]);
243 cosTeta = 1. - onecost;
244 sinTeta = std::sqrt(sint2);
252 Phi =
twopi * rndm[0];
254 PrintWarning(aDynamicGamma, nloop, greject, onecost, Phi,
255 "too many iterations");
259 G4double gdiced = 1./epsilon + epsilon - sint2
260 + std::abs(theBeamPolarization.
p3())*
261 ( std::abs((1./epsilon-epsilon)*cosTeta*theTargetPolarization.
p3())
262 +(1.-epsilon)*sinTeta*(std::sqrt(
sqr(theTargetPolarization.
p1())
263 +
sqr(theTargetPolarization.
p2()))))
264 +sint2*(std::sqrt(
sqr(theBeamPolarization.
p1()) +
265 sqr(theBeamPolarization.
p2())));
267 G4double gdist = 1./epsilon + epsilon - sint2
268 + theBeamPolarization.
p3()*
269 ((1./epsilon-
epsilon)*cosTeta*theTargetPolarization.
p3()
270 +(1.-
epsilon)*sinTeta*(std::cos(Phi)*theTargetPolarization.
p1()+
271 std::sin(Phi)*theTargetPolarization.
p2()))
272 -sint2*(std::cos(2.*Phi)*theBeamPolarization.
p1()
273 +std::sin(2.*Phi)*theBeamPolarization.
p2());
274 greject = gdist/gdiced;
277 PrintWarning(aDynamicGamma, nloop, greject, onecost, Phi,
278 "phi majoranta wrong");
281 if(greject < 1.e-3) {
282 PrintWarning(aDynamicGamma, nloop, greject, onecost, Phi,
283 "phi loop ineffective");
290 }
while (greject < rndm[1]);
292 G4double dirx = sinTeta*std::cos(Phi), diry = sinTeta*std::sin(Phi),
300 gamDirection1.rotateUz(gamDirection0);
301 G4double gamEnergy1 = epsilon*gamEnergy0;
321 if (verboseLevel>=1) {
322 G4cout <<
"========================================\n";
323 G4cout <<
" nInteractionFrame = " <<nInteractionFrame<<
"\n";
324 G4cout <<
" GammaDirection0 = " <<gamDirection0<<
"\n";
325 G4cout <<
" gammaPolarization = " <<theBeamPolarization<<
"\n";
326 G4cout <<
" electronPolarization = " <<theTargetPolarization<<
"\n";
329 theBeamPolarization.
InvRotateAz(nInteractionFrame,gamDirection0);
330 theTargetPolarization.
InvRotateAz(nInteractionFrame,gamDirection0);
332 if (verboseLevel>=1) {
333 G4cout <<
"----------------------------------------\n";
334 G4cout <<
" gammaPolarization = " <<theBeamPolarization<<
"\n";
335 G4cout <<
" electronPolarization = " <<theTargetPolarization<<
"\n";
336 G4cout <<
"----------------------------------------\n";
340 crossSectionCalculator->
Initialize(epsilon,E0_m,0.,
342 theTargetPolarization,2);
348 finalGammaPolarization = crossSectionCalculator->
GetPol2();
349 if (verboseLevel>=1) {
350 G4cout <<
" gammaPolarization1 = " <<finalGammaPolarization<<
"\n";
355 finalGammaPolarization.
RotateAz(nInteractionFrame,gamDirection1);
356 if (finalGammaPolarization.
mag() > 1.+1.e-8){
357 G4cout<<
"ERROR in Polarizaed Compton Scattering !"<<
G4endl;
358 G4cout<<
"Polarization of final photon more than 100%"<<
G4endl;
359 G4cout<<finalGammaPolarization<<
" mag = "
364 if (verboseLevel>=1) {
365 G4cout <<
" gammaPolarization1 = " <<finalGammaPolarization<<
"\n";
366 G4cout <<
" GammaDirection1 = " <<gamDirection1<<
"\n";
373 G4double eKinEnergy = gamEnergy0 - gamEnergy1;
378 gamEnergy0*gamDirection0 - gamEnergy1*gamDirection1;
379 eDirection = eDirection.
unit();
381 finalElectronPolarization = crossSectionCalculator->
GetPol3();
382 if (verboseLevel>=1) {
383 G4cout <<
" electronPolarization1 = "
384 <<finalElectronPolarization<<
"\n";
387 finalElectronPolarization.
RotateAz(nInteractionFrame,eDirection);
388 if (verboseLevel>=1) {
389 G4cout <<
" electronPolarization1 = "
390 <<finalElectronPolarization<<
"\n";
391 G4cout <<
" ElecDirection = " <<eDirection<<
"\n";
398 if (finalElectronPolarization.
mag() > 1.+1.e-8){
399 G4cout<<
"ERROR in Polarizaed Compton Scattering !"<<
G4endl;
400 G4cout<<
"Polarization of final electron more than 100%"<<
G4endl;
401 G4cout<<finalElectronPolarization<<
" mag = "
402 <<finalElectronPolarization.
mag()<<
G4endl;
405 finalElectronPolarization.
p2(),
406 finalElectronPolarization.
p3());
407 fvect->push_back(aElectron);
G4double LowEnergyLimit() const
G4ParticleChangeForGamma * fParticleChange
G4double GetKineticEnergy() const
G4StokesVector GetPol3() override
static const G4int nlooplim
static G4PolarizationManager * GetInstance()
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static constexpr double twopi
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
void ProposePolarization(const G4ThreeVector &dir)
static G4ThreeVector GetFrame(const G4ThreeVector &, const G4ThreeVector &)
G4double lowestSecondaryEnergy
void SetPolarization(G4double polX, G4double polY, G4double polZ)
G4StokesVector GetPol2() override
G4double G4Log(G4double x)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
void InvRotateAz(G4ThreeVector nInteractionFrame, G4ThreeVector particleDirection)
G4LogicalVolume * GetLogicalVolume() const
bool IsPolarized(G4LogicalVolume *lVol) const
const G4ThreeVector & GetPolarization() const
G4VPhysicalVolume * GetVolume() const
const G4Track * GetCurrentTrack() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
virtual void Initialize(G4double eps, G4double X, G4double phi, const G4StokesVector &p0, const G4StokesVector &p1, G4int flag=0) override
G4ParticleDefinition * theElectron
const G4String & GetName() const
void ProposeTrackStatus(G4TrackStatus status)
virtual void flatArray(const int size, double *vect)=0
const G4ThreeVector & GetVolumePolarization(G4LogicalVolume *lVol) const
double epsilon(double density, double temperature)
void RotateAz(G4ThreeVector nInteractionFrame, G4ThreeVector particleDirection)