Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 97175 2016-05-27 12:42:05Z 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 
66 #include "G4PhysicsTableHelper.hh"
67 #include "G4ProductionCutsTable.hh"
68 #include "G4PolarizationManager.hh"
69 #include "G4PolarizationHelper.hh"
70 #include "G4StokesVector.hh"
71 
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73 
75  : G4eplusAnnihilation(name), isInitialised(false),
76  theAsymmetryTable(nullptr),
77  theTransverseAsymmetryTable(nullptr)
78 {
79  emModel = new G4PolarizedAnnihilationModel();
80  SetEmModel(emModel, 1);
81 }
82 
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84 
86 {
87  CleanTables();
88 }
89 
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91 
92 void G4eplusPolarizedAnnihilation::CleanTables()
93 {
94  if(theAsymmetryTable) {
95  theAsymmetryTable->clearAndDestroy();
96  delete theAsymmetryTable;
97  theAsymmetryTable = nullptr;
98  }
99  if(theTransverseAsymmetryTable) {
100  theTransverseAsymmetryTable->clearAndDestroy();
101  delete theTransverseAsymmetryTable;
102  theTransverseAsymmetryTable = nullptr;
103  }
104 }
105 
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107 
109  G4double previousStepSize,
111 {
112  G4double mfp = G4VEmProcess::GetMeanFreePath(track, previousStepSize, condition);
113 
114  if(theAsymmetryTable && theTransverseAsymmetryTable && mfp < DBL_MAX) {
115  mfp *= ComputeSaturationFactor(track);
116  }
117  if (verboseLevel>=2) {
118  G4cout << "G4eplusPolarizedAnnihilation::MeanFreePath: "
119  << mfp / mm << " mm " << G4endl;
120  }
121  return mfp;
122 }
123 
124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
125 
127  const G4Track& track,
128  G4double previousStepSize,
130 {
131  // save previous values
134 
135  // *** compute unpolarized step limit ***
136  // this changes theNumberOfInteractionLengthLeft and currentInteractionLength
138  previousStepSize,
139  condition);
140  G4double x0 = x;
141  G4double satFact = 1.0;
142 
143  // *** add corrections on polarisation ***
144  if(theAsymmetryTable && theTransverseAsymmetryTable && x < DBL_MAX) {
145  satFact = ComputeSaturationFactor(track);
146  G4double curLength = currentInteractionLength*satFact;
147  G4double prvLength = iLength*satFact;
148  if(nLength > 0.0) {
150  std::max(nLength - previousStepSize/prvLength, 0.0);
151  }
152  x = theNumberOfInteractionLengthLeft * curLength;
153  }
154  if (verboseLevel>=2) {
155  G4cout << "G4eplusPolarizedAnnihilation::PostStepGPIL: "
156  << std::setprecision(8) << x/mm << " mm;" << G4endl
157  << " unpolarized value: "
158  << std::setprecision(8) << x0/mm << " mm." << G4endl;
159  }
160  return x;
161 }
162 
163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
164 
165 G4double
166 G4eplusPolarizedAnnihilation::ComputeSaturationFactor(const G4Track& track)
167 {
168  G4Material* aMaterial = track.GetMaterial();
169  G4VPhysicalVolume* aPVolume = track.GetVolume();
170  G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
171 
173 
174  const G4bool volumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
175  G4StokesVector electronPolarization = polarizationManger->GetVolumePolarization(aLVolume);
176 
177  G4double factor = 1.0;
178 
179  if (volumeIsPolarized) {
180 
181  // *** get asymmetry, if target is polarized ***
182  const G4DynamicParticle* aDynamicPositron = track.GetDynamicParticle();
183  const G4double positronEnergy = aDynamicPositron->GetKineticEnergy();
184  const G4StokesVector positronPolarization = track.GetPolarization();
185  const G4ParticleMomentum positronDirection0 = aDynamicPositron->GetMomentumDirection();
186 
187  if (verboseLevel>=2) {
188  G4cout << "G4eplusPolarizedAnnihilation::ComputeSaturationFactor: " << G4endl;
189  G4cout << " Mom " << positronDirection0 << G4endl;
190  G4cout << " Polarization " << positronPolarization << G4endl;
191  G4cout << " MaterialPol. " << electronPolarization << G4endl;
192  G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl;
193  G4cout << " Log. Volume " << aLVolume->GetName() << G4endl;
194  G4cout << " Material " << aMaterial << G4endl;
195  }
196 
197  size_t midx = CurrentMaterialCutsCoupleIndex();
198  const G4PhysicsVector* aVector = nullptr;
199  const G4PhysicsVector* bVector = nullptr;
200  if(midx < theAsymmetryTable->size()) {
201  aVector = (*theAsymmetryTable)(midx);
202  }
203  if(midx < theTransverseAsymmetryTable->size()) {
204  bVector = (*theTransverseAsymmetryTable)(midx);
205  }
206  if(aVector && bVector) {
207  G4double lAsymmetry = aVector->Value(positronEnergy);
208  G4double tAsymmetry = bVector->Value(positronEnergy);
209  G4double polZZ = positronPolarization.z()*
210  (electronPolarization*positronDirection0);
211  G4double polXX = positronPolarization.x()*
212  (electronPolarization*G4PolarizationHelper::GetParticleFrameX(positronDirection0));
213  G4double polYY = positronPolarization.y()*
214  (electronPolarization*G4PolarizationHelper::GetParticleFrameY(positronDirection0));
215 
216  factor /= (1. + polZZ*lAsymmetry + (polXX + polYY)*tAsymmetry);
217 
218  if (verboseLevel>=2) {
219  G4cout << " Asymmetry: " << lAsymmetry << ", " << tAsymmetry << G4endl;
220  G4cout << " PolProduct: " << polXX << ", " << polYY << ", " << polZZ << G4endl;
221  G4cout << " Factor: " << factor << G4endl;
222  }
223  } else {
225  ed << "Problem with asymmetry tables: material index " << midx
226  << " is out of range or tables are not filled";
227  G4Exception("G4eplusPolarizedAnnihilation::ComputeSaturationFactor","em0048",
228  JustWarning, ed, "");
229  }
230  }
231  return factor;
232 }
233 
234 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
235 
237  const G4ParticleDefinition& part)
238 {
240  G4bool isMaster = true;
241  const G4eplusPolarizedAnnihilation* masterProcess =
242  static_cast<const G4eplusPolarizedAnnihilation*>(GetMasterProcess());
243  if(masterProcess && masterProcess != this) { isMaster = false; }
244  if(isMaster) { BuildAsymmetryTables(part); }
245 }
246 
247 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
248 
249 void G4eplusPolarizedAnnihilation::BuildAsymmetryTables(
250  const G4ParticleDefinition& part)
251 {
252  // cleanup old, initialise new table
253  CleanTables();
254  theAsymmetryTable =
256  theTransverseAsymmetryTable =
257  G4PhysicsTableHelper::PreparePhysicsTable(theTransverseAsymmetryTable);
258 
259  // Access to materials
260  const G4ProductionCutsTable* theCoupleTable=
262  size_t numOfCouples = theCoupleTable->GetTableSize();
263  //G4cout<<" annih-numOfCouples="<<numOfCouples<<"\n";
264  for(size_t i=0; i<numOfCouples; ++i) {
265  //G4cout<<"annih- "<<i<<"/"<<numOfCouples<<"\n";
266  if (!theAsymmetryTable) break;
267  //G4cout<<"annih- "<<theAsymmetryTable->GetFlag(i)<<"\n";
268  if (theAsymmetryTable->GetFlag(i)) {
269  //G4cout<<" building pol-annih ... \n";
270 
271  // create physics vector and fill it
272  const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
273 
274  // use same parameters as for lambda
275  G4PhysicsVector* aVector = LambdaPhysicsVector(couple);
276  G4PhysicsVector* tVector = LambdaPhysicsVector(couple);
277 
278  for (G4int j = 0 ; j < LambdaBinning() ; ++j ) {
279  G4double lowEdgeEnergy = aVector->GetLowEdgeEnergy(j);
280  G4double tasm=0.;
281  G4double asym = ComputeAsymmetry(lowEdgeEnergy, couple, part, 0., tasm);
282  aVector->PutValue(j,asym);
283  tVector->PutValue(j,tasm);
284  }
285  G4PhysicsTableHelper::SetPhysicsVector(theAsymmetryTable, i, aVector);
286  G4PhysicsTableHelper::SetPhysicsVector(theTransverseAsymmetryTable, i, tVector);
287  }
288  }
289 }
290 
291 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
292 
293 G4double G4eplusPolarizedAnnihilation::ComputeAsymmetry(G4double energy,
294  const G4MaterialCutsCouple* couple,
295  const G4ParticleDefinition& aParticle,
296  G4double cut,
297  G4double &tAsymmetry)
298 {
299  G4double lAsymmetry = 0.0;
300  tAsymmetry = 0.0;
301 
302  // calculate polarized cross section
303  theTargetPolarization=G4ThreeVector(0.,0.,1.);
304  emModel->SetTargetPolarization(theTargetPolarization);
305  emModel->SetBeamPolarization(theTargetPolarization);
306  G4double sigma2=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
307 
308  // calculate transversely polarized cross section
309  theTargetPolarization=G4ThreeVector(1.,0.,0.);
310  emModel->SetTargetPolarization(theTargetPolarization);
311  emModel->SetBeamPolarization(theTargetPolarization);
312  G4double sigma3=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
313 
314  // calculate unpolarized cross section
315  theTargetPolarization=G4ThreeVector();
316  emModel->SetTargetPolarization(theTargetPolarization);
317  emModel->SetBeamPolarization(theTargetPolarization);
318  G4double sigma0=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
319 
320  // determine assymmetries
321  if (sigma0>0.) {
322  lAsymmetry=sigma2/sigma0-1.;
323  tAsymmetry=sigma3/sigma0-1.;
324  }
325  return lAsymmetry;
326 }
327 
328 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
329 
331 {
332  G4cout << " Polarized model for annihilation into 2 photons"
333  << G4endl;
334 }
335 
336 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4double condition(const G4ErrorSymMatrix &m)
const G4VProcess * GetMasterProcess() const
Definition: G4VProcess.hh:538
void SetBeamPolarization(const G4ThreeVector &pBeam)
const XML_Char * name
Definition: expat.h:151
const G4ThreeVector & GetPolarization() const
static void SetPhysicsVector(G4PhysicsTable *physTable, size_t idx, G4PhysicsVector *vec)
virtual void BuildPhysicsTable(const G4ParticleDefinition &) override
static constexpr double mm
Definition: G4SIunits.hh:115
G4int verboseLevel
Definition: G4VProcess.hh:368
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
double x() const
const G4DynamicParticle * GetDynamicParticle() const
static G4PolarizationManager * GetInstance()
G4double theNumberOfInteractionLengthLeft
Definition: G4VProcess.hh:293
virtual G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition) override
tuple x
Definition: test.py:50
G4double GetLowEdgeEnergy(size_t binNumber) const
virtual G4double GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition) override
int G4int
Definition: G4Types.hh:78
static G4PhysicsTable * PreparePhysicsTable(G4PhysicsTable *physTable)
double z() const
G4int LambdaBinning() const
void SetEmModel(G4VEmModel *, G4int index=1)
size_t CurrentMaterialCutsCoupleIndex() const
G4GLOB_DLL std::ostream G4cout
const G4String & GetName() const
void SetTargetPolarization(const G4ThreeVector &pTarget)
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
virtual void BuildPhysicsTable(const G4ParticleDefinition &) override
G4double CrossSection(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:503
G4double currentInteractionLength
Definition: G4VProcess.hh:297
void PutValue(size_t index, G4double theValue)
G4double Value(G4double theEnergy, size_t &lastidx) const
G4Material * GetMaterial() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4ProductionCutsTable * GetProductionCutsTable()
G4PhysicsVector * LambdaPhysicsVector(const G4MaterialCutsCouple *)
G4LogicalVolume * GetLogicalVolume() const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
virtual G4double GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition) override
double y() const
bool IsPolarized(G4LogicalVolume *lVol) const
G4VPhysicalVolume * GetVolume() const
static G4ThreeVector GetParticleFrameY(const G4ThreeVector &)
#define G4endl
Definition: G4ios.hh:61
G4bool GetFlag(size_t i) const
const G4String & GetName() const
double G4double
Definition: G4Types.hh:76
G4ForceCondition
#define DBL_MAX
Definition: templates.hh:83
static G4ThreeVector GetParticleFrameX(const G4ThreeVector &)
const G4ThreeVector & GetVolumePolarization(G4LogicalVolume *lVol) const
void clearAndDestroy()
G4eplusPolarizedAnnihilation(const G4String &name="pol-annihil")
virtual G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition) override