G4PolarizedComptonModel Class Reference

#include <G4PolarizedComptonModel.hh>

Inheritance diagram for G4PolarizedComptonModel:

G4KleinNishinaCompton G4VEmModel

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 G4ThreeVectorGetTargetPolarization () const
const G4ThreeVectorGetBeamPolarization () const
const G4ThreeVectorGetFinalGammaPolarization () const
const G4ThreeVectorGetFinalElectronPolarization () 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.

00069   : G4KleinNishinaCompton(0,nam),
00070     verboseLevel(0)
00071 {
00072   crossSectionCalculator=new G4PolarizedComptonCrossSection();
00073 }

G4PolarizedComptonModel::~G4PolarizedComptonModel (  )  [virtual]

Definition at line 77 of file G4PolarizedComptonModel.cc.

00078 {
00079   if (crossSectionCalculator) delete crossSectionCalculator;
00080 }


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().

00087 {
00088  G4double asymmetry = 0.0 ;
00089 
00090  G4double k0 = gammaEnergy / electron_mass_c2 ;
00091  G4double k1 = 1 + 2*k0 ;
00092 
00093  asymmetry = -k0;
00094  asymmetry *= (k0 + 1.)*sqr(k1)*std::log(k1) - 2.*k0*(5.*sqr(k0) + 4.*k0 + 1.);
00095  asymmetry /= ((k0 - 2.)*k0  -2.)*sqr(k1)*std::log(k1) + 2.*k0*(k0*(k0 + 1.)*(k0 + 8.) + 2.);           
00096 
00097  // G4cout<<"energy = "<<GammaEnergy<<"  asymmetry = "<<asymmetry<<"\t\t GAM = "<<k0<<G4endl;
00098  if (asymmetry>1.) G4cout<<"ERROR in G4PolarizedComptonModel::ComputeAsymmetryPerAtom"<<G4endl;
00099 
00100  return asymmetry;
00101 }

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().

00111 {
00112   double xs = 
00113     G4KleinNishinaCompton::ComputeCrossSectionPerAtom(pd,kinEnergy,
00114                                                       Z,A,cut,emax);
00115   G4double polzz = theBeamPolarization.p3()*theTargetPolarization.z();
00116   if (polzz!=0) {
00117     G4double asym=ComputeAsymmetryPerAtom(kinEnergy, Z);  
00118     xs*=(1.+polzz*asym);
00119   }
00120   return xs;
00121 }

const G4ThreeVector & G4PolarizedComptonModel::GetBeamPolarization (  )  const [inline]

Definition at line 123 of file G4PolarizedComptonModel.hh.

00124 {
00125   return theBeamPolarization;
00126 }

const G4ThreeVector & G4PolarizedComptonModel::GetFinalElectronPolarization (  )  const [inline]

Definition at line 131 of file G4PolarizedComptonModel.hh.

00132 {
00133   return finalElectronPolarization;
00134 }

const G4ThreeVector & G4PolarizedComptonModel::GetFinalGammaPolarization (  )  const [inline]

Definition at line 127 of file G4PolarizedComptonModel.hh.

00128 {
00129   return finalGammaPolarization;
00130 }

const G4ThreeVector & G4PolarizedComptonModel::GetTargetPolarization (  )  const [inline]

Definition at line 119 of file G4PolarizedComptonModel.hh.

00120 {
00121   return theTargetPolarization;
00122 }

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.

00131 {
00132   const G4Track * aTrack = fParticleChange->GetCurrentTrack();
00133   G4VPhysicalVolume*  aPVolume  = aTrack->GetVolume();
00134   G4LogicalVolume*    aLVolume  = aPVolume->GetLogicalVolume();
00135 
00136   if (verboseLevel>=1) 
00137     G4cout<<"G4PolarizedComptonModel::SampleSecondaries in "
00138           <<  aLVolume->GetName() <<G4endl;
00139 
00140   G4PolarizationManager * polarizationManager = G4PolarizationManager::GetInstance();
00141 
00142   // obtain polarization of the beam
00143   theBeamPolarization =  aDynamicGamma->GetPolarization();
00144   theBeamPolarization.SetPhoton();
00145 
00146   // obtain polarization of the media
00147   const G4bool targetIsPolarized = polarizationManager->IsPolarized(aLVolume);
00148   theTargetPolarization = polarizationManager->GetVolumePolarization(aLVolume);
00149 
00150   // if beam is linear polarized or target is transversely polarized 
00151   // determine the angle to x-axis
00152   // (assumes same PRF as in the polarization definition)
00153 
00154   G4ThreeVector gamDirection0 = aDynamicGamma->GetMomentumDirection();
00155 
00156   // transfere theTargetPolarization 
00157   // into the gamma frame (problem electron is at rest)
00158   if (targetIsPolarized)
00159     theTargetPolarization.rotateUz(gamDirection0);
00160 
00161   // The scattered gamma energy is sampled according to Klein - Nishina formula.
00162   // The random number techniques of Butcher & Messel are used 
00163   // (Nuc Phys 20(1960),15).
00164   // Note : Effects due to binding of atomic electrons are negliged.
00165  
00166   G4double gamEnergy0 = aDynamicGamma->GetKineticEnergy();
00167   G4double E0_m = gamEnergy0 / electron_mass_c2 ;
00168 
00169 
00170   //
00171   // sample the energy rate of the scattered gamma 
00172   //
00173 
00174   G4double epsilon, epsilonsq, onecost, sint2, greject ;
00175 
00176   G4double eps0       = 1./(1. + 2.*E0_m);
00177   G4double epsilon0sq = eps0*eps0;
00178   G4double alpha1     = - std::log(eps0);
00179   G4double alpha2     = 0.5*(1.- epsilon0sq);
00180 
00181   G4double polarization = theBeamPolarization.p3()*theTargetPolarization.p3();
00182   do {
00183     if ( alpha1/(alpha1+alpha2) > G4UniformRand() ) {
00184       epsilon   = std::exp(-alpha1*G4UniformRand());   // epsilon0**r
00185       epsilonsq = epsilon*epsilon; 
00186 
00187     } else {
00188       epsilonsq = epsilon0sq + (1.- epsilon0sq)*G4UniformRand();
00189       epsilon   = std::sqrt(epsilonsq);
00190     };
00191 
00192     onecost = (1.- epsilon)/(epsilon*E0_m);
00193     sint2   = onecost*(2.-onecost);
00194 
00195 
00196     G4double gdiced = 2.*(1./epsilon+epsilon);
00197     G4double gdist  = 1./epsilon + epsilon - sint2 
00198       - polarization*(1./epsilon-epsilon)*(1.-onecost);
00199 
00200     greject = gdist/gdiced;
00201 
00202     if (greject>1) G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n"
00203                          <<" costh rejection does not work properly: "<<greject<<G4endl;
00204 
00205   } while (greject < G4UniformRand());
00206  
00207   //
00208   // scattered gamma angles. ( Z - axis along the parent gamma)
00209   //
00210 
00211   G4double cosTeta = 1. - onecost; 
00212   G4double sinTeta = std::sqrt (sint2);
00213   G4double Phi;
00214   do {
00215     Phi     = twopi * G4UniformRand();
00216      G4double gdiced = 1./epsilon + epsilon - sint2 
00217        + std::abs(theBeamPolarization.p3())*
00218        ( std::abs((1./epsilon-epsilon)*cosTeta*theTargetPolarization.p3())
00219         +(1.-epsilon)*sinTeta*(std::sqrt(sqr(theTargetPolarization.p1()) 
00220                                     + sqr(theTargetPolarization.p2()))))
00221        +sint2*(std::sqrt(sqr(theBeamPolarization.p1()) + sqr(theBeamPolarization.p2())));
00222 
00223      G4double gdist = 1./epsilon + epsilon - sint2 
00224        + theBeamPolarization.p3()*
00225        ((1./epsilon-epsilon)*cosTeta*theTargetPolarization.p3()
00226         +(1.-epsilon)*sinTeta*(std::cos(Phi)*theTargetPolarization.p1()+
00227                                std::sin(Phi)*theTargetPolarization.p2()))
00228        -sint2*(std::cos(2.*Phi)*theBeamPolarization.p1()
00229                +std::sin(2.*Phi)*theBeamPolarization.p2());
00230      greject = gdist/gdiced;
00231 
00232     if (greject>1.+1.e-10 || greject<0) G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n"
00233                                       <<" phi rejection does not work properly: "<<greject<<G4endl;
00234 
00235     if (greject<1.e-3) {
00236       G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n"
00237             <<" phi rejection does not work properly: "<<greject<<"\n";
00238       G4cout<<" greject="<<greject<<"  phi="<<Phi<<"   cost="<<cosTeta<<"\n";
00239       G4cout<<" gdiced="<<gdiced<<"   gdist="<<gdist<<"\n";
00240       G4cout<<" eps="<<epsilon<<"    1/eps="<<1./epsilon<<"\n";
00241     }
00242      
00243   } while (greject < G4UniformRand());
00244   G4double dirx = sinTeta*std::cos(Phi), diry = sinTeta*std::sin(Phi), dirz = cosTeta;
00245 
00246   //
00247   // update G4VParticleChange for the scattered gamma
00248   //
00249    
00250   G4ThreeVector gamDirection1 ( dirx,diry,dirz );
00251   gamDirection1.rotateUz(gamDirection0);
00252   G4double gamEnergy1 = epsilon*gamEnergy0;
00253   fParticleChange->SetProposedKineticEnergy(gamEnergy1);
00254 
00255 
00256   if(gamEnergy1 > lowestGammaEnergy) {
00257     fParticleChange->ProposeMomentumDirection(gamDirection1);
00258   } else { 
00259     fParticleChange->ProposeTrackStatus(fStopAndKill);
00260     gamEnergy1 += fParticleChange->GetLocalEnergyDeposit();
00261     fParticleChange->ProposeLocalEnergyDeposit(gamEnergy1);
00262   }
00263  
00264   //
00265   // kinematic of the scattered electron
00266   //
00267 
00268   G4double eKinEnergy = gamEnergy0 - gamEnergy1;
00269   G4ThreeVector eDirection = gamEnergy0*gamDirection0 - gamEnergy1*gamDirection1;
00270   eDirection = eDirection.unit();
00271 
00272   // 
00273   // calculate Stokesvector of final state photon and electron
00274   //
00275   G4ThreeVector  nInteractionFrame;
00276   if((gamEnergy1 > lowestGammaEnergy) ||
00277      (eKinEnergy > DBL_MIN)) {
00278 
00279     // determine interaction plane
00280 //     nInteractionFrame = 
00281 //       G4PolarizationHelper::GetFrame(gamDirection1,eDirection);
00282     if (gamEnergy1 > lowestGammaEnergy) 
00283       nInteractionFrame = G4PolarizationHelper::GetFrame(gamDirection1,gamDirection0);
00284     else 
00285       nInteractionFrame = G4PolarizationHelper::GetFrame(gamDirection0, eDirection);
00286 
00287     // transfere theBeamPolarization and theTargetPolarization 
00288     // into the interaction frame (note electron is in gamma frame)
00289     if (verboseLevel>=1) {
00290       G4cout << "========================================\n";
00291       G4cout << " nInteractionFrame = " <<nInteractionFrame<<"\n";
00292       G4cout << " GammaDirection0 = " <<gamDirection0<<"\n";
00293       G4cout << " gammaPolarization = " <<theBeamPolarization<<"\n";
00294       G4cout << " electronPolarization = " <<theTargetPolarization<<"\n";
00295     }
00296 
00297     theBeamPolarization.InvRotateAz(nInteractionFrame,gamDirection0);
00298     theTargetPolarization.InvRotateAz(nInteractionFrame,gamDirection0);
00299 
00300     if (verboseLevel>=1) {
00301       G4cout << "----------------------------------------\n";
00302       G4cout << " gammaPolarization = " <<theBeamPolarization<<"\n";
00303       G4cout << " electronPolarization = " <<theTargetPolarization<<"\n";
00304       G4cout << "----------------------------------------\n";
00305     }
00306 
00307     // initialize the polarization transfer matrix
00308     crossSectionCalculator->Initialize(epsilon,E0_m,0.,
00309                                        theBeamPolarization,
00310                                        theTargetPolarization,2);
00311   }
00312 
00313   //  if(eKinEnergy > DBL_MIN)
00314   {
00315     // in interaction frame
00316     // calculate polarization transfer to the photon (in interaction plane)
00317     finalGammaPolarization = crossSectionCalculator->GetPol2();
00318     if (verboseLevel>=1) G4cout << " gammaPolarization1 = " <<finalGammaPolarization<<"\n";
00319     finalGammaPolarization.SetPhoton();
00320 
00321     // translate polarization into particle reference frame
00322     finalGammaPolarization.RotateAz(nInteractionFrame,gamDirection1);
00323     //store polarization vector
00324     fParticleChange->ProposePolarization(finalGammaPolarization);
00325     if (finalGammaPolarization.mag() > 1.+1.e-8){
00326       G4cout<<"ERROR in Polarizaed Compton Scattering !"<<G4endl;
00327       G4cout<<"Polarization of final photon more than 100%"<<G4endl;
00328       G4cout<<finalGammaPolarization<<" mag = "<<finalGammaPolarization.mag()<<G4endl;
00329     }
00330     if (verboseLevel>=1) {
00331       G4cout << " gammaPolarization1 = " <<finalGammaPolarization<<"\n";
00332       G4cout << " GammaDirection1 = " <<gamDirection1<<"\n";
00333     }
00334   }
00335 
00336   //    if (ElecKineEnergy > fminimalEnergy) {
00337   {
00338     finalElectronPolarization = crossSectionCalculator->GetPol3();
00339     if (verboseLevel>=1) 
00340       G4cout << " electronPolarization1 = " <<finalElectronPolarization<<"\n";
00341 
00342     // transfer into particle reference frame
00343     finalElectronPolarization.RotateAz(nInteractionFrame,eDirection);
00344     if (verboseLevel>=1) {
00345       G4cout << " electronPolarization1 = " <<finalElectronPolarization<<"\n";
00346       G4cout << " ElecDirection = " <<eDirection<<"\n";
00347     }
00348   }
00349   if (verboseLevel>=1)
00350     G4cout << "========================================\n";
00351        
00352 
00353   if(eKinEnergy > DBL_MIN) {
00354 
00355     // create G4DynamicParticle object for the electron.
00356     G4DynamicParticle* aElectron = new G4DynamicParticle(theElectron,eDirection,eKinEnergy);
00357     //store polarization vector
00358     if (finalElectronPolarization.mag() > 1.+1.e-8){
00359       G4cout<<"ERROR in Polarizaed Compton Scattering !"<<G4endl;
00360       G4cout<<"Polarization of final electron more than 100%"<<G4endl;
00361       G4cout<<finalElectronPolarization<<" mag = "<<finalElectronPolarization.mag()<<G4endl;
00362     }
00363     aElectron->SetPolarization(finalElectronPolarization.p1(),
00364                                finalElectronPolarization.p2(),
00365                                finalElectronPolarization.p3());
00366     fvect->push_back(aElectron);
00367   }
00368 }

void G4PolarizedComptonModel::SetBeamPolarization ( const G4ThreeVector pBeam  )  [inline]

Definition at line 115 of file G4PolarizedComptonModel.hh.

Referenced by G4PolarizedCompton::ComputeAsymmetry().

00116 {
00117   theBeamPolarization = pBeam;
00118 }

void G4PolarizedComptonModel::SetTargetPolarization ( const G4ThreeVector pTarget  )  [inline]

Definition at line 111 of file G4PolarizedComptonModel.hh.

Referenced by G4PolarizedCompton::ComputeAsymmetry().

00112 {
00113   theTargetPolarization = pTarget;
00114 }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:57 2013 for Geant4 by  doxygen 1.4.7