Geant4  10.02.p03
Public Member Functions | Private Member Functions | Private Attributes | List of all members
G4PolarizedAnnihilationModel Class Reference

#include <G4PolarizedAnnihilationModel.hh>

Inheritance diagram for G4PolarizedAnnihilationModel:
Collaboration diagram for G4PolarizedAnnihilationModel:

Public Member Functions

 G4PolarizedAnnihilationModel (const G4ParticleDefinition *p=0, const G4String &nam="Polarized-Annihilation")
 
virtual ~G4PolarizedAnnihilationModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)
 
virtual G4double ComputeCrossSectionPerElectron (const G4ParticleDefinition *, G4double kinEnergy, G4double cut, G4double emax)
 
void ComputeAsymmetriesPerElectron (G4double gammaEnergy, G4double &valueX, G4double &valueA, G4double &valueT)
 
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 & GetFinalGamma1Polarization () const
 
const G4ThreeVector & GetFinalGamma2Polarization () const
 
- Public Member Functions inherited from G4eeToTwoGammaModel
 G4eeToTwoGammaModel (const G4ParticleDefinition *p=0, const G4String &nam="eplus2gg")
 
virtual ~G4eeToTwoGammaModel ()
 
virtual G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0., G4double maxEnergy=DBL_MAX)
 
virtual G4double CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
 
- Public Member Functions inherited from G4VEmModel
 G4VEmModel (const G4String &nam)
 
virtual ~G4VEmModel ()
 
virtual void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel)
 
virtual void InitialiseForMaterial (const G4ParticleDefinition *, const G4Material *)
 
virtual void InitialiseForElement (const G4ParticleDefinition *, G4int Z)
 
virtual G4double ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
virtual G4double GetPartialCrossSection (const G4Material *, G4int, const G4ParticleDefinition *, G4double)
 
virtual G4double ComputeCrossSectionPerShell (const G4ParticleDefinition *, G4int Z, G4int shellIdx, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double ChargeSquareRatio (const G4Track &)
 
virtual G4double GetChargeSquareRatio (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual G4double GetParticleCharge (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void StartTracking (G4Track *)
 
virtual void CorrectionsAlongStep (const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length)
 
virtual G4double Value (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
 
virtual G4double MinPrimaryEnergy (const G4Material *, const G4ParticleDefinition *, G4double cut=0.0)
 
virtual G4double MinEnergyCut (const G4ParticleDefinition *, const G4MaterialCutsCouple *)
 
virtual void SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void DefineForRegion (const G4Region *)
 
virtual void ModelDescription (std::ostream &outFile) const
 
void InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
 
std::vector< G4EmElementSelector * > * GetElementSelectors ()
 
void SetElementSelectors (std::vector< G4EmElementSelector *> *)
 
G4double ComputeDEDX (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
G4double CrossSection (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeMeanFreePath (const G4ParticleDefinition *, G4double kineticEnergy, const G4Material *, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, const G4Element *, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4int SelectIsotopeNumber (const G4Element *)
 
const G4Element * SelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4Element * SelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4int SelectRandomAtomNumber (const G4Material *)
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=0)
 
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
 
G4ElementData * GetElementData ()
 
G4PhysicsTable * GetCrossSectionTable ()
 
G4VEmFluctuationModel * GetModelOfFluctuations ()
 
G4VEmAngularDistribution * GetAngularDistribution ()
 
void SetAngularDistribution (G4VEmAngularDistribution *)
 
G4double HighEnergyLimit () const
 
G4double LowEnergyLimit () const
 
G4double HighEnergyActivationLimit () const
 
G4double LowEnergyActivationLimit () const
 
G4double PolarAngleLimit () const
 
G4double SecondaryThreshold () const
 
G4bool LPMFlag () const
 
G4bool DeexcitationFlag () const
 
G4bool ForceBuildTableFlag () const
 
G4bool UseAngularGeneratorFlag () const
 
void SetAngularGeneratorFlag (G4bool)
 
void SetHighEnergyLimit (G4double)
 
void SetLowEnergyLimit (G4double)
 
void SetActivationHighEnergyLimit (G4double)
 
void SetActivationLowEnergyLimit (G4double)
 
G4bool IsActive (G4double kinEnergy)
 
void SetPolarAngleLimit (G4double)
 
void SetSecondaryThreshold (G4double)
 
void SetLPMFlag (G4bool val)
 
void SetDeexcitationFlag (G4bool val)
 
void SetForceBuildTable (G4bool val)
 
void SetMasterThread (G4bool val)
 
G4bool IsMaster () const
 
G4double MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
 
const G4String & GetName () const
 
void SetCurrentCouple (const G4MaterialCutsCouple *)
 
const G4Element * GetCurrentElement () const
 
const G4Isotope * GetCurrentIsotope () const
 
G4bool IsLocked () const
 
void SetLocked (G4bool)
 

Private Member Functions

G4PolarizedAnnihilationModel & operator= (const G4PolarizedAnnihilationModel &right)
 
 G4PolarizedAnnihilationModel (const G4PolarizedAnnihilationModel &)
 

Private Attributes

G4PolarizedAnnihilationCrossSection * crossSectionCalculator
 
G4StokesVector theBeamPolarization
 
G4StokesVector theTargetPolarization
 
G4StokesVector finalGamma1Polarization
 
G4StokesVector finalGamma2Polarization
 
G4int verboseLevel
 
G4ParticleChangeForGamma * gParticleChange
 
G4bool gIsInitialised
 

Additional Inherited Members

- Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLoss * GetParticleChangeForLoss ()
 
G4ParticleChangeForGamma * GetParticleChangeForGamma ()
 
virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
 
const G4MaterialCutsCouple * CurrentCouple () const
 
void SetCurrentElement (const G4Element *)
 
- Protected Attributes inherited from G4VEmModel
G4ElementData * fElementData
 
G4VParticleChange * pParticleChange
 
G4PhysicsTable * xSectionTable
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 
size_t idxTable
 
- Static Protected Attributes inherited from G4VEmModel
static const G4double inveplus = 1.0/CLHEP::eplus
 

Detailed Description

Definition at line 63 of file G4PolarizedAnnihilationModel.hh.

Constructor & Destructor Documentation

◆ G4PolarizedAnnihilationModel() [1/2]

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

Definition at line 64 of file G4PolarizedAnnihilationModel.cc.

66  : G4eeToTwoGammaModel(p,nam),crossSectionCalculator(0),verboseLevel(0),gParticleChange(0),
67  gIsInitialised(false)
68 {
69  crossSectionCalculator=new G4PolarizedAnnihilationCrossSection();
70 }
G4PolarizedAnnihilationModel::gParticleChange
G4ParticleChangeForGamma * gParticleChange
Definition: G4PolarizedAnnihilationModel.hh:114
G4PolarizedAnnihilationModel::crossSectionCalculator
G4PolarizedAnnihilationCrossSection * crossSectionCalculator
Definition: G4PolarizedAnnihilationModel.hh:104
G4eeToTwoGammaModel::G4eeToTwoGammaModel
G4eeToTwoGammaModel(const G4ParticleDefinition *p=0, const G4String &nam="eplus2gg")
Definition: G4eeToTwoGammaModel.cc:87

◆ ~G4PolarizedAnnihilationModel()

G4PolarizedAnnihilationModel::~G4PolarizedAnnihilationModel ( )
virtual

Definition at line 72 of file G4PolarizedAnnihilationModel.cc.

73 {
74  if (crossSectionCalculator) delete crossSectionCalculator;
75 }
G4PolarizedAnnihilationModel::crossSectionCalculator
G4PolarizedAnnihilationCrossSection * crossSectionCalculator
Definition: G4PolarizedAnnihilationModel.hh:104

◆ G4PolarizedAnnihilationModel() [2/2]

G4PolarizedAnnihilationModel::G4PolarizedAnnihilationModel ( const G4PolarizedAnnihilationModel &  )
private

Member Function Documentation

◆ ComputeAsymmetriesPerElectron()

void G4PolarizedAnnihilationModel::ComputeAsymmetriesPerElectron ( G4double  gammaEnergy,
G4double &  valueX,
G4double &  valueA,
G4double &  valueT 
)

Definition at line 110 of file G4PolarizedAnnihilationModel.cc.

114 {
115  // *** calculate asymmetries
116  G4double gam = 1. + ene/electron_mass_c2;
117  G4double xs0=crossSectionCalculator->TotalXSection(0.,1.,gam,
118  G4StokesVector::ZERO,
119  G4StokesVector::ZERO);
120  G4double xsA=crossSectionCalculator->TotalXSection(0.,1.,gam,
121  G4StokesVector::P3,
122  G4StokesVector::P3);
123  G4double xsT1=crossSectionCalculator->TotalXSection(0.,1.,gam,
124  G4StokesVector::P1,
125  G4StokesVector::P1);
126  G4double xsT2=crossSectionCalculator->TotalXSection(0.,1.,gam,
127  G4StokesVector::P2,
128  G4StokesVector::P2);
129  G4double xsT=0.5*(xsT1+xsT2);
130 
131  valueX=xs0;
132  valueA=xsA/xs0-1.;
133  valueT=xsT/xs0-1.;
134  // G4cout<<valueX<<"\t"<<valueA<<"\t"<<valueT<<" energy = "<<gam<<G4endl;
135  if ( (valueA < -1) || (1 < valueA)) {
136  G4cout<< " ERROR PolarizedAnnihilationPS::ComputeAsymmetries \n";
137  G4cout<< " something wrong in total cross section calculation (valueA)\n";
138  G4cout<<"*********** LONG "<<valueX<<"\t"<<valueA<<"\t"<<valueT<<" energy = "<<gam<<G4endl;
139  }
140  if ( (valueT < -1) || (1 < valueT)) {
141  G4cout<< " ERROR PolarizedAnnihilationPS::ComputeAsymmetries \n";
142  G4cout<< " something wrong in total cross section calculation (valueT)\n";
143  G4cout<<"****** TRAN "<<valueX<<"\t"<<valueA<<"\t"<<valueT<<" energy = "<<gam<<G4endl;
144  }
145 }
G4PolarizedAnnihilationCrossSection::TotalXSection
virtual G4double TotalXSection(G4double xmin, G4double xmax, G4double y, const G4StokesVector &pol0, const G4StokesVector &pol1)
Definition: G4PolarizedAnnihilationCrossSection.cc:238
G4StokesVector::P3
static const G4StokesVector P3
Definition: G4StokesVector.hh:65
G4cout
G4GLOB_DLL std::ostream G4cout
G4StokesVector::P2
static const G4StokesVector P2
Definition: G4StokesVector.hh:64
G4PolarizedAnnihilationModel::crossSectionCalculator
G4PolarizedAnnihilationCrossSection * crossSectionCalculator
Definition: G4PolarizedAnnihilationModel.hh:104
G4StokesVector::P1
static const G4StokesVector P1
Definition: G4StokesVector.hh:63
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4StokesVector::ZERO
static const G4StokesVector ZERO
Definition: G4StokesVector.hh:62
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ComputeCrossSectionPerElectron()

G4double G4PolarizedAnnihilationModel::ComputeCrossSectionPerElectron ( const G4ParticleDefinition *  pd,
G4double  kinEnergy,
G4double  cut,
G4double  emax 
)
virtual

Reimplemented from G4eeToTwoGammaModel.

Definition at line 89 of file G4PolarizedAnnihilationModel.cc.

94 {
95  G4double xs = G4eeToTwoGammaModel::ComputeCrossSectionPerElectron(pd,kinEnergy,
96  cut,emax);
97 
98  G4double polzz = theBeamPolarization.z()*theTargetPolarization.z();
99  G4double poltt = theBeamPolarization.x()*theTargetPolarization.x()
100  + theBeamPolarization.y()*theTargetPolarization.y();
101  if (polzz!=0 || poltt!=0) {
102  G4double xval,lasym,tasym;
103  ComputeAsymmetriesPerElectron(kinEnergy,xval,lasym,tasym);
104  xs*=(1.+polzz*lasym+poltt*tasym);
105  }
106 
107  return xs;
108 }
G4eeToTwoGammaModel::ComputeCrossSectionPerElectron
virtual G4double ComputeCrossSectionPerElectron(const G4ParticleDefinition *, G4double kinEnergy, G4double cutEnergy=0., G4double maxEnergy=DBL_MAX)
Definition: G4eeToTwoGammaModel.cc:114
CLHEP::Hep3Vector::x
double x() const
G4PolarizedAnnihilationModel::ComputeAsymmetriesPerElectron
void ComputeAsymmetriesPerElectron(G4double gammaEnergy, G4double &valueX, G4double &valueA, G4double &valueT)
Definition: G4PolarizedAnnihilationModel.cc:110
emax
static const G4double emax
Definition: G4ChatterjeeCrossSection.hh:39
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::z
double z() const
G4double
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:

◆ GetBeamPolarization()

const G4ThreeVector & G4PolarizedAnnihilationModel::GetBeamPolarization ( ) const
inline

Definition at line 132 of file G4PolarizedAnnihilationModel.hh.

133 {
134  return theBeamPolarization;
135 }

◆ GetFinalGamma1Polarization()

const G4ThreeVector & G4PolarizedAnnihilationModel::GetFinalGamma1Polarization ( ) const
inline

Definition at line 136 of file G4PolarizedAnnihilationModel.hh.

137 {
138  return finalGamma1Polarization;
139 }

◆ GetFinalGamma2Polarization()

const G4ThreeVector & G4PolarizedAnnihilationModel::GetFinalGamma2Polarization ( ) const
inline

Definition at line 140 of file G4PolarizedAnnihilationModel.hh.

141 {
142  return finalGamma2Polarization;
143 }

◆ GetTargetPolarization()

const G4ThreeVector & G4PolarizedAnnihilationModel::GetTargetPolarization ( ) const
inline

Definition at line 128 of file G4PolarizedAnnihilationModel.hh.

129 {
130  return theTargetPolarization;
131 }

◆ Initialise()

void G4PolarizedAnnihilationModel::Initialise ( const G4ParticleDefinition *  ,
const G4DataVector &   
)
virtual

Reimplemented from G4eeToTwoGammaModel.

Definition at line 80 of file G4PolarizedAnnihilationModel.cc.

82 {
83  // G4eeToTwoGammaModel::Initialise(part,dv);
84  if(gIsInitialised) return;
85  gParticleChange = GetParticleChangeForGamma();
86  gIsInitialised = true;
87 }
G4PolarizedAnnihilationModel::gParticleChange
G4ParticleChangeForGamma * gParticleChange
Definition: G4PolarizedAnnihilationModel.hh:114
G4VEmModel::GetParticleChangeForGamma
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:134
Here is the call graph for this function:

◆ operator=()

G4PolarizedAnnihilationModel& G4PolarizedAnnihilationModel::operator= ( const G4PolarizedAnnihilationModel &  right)
private

◆ SampleSecondaries()

void G4PolarizedAnnihilationModel::SampleSecondaries ( std::vector< G4DynamicParticle *> *  fvect,
const G4MaterialCutsCouple *  ,
const G4DynamicParticle *  dp,
G4double  tmin,
G4double  maxEnergy 
)
virtual

Reimplemented from G4eeToTwoGammaModel.

Definition at line 148 of file G4PolarizedAnnihilationModel.cc.

153 {
154 // G4ParticleChangeForGamma* gParticleChange
155 // = dynamic_cast<G4ParticleChangeForGamma*>(pParticleChange);
156  const G4Track * aTrack = gParticleChange->GetCurrentTrack();
157 
158  // kill primary
159  gParticleChange->SetProposedKineticEnergy(0.);
160  gParticleChange->ProposeTrackStatus(fStopAndKill);
161 
162  // V.Ivanchenko add protection against zero kin energy
163  G4double PositKinEnergy = dp->GetKineticEnergy();
164 
165  if(PositKinEnergy < DBL_MIN) {
166 
167  G4double cosTeta = 2.*G4UniformRand()-1.;
168  G4double sinTeta = std::sqrt((1.0 - cosTeta)*(1.0 + cosTeta));
169  G4double phi = twopi * G4UniformRand();
170  G4ThreeVector dir(sinTeta*std::cos(phi), sinTeta*std::sin(phi), cosTeta);
171  fvect->push_back( new G4DynamicParticle(G4Gamma::Gamma(), dir, electron_mass_c2));
172  fvect->push_back( new G4DynamicParticle(G4Gamma::Gamma(),-dir, electron_mass_c2));
173  return;
174  }
175 
176  // *** obtain and save target and beam polarization ***
177  G4PolarizationManager * polarizationManager = G4PolarizationManager::GetInstance();
178 
179  // obtain polarization of the beam
180  theBeamPolarization = aTrack->GetPolarization();
181 
182  // obtain polarization of the media
183  G4VPhysicalVolume* aPVolume = aTrack->GetVolume();
184  G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
185  const G4bool targetIsPolarized = polarizationManager->IsPolarized(aLVolume);
186  theTargetPolarization = polarizationManager->GetVolumePolarization(aLVolume);
187 
188  // transfer target electron polarization in frame of positron
189  if (targetIsPolarized)
190  theTargetPolarization.rotateUz(dp->GetMomentumDirection());
191 
192  G4ParticleMomentum PositDirection = dp->GetMomentumDirection();
193 
194  // polar asymmetry:
195  G4double polarization = theBeamPolarization.p3()*theTargetPolarization.p3();
196 
197  G4double gamam1 = PositKinEnergy/electron_mass_c2;
198  G4double gama = gamam1+1. , gamap1 = gamam1+2.;
199  G4double sqgrate = std::sqrt(gamam1/gamap1)/2. , sqg2m1 = std::sqrt(gamam1*gamap1);
200 
201  // limits of the energy sampling
202  G4double epsilmin = 0.5 - sqgrate , epsilmax = 0.5 + sqgrate;
203  G4double epsilqot = epsilmax/epsilmin;
204 
205  //
206  // sample the energy rate of the created gammas
207  // note: for polarized partices, the actual dicing strategy
208  // will depend on the energy, and the degree of polarization !!
209  //
210  G4double epsil;
211  G4double gmax=1. + std::fabs(polarization); // crude estimate
212 
213  //G4bool check_range=true;
214 
215  crossSectionCalculator->Initialize(epsilmin, gama, 0., theBeamPolarization, theTargetPolarization);
216  if (crossSectionCalculator->DiceEpsilon()<0) {
217  G4cout<<"ERROR in PolarizedAnnihilationPS::PostStepDoIt\n"
218  <<"epsilmin DiceRoutine not appropriate ! "<<crossSectionCalculator->DiceEpsilon()<<G4endl;
219  //check_range=false;
220  }
221 
222  crossSectionCalculator->Initialize(epsilmax, gama, 0., theBeamPolarization, theTargetPolarization);
223  if (crossSectionCalculator->DiceEpsilon()<0) {
224  G4cout<<"ERROR in PolarizedAnnihilationPS::PostStepDoIt\n"
225  <<"epsilmax DiceRoutine not appropriate ! "<<crossSectionCalculator->DiceEpsilon()<<G4endl;
226  //check_range=false;
227  }
228 
229  G4int ncount=0;
230  G4double trejectmax=0.;
231  G4double treject;
232 
233 
234  do {
235  //
236  epsil = epsilmin*std::pow(epsilqot,G4UniformRand());
237 
238  crossSectionCalculator->Initialize(epsil, gama, 0., theBeamPolarization, theTargetPolarization,1);
239 
240  treject = crossSectionCalculator->DiceEpsilon();
241  treject*=epsil;
242 
243  if (treject>gmax || treject<0.)
244  G4cout<<"ERROR in PolarizedAnnihilationPS::PostStepDoIt\n"
245  <<" eps ("<<epsil<<") rejection does not work properly: "<<treject<<G4endl;
246  ++ncount;
247  if (treject>trejectmax) trejectmax=treject;
248  if (ncount>1000) {
249  G4cout<<"WARNING in PolarizedAnnihilationPS::PostStepDoIt\n"
250  <<"eps dicing very inefficient ="<<trejectmax/gmax
251  <<", "<<treject/gmax<<". For secondary energy = "<<epsil<<" "<<ncount<<G4endl;
252  break;
253  }
254 
255  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
256  } while( treject < gmax*G4UniformRand() );
257 
258  //
259  // scattered Gamma angles. ( Z - axis along the parent positron)
260  //
261 
262  G4double cost = (epsil*gamap1-1.)/(epsil*sqg2m1);
263  G4double sint = std::sqrt((1.+cost)*(1.-cost));
264  G4double phi = 0.;
265  G4double beamTrans = std::sqrt(sqr(theBeamPolarization.p1()) + sqr(theBeamPolarization.p2()));
266  G4double targetTrans = std::sqrt(sqr(theTargetPolarization.p1()) + sqr(theTargetPolarization.p2()));
267 
268  // G4cout<<"phi dicing START"<<G4endl;
269  do{
270  phi = twopi * G4UniformRand();
271  crossSectionCalculator->Initialize(epsil, gama, 0., theBeamPolarization, theTargetPolarization,2);
272 
273  G4double gdiced =crossSectionCalculator->getVar(0);
274  gdiced += crossSectionCalculator->getVar(3)*theBeamPolarization.p3()*theTargetPolarization.p3();
275  gdiced += 1.*(std::fabs(crossSectionCalculator->getVar(1))
276  + std::fabs(crossSectionCalculator->getVar(2)))*beamTrans*targetTrans;
277  gdiced += 1.*std::fabs(crossSectionCalculator->getVar(4))
278  *(std::fabs(theBeamPolarization.p3())*targetTrans + std::fabs(theTargetPolarization.p3())*beamTrans);
279 
280  G4double gdist = crossSectionCalculator->getVar(0);
281  gdist += crossSectionCalculator->getVar(3)*theBeamPolarization.p3()*theTargetPolarization.p3();
282  gdist += crossSectionCalculator->getVar(1)*(std::cos(phi)*theBeamPolarization.p1()
283  + std::sin(phi)*theBeamPolarization.p2())
284  *(std::cos(phi)*theTargetPolarization.p1()
285  + std::sin(phi)*theTargetPolarization.p2());
286  gdist += crossSectionCalculator->getVar(2)*(std::cos(phi)*theBeamPolarization.p2()
287  - std::sin(phi)*theBeamPolarization.p1())
288  *(std::cos(phi)*theTargetPolarization.p2()
289  - std::sin(phi)*theTargetPolarization.p1());
290  gdist += crossSectionCalculator->getVar(4)
291  *(std::cos(phi)*theBeamPolarization.p3()*theTargetPolarization.p1()
292  + std::cos(phi)*theBeamPolarization.p1()*theTargetPolarization.p3()
293  + std::sin(phi)*theBeamPolarization.p3()*theTargetPolarization.p2()
294  + std::sin(phi)*theBeamPolarization.p2()*theTargetPolarization.p3());
295 
296  treject = gdist/gdiced;
297  //G4cout<<" treject = "<<treject<<" at phi = "<<phi<<G4endl;
298  if (treject>1.+1.e-10 || treject<0){
299  G4cout<<"!!!ERROR in PolarizedAnnihilationPS::PostStepDoIt\n"
300  <<" phi rejection does not work properly: "<<treject<<G4endl;
301  G4cout<<" gdiced = "<<gdiced<<G4endl;
302  G4cout<<" gdist = "<<gdist<<G4endl;
303  G4cout<<" epsil = "<<epsil<<G4endl;
304  }
305 
306  if (treject<1.e-3) {
307  G4cout<<"!!!ERROR in PolarizedAnnihilationPS::PostStepDoIt\n"
308  <<" phi rejection does not work properly: "<<treject<<"\n";
309  G4cout<<" gdiced="<<gdiced<<" gdist="<<gdist<<"\n";
310  G4cout<<" epsil = "<<epsil<<G4endl;
311  }
312 
313  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
314  } while( treject < G4UniformRand() );
315  // G4cout<<"phi dicing END"<<G4endl;
316 
317  G4double dirx = sint*std::cos(phi) , diry = sint*std::sin(phi) , dirz = cost;
318 
319  //
320  // kinematic of the created pair
321  //
322  G4double TotalAvailableEnergy = PositKinEnergy + 2*electron_mass_c2;
323  G4double Phot1Energy = epsil*TotalAvailableEnergy;
324  G4double Phot2Energy =(1.-epsil)*TotalAvailableEnergy;
325 
326  // *** prepare calculation of polarization transfer ***
327  G4ThreeVector Phot1Direction (dirx, diry, dirz);
328 
329  // get interaction frame
330  G4ThreeVector nInteractionFrame =
331  G4PolarizationHelper::GetFrame(PositDirection,Phot1Direction);
332 
333  // define proper in-plane and out-of-plane component of initial spins
334  theBeamPolarization.InvRotateAz(nInteractionFrame,PositDirection);
335  theTargetPolarization.InvRotateAz(nInteractionFrame,PositDirection);
336 
337  // calculate spin transfere matrix
338 
339  crossSectionCalculator->Initialize(epsil,gama,phi,theBeamPolarization,theTargetPolarization,2);
340 
341  // **********************************************************************
342 
343  Phot1Direction.rotateUz(PositDirection);
344  // create G4DynamicParticle object for the particle1
345  G4DynamicParticle* aParticle1= new G4DynamicParticle (G4Gamma::Gamma(),
346  Phot1Direction, Phot1Energy);
347  finalGamma1Polarization=crossSectionCalculator->GetPol2();
348  G4double n1=finalGamma1Polarization.mag2();
349  if (n1>1) {
350  G4cout<<"ERROR: PolarizedAnnihilation Polarization Vector at epsil = "
351  <<epsil<<" is too large!!! \n"
352  <<"annihi pol1= "<<finalGamma1Polarization<<", ("<<n1<<")\n";
353  finalGamma1Polarization*=1./std::sqrt(n1);
354  }
355 
356  // define polarization of first final state photon
357  finalGamma1Polarization.SetPhoton();
358  finalGamma1Polarization.RotateAz(nInteractionFrame,Phot1Direction);
359  aParticle1->SetPolarization(finalGamma1Polarization.p1(),
360  finalGamma1Polarization.p2(),
361  finalGamma1Polarization.p3());
362 
363  fvect->push_back(aParticle1);
364 
365 
366  // **********************************************************************
367 
368  G4double Eratio= Phot1Energy/Phot2Energy;
369  G4double PositP= std::sqrt(PositKinEnergy*(PositKinEnergy+2.*electron_mass_c2));
370  G4ThreeVector Phot2Direction (-dirx*Eratio, -diry*Eratio,
371  (PositP-dirz*Phot1Energy)/Phot2Energy);
372  Phot2Direction.rotateUz(PositDirection);
373  // create G4DynamicParticle object for the particle2
374  G4DynamicParticle* aParticle2= new G4DynamicParticle (G4Gamma::Gamma(),
375  Phot2Direction, Phot2Energy);
376 
377  // define polarization of second final state photon
378  finalGamma2Polarization=crossSectionCalculator->GetPol3();
379  G4double n2=finalGamma2Polarization.mag2();
380  if (n2>1) {
381  G4cout<<"ERROR: PolarizedAnnihilation Polarization Vector at epsil = "<<epsil<<" is too large!!! \n";
382  G4cout<<"annihi pol2= "<<finalGamma2Polarization<<", ("<<n2<<")\n";
383 
384  finalGamma2Polarization*=1./std::sqrt(n2);
385  }
386  finalGamma2Polarization.SetPhoton();
387  finalGamma2Polarization.RotateAz(nInteractionFrame,Phot2Direction);
388  aParticle2->SetPolarization(finalGamma2Polarization.p1(),
389  finalGamma2Polarization.p2(),
390  finalGamma2Polarization.p3());
391 
392  fvect->push_back(aParticle2);
393 }
G4PolarizationManager::GetVolumePolarization
const G4ThreeVector & GetVolumePolarization(G4LogicalVolume *lVol) const
Definition: G4PolarizationManager.hh:98
G4PolarizationManager::GetInstance
static G4PolarizationManager * GetInstance()
Definition: G4PolarizationManager.cc:52
G4StokesVector::p1
G4double p1() const
Definition: G4StokesVector.hh:74
dir
TDirectory * dir
Definition: human_phantom/macro.C:5
G4int
int G4int
Definition: G4Types.hh:78
CLHEP::Hep3Vector::mag2
double mag2() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4PolarizationManager::IsPolarized
bool IsPolarized(G4LogicalVolume *lVol) const
Definition: G4PolarizationManager.hh:108
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
G4PolarizedAnnihilationModel::gParticleChange
G4ParticleChangeForGamma * gParticleChange
Definition: G4PolarizedAnnihilationModel.hh:114
G4bool
bool G4bool
Definition: G4Types.hh:79
G4StokesVector::p3
G4double p3() const
Definition: G4StokesVector.hh:76
G4PolarizedAnnihilationModel::crossSectionCalculator
G4PolarizedAnnihilationCrossSection * crossSectionCalculator
Definition: G4PolarizedAnnihilationModel.hh:104
CLHEP::Hep3Vector::rotateUz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
twopi
static const double twopi
Definition: G4SIunits.hh:75
G4PolarizationHelper::GetFrame
static G4ThreeVector GetFrame(const G4ThreeVector &, const G4ThreeVector &)
Definition: G4PolarizationHelper.cc:48
G4StokesVector::p2
G4double p2() const
Definition: G4StokesVector.hh:75
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4DynamicParticle::SetPolarization
void SetPolarization(G4double polX, G4double polY, G4double polZ)
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4PolarizedAnnihilationCrossSection::Initialize
virtual void Initialize(G4double eps, G4double gamma, G4double phi, const G4StokesVector &p0, const G4StokesVector &p1, G4int flag=0)
Definition: G4PolarizedAnnihilationCrossSection.cc:90
G4StokesVector::InvRotateAz
void InvRotateAz(G4ThreeVector nInteractionFrame, G4ThreeVector particleDirection)
Definition: G4StokesVector.cc:108
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
DBL_MIN
#define DBL_MIN
Definition: templates.hh:75
G4endl
#define G4endl
Definition: G4ios.hh:61
sqr
T sqr(const T &x)
Definition: templates.hh:145
G4double
double G4double
Definition: G4Types.hh:76
G4VPhysicalVolume::GetLogicalVolume
G4LogicalVolume * GetLogicalVolume() const
G4StokesVector::RotateAz
void RotateAz(G4ThreeVector nInteractionFrame, G4ThreeVector particleDirection)
Definition: G4StokesVector.cc:71
Here is the call graph for this function:

◆ SetBeamPolarization()

void G4PolarizedAnnihilationModel::SetBeamPolarization ( const G4ThreeVector &  pBeam)
inline

Definition at line 124 of file G4PolarizedAnnihilationModel.hh.

125 {
126  theBeamPolarization = pBeam;
127 }
Here is the caller graph for this function:

◆ SetTargetPolarization()

void G4PolarizedAnnihilationModel::SetTargetPolarization ( const G4ThreeVector &  pTarget)
inline

Definition at line 120 of file G4PolarizedAnnihilationModel.hh.

121 {
122  theTargetPolarization = pTarget;
123 }
Here is the caller graph for this function:

Member Data Documentation

◆ crossSectionCalculator

G4PolarizedAnnihilationCrossSection* G4PolarizedAnnihilationModel::crossSectionCalculator
private

Definition at line 104 of file G4PolarizedAnnihilationModel.hh.

◆ finalGamma1Polarization

G4StokesVector G4PolarizedAnnihilationModel::finalGamma1Polarization
private

Definition at line 109 of file G4PolarizedAnnihilationModel.hh.

◆ finalGamma2Polarization

G4StokesVector G4PolarizedAnnihilationModel::finalGamma2Polarization
private

Definition at line 110 of file G4PolarizedAnnihilationModel.hh.

◆ gIsInitialised

G4bool G4PolarizedAnnihilationModel::gIsInitialised
private

Definition at line 115 of file G4PolarizedAnnihilationModel.hh.

◆ gParticleChange

G4ParticleChangeForGamma* G4PolarizedAnnihilationModel::gParticleChange
private

Definition at line 114 of file G4PolarizedAnnihilationModel.hh.

◆ theBeamPolarization

G4StokesVector G4PolarizedAnnihilationModel::theBeamPolarization
private

Definition at line 106 of file G4PolarizedAnnihilationModel.hh.

◆ theTargetPolarization

G4StokesVector G4PolarizedAnnihilationModel::theTargetPolarization
private

Definition at line 107 of file G4PolarizedAnnihilationModel.hh.

◆ verboseLevel

G4int G4PolarizedAnnihilationModel::verboseLevel
private

Definition at line 112 of file G4PolarizedAnnihilationModel.hh.

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