G4PolarizedComptonModel Class Reference

#include <G4PolarizedComptonModel.hh>

Inheritance diagram for G4PolarizedComptonModel:

Public Member Functions
	G4PolarizedComptonModel (const G4ParticleDefinition *p=0, const G4String &nam="Polarized-Compton")
virtual G4double	ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A, G4double cut, G4double emax)
G4double	ComputeAsymmetryPerAtom (G4double gammaEnergy, G4double Z)
virtual	~G4PolarizedComptonModel ()
virtual void	SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
void	SetTargetPolarization (const G4ThreeVector &pTarget)
void	SetBeamPolarization (const G4ThreeVector &pBeam)
const G4ThreeVector &	GetTargetPolarization () const
const G4ThreeVector &	GetBeamPolarization () const
const G4ThreeVector &	GetFinalGammaPolarization () const
const G4ThreeVector &	GetFinalElectronPolarization () const

---

Detailed Description

Definition at line 61 of file G4PolarizedComptonModel.hh.

---

Constructor & Destructor Documentation

G4PolarizedComptonModel::G4PolarizedComptonModel	(	const G4ParticleDefinition *	p = 0,
		const G4String &	nam = "Polarized-Compton"
	)

Definition at line 67 of file G4PolarizedComptonModel.cc.

  : G4KleinNishinaCompton(0,nam),
    verboseLevel(0)
{
  crossSectionCalculator=new G4PolarizedComptonCrossSection();
}

G4PolarizedComptonModel::~G4PolarizedComptonModel

(

)

[virtual]

Definition at line 77 of file G4PolarizedComptonModel.cc.

{
  if (crossSectionCalculator) delete crossSectionCalculator;
}

---

Member Function Documentation

G4double G4PolarizedComptonModel::ComputeAsymmetryPerAtom	(	G4double	gammaEnergy,
		G4double	Z
	)

Definition at line 85 of file G4PolarizedComptonModel.cc.

References G4cout, G4endl, G4InuclParticleNames::k0, and sqr().

Referenced by ComputeCrossSectionPerAtom().

{
 G4double asymmetry = 0.0 ;

 G4double k0 = gammaEnergy / electron_mass_c2 ;
 G4double k1 = 1 + 2*k0 ;

 asymmetry = -k0;
 asymmetry *= (k0 + 1.)*sqr(k1)*std::log(k1) - 2.*k0*(5.*sqr(k0) + 4.*k0 + 1.);
 asymmetry /= ((k0 - 2.)*k0  -2.)*sqr(k1)*std::log(k1) + 2.*k0*(k0*(k0 + 1.)*(k0 + 8.) + 2.);           

 // G4cout<<"energy = "<<GammaEnergy<<"  asymmetry = "<<asymmetry<<"\t\t GAM = "<<k0<<G4endl;
 if (asymmetry>1.) G4cout<<"ERROR in G4PolarizedComptonModel::ComputeAsymmetryPerAtom"<<G4endl;

 return asymmetry;
}

G4double G4PolarizedComptonModel::ComputeCrossSectionPerAtom	(	const G4ParticleDefinition *	,
		G4double	kinEnergy,
		G4double	Z,
		G4double	A,
		G4double	cut,
		G4double	emax
	)			[virtual]

Reimplemented from G4KleinNishinaCompton.

Definition at line 104 of file G4PolarizedComptonModel.cc.

References ComputeAsymmetryPerAtom(), G4KleinNishinaCompton::ComputeCrossSectionPerAtom(), and G4StokesVector::p3().

{
  double xs = 
    G4KleinNishinaCompton::ComputeCrossSectionPerAtom(pd,kinEnergy,
                                                      Z,A,cut,emax);
  G4double polzz = theBeamPolarization.p3()*theTargetPolarization.z();
  if (polzz!=0) {
    G4double asym=ComputeAsymmetryPerAtom(kinEnergy, Z);  
    xs*=(1.+polzz*asym);
  }
  return xs;
}

const G4ThreeVector & G4PolarizedComptonModel::GetBeamPolarization

(

)

const [inline]

Definition at line 123 of file G4PolarizedComptonModel.hh.

{
  return theBeamPolarization;
}

const G4ThreeVector & G4PolarizedComptonModel::GetFinalElectronPolarization

(

)

const [inline]

Definition at line 131 of file G4PolarizedComptonModel.hh.

{
  return finalElectronPolarization;
}

const G4ThreeVector & G4PolarizedComptonModel::GetFinalGammaPolarization

(

)

const [inline]

Definition at line 127 of file G4PolarizedComptonModel.hh.

{
  return finalGammaPolarization;
}

const G4ThreeVector & G4PolarizedComptonModel::GetTargetPolarization

(

)

const [inline]

Definition at line 119 of file G4PolarizedComptonModel.hh.

{
  return theTargetPolarization;
}

void G4PolarizedComptonModel::SampleSecondaries	(	std::vector< G4DynamicParticle * > *	,
		const G4MaterialCutsCouple *	,
		const G4DynamicParticle *	,
		G4double	tmin,
		G4double	maxEnergy
	)			[virtual]

Reimplemented from G4KleinNishinaCompton.

Definition at line 126 of file G4PolarizedComptonModel.cc.

References DBL_MIN, G4KleinNishinaCompton::fParticleChange, fStopAndKill, G4cout, G4endl, G4UniformRand, G4ParticleChangeForGamma::GetCurrentTrack(), G4PolarizationHelper::GetFrame(), G4PolarizationManager::GetInstance(), G4DynamicParticle::GetKineticEnergy(), G4VParticleChange::GetLocalEnergyDeposit(), G4VPhysicalVolume::GetLogicalVolume(), G4DynamicParticle::GetMomentumDirection(), G4LogicalVolume::GetName(), G4PolarizedComptonCrossSection::GetPol2(), G4PolarizedComptonCrossSection::GetPol3(), G4DynamicParticle::GetPolarization(), G4Track::GetVolume(), G4PolarizationManager::GetVolumePolarization(), G4PolarizedComptonCrossSection::Initialize(), G4StokesVector::InvRotateAz(), G4PolarizationManager::IsPolarized(), G4KleinNishinaCompton::lowestGammaEnergy, G4StokesVector::p1(), G4StokesVector::p2(), G4StokesVector::p3(), G4VParticleChange::ProposeLocalEnergyDeposit(), G4ParticleChangeForGamma::ProposeMomentumDirection(), G4ParticleChangeForGamma::ProposePolarization(), G4VParticleChange::ProposeTrackStatus(), G4StokesVector::RotateAz(), G4StokesVector::SetPhoton(), G4DynamicParticle::SetPolarization(), G4ParticleChangeForGamma::SetProposedKineticEnergy(), sqr(), and G4KleinNishinaCompton::theElectron.

{
  const G4Track * aTrack = fParticleChange->GetCurrentTrack();
  G4VPhysicalVolume*  aPVolume  = aTrack->GetVolume();
  G4LogicalVolume*    aLVolume  = aPVolume->GetLogicalVolume();

  if (verboseLevel>=1) 
    G4cout<<"G4PolarizedComptonModel::SampleSecondaries in "
          <<  aLVolume->GetName() <<G4endl;

  G4PolarizationManager * polarizationManager = G4PolarizationManager::GetInstance();

  // obtain polarization of the beam
  theBeamPolarization =  aDynamicGamma->GetPolarization();
  theBeamPolarization.SetPhoton();

  // obtain polarization of the media
  const G4bool targetIsPolarized = polarizationManager->IsPolarized(aLVolume);
  theTargetPolarization = polarizationManager->GetVolumePolarization(aLVolume);

  // if beam is linear polarized or target is transversely polarized 
  // determine the angle to x-axis
  // (assumes same PRF as in the polarization definition)

  G4ThreeVector gamDirection0 = aDynamicGamma->GetMomentumDirection();

  // transfere theTargetPolarization 
  // into the gamma frame (problem electron is at rest)
  if (targetIsPolarized)
    theTargetPolarization.rotateUz(gamDirection0);

  // The scattered gamma energy is sampled according to Klein - Nishina formula.
  // The random number techniques of Butcher & Messel are used 
  // (Nuc Phys 20(1960),15).
  // Note : Effects due to binding of atomic electrons are negliged.
 
  G4double gamEnergy0 = aDynamicGamma->GetKineticEnergy();
  G4double E0_m = gamEnergy0 / electron_mass_c2 ;


  //
  // sample the energy rate of the scattered gamma 
  //

  G4double epsilon, epsilonsq, onecost, sint2, greject ;

  G4double eps0       = 1./(1. + 2.*E0_m);
  G4double epsilon0sq = eps0*eps0;
  G4double alpha1     = - std::log(eps0);
  G4double alpha2     = 0.5*(1.- epsilon0sq);

  G4double polarization = theBeamPolarization.p3()*theTargetPolarization.p3();
  do {
    if ( alpha1/(alpha1+alpha2) > G4UniformRand() ) {
      epsilon   = std::exp(-alpha1*G4UniformRand());   // epsilon0**r
      epsilonsq = epsilon*epsilon; 

    } else {
      epsilonsq = epsilon0sq + (1.- epsilon0sq)*G4UniformRand();
      epsilon   = std::sqrt(epsilonsq);
    };

    onecost = (1.- epsilon)/(epsilon*E0_m);
    sint2   = onecost*(2.-onecost);


    G4double gdiced = 2.*(1./epsilon+epsilon);
    G4double gdist  = 1./epsilon + epsilon - sint2 
      - polarization*(1./epsilon-epsilon)*(1.-onecost);

    greject = gdist/gdiced;

    if (greject>1) G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n"
                         <<" costh rejection does not work properly: "<<greject<<G4endl;

  } while (greject < G4UniformRand());
 
  //
  // scattered gamma angles. ( Z - axis along the parent gamma)
  //

  G4double cosTeta = 1. - onecost; 
  G4double sinTeta = std::sqrt (sint2);
  G4double Phi;
  do {
    Phi     = twopi * G4UniformRand();
     G4double gdiced = 1./epsilon + epsilon - sint2 
       + std::abs(theBeamPolarization.p3())*
       ( std::abs((1./epsilon-epsilon)*cosTeta*theTargetPolarization.p3())
        +(1.-epsilon)*sinTeta*(std::sqrt(sqr(theTargetPolarization.p1()) 
                                    + sqr(theTargetPolarization.p2()))))
       +sint2*(std::sqrt(sqr(theBeamPolarization.p1()) + sqr(theBeamPolarization.p2())));

     G4double gdist = 1./epsilon + epsilon - sint2 
       + theBeamPolarization.p3()*
       ((1./epsilon-epsilon)*cosTeta*theTargetPolarization.p3()
        +(1.-epsilon)*sinTeta*(std::cos(Phi)*theTargetPolarization.p1()+
                               std::sin(Phi)*theTargetPolarization.p2()))
       -sint2*(std::cos(2.*Phi)*theBeamPolarization.p1()
               +std::sin(2.*Phi)*theBeamPolarization.p2());
     greject = gdist/gdiced;

    if (greject>1.+1.e-10 || greject<0) G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n"
                                      <<" phi rejection does not work properly: "<<greject<<G4endl;

    if (greject<1.e-3) {
      G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n"
            <<" phi rejection does not work properly: "<<greject<<"\n";
      G4cout<<" greject="<<greject<<"  phi="<<Phi<<"   cost="<<cosTeta<<"\n";
      G4cout<<" gdiced="<<gdiced<<"   gdist="<<gdist<<"\n";
      G4cout<<" eps="<<epsilon<<"    1/eps="<<1./epsilon<<"\n";
    }
     
  } while (greject < G4UniformRand());
  G4double dirx = sinTeta*std::cos(Phi), diry = sinTeta*std::sin(Phi), dirz = cosTeta;

  //
  // update G4VParticleChange for the scattered gamma
  //
   
  G4ThreeVector gamDirection1 ( dirx,diry,dirz );
  gamDirection1.rotateUz(gamDirection0);
  G4double gamEnergy1 = epsilon*gamEnergy0;
  fParticleChange->SetProposedKineticEnergy(gamEnergy1);


  if(gamEnergy1 > lowestGammaEnergy) {
    fParticleChange->ProposeMomentumDirection(gamDirection1);
  } else { 
    fParticleChange->ProposeTrackStatus(fStopAndKill);
    gamEnergy1 += fParticleChange->GetLocalEnergyDeposit();
    fParticleChange->ProposeLocalEnergyDeposit(gamEnergy1);
  }
 
  //
  // kinematic of the scattered electron
  //

  G4double eKinEnergy = gamEnergy0 - gamEnergy1;
  G4ThreeVector eDirection = gamEnergy0*gamDirection0 - gamEnergy1*gamDirection1;
  eDirection = eDirection.unit();

  // 
  // calculate Stokesvector of final state photon and electron
  //
  G4ThreeVector  nInteractionFrame;
  if((gamEnergy1 > lowestGammaEnergy) ||
     (eKinEnergy > DBL_MIN)) {

    // determine interaction plane
//     nInteractionFrame = 
//       G4PolarizationHelper::GetFrame(gamDirection1,eDirection);
    if (gamEnergy1 > lowestGammaEnergy) 
      nInteractionFrame = G4PolarizationHelper::GetFrame(gamDirection1,gamDirection0);
    else 
      nInteractionFrame = G4PolarizationHelper::GetFrame(gamDirection0, eDirection);

    // transfere theBeamPolarization and theTargetPolarization 
    // into the interaction frame (note electron is in gamma frame)
    if (verboseLevel>=1) {
      G4cout << "========================================\n";
      G4cout << " nInteractionFrame = " <<nInteractionFrame<<"\n";
      G4cout << " GammaDirection0 = " <<gamDirection0<<"\n";
      G4cout << " gammaPolarization = " <<theBeamPolarization<<"\n";
      G4cout << " electronPolarization = " <<theTargetPolarization<<"\n";
    }

    theBeamPolarization.InvRotateAz(nInteractionFrame,gamDirection0);
    theTargetPolarization.InvRotateAz(nInteractionFrame,gamDirection0);

    if (verboseLevel>=1) {
      G4cout << "----------------------------------------\n";
      G4cout << " gammaPolarization = " <<theBeamPolarization<<"\n";
      G4cout << " electronPolarization = " <<theTargetPolarization<<"\n";
      G4cout << "----------------------------------------\n";
    }

    // initialize the polarization transfer matrix
    crossSectionCalculator->Initialize(epsilon,E0_m,0.,
                                       theBeamPolarization,
                                       theTargetPolarization,2);
  }

  //  if(eKinEnergy > DBL_MIN)
  {
    // in interaction frame
    // calculate polarization transfer to the photon (in interaction plane)
    finalGammaPolarization = crossSectionCalculator->GetPol2();
    if (verboseLevel>=1) G4cout << " gammaPolarization1 = " <<finalGammaPolarization<<"\n";
    finalGammaPolarization.SetPhoton();

    // translate polarization into particle reference frame
    finalGammaPolarization.RotateAz(nInteractionFrame,gamDirection1);
    //store polarization vector
    fParticleChange->ProposePolarization(finalGammaPolarization);
    if (finalGammaPolarization.mag() > 1.+1.e-8){
      G4cout<<"ERROR in Polarizaed Compton Scattering !"<<G4endl;
      G4cout<<"Polarization of final photon more than 100%"<<G4endl;
      G4cout<<finalGammaPolarization<<" mag = "<<finalGammaPolarization.mag()<<G4endl;
    }
    if (verboseLevel>=1) {
      G4cout << " gammaPolarization1 = " <<finalGammaPolarization<<"\n";
      G4cout << " GammaDirection1 = " <<gamDirection1<<"\n";
    }
  }

  //    if (ElecKineEnergy > fminimalEnergy) {
  {
    finalElectronPolarization = crossSectionCalculator->GetPol3();
    if (verboseLevel>=1) 
      G4cout << " electronPolarization1 = " <<finalElectronPolarization<<"\n";

    // transfer into particle reference frame
    finalElectronPolarization.RotateAz(nInteractionFrame,eDirection);
    if (verboseLevel>=1) {
      G4cout << " electronPolarization1 = " <<finalElectronPolarization<<"\n";
      G4cout << " ElecDirection = " <<eDirection<<"\n";
    }
  }
  if (verboseLevel>=1)
    G4cout << "========================================\n";
       

  if(eKinEnergy > DBL_MIN) {

    // create G4DynamicParticle object for the electron.
    G4DynamicParticle* aElectron = new G4DynamicParticle(theElectron,eDirection,eKinEnergy);
    //store polarization vector
    if (finalElectronPolarization.mag() > 1.+1.e-8){
      G4cout<<"ERROR in Polarizaed Compton Scattering !"<<G4endl;
      G4cout<<"Polarization of final electron more than 100%"<<G4endl;
      G4cout<<finalElectronPolarization<<" mag = "<<finalElectronPolarization.mag()<<G4endl;
    }
    aElectron->SetPolarization(finalElectronPolarization.p1(),
                               finalElectronPolarization.p2(),
                               finalElectronPolarization.p3());
    fvect->push_back(aElectron);
  }
}

void G4PolarizedComptonModel::SetBeamPolarization

(

const G4ThreeVector &

pBeam

)

[inline]

Definition at line 115 of file G4PolarizedComptonModel.hh.

Referenced by G4PolarizedCompton::ComputeAsymmetry().

{
  theBeamPolarization = pBeam;
}

void G4PolarizedComptonModel::SetTargetPolarization

(

const G4ThreeVector &

pTarget

)

[inline]

Definition at line 111 of file G4PolarizedComptonModel.hh.

Referenced by G4PolarizedCompton::ComputeAsymmetry().

{
  theTargetPolarization = pTarget;
}

---

The documentation for this class was generated from the following files:

• G4PolarizedComptonModel.hh
• G4PolarizedComptonModel.cc

---