Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | List of all members
G4PolarizedAnnihilationModel Class Reference

#include <G4PolarizedAnnihilationModel.hh>

Inheritance diagram for G4PolarizedAnnihilationModel:
Inheritance graph
[legend]
Collaboration diagram for G4PolarizedAnnihilationModel:
Collaboration graph
[legend]

Public Member Functions

 G4PolarizedAnnihilationModel (const G4ParticleDefinition *p=nullptr, const G4String &nam="Polarized-Annihilation")
 
virtual ~G4PolarizedAnnihilationModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &) override
 
virtual G4double ComputeCrossSectionPerElectron (const G4ParticleDefinition *, G4double kinEnergy, G4double cut, G4double emax) override
 
void ComputeAsymmetriesPerElectron (G4double gammaEnergy, G4double &valueX, G4double &valueA, G4double &valueT)
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
 
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) final
 
virtual G4double CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy) final
 
- 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 level, const G4ParticleDefinition *, G4double kineticEnergy)
 
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 * > *)
 
virtual 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=nullptr)
 
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)
 

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::G4PolarizedAnnihilationModel ( const G4ParticleDefinition *  p = nullptr,
const G4String &  nam = "Polarized-Annihilation" 
)
explicit

Definition at line 64 of file G4PolarizedAnnihilationModel.cc.

66  : G4eeToTwoGammaModel(p,nam),
67  crossSectionCalculator(nullptr),
68  verboseLevel(0),
69  gParticleChange(nullptr),
70  gIsInitialised(false)
71 {
72  crossSectionCalculator=new G4PolarizedAnnihilationCrossSection();
73 }
G4eeToTwoGammaModel::G4eeToTwoGammaModel
G4eeToTwoGammaModel(const G4ParticleDefinition *p=0, const G4String &nam="eplus2gg")
Definition: G4eeToTwoGammaModel.cc:87
G4PolarizedAnnihilationModel::~G4PolarizedAnnihilationModel ( )
virtual

Definition at line 75 of file G4PolarizedAnnihilationModel.cc.

76 {
77  if (crossSectionCalculator) delete crossSectionCalculator;
78 }

Member Function Documentation

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

Definition at line 113 of file G4PolarizedAnnihilationModel.cc.

117 {
118  // *** calculate asymmetries
119  G4double gam = 1. + ene/electron_mass_c2;
120  G4double xs0=crossSectionCalculator->TotalXSection(0.,1.,gam,
121  G4StokesVector::ZERO,
122  G4StokesVector::ZERO);
123  G4double xsA=crossSectionCalculator->TotalXSection(0.,1.,gam,
124  G4StokesVector::P3,
125  G4StokesVector::P3);
126  G4double xsT1=crossSectionCalculator->TotalXSection(0.,1.,gam,
127  G4StokesVector::P1,
128  G4StokesVector::P1);
129  G4double xsT2=crossSectionCalculator->TotalXSection(0.,1.,gam,
130  G4StokesVector::P2,
131  G4StokesVector::P2);
132  G4double xsT=0.5*(xsT1+xsT2);
133 
134  valueX=xs0;
135  valueA=xsA/xs0-1.;
136  valueT=xsT/xs0-1.;
137  // G4cout<<valueX<<"\t"<<valueA<<"\t"<<valueT<<" energy = "<<gam<<G4endl;
138  if ( (valueA < -1) || (1 < valueA)) {
139  G4cout<< " ERROR PolarizedAnnihilationPS::ComputeAsymmetries \n";
140  G4cout<< " something wrong in total cross section calculation (valueA)\n";
141  G4cout<<"*********** LONG "<<valueX<<"\t"<<valueA<<"\t"<<valueT<<" energy = "<<gam<<G4endl;
142  }
143  if ( (valueT < -1) || (1 < valueT)) {
144  G4cout<< " ERROR PolarizedAnnihilationPS::ComputeAsymmetries \n";
145  G4cout<< " something wrong in total cross section calculation (valueT)\n";
146  G4cout<<"****** TRAN "<<valueX<<"\t"<<valueA<<"\t"<<valueT<<" energy = "<<gam<<G4endl;
147  }
148 }
CLHEP::electron_mass_c2
static constexpr double electron_mass_c2
Definition: PhysicalConstants.h:79
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
G4StokesVector::P1
static const G4StokesVector P1
Definition: G4StokesVector.hh:63
G4PolarizedAnnihilationCrossSection::TotalXSection
virtual G4double TotalXSection(G4double xmin, G4double xmax, G4double y, const G4StokesVector &pol0, const G4StokesVector &pol1) override
Definition: G4PolarizedAnnihilationCrossSection.cc:238
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:

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

Reimplemented from G4eeToTwoGammaModel.

Definition at line 92 of file G4PolarizedAnnihilationModel.cc.

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

Here is the call graph for this function:

const G4ThreeVector & G4PolarizedAnnihilationModel::GetBeamPolarization ( ) const
inline

Definition at line 133 of file G4PolarizedAnnihilationModel.hh.

134 {
135  return theBeamPolarization;
136 }
const G4ThreeVector & G4PolarizedAnnihilationModel::GetFinalGamma1Polarization ( ) const
inline

Definition at line 137 of file G4PolarizedAnnihilationModel.hh.

138 {
139  return finalGamma1Polarization;
140 }
const G4ThreeVector & G4PolarizedAnnihilationModel::GetFinalGamma2Polarization ( ) const
inline

Definition at line 141 of file G4PolarizedAnnihilationModel.hh.

142 {
143  return finalGamma2Polarization;
144 }
const G4ThreeVector & G4PolarizedAnnihilationModel::GetTargetPolarization ( ) const
inline

Definition at line 129 of file G4PolarizedAnnihilationModel.hh.

130 {
131  return theTargetPolarization;
132 }
void G4PolarizedAnnihilationModel::Initialise ( const G4ParticleDefinition *  ,
const G4DataVector &   
)
overridevirtual

Reimplemented from G4eeToTwoGammaModel.

Definition at line 83 of file G4PolarizedAnnihilationModel.cc.

85 {
86  // G4eeToTwoGammaModel::Initialise(part,dv);
87  if(gIsInitialised) return;
88  gParticleChange = GetParticleChangeForGamma();
89  gIsInitialised = true;
90 }
G4VEmModel::GetParticleChangeForGamma
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:132

Here is the call graph for this function:

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

Reimplemented from G4eeToTwoGammaModel.

Definition at line 151 of file G4PolarizedAnnihilationModel.cc.

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

Here is the call graph for this function:

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

Definition at line 125 of file G4PolarizedAnnihilationModel.hh.

126 {
127  theBeamPolarization = pBeam;
128 }
void G4PolarizedAnnihilationModel::SetTargetPolarization ( const G4ThreeVector &  pTarget)
inline

Definition at line 121 of file G4PolarizedAnnihilationModel.hh.

122 {
123  theTargetPolarization = pTarget;
124 }
The documentation for this class was generated from the following files: