Geant4  10.01
G4eplusPolarizedAnnihilation.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4eplusPolarizedAnnihilation.cc 76472 2013-11-11 10:34:07Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4eplusPolarizedAnnihilation
34 //
35 // Author: A. Schaelicke on base of Vladimir Ivanchenko / Michel Maire code
36 //
37 // Creation date: 02.07.2006
38 //
39 // Modifications:
40 // 26-07-06 modified cross section (P. Starovoitov)
41 // 21-08-06 interface updated (A. Schaelicke)
42 // 11-06-07, add PostStepGetPhysicalInteractionLength (A.Schalicke)
43 // 02-10-07, enable AtRest (V.Ivanchenko)
44 //
45 //
46 // Class Description:
47 //
48 // Polarized process of e+ annihilation into 2 gammas
49 //
50 
51 //
52 // -------------------------------------------------------------------
53 //
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
56 
58 #include "G4PhysicalConstants.hh"
59 #include "G4SystemOfUnits.hh"
60 #include "G4MaterialCutsCouple.hh"
61 #include "G4Gamma.hh"
62 #include "G4PhysicsVector.hh"
63 #include "G4PhysicsLogVector.hh"
64 
65 
67 #include "G4PhysicsTableHelper.hh"
68 #include "G4ProductionCutsTable.hh"
69 #include "G4PolarizationManager.hh"
70 #include "G4PolarizationHelper.hh"
71 #include "G4StokesVector.hh"
72 
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
74 
76  : G4VEmProcess(name), isInitialised(false),
77  theAsymmetryTable(NULL),
78  theTransverseAsymmetryTable(NULL)
79 {
80  enableAtRestDoIt = true;
82  emModel = 0;
83 }
84 
85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
86 
88 {
89  delete theAsymmetryTable;
91 }
92 
93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
94 
96 {
97  if(!isInitialised) {
98  isInitialised = true;
99  // SetVerboseLevel(3);
100  SetBuildTableFlag(true);
101  SetStartFromNullFlag(false);
103  G4double emin = 0.1*keV;
104  G4double emax = 100.*TeV;
105  SetLambdaBinning(120);
106  SetMinKinEnergy(emin);
107  SetMaxKinEnergy(emax);
109  emModel->SetLowEnergyLimit(emin);
111  AddEmModel(1, emModel);
112  }
113 }
114 
115 
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
117 
118  // for polarization
119 
121  G4double previousStepSize,
123 {
124  G4double mfp = G4VEmProcess::GetMeanFreePath(track, previousStepSize, condition);
125 
126  if (theAsymmetryTable) {
127 
128  G4Material* aMaterial = track.GetMaterial();
129  G4VPhysicalVolume* aPVolume = track.GetVolume();
130  G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
131 
132  // G4Material* bMaterial = aLVolume->GetMaterial();
134 
135  const G4bool volumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
136  G4StokesVector electronPolarization = polarizationManger->GetVolumePolarization(aLVolume);
137 
138  if (!volumeIsPolarized || mfp == DBL_MAX) return mfp;
139 
140  // *** get asymmetry, if target is polarized ***
141  const G4DynamicParticle* aDynamicPositron = track.GetDynamicParticle();
142  const G4double positronEnergy = aDynamicPositron->GetKineticEnergy();
143  const G4StokesVector positronPolarization = track.GetPolarization();
144  const G4ParticleMomentum positronDirection0 = aDynamicPositron->GetMomentumDirection();
145 
146  if (verboseLevel>=2) {
147 
148  G4cout << " Mom " << positronDirection0 << G4endl;
149  G4cout << " Polarization " << positronPolarization << G4endl;
150  G4cout << " MaterialPol. " << electronPolarization << G4endl;
151  G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl;
152  G4cout << " Log. Volume " << aLVolume->GetName() << G4endl;
153  G4cout << " Material " << aMaterial << G4endl;
154  }
155 
156  G4bool isOutRange;
158  G4double lAsymmetry = (*theAsymmetryTable)(idx)->
159  GetValue(positronEnergy, isOutRange);
160  G4double tAsymmetry = (*theTransverseAsymmetryTable)(idx)->
161  GetValue(positronEnergy, isOutRange);
162 
163  G4double polZZ = positronPolarization.z()*
164  electronPolarization*positronDirection0;
165  G4double polXX = positronPolarization.x()*
166  electronPolarization*G4PolarizationHelper::GetParticleFrameX(positronDirection0);
167  G4double polYY = positronPolarization.y()*
168  electronPolarization*G4PolarizationHelper::GetParticleFrameY(positronDirection0);
169 
170  G4double impact = 1. + polZZ*lAsymmetry + (polXX + polYY)*tAsymmetry;
171 
172  mfp *= 1. / impact;
173 
174  if (verboseLevel>=2) {
175  G4cout << " MeanFreePath: " << mfp / mm << " mm " << G4endl;
176  G4cout << " Asymmetry: " << lAsymmetry << ", " << tAsymmetry << G4endl;
177  G4cout << " PolProduct: " << polXX << ", " << polYY << ", " << polZZ << G4endl;
178  }
179  }
180 
181  return mfp;
182 }
183 
184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
185 
187  const G4Track& track,
188  G4double previousStepSize,
190 {
191  G4double mfp = G4VEmProcess::PostStepGetPhysicalInteractionLength(track, previousStepSize, condition);
192 
193  if (theAsymmetryTable) {
194 
195  G4Material* aMaterial = track.GetMaterial();
196  G4VPhysicalVolume* aPVolume = track.GetVolume();
197  G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
198 
199  // G4Material* bMaterial = aLVolume->GetMaterial();
201 
202  const G4bool volumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
203  G4StokesVector electronPolarization = polarizationManger->GetVolumePolarization(aLVolume);
204 
205  if (!volumeIsPolarized || mfp == DBL_MAX) return mfp;
206 
207  // *** get asymmetry, if target is polarized ***
208  const G4DynamicParticle* aDynamicPositron = track.GetDynamicParticle();
209  const G4double positronEnergy = aDynamicPositron->GetKineticEnergy();
210  const G4StokesVector positronPolarization = track.GetPolarization();
211  const G4ParticleMomentum positronDirection0 = aDynamicPositron->GetMomentumDirection();
212 
213  if (verboseLevel>=2) {
214 
215  G4cout << " Mom " << positronDirection0 << G4endl;
216  G4cout << " Polarization " << positronPolarization << G4endl;
217  G4cout << " MaterialPol. " << electronPolarization << G4endl;
218  G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl;
219  G4cout << " Log. Volume " << aLVolume->GetName() << G4endl;
220  G4cout << " Material " << aMaterial << G4endl;
221  }
222 
223  G4bool isOutRange;
225  G4double lAsymmetry = (*theAsymmetryTable)(idx)->
226  GetValue(positronEnergy, isOutRange);
227  G4double tAsymmetry = (*theTransverseAsymmetryTable)(idx)->
228  GetValue(positronEnergy, isOutRange);
229 
230  G4double polZZ = positronPolarization.z()*
231  electronPolarization*positronDirection0;
232  G4double polXX = positronPolarization.x()*
233  electronPolarization*G4PolarizationHelper::GetParticleFrameX(positronDirection0);
234  G4double polYY = positronPolarization.y()*
235  electronPolarization*G4PolarizationHelper::GetParticleFrameY(positronDirection0);
236 
237  G4double impact = 1. + polZZ*lAsymmetry + (polXX + polYY)*tAsymmetry;
238 
239  mfp *= 1. / impact;
240 
241  if (verboseLevel>=2) {
242  G4cout << " MeanFreePath: " << mfp / mm << " mm " << G4endl;
243  G4cout << " Asymmetry: " << lAsymmetry << ", " << tAsymmetry << G4endl;
244  G4cout << " PolProduct: " << polXX << ", " << polYY << ", " << polZZ << G4endl;
245  }
246  }
247 
248  return mfp;
249 }
250 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
251 
253 {
255  BuildAsymmetryTable(pd);
256 }
257 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
258 
260 {
264 }
265 
267 {
268  // Access to materials
269  const G4ProductionCutsTable* theCoupleTable=
271  size_t numOfCouples = theCoupleTable->GetTableSize();
272  G4cout<<" annih-numOfCouples="<<numOfCouples<<"\n";
273  for(size_t i=0; i<numOfCouples; ++i) {
274  G4cout<<"annih- "<<i<<"/"<<numOfCouples<<"\n";
275  if (!theAsymmetryTable) break;
276  G4cout<<"annih- "<<theAsymmetryTable->GetFlag(i)<<"\n";
277  if (theAsymmetryTable->GetFlag(i)) {
278  G4cout<<" building pol-annih ... \n";
279 
280  // create physics vector and fill it
281  const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
282 
283  // use same parameters as for lambda
284  G4PhysicsVector* aVector = LambdaPhysicsVector(couple);
285  G4PhysicsVector* tVector = LambdaPhysicsVector(couple);
286 
287  for (G4int j = 0 ; j < LambdaBinning() ; ++j ) {
288  G4double lowEdgeEnergy = aVector->GetLowEdgeEnergy(j);
289  G4double tasm=0.;
290  G4double asym = ComputeAsymmetry(lowEdgeEnergy, couple, part, 0., tasm);
291  aVector->PutValue(j,asym);
292  tVector->PutValue(j,tasm);
293  }
294 
297  }
298  }
299 
300 }
301 
302 
303 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
304 
306  const G4MaterialCutsCouple* couple,
307  const G4ParticleDefinition& aParticle,
308  G4double cut,
309  G4double &tAsymmetry)
310 {
311  G4double lAsymmetry = 0.0;
312  tAsymmetry = 0.0;
313 
314  // calculate polarized cross section
318  G4double sigma2=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
319 
320  // calculate transversely polarized cross section
324  G4double sigma3=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
325 
326  // calculate unpolarized cross section
330  G4double sigma0=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
331 
332  // determine assymmetries
333  if (sigma0>0.) {
334  lAsymmetry=sigma2/sigma0-1.;
335  tAsymmetry=sigma3/sigma0-1.;
336  }
337  return lAsymmetry;
338 
339 }
340 
341 
342 
343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
344 
346 {
347  G4cout << " Polarized model for annihilation into 2 photons"
348  << G4endl;
349 }
350 
351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
352 
354  const G4Step& )
355 //
356 // Performs the e+ e- annihilation when both particles are assumed at rest.
357 // It generates two back to back photons with energy = electron_mass.
358 // The angular distribution is isotropic.
359 // GEANT4 internal units
360 //
361 // Note : Effects due to binding of atomic electrons are negliged.
362 {
364 
366 
367  G4double cosTeta = 2.*G4UniformRand()-1. , sinTeta = std::sqrt(1.-cosTeta*cosTeta);
368  G4double phi = twopi * G4UniformRand();
369  G4ThreeVector direction (sinTeta*std::cos(phi), sinTeta*std::sin(phi), cosTeta);
371  direction, electron_mass_c2) );
373  -direction, electron_mass_c2) );
374  // Kill the incident positron
375  //
377  return &fParticleChange;
378 }
379 
380 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4double condition(const G4ErrorSymMatrix &m)
void SetBeamPolarization(const G4ThreeVector &pBeam)
const G4ThreeVector & GetPolarization() const
static void SetPhysicsVector(G4PhysicsTable *physTable, size_t idx, G4PhysicsVector *vec)
G4int verboseLevel
Definition: G4VProcess.hh:368
G4String GetName() const
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
const G4DynamicParticle * GetDynamicParticle() const
void SetBuildTableFlag(G4bool val)
G4String name
Definition: TRTMaterials.hh:40
G4ParticleChangeForGamma fParticleChange
static G4PolarizationManager * GetInstance()
virtual G4VParticleChange * AtRestDoIt(const G4Track &track, const G4Step &stepData)
void BuildPhysicsTable(const G4ParticleDefinition &)
virtual void InitialiseProcess(const G4ParticleDefinition *)
G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
G4double GetLowEdgeEnergy(size_t binNumber) const
void SetStartFromNullFlag(G4bool val)
G4bool enableAtRestDoIt
Definition: G4VProcess.hh:350
int G4int
Definition: G4Types.hh:78
static G4PhysicsTable * PreparePhysicsTable(G4PhysicsTable *physTable)
G4int LambdaBinning() const
G4double ComputeAsymmetry(G4double energy, const G4MaterialCutsCouple *couple, const G4ParticleDefinition &particle, G4double cut, G4double &tasm)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:706
void AddSecondary(G4DynamicParticle *aParticle)
void SetLambdaBinning(G4int nbins)
size_t CurrentMaterialCutsCoupleIndex() const
#define G4UniformRand()
Definition: Randomize.hh:95
G4GLOB_DLL std::ostream G4cout
const G4String & GetName() const
void SetTargetPolarization(const G4ThreeVector &pTarget)
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
G4double CrossSection(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:482
void PutValue(size_t index, G4double theValue)
G4double GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
Definition: G4Step.hh:76
void PreparePhysicsTable(const G4ParticleDefinition &)
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
G4Material * GetMaterial() const
virtual void PreparePhysicsTable(const G4ParticleDefinition &)
G4double GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
static G4ProductionCutsTable * GetProductionCutsTable()
G4PhysicsVector * LambdaPhysicsVector(const G4MaterialCutsCouple *)
void SetMaxKinEnergy(G4double e)
G4LogicalVolume * GetLogicalVolume() const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
void AddEmModel(G4int, G4VEmModel *, const G4Region *region=0)
G4double energy(const ThreeVector &p, const G4double m)
void SetNumberOfSecondaries(G4int totSecondaries)
void SetSecondaryParticle(const G4ParticleDefinition *p)
G4PolarizedAnnihilationModel * emModel
bool IsPolarized(G4LogicalVolume *lVol) const
G4VPhysicalVolume * GetVolume() const
static G4ThreeVector GetParticleFrameY(const G4ThreeVector &)
#define G4endl
Definition: G4ios.hh:61
static const double TeV
Definition: G4SIunits.hh:197
static const double keV
Definition: G4SIunits.hh:195
void SetMinKinEnergy(G4double e)
G4bool GetFlag(size_t i) const
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:713
void InitializeForPostStep(const G4Track &)
G4ForceCondition
G4ThreeVector G4ParticleMomentum
#define DBL_MAX
Definition: templates.hh:83
static G4ThreeVector GetParticleFrameX(const G4ThreeVector &)
static const double mm
Definition: G4SIunits.hh:102
const G4ThreeVector & GetVolumePolarization(G4LogicalVolume *lVol) const
G4eplusPolarizedAnnihilation(const G4String &name="pol-annihil")
void BuildAsymmetryTable(const G4ParticleDefinition &part)