G4PolarizedMollerBhabhaModel Class Reference

#include <G4PolarizedMollerBhabhaModel.hh>

Inheritance diagram for G4PolarizedMollerBhabhaModel:

G4MollerBhabhaModel G4VEmModel

Public Member Functions

 G4PolarizedMollerBhabhaModel (const G4ParticleDefinition *p=0, const G4String &nam="PolarizedMollerBhabha")
virtual ~G4PolarizedMollerBhabhaModel ()
virtual G4double ComputeCrossSectionPerElectron (const G4ParticleDefinition *, G4double kinEnergy, G4double cut, G4double emax)
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 G4StokesVector & GetTargetPolarization ()
const G4StokesVector & GetBeamPolarization ()
const G4StokesVector & GetFinalElectronPolarization ()
const G4StokesVector & GetFinalPositronPolarization ()

Detailed Description

Definition at line 59 of file G4PolarizedMollerBhabhaModel.hh.

Constructor & Destructor Documentation

G4PolarizedMollerBhabhaModel::G4PolarizedMollerBhabhaModel ( const G4ParticleDefinition *  p = 0,
const G4String &  nam = "PolarizedMollerBhabha" 
)

Definition at line 65 of file G4PolarizedMollerBhabhaModel.cc.

References G4cout, G4endl, G4MollerBhabhaModel::isElectron, and G4MollerBhabhaModel::theElectron.

00067   : G4MollerBhabhaModel(p,nam)
00068 {
00069 
00070   //   G4cout<<" particle==electron "<<(p==theElectron)<<G4endl;
00071   isElectron=(p==theElectron);  // necessary due to wrong order in G4MollerBhabhaModel constructor!
00072 
00073   if (p==0) { 
00074     
00075   }
00076   if (!isElectron) {
00077     G4cout<<" buildBhabha cross section "<<isElectron<<G4endl;
00078     crossSectionCalculator = new G4PolarizedBhabhaCrossSection();
00079   }  else {
00080     G4cout<<" buildMoller cross section "<<isElectron<<G4endl;
00081     crossSectionCalculator = new G4PolarizedMollerCrossSection();
00082   }
00083 }

G4PolarizedMollerBhabhaModel::~G4PolarizedMollerBhabhaModel (  )  [virtual]

Definition at line 87 of file G4PolarizedMollerBhabhaModel.cc.

00088 {
00089   if (crossSectionCalculator) {
00090     delete crossSectionCalculator;
00091   }
00092 }

Member Function Documentation

G4double G4PolarizedMollerBhabhaModel::ComputeCrossSectionPerElectron ( const G4ParticleDefinition *  ,
G4double  kinEnergy,
G4double  cut,
G4double  emax 
) [virtual]

Reimplemented from G4MollerBhabhaModel.

Definition at line 97 of file G4PolarizedMollerBhabhaModel.cc.

References G4MollerBhabhaModel::ComputeCrossSectionPerElectron(), G4InuclParticleNames::gam, G4MollerBhabhaModel::MaxSecondaryEnergy(), and G4StokesVector::ZERO.

00102 {
00103   G4double xs = 
00104     G4MollerBhabhaModel::ComputeCrossSectionPerElectron(pd,kinEnergy,
00105                                                         cut,emax);
00106 //   G4cout<<"calc eIoni xsec "<<xs<<G4endl;
00107 //   G4cout<<" "<<kinEnergy<<" "<<cut<<" "<<emax<<G4endl;
00108   G4double factor=1.;
00109   if (xs!=0) {
00110     //    G4cout<<"calc asym"<<G4endl;
00111     G4double tmax = MaxSecondaryEnergy(pd, kinEnergy);
00112     tmax = std::min(emax, tmax);
00113 
00114     if (std::fabs(cut/emax-1.)<1.e-10) return xs;
00115 
00116     if(cut < tmax) {
00117 
00118       G4double xmin  = cut/kinEnergy;
00119       G4double xmax  = tmax/kinEnergy;
00120 //       G4cout<<"calc asym "<<xmin<<","<<xmax<<G4endl;
00121       G4double gam   = kinEnergy/electron_mass_c2 + 1.0;
00122 
00123       G4double crossPol=crossSectionCalculator->
00124         TotalXSection(xmin,xmax,gam,
00125                       theBeamPolarization,
00126                       theTargetPolarization);
00127       G4double crossUnpol=crossSectionCalculator->
00128         TotalXSection(xmin,xmax,gam,
00129                       G4StokesVector::ZERO,
00130                       G4StokesVector::ZERO);
00131       if (crossUnpol>0.)  factor=crossPol/crossUnpol;
00132       //     G4cout<<" factor="<<factor<<G4endl;
00133     }
00134   }
00135   return xs*factor;
00136 }

const G4StokesVector& G4PolarizedMollerBhabhaModel::GetBeamPolarization (  )  [inline]

Definition at line 96 of file G4PolarizedMollerBhabhaModel.hh.

00097   {
00098     return theBeamPolarization;
00099   }

const G4StokesVector& G4PolarizedMollerBhabhaModel::GetFinalElectronPolarization (  )  [inline]

Definition at line 100 of file G4PolarizedMollerBhabhaModel.hh.

00101   {
00102     return fElectronPolarization;
00103   }

const G4StokesVector& G4PolarizedMollerBhabhaModel::GetFinalPositronPolarization (  )  [inline]

Definition at line 104 of file G4PolarizedMollerBhabhaModel.hh.

00105   {
00106     return fPositronPolarization;
00107   }

const G4StokesVector& G4PolarizedMollerBhabhaModel::GetTargetPolarization (  )  [inline]

Definition at line 92 of file G4PolarizedMollerBhabhaModel.hh.

00093   {
00094     return theTargetPolarization;
00095   }

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

Reimplemented from G4MollerBhabhaModel.

Definition at line 140 of file G4PolarizedMollerBhabhaModel.cc.

References DBL_MIN, G4MollerBhabhaModel::fParticleChange, G4cout, G4endl, G4UniformRand, G4InuclParticleNames::gam, G4ParticleChangeForLoss::GetCurrentTrack(), G4PolarizationHelper::GetFrame(), G4PolarizationManager::GetInstance(), G4DynamicParticle::GetKineticEnergy(), G4VPhysicalVolume::GetLogicalVolume(), G4DynamicParticle::GetMomentumDirection(), G4VPolarizedCrossSection::GetPol2(), G4VPolarizedCrossSection::GetPol3(), G4DynamicParticle::GetPolarization(), G4Track::GetVolume(), G4PolarizationManager::GetVolumePolarization(), G4VPolarizedCrossSection::Initialize(), G4StokesVector::InvRotateAz(), G4MollerBhabhaModel::isElectron, G4PolarizationManager::IsPolarized(), G4VEmModel::MaxSecondaryKinEnergy(), G4ParticleChangeForLoss::ProposePolarization(), G4StokesVector::RotateAz(), G4DynamicParticle::SetPolarization(), G4ParticleChangeForLoss::SetProposedKineticEnergy(), G4ParticleChangeForLoss::SetProposedMomentumDirection(), sqr(), G4MollerBhabhaModel::theElectron, G4VPolarizedCrossSection::XSection(), and G4StokesVector::ZERO.

00145 {
00146   // *** obtain and save target and beam polarization ***
00147   G4PolarizationManager * polarizationManger = G4PolarizationManager::GetInstance();
00148 
00149   const G4Track * aTrack = fParticleChange->GetCurrentTrack();
00150 
00151   // obtain polarization of the beam
00152   theBeamPolarization = dp->GetPolarization();
00153 
00154   // obtain polarization of the media
00155   G4VPhysicalVolume*  aPVolume  = aTrack->GetVolume();
00156   G4LogicalVolume*    aLVolume  = aPVolume->GetLogicalVolume();
00157   const G4bool targetIsPolarized = polarizationManger->IsPolarized(aLVolume);
00158   theTargetPolarization = polarizationManger->GetVolumePolarization(aLVolume);
00159 
00160   // transfer target polarization in interaction frame
00161   if (targetIsPolarized)
00162       theTargetPolarization.rotateUz(dp->GetMomentumDirection());
00163 
00164 
00165 
00166 
00167   G4double tmax = std::min(maxEnergy, MaxSecondaryKinEnergy(dp));
00168   if(tmin >= tmax) return;
00169   //  if(tmin > tmax) tmin = tmax;
00170 
00171   G4double polL = theBeamPolarization.z()*theTargetPolarization.z();
00172                 polL=std::fabs(polL);
00173   G4double polT = theBeamPolarization.x()*theTargetPolarization.x() +
00174                   theBeamPolarization.y()*theTargetPolarization.y();
00175                 polT=std::fabs(polT);
00176 
00177   G4double kineticEnergy = dp->GetKineticEnergy();
00178   G4double energy = kineticEnergy + electron_mass_c2;
00179   G4double totalMomentum = std::sqrt(kineticEnergy*(energy + electron_mass_c2));
00180   G4double xmin   = tmin/kineticEnergy;
00181   G4double xmax   = tmax/kineticEnergy;
00182   G4double gam    = energy/electron_mass_c2;
00183   G4double gamma2 = gam*gam;
00184     G4double gmo  = gam - 1.;
00185     G4double gmo2 = gmo*gmo;
00186     G4double gmo3 = gmo2*gmo;
00187     G4double gpo  = gam + 1.;
00188     G4double gpo2 = gpo*gpo;
00189     G4double gpo3 = gpo2*gpo;
00190   G4double x, y, q, grej, grej2;
00191   G4double z = 0.;
00192   G4double xs = 0., phi =0.;
00193   G4ThreeVector direction = dp->GetMomentumDirection();
00194 
00195   //(Polarized) Moller (e-e-) scattering
00196   if (isElectron) {
00197     // *** dice according to polarized cross section
00198     G4double G = ((2.0*gam - 1.0)/gamma2)*(1. - polT - polL*gam);
00199     G4double H =  (sqr(gam - 1.0)/gamma2)*(1. + polT + polL*((gam + 3.)/(gam - 1.)));
00200 
00201     y  = 1.0 - xmax;
00202     grej  = 1.0 - G*xmax + xmax*xmax*(H + (1.0 - G*y)/(y*y));
00203     grej2 = 1.0 - G*xmin + xmin*xmin*(H + (1.0 - G*y)/(y*y));
00204     if (grej2 > grej) grej = grej2;
00205     G4double prefM = gamma2*classic_electr_radius*classic_electr_radius/(gmo2*(gam + 1.0));
00206     grej *= prefM;
00207     do {
00208       q = G4UniformRand();
00209       x = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
00210       if (crossSectionCalculator) {
00211         crossSectionCalculator->Initialize(x,gam,phi,theBeamPolarization,
00212                                            theTargetPolarization,1);
00213         xs=crossSectionCalculator->XSection(G4StokesVector::ZERO,
00214                                                      G4StokesVector::ZERO);
00215         z=xs*sqr(x)*4.;
00216         if (grej < z) {
00217           G4cout<<"WARNING : error in Moller rejection routine! \n"
00218                 <<" z = "<<z<<" grej="<<grej<<"\n";
00219         }
00220       } else {
00221         G4cout<<"No calculator in Moller scattering"<<G4endl;
00222       }
00223        } while(grej * G4UniformRand() > z);
00224     //Bhabha (e+e-) scattering
00225   } else {
00226     // *** dice according to polarized cross section
00227     y     = xmax*xmax;
00228     grej = 0.;
00229     grej += y*y*gmo3*(1. + (polL + polT)*(gam + 3.)/gmo);
00230     grej += -2.*xmin*xmin*xmin*gam*gmo2*(1. - (polL + polT)*(gam + 3.)/gmo);
00231     grej += y*y*gmo*(3.*gamma2 + 6.*gam + 4.)*(1. + (polL*(3.*gam + 1.)*(gamma2 + gam + 1.) + polT*((gam + 2.)*gamma2 + 1.))/(gmo*(3.*gam*(gam + 2.) + 4.)));
00232     grej /= gpo3;
00233     grej += -xmin*(2.*gamma2 + 4.*gam + 1.)*(1. - gam*(polL*(2.*gam + 1.) + polT)/(2.*gam*(gam + 2.) + 1.))/gpo2;
00234     grej += gamma2/(gamma2 - 1.);
00235     G4double prefB = classic_electr_radius*classic_electr_radius/(gam - 1.0);
00236     grej *= prefB;
00237 
00238     do {
00239       q  = G4UniformRand();
00240       x  = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
00241       if (crossSectionCalculator) {
00242         crossSectionCalculator->Initialize(x,gam,phi,theBeamPolarization,
00243                                            theTargetPolarization,1);
00244         xs=crossSectionCalculator->XSection(G4StokesVector::ZERO,
00245                                                      G4StokesVector::ZERO);
00246         z=xs*sqr(x)*4.;
00247       } else {
00248         G4cout<<"No calculator in Bhabha scattering"<<G4endl;
00249       }
00250 
00251       if(z > grej) {
00252         G4cout<<"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&"<<G4endl;
00253         G4cout << "G4PolarizedMollerBhabhaModel::SampleSecondaries Warning! "<<G4endl
00254                << "Majorant " << grej << " < "
00255                << z << " for x= " << x<<G4endl
00256                << " e+e- (Bhabha) scattering"<<" at KinEnergy "<<kineticEnergy<<G4endl;
00257         G4cout<<"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&"<<G4endl;
00258       }
00259     } while(grej * G4UniformRand() > z);
00260   }
00261   //
00262   //
00263   // polar asymmetries (due to transverse polarizations)
00264   //
00265   //
00266   if (crossSectionCalculator) {
00267    // grej*=1./(sqr(x)*sqr(gamma2-1))*sqr(gam*(1+gam));
00268     grej=xs*2.;
00269     do {
00270       phi = twopi * G4UniformRand() ;
00271       crossSectionCalculator->Initialize(x,gam,phi,theBeamPolarization,
00272                                                    theTargetPolarization,1);
00273       xs=crossSectionCalculator->XSection(G4StokesVector::ZERO,
00274                                           G4StokesVector::ZERO);
00275       if(xs > grej) {
00276         if (isElectron){
00277           G4cout << "G4PolarizedMollerBhabhaModel::SampleSecondaries Warning! "<<G4endl
00278                  << "Majorant " << grej << " < "
00279                  << xs << " for phi= " << phi<<G4endl
00280                  << " e-e- (Moller) scattering"<< G4endl
00281                  <<"PHI DICING"<<G4endl;
00282         } else {
00283           G4cout << "G4PolarizedMollerBhabhaModel::SampleSecondaries Warning! "<<G4endl
00284                  << "Majorant " << grej << " < "
00285                  << xs << " for phi= " << phi<<G4endl
00286                  << " e+e- (Bhabha) scattering"<< G4endl
00287                  <<"PHI DICING"<<G4endl;
00288         }
00289       }
00290     } while(grej * G4UniformRand() > xs);
00291   }
00292 
00293   // fix kinematics of delta electron
00294   G4double deltaKinEnergy = x * kineticEnergy;
00295   G4double deltaMomentum =
00296            std::sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
00297   G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
00298                                    (deltaMomentum * totalMomentum);
00299   G4double sint = 1.0 - cost*cost;
00300   if(sint > 0.0) sint = std::sqrt(sint);
00301 
00302 
00303   G4ThreeVector deltaDirection(-sint*std::cos(phi),-sint*std::sin(phi), cost) ;
00304   deltaDirection.rotateUz(direction);
00305 
00306   // primary change
00307   kineticEnergy -= deltaKinEnergy;
00308   fParticleChange->SetProposedKineticEnergy(kineticEnergy);
00309 
00310   if(kineticEnergy > DBL_MIN) {
00311     G4ThreeVector dir = totalMomentum*direction - deltaMomentum*deltaDirection;
00312     direction = dir.unit();
00313     fParticleChange->SetProposedMomentumDirection(direction);
00314   }
00315 
00316   // create G4DynamicParticle object for delta ray
00317   G4DynamicParticle* delta = new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
00318   vdp->push_back(delta);
00319 
00320   // get interaction frame
00321   G4ThreeVector  nInteractionFrame = 
00322     G4PolarizationHelper::GetFrame(direction,deltaDirection);
00323 
00324   if (crossSectionCalculator) {
00325     // calculate mean final state polarizations
00326 
00327     theBeamPolarization.InvRotateAz(nInteractionFrame,direction);
00328     theTargetPolarization.InvRotateAz(nInteractionFrame,direction);
00329     crossSectionCalculator->Initialize(x,gam,phi,theBeamPolarization,
00330                                        theTargetPolarization,2);
00331 
00332     // electron/positron
00333     fPositronPolarization=crossSectionCalculator->GetPol2();
00334     fPositronPolarization.RotateAz(nInteractionFrame,direction);
00335 
00336     fParticleChange->ProposePolarization(fPositronPolarization);
00337 
00338     // electron
00339     fElectronPolarization=crossSectionCalculator->GetPol3();
00340     fElectronPolarization.RotateAz(nInteractionFrame,deltaDirection);
00341     delta->SetPolarization(fElectronPolarization.x(),
00342                            fElectronPolarization.y(),
00343                            fElectronPolarization.z());
00344   }
00345   else {
00346     fPositronPolarization=G4ThreeVector();
00347     fElectronPolarization=G4ThreeVector();
00348   }
00349 }

void G4PolarizedMollerBhabhaModel::SetBeamPolarization ( const G4ThreeVector &  pBeam  )  [inline]

Definition at line 88 of file G4PolarizedMollerBhabhaModel.hh.

Referenced by G4ePolarizedIonisation::ComputeAsymmetry().

00089   {
00090     theBeamPolarization = pBeam;
00091   }

void G4PolarizedMollerBhabhaModel::SetTargetPolarization ( const G4ThreeVector &  pTarget  )  [inline]

Definition at line 84 of file G4PolarizedMollerBhabhaModel.hh.

Referenced by G4ePolarizedIonisation::ComputeAsymmetry().

00085   {
00086     theTargetPolarization = pTarget;
00087   }

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