Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/polarisation/src/G4PolarizedAnnihilation.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 // -------------------------------------------------------------------
 27 //
 28 // Geant4 Class file
 29 //
 30 // File name:     G4PolarizedAnnihilation
 31 //
 32 // Author:  A. Schaelicke on base of Vladimir Ivanchenko / Michel Maire code
 33 //
 34 // Class Description:
 35 //   Polarized process of e+ annihilation into 2 gammas
 36 
 37 #include "G4PolarizedAnnihilation.hh"
 38 
 39 #include "G4DynamicParticle.hh"
 40 #include "G4MaterialCutsCouple.hh"
 41 #include "G4PhysicsTableHelper.hh"
 42 #include "G4PhysicsVector.hh"
 43 #include "G4PolarizationHelper.hh"
 44 #include "G4PolarizationManager.hh"
 45 #include "G4PolarizedAnnihilationModel.hh"
 46 #include "G4ProductionCutsTable.hh"
 47 #include "G4SystemOfUnits.hh"
 48 #include "G4StokesVector.hh"
 49 
 50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 51 G4PolarizedAnnihilation::G4PolarizedAnnihilation(const G4String& name)
 52   : G4eplusAnnihilation(name)
 53   , fAsymmetryTable(nullptr)
 54   , fTransverseAsymmetryTable(nullptr)
 55 {
 56   fEmModel = new G4PolarizedAnnihilationModel();
 57   SetEmModel(fEmModel);
 58 }
 59 
 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 61 G4PolarizedAnnihilation::~G4PolarizedAnnihilation() { CleanTables(); }
 62 
 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 64 void G4PolarizedAnnihilation::CleanTables()
 65 {
 66   if(fAsymmetryTable)
 67   {
 68     fAsymmetryTable->clearAndDestroy();
 69     delete fAsymmetryTable;
 70     fAsymmetryTable = nullptr;
 71   }
 72   if(fTransverseAsymmetryTable)
 73   {
 74     fTransverseAsymmetryTable->clearAndDestroy();
 75     delete fTransverseAsymmetryTable;
 76     fTransverseAsymmetryTable = nullptr;
 77   }
 78 }
 79 
 80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 81 G4double G4PolarizedAnnihilation::GetMeanFreePath(const G4Track& track,
 82                                                   G4double previousStepSize,
 83                                                   G4ForceCondition* condition)
 84 {
 85   G4double mfp =
 86     G4VEmProcess::GetMeanFreePath(track, previousStepSize, condition);
 87 
 88   if(nullptr != fAsymmetryTable && nullptr != fTransverseAsymmetryTable && mfp < DBL_MAX)
 89   {
 90     mfp *= ComputeSaturationFactor(track);
 91   }
 92   if(verboseLevel >= 2)
 93   {
 94     G4cout << "G4PolarizedAnnihilation::MeanFreePath:  " << mfp / mm << " mm "
 95            << G4endl;
 96   }
 97   return mfp;
 98 }
 99 
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101 G4double G4PolarizedAnnihilation::PostStepGetPhysicalInteractionLength(
102   const G4Track& track, G4double previousStepSize, G4ForceCondition* condition)
103 {
104   // save previous values
105   G4double nLength = theNumberOfInteractionLengthLeft;
106   G4double iLength = currentInteractionLength;
107 
108   // *** compute unpolarized step limit ***
109   // this changes theNumberOfInteractionLengthLeft and currentInteractionLength
110   G4double x = G4VEmProcess::PostStepGetPhysicalInteractionLength(
111     track, previousStepSize, condition);
112   G4double x0      = x;
113   G4double satFact = 1.0;
114 
115   // *** add corrections on polarisation ***
116   if(nullptr != fAsymmetryTable && nullptr != fTransverseAsymmetryTable && x < DBL_MAX)
117   {
118     satFact            = ComputeSaturationFactor(track);
119     G4double curLength = currentInteractionLength * satFact;
120     G4double prvLength = iLength * satFact;
121     if(nLength > 0.0)
122     {
123       theNumberOfInteractionLengthLeft =
124         std::max(nLength - previousStepSize / prvLength, 0.0);
125     }
126     x = theNumberOfInteractionLengthLeft * curLength;
127   }
128   if(verboseLevel >= 2)
129   {
130     G4cout << "G4PolarizedAnnihilation::PostStepGPIL: " << std::setprecision(8)
131            << x / mm << " mm;" << G4endl
132            << "                         unpolarized value: "
133            << std::setprecision(8) << x0 / mm << " mm." << G4endl;
134   }
135   return x;
136 }
137 
138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
139 G4double G4PolarizedAnnihilation::ComputeSaturationFactor(const G4Track& track)
140 {
141   const G4Material* aMaterial = track.GetMaterial();
142   G4VPhysicalVolume* aPVolume = track.GetVolume();
143   G4LogicalVolume* aLVolume   = aPVolume->GetLogicalVolume();
144 
145   G4PolarizationManager* polarizationManager =
146     G4PolarizationManager::GetInstance();
147 
148   const G4bool volumeIsPolarized = polarizationManager->IsPolarized(aLVolume);
149   G4StokesVector electronPolarization =
150     polarizationManager->GetVolumePolarization(aLVolume);
151 
152   G4double factor = 1.0;
153 
154   if(volumeIsPolarized)
155   {
156     // *** get asymmetry, if target is polarized ***
157     const G4DynamicParticle* aDynamicPositron = track.GetDynamicParticle();
158     const G4double positronEnergy = aDynamicPositron->GetKineticEnergy();
159     const G4StokesVector positronPolarization =
160       G4StokesVector(track.GetPolarization());
161     const G4ParticleMomentum positronDirection0 =
162       aDynamicPositron->GetMomentumDirection();
163 
164     if(verboseLevel >= 2)
165     {
166       G4cout << "G4PolarizedAnnihilation::ComputeSaturationFactor: " << G4endl;
167       G4cout << " Mom " << positronDirection0 << G4endl;
168       G4cout << " Polarization " << positronPolarization << G4endl;
169       G4cout << " MaterialPol. " << electronPolarization << G4endl;
170       G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl;
171       G4cout << " Log. Volume  " << aLVolume->GetName() << G4endl;
172       G4cout << " Material     " << aMaterial << G4endl;
173     }
174 
175     std::size_t midx               = CurrentMaterialCutsCoupleIndex();
176     const G4PhysicsVector* aVector = nullptr;
177     const G4PhysicsVector* bVector = nullptr;
178     if(midx < fAsymmetryTable->size())
179     {
180       aVector = (*fAsymmetryTable)(midx);
181     }
182     if(midx < fTransverseAsymmetryTable->size())
183     {
184       bVector = (*fTransverseAsymmetryTable)(midx);
185     }
186     if(aVector && bVector)
187     {
188       G4double lAsymmetry = aVector->Value(positronEnergy);
189       G4double tAsymmetry = bVector->Value(positronEnergy);
190       G4double polZZ =
191         positronPolarization.z() * (electronPolarization * positronDirection0);
192       G4double polXX =
193         positronPolarization.x() *
194         (electronPolarization *
195          G4PolarizationHelper::GetParticleFrameX(positronDirection0));
196       G4double polYY =
197         positronPolarization.y() *
198         (electronPolarization *
199          G4PolarizationHelper::GetParticleFrameY(positronDirection0));
200 
201       factor /= (1. + polZZ * lAsymmetry + (polXX + polYY) * tAsymmetry);
202 
203       if(verboseLevel >= 2)
204       {
205         G4cout << " Asymmetry:     " << lAsymmetry << ", " << tAsymmetry
206                << G4endl;
207         G4cout << " PolProduct:    " << polXX << ", " << polYY << ", " << polZZ
208                << G4endl;
209         G4cout << " Factor:        " << factor << G4endl;
210       }
211     }
212     else
213     {
214       G4ExceptionDescription ed;
215       ed << "Problem with asymmetry tables: material index " << midx
216          << " is out of range or tables are not filled";
217       G4Exception("G4PolarizedAnnihilation::ComputeSaturationFactor", "em0048",
218                   JustWarning, ed, "");
219     }
220   }
221   return factor;
222 }
223 
224 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
225 void G4PolarizedAnnihilation::BuildPhysicsTable(
226   const G4ParticleDefinition& part)
227 {
228   G4VEmProcess::BuildPhysicsTable(part);
229   if(isTheMaster)
230   {
231     BuildAsymmetryTables(part);
232   }
233 }
234 
235 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
236 void G4PolarizedAnnihilation::BuildAsymmetryTables(
237   const G4ParticleDefinition& part)
238 {
239   // cleanup old, initialise new table
240   CleanTables();
241   fAsymmetryTable = G4PhysicsTableHelper::PreparePhysicsTable(fAsymmetryTable);
242   fTransverseAsymmetryTable =
243     G4PhysicsTableHelper::PreparePhysicsTable(fTransverseAsymmetryTable);
244   if(nullptr == fAsymmetryTable) return;
245 
246   // Access to materials
247   const G4ProductionCutsTable* theCoupleTable =
248     G4ProductionCutsTable::GetProductionCutsTable();
249   G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
250   for(G4int i = 0; i < numOfCouples; ++i)
251   {
252     if(fAsymmetryTable->GetFlag(i))
253     {
254       // create physics vector and fill it
255       const G4MaterialCutsCouple* couple =
256         theCoupleTable->GetMaterialCutsCouple(i);
257 
258       // use same parameters as for lambda
259       G4PhysicsVector* aVector = LambdaPhysicsVector(couple);
260       G4PhysicsVector* tVector = LambdaPhysicsVector(couple);
261       G4int nn = (G4int)aVector->GetVectorLength();
262       for(G4int j = 0; j < nn; ++j)
263       {
264         G4double energy = aVector->Energy(j);
265         G4double tasm   = 0.;
266         G4double asym = ComputeAsymmetry(energy, couple, part, 0., tasm);
267         aVector->PutValue(j, asym);
268         tVector->PutValue(j, tasm);
269       }
270       if(aVector->GetSpline()) {
271   aVector->FillSecondDerivatives();
272   tVector->FillSecondDerivatives();
273       }
274       G4PhysicsTableHelper::SetPhysicsVector(fAsymmetryTable, i, aVector);
275       G4PhysicsTableHelper::SetPhysicsVector(fTransverseAsymmetryTable, i,
276                                              tVector);
277     }
278   }
279 }
280 
281 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
282 G4double G4PolarizedAnnihilation::ComputeAsymmetry(
283   G4double energy, const G4MaterialCutsCouple* couple,
284   const G4ParticleDefinition& aParticle, G4double cut, G4double& tAsymmetry)
285 {
286   G4double lAsymmetry = 0.0;
287   tAsymmetry          = 0.0;
288 
289   // calculate polarized cross section
290   G4ThreeVector targetPolarization = G4ThreeVector(0., 0., 1.);
291   fEmModel->SetTargetPolarization(targetPolarization);
292   fEmModel->SetBeamPolarization(targetPolarization);
293   G4double sigma2 =
294     fEmModel->CrossSection(couple, &aParticle, energy, cut, energy);
295 
296   // calculate transversely polarized cross section
297   targetPolarization = G4ThreeVector(1., 0., 0.);
298   fEmModel->SetTargetPolarization(targetPolarization);
299   fEmModel->SetBeamPolarization(targetPolarization);
300   G4double sigma3 =
301     fEmModel->CrossSection(couple, &aParticle, energy, cut, energy);
302 
303   // calculate unpolarized cross section
304   targetPolarization = G4ThreeVector();
305   fEmModel->SetTargetPolarization(targetPolarization);
306   fEmModel->SetBeamPolarization(targetPolarization);
307   G4double sigma0 =
308     fEmModel->CrossSection(couple, &aParticle, energy, cut, energy);
309 
310   // determine asymmetries
311   if(sigma0 > 0.)
312   {
313     lAsymmetry = sigma2 / sigma0 - 1.;
314     tAsymmetry = sigma3 / sigma0 - 1.;
315   }
316   return lAsymmetry;
317 }
318 
319 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
320 void G4PolarizedAnnihilation::ProcessDescription(std::ostream& out) const
321 {
322   out << "Polarized model for positron annihilation into 2 photons.\n";
323 }
324