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 ]

Diff markup

Differences between /processes/electromagnetic/polarisation/src/G4PolarizedAnnihilation.cc (Version 11.3.0) and /processes/electromagnetic/polarisation/src/G4PolarizedAnnihilation.cc (Version 11.1.1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // -------------------------------------------     26 // -------------------------------------------------------------------
 27 //                                                 27 //
 28 // Geant4 Class file                               28 // Geant4 Class file
 29 //                                                 29 //
 30 // File name:     G4PolarizedAnnihilation          30 // File name:     G4PolarizedAnnihilation
 31 //                                                 31 //
 32 // Author:  A. Schaelicke on base of Vladimir      32 // Author:  A. Schaelicke on base of Vladimir Ivanchenko / Michel Maire code
 33 //                                                 33 //
 34 // Class Description:                              34 // Class Description:
 35 //   Polarized process of e+ annihilation into     35 //   Polarized process of e+ annihilation into 2 gammas
 36                                                    36 
 37 #include "G4PolarizedAnnihilation.hh"              37 #include "G4PolarizedAnnihilation.hh"
 38                                                    38 
 39 #include "G4DynamicParticle.hh"                    39 #include "G4DynamicParticle.hh"
 40 #include "G4MaterialCutsCouple.hh"                 40 #include "G4MaterialCutsCouple.hh"
 41 #include "G4PhysicsTableHelper.hh"                 41 #include "G4PhysicsTableHelper.hh"
 42 #include "G4PhysicsVector.hh"                      42 #include "G4PhysicsVector.hh"
 43 #include "G4PolarizationHelper.hh"                 43 #include "G4PolarizationHelper.hh"
 44 #include "G4PolarizationManager.hh"                44 #include "G4PolarizationManager.hh"
 45 #include "G4PolarizedAnnihilationModel.hh"         45 #include "G4PolarizedAnnihilationModel.hh"
 46 #include "G4ProductionCutsTable.hh"                46 #include "G4ProductionCutsTable.hh"
 47 #include "G4SystemOfUnits.hh"                      47 #include "G4SystemOfUnits.hh"
 48 #include "G4StokesVector.hh"                       48 #include "G4StokesVector.hh"
 49                                                    49 
 50 //....oooOO0OOooo........oooOO0OOooo........oo     50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 51 G4PolarizedAnnihilation::G4PolarizedAnnihilati     51 G4PolarizedAnnihilation::G4PolarizedAnnihilation(const G4String& name)
 52   : G4eplusAnnihilation(name)                      52   : G4eplusAnnihilation(name)
 53   , fAsymmetryTable(nullptr)                       53   , fAsymmetryTable(nullptr)
 54   , fTransverseAsymmetryTable(nullptr)             54   , fTransverseAsymmetryTable(nullptr)
 55 {                                                  55 {
 56   fEmModel = new G4PolarizedAnnihilationModel(     56   fEmModel = new G4PolarizedAnnihilationModel();
 57   SetEmModel(fEmModel);                            57   SetEmModel(fEmModel);
 58 }                                                  58 }
 59                                                    59 
 60 //....oooOO0OOooo........oooOO0OOooo........oo     60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 61 G4PolarizedAnnihilation::~G4PolarizedAnnihilat     61 G4PolarizedAnnihilation::~G4PolarizedAnnihilation() { CleanTables(); }
 62                                                    62 
 63 //....oooOO0OOooo........oooOO0OOooo........oo     63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 64 void G4PolarizedAnnihilation::CleanTables()        64 void G4PolarizedAnnihilation::CleanTables()
 65 {                                                  65 {
 66   if(fAsymmetryTable)                              66   if(fAsymmetryTable)
 67   {                                                67   {
 68     fAsymmetryTable->clearAndDestroy();            68     fAsymmetryTable->clearAndDestroy();
 69     delete fAsymmetryTable;                        69     delete fAsymmetryTable;
 70     fAsymmetryTable = nullptr;                     70     fAsymmetryTable = nullptr;
 71   }                                                71   }
 72   if(fTransverseAsymmetryTable)                    72   if(fTransverseAsymmetryTable)
 73   {                                                73   {
 74     fTransverseAsymmetryTable->clearAndDestroy     74     fTransverseAsymmetryTable->clearAndDestroy();
 75     delete fTransverseAsymmetryTable;              75     delete fTransverseAsymmetryTable;
 76     fTransverseAsymmetryTable = nullptr;           76     fTransverseAsymmetryTable = nullptr;
 77   }                                                77   }
 78 }                                                  78 }
 79                                                    79 
 80 //....oooOO0OOooo........oooOO0OOooo........oo     80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 81 G4double G4PolarizedAnnihilation::GetMeanFreeP     81 G4double G4PolarizedAnnihilation::GetMeanFreePath(const G4Track& track,
 82                                                    82                                                   G4double previousStepSize,
 83                                                    83                                                   G4ForceCondition* condition)
 84 {                                                  84 {
 85   G4double mfp =                                   85   G4double mfp =
 86     G4VEmProcess::GetMeanFreePath(track, previ     86     G4VEmProcess::GetMeanFreePath(track, previousStepSize, condition);
 87                                                    87 
 88   if(nullptr != fAsymmetryTable && nullptr !=      88   if(nullptr != fAsymmetryTable && nullptr != fTransverseAsymmetryTable && mfp < DBL_MAX)
 89   {                                                89   {
 90     mfp *= ComputeSaturationFactor(track);         90     mfp *= ComputeSaturationFactor(track);
 91   }                                                91   }
 92   if(verboseLevel >= 2)                            92   if(verboseLevel >= 2)
 93   {                                                93   {
 94     G4cout << "G4PolarizedAnnihilation::MeanFr     94     G4cout << "G4PolarizedAnnihilation::MeanFreePath:  " << mfp / mm << " mm "
 95            << G4endl;                              95            << G4endl;
 96   }                                                96   }
 97   return mfp;                                      97   return mfp;
 98 }                                                  98 }
 99                                                    99 
100 //....oooOO0OOooo........oooOO0OOooo........oo    100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101 G4double G4PolarizedAnnihilation::PostStepGetP    101 G4double G4PolarizedAnnihilation::PostStepGetPhysicalInteractionLength(
102   const G4Track& track, G4double previousStepS    102   const G4Track& track, G4double previousStepSize, G4ForceCondition* condition)
103 {                                                 103 {
104   // save previous values                         104   // save previous values
105   G4double nLength = theNumberOfInteractionLen    105   G4double nLength = theNumberOfInteractionLengthLeft;
106   G4double iLength = currentInteractionLength;    106   G4double iLength = currentInteractionLength;
107                                                   107 
108   // *** compute unpolarized step limit ***       108   // *** compute unpolarized step limit ***
109   // this changes theNumberOfInteractionLength    109   // this changes theNumberOfInteractionLengthLeft and currentInteractionLength
110   G4double x = G4VEmProcess::PostStepGetPhysic    110   G4double x = G4VEmProcess::PostStepGetPhysicalInteractionLength(
111     track, previousStepSize, condition);          111     track, previousStepSize, condition);
112   G4double x0      = x;                           112   G4double x0      = x;
113   G4double satFact = 1.0;                         113   G4double satFact = 1.0;
114                                                   114 
115   // *** add corrections on polarisation ***      115   // *** add corrections on polarisation ***
116   if(nullptr != fAsymmetryTable && nullptr !=     116   if(nullptr != fAsymmetryTable && nullptr != fTransverseAsymmetryTable && x < DBL_MAX)
117   {                                               117   {
118     satFact            = ComputeSaturationFact    118     satFact            = ComputeSaturationFactor(track);
119     G4double curLength = currentInteractionLen    119     G4double curLength = currentInteractionLength * satFact;
120     G4double prvLength = iLength * satFact;       120     G4double prvLength = iLength * satFact;
121     if(nLength > 0.0)                             121     if(nLength > 0.0)
122     {                                             122     {
123       theNumberOfInteractionLengthLeft =          123       theNumberOfInteractionLengthLeft =
124         std::max(nLength - previousStepSize /     124         std::max(nLength - previousStepSize / prvLength, 0.0);
125     }                                             125     }
126     x = theNumberOfInteractionLengthLeft * cur    126     x = theNumberOfInteractionLengthLeft * curLength;
127   }                                               127   }
128   if(verboseLevel >= 2)                           128   if(verboseLevel >= 2)
129   {                                               129   {
130     G4cout << "G4PolarizedAnnihilation::PostSt    130     G4cout << "G4PolarizedAnnihilation::PostStepGPIL: " << std::setprecision(8)
131            << x / mm << " mm;" << G4endl          131            << x / mm << " mm;" << G4endl
132            << "                         unpola    132            << "                         unpolarized value: "
133            << std::setprecision(8) << x0 / mm     133            << std::setprecision(8) << x0 / mm << " mm." << G4endl;
134   }                                               134   }
135   return x;                                       135   return x;
136 }                                                 136 }
137                                                   137 
138 //....oooOO0OOooo........oooOO0OOooo........oo    138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
139 G4double G4PolarizedAnnihilation::ComputeSatur    139 G4double G4PolarizedAnnihilation::ComputeSaturationFactor(const G4Track& track)
140 {                                                 140 {
141   const G4Material* aMaterial = track.GetMater    141   const G4Material* aMaterial = track.GetMaterial();
142   G4VPhysicalVolume* aPVolume = track.GetVolum    142   G4VPhysicalVolume* aPVolume = track.GetVolume();
143   G4LogicalVolume* aLVolume   = aPVolume->GetL    143   G4LogicalVolume* aLVolume   = aPVolume->GetLogicalVolume();
144                                                   144 
145   G4PolarizationManager* polarizationManager =    145   G4PolarizationManager* polarizationManager =
146     G4PolarizationManager::GetInstance();         146     G4PolarizationManager::GetInstance();
147                                                   147 
148   const G4bool volumeIsPolarized = polarizatio    148   const G4bool volumeIsPolarized = polarizationManager->IsPolarized(aLVolume);
149   G4StokesVector electronPolarization =           149   G4StokesVector electronPolarization =
150     polarizationManager->GetVolumePolarization    150     polarizationManager->GetVolumePolarization(aLVolume);
151                                                   151 
152   G4double factor = 1.0;                          152   G4double factor = 1.0;
153                                                   153 
154   if(volumeIsPolarized)                           154   if(volumeIsPolarized)
155   {                                               155   {
156     // *** get asymmetry, if target is polariz    156     // *** get asymmetry, if target is polarized ***
157     const G4DynamicParticle* aDynamicPositron     157     const G4DynamicParticle* aDynamicPositron = track.GetDynamicParticle();
158     const G4double positronEnergy = aDynamicPo    158     const G4double positronEnergy = aDynamicPositron->GetKineticEnergy();
159     const G4StokesVector positronPolarization     159     const G4StokesVector positronPolarization =
160       G4StokesVector(track.GetPolarization());    160       G4StokesVector(track.GetPolarization());
161     const G4ParticleMomentum positronDirection    161     const G4ParticleMomentum positronDirection0 =
162       aDynamicPositron->GetMomentumDirection()    162       aDynamicPositron->GetMomentumDirection();
163                                                   163 
164     if(verboseLevel >= 2)                         164     if(verboseLevel >= 2)
165     {                                             165     {
166       G4cout << "G4PolarizedAnnihilation::Comp    166       G4cout << "G4PolarizedAnnihilation::ComputeSaturationFactor: " << G4endl;
167       G4cout << " Mom " << positronDirection0     167       G4cout << " Mom " << positronDirection0 << G4endl;
168       G4cout << " Polarization " << positronPo    168       G4cout << " Polarization " << positronPolarization << G4endl;
169       G4cout << " MaterialPol. " << electronPo    169       G4cout << " MaterialPol. " << electronPolarization << G4endl;
170       G4cout << " Phys. Volume " << aPVolume->    170       G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl;
171       G4cout << " Log. Volume  " << aLVolume->    171       G4cout << " Log. Volume  " << aLVolume->GetName() << G4endl;
172       G4cout << " Material     " << aMaterial     172       G4cout << " Material     " << aMaterial << G4endl;
173     }                                             173     }
174                                                   174 
175     std::size_t midx               = CurrentMa    175     std::size_t midx               = CurrentMaterialCutsCoupleIndex();
176     const G4PhysicsVector* aVector = nullptr;     176     const G4PhysicsVector* aVector = nullptr;
177     const G4PhysicsVector* bVector = nullptr;     177     const G4PhysicsVector* bVector = nullptr;
178     if(midx < fAsymmetryTable->size())            178     if(midx < fAsymmetryTable->size())
179     {                                             179     {
180       aVector = (*fAsymmetryTable)(midx);         180       aVector = (*fAsymmetryTable)(midx);
181     }                                             181     }
182     if(midx < fTransverseAsymmetryTable->size(    182     if(midx < fTransverseAsymmetryTable->size())
183     {                                             183     {
184       bVector = (*fTransverseAsymmetryTable)(m    184       bVector = (*fTransverseAsymmetryTable)(midx);
185     }                                             185     }
186     if(aVector && bVector)                        186     if(aVector && bVector)
187     {                                             187     {
188       G4double lAsymmetry = aVector->Value(pos    188       G4double lAsymmetry = aVector->Value(positronEnergy);
189       G4double tAsymmetry = bVector->Value(pos    189       G4double tAsymmetry = bVector->Value(positronEnergy);
190       G4double polZZ =                            190       G4double polZZ =
191         positronPolarization.z() * (electronPo    191         positronPolarization.z() * (electronPolarization * positronDirection0);
192       G4double polXX =                            192       G4double polXX =
193         positronPolarization.x() *                193         positronPolarization.x() *
194         (electronPolarization *                   194         (electronPolarization *
195          G4PolarizationHelper::GetParticleFram    195          G4PolarizationHelper::GetParticleFrameX(positronDirection0));
196       G4double polYY =                            196       G4double polYY =
197         positronPolarization.y() *                197         positronPolarization.y() *
198         (electronPolarization *                   198         (electronPolarization *
199          G4PolarizationHelper::GetParticleFram    199          G4PolarizationHelper::GetParticleFrameY(positronDirection0));
200                                                   200 
201       factor /= (1. + polZZ * lAsymmetry + (po    201       factor /= (1. + polZZ * lAsymmetry + (polXX + polYY) * tAsymmetry);
202                                                   202 
203       if(verboseLevel >= 2)                       203       if(verboseLevel >= 2)
204       {                                           204       {
205         G4cout << " Asymmetry:     " << lAsymm    205         G4cout << " Asymmetry:     " << lAsymmetry << ", " << tAsymmetry
206                << G4endl;                         206                << G4endl;
207         G4cout << " PolProduct:    " << polXX     207         G4cout << " PolProduct:    " << polXX << ", " << polYY << ", " << polZZ
208                << G4endl;                         208                << G4endl;
209         G4cout << " Factor:        " << factor    209         G4cout << " Factor:        " << factor << G4endl;
210       }                                           210       }
211     }                                             211     }
212     else                                          212     else
213     {                                             213     {
214       G4ExceptionDescription ed;                  214       G4ExceptionDescription ed;
215       ed << "Problem with asymmetry tables: ma    215       ed << "Problem with asymmetry tables: material index " << midx
216          << " is out of range or tables are no    216          << " is out of range or tables are not filled";
217       G4Exception("G4PolarizedAnnihilation::Co    217       G4Exception("G4PolarizedAnnihilation::ComputeSaturationFactor", "em0048",
218                   JustWarning, ed, "");           218                   JustWarning, ed, "");
219     }                                             219     }
220   }                                               220   }
221   return factor;                                  221   return factor;
222 }                                                 222 }
223                                                   223 
224 //....oooOO0OOooo........oooOO0OOooo........oo    224 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
225 void G4PolarizedAnnihilation::BuildPhysicsTabl    225 void G4PolarizedAnnihilation::BuildPhysicsTable(
226   const G4ParticleDefinition& part)               226   const G4ParticleDefinition& part)
227 {                                                 227 {
228   G4VEmProcess::BuildPhysicsTable(part);          228   G4VEmProcess::BuildPhysicsTable(part);
229   if(isTheMaster)                                 229   if(isTheMaster)
230   {                                               230   {
231     BuildAsymmetryTables(part);                   231     BuildAsymmetryTables(part);
232   }                                               232   }
233 }                                                 233 }
234                                                   234 
235 //....oooOO0OOooo........oooOO0OOooo........oo    235 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
236 void G4PolarizedAnnihilation::BuildAsymmetryTa    236 void G4PolarizedAnnihilation::BuildAsymmetryTables(
237   const G4ParticleDefinition& part)               237   const G4ParticleDefinition& part)
238 {                                                 238 {
239   // cleanup old, initialise new table            239   // cleanup old, initialise new table
240   CleanTables();                                  240   CleanTables();
241   fAsymmetryTable = G4PhysicsTableHelper::Prep    241   fAsymmetryTable = G4PhysicsTableHelper::PreparePhysicsTable(fAsymmetryTable);
242   fTransverseAsymmetryTable =                     242   fTransverseAsymmetryTable =
243     G4PhysicsTableHelper::PreparePhysicsTable(    243     G4PhysicsTableHelper::PreparePhysicsTable(fTransverseAsymmetryTable);
244   if(nullptr == fAsymmetryTable) return;          244   if(nullptr == fAsymmetryTable) return;
245                                                   245 
246   // Access to materials                          246   // Access to materials
247   const G4ProductionCutsTable* theCoupleTable     247   const G4ProductionCutsTable* theCoupleTable =
248     G4ProductionCutsTable::GetProductionCutsTa    248     G4ProductionCutsTable::GetProductionCutsTable();
249   G4int numOfCouples = (G4int)theCoupleTable->    249   G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
250   for(G4int i = 0; i < numOfCouples; ++i)         250   for(G4int i = 0; i < numOfCouples; ++i)
251   {                                               251   {
252     if(fAsymmetryTable->GetFlag(i))               252     if(fAsymmetryTable->GetFlag(i))
253     {                                             253     {
254       // create physics vector and fill it        254       // create physics vector and fill it
255       const G4MaterialCutsCouple* couple =        255       const G4MaterialCutsCouple* couple =
256         theCoupleTable->GetMaterialCutsCouple(    256         theCoupleTable->GetMaterialCutsCouple(i);
257                                                   257 
258       // use same parameters as for lambda        258       // use same parameters as for lambda
259       G4PhysicsVector* aVector = LambdaPhysics    259       G4PhysicsVector* aVector = LambdaPhysicsVector(couple);
260       G4PhysicsVector* tVector = LambdaPhysics    260       G4PhysicsVector* tVector = LambdaPhysicsVector(couple);
261       G4int nn = (G4int)aVector->GetVectorLeng    261       G4int nn = (G4int)aVector->GetVectorLength();
262       for(G4int j = 0; j < nn; ++j)               262       for(G4int j = 0; j < nn; ++j)
263       {                                           263       {
264         G4double energy = aVector->Energy(j);     264         G4double energy = aVector->Energy(j);
265         G4double tasm   = 0.;                     265         G4double tasm   = 0.;
266         G4double asym = ComputeAsymmetry(energ    266         G4double asym = ComputeAsymmetry(energy, couple, part, 0., tasm);
267         aVector->PutValue(j, asym);               267         aVector->PutValue(j, asym);
268         tVector->PutValue(j, tasm);               268         tVector->PutValue(j, tasm);
269       }                                           269       }
270       if(aVector->GetSpline()) {                  270       if(aVector->GetSpline()) {
271   aVector->FillSecondDerivatives();               271   aVector->FillSecondDerivatives();
272   tVector->FillSecondDerivatives();               272   tVector->FillSecondDerivatives();
273       }                                           273       }
274       G4PhysicsTableHelper::SetPhysicsVector(f    274       G4PhysicsTableHelper::SetPhysicsVector(fAsymmetryTable, i, aVector);
275       G4PhysicsTableHelper::SetPhysicsVector(f    275       G4PhysicsTableHelper::SetPhysicsVector(fTransverseAsymmetryTable, i,
276                                              t    276                                              tVector);
277     }                                             277     }
278   }                                               278   }
279 }                                                 279 }
280                                                   280 
281 //....oooOO0OOooo........oooOO0OOooo........oo    281 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
282 G4double G4PolarizedAnnihilation::ComputeAsymm    282 G4double G4PolarizedAnnihilation::ComputeAsymmetry(
283   G4double energy, const G4MaterialCutsCouple*    283   G4double energy, const G4MaterialCutsCouple* couple,
284   const G4ParticleDefinition& aParticle, G4dou    284   const G4ParticleDefinition& aParticle, G4double cut, G4double& tAsymmetry)
285 {                                                 285 {
286   G4double lAsymmetry = 0.0;                      286   G4double lAsymmetry = 0.0;
287   tAsymmetry          = 0.0;                      287   tAsymmetry          = 0.0;
288                                                   288 
289   // calculate polarized cross section            289   // calculate polarized cross section
290   G4ThreeVector targetPolarization = G4ThreeVe    290   G4ThreeVector targetPolarization = G4ThreeVector(0., 0., 1.);
291   fEmModel->SetTargetPolarization(targetPolari    291   fEmModel->SetTargetPolarization(targetPolarization);
292   fEmModel->SetBeamPolarization(targetPolariza    292   fEmModel->SetBeamPolarization(targetPolarization);
293   G4double sigma2 =                               293   G4double sigma2 =
294     fEmModel->CrossSection(couple, &aParticle,    294     fEmModel->CrossSection(couple, &aParticle, energy, cut, energy);
295                                                   295 
296   // calculate transversely polarized cross se    296   // calculate transversely polarized cross section
297   targetPolarization = G4ThreeVector(1., 0., 0    297   targetPolarization = G4ThreeVector(1., 0., 0.);
298   fEmModel->SetTargetPolarization(targetPolari    298   fEmModel->SetTargetPolarization(targetPolarization);
299   fEmModel->SetBeamPolarization(targetPolariza    299   fEmModel->SetBeamPolarization(targetPolarization);
300   G4double sigma3 =                               300   G4double sigma3 =
301     fEmModel->CrossSection(couple, &aParticle,    301     fEmModel->CrossSection(couple, &aParticle, energy, cut, energy);
302                                                   302 
303   // calculate unpolarized cross section          303   // calculate unpolarized cross section
304   targetPolarization = G4ThreeVector();           304   targetPolarization = G4ThreeVector();
305   fEmModel->SetTargetPolarization(targetPolari    305   fEmModel->SetTargetPolarization(targetPolarization);
306   fEmModel->SetBeamPolarization(targetPolariza    306   fEmModel->SetBeamPolarization(targetPolarization);
307   G4double sigma0 =                               307   G4double sigma0 =
308     fEmModel->CrossSection(couple, &aParticle,    308     fEmModel->CrossSection(couple, &aParticle, energy, cut, energy);
309                                                   309 
310   // determine asymmetries                        310   // determine asymmetries
311   if(sigma0 > 0.)                                 311   if(sigma0 > 0.)
312   {                                               312   {
313     lAsymmetry = sigma2 / sigma0 - 1.;            313     lAsymmetry = sigma2 / sigma0 - 1.;
314     tAsymmetry = sigma3 / sigma0 - 1.;            314     tAsymmetry = sigma3 / sigma0 - 1.;
315   }                                               315   }
316   return lAsymmetry;                              316   return lAsymmetry;
317 }                                                 317 }
318                                                   318 
319 //....oooOO0OOooo........oooOO0OOooo........oo    319 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
320 void G4PolarizedAnnihilation::ProcessDescripti    320 void G4PolarizedAnnihilation::ProcessDescription(std::ostream& out) const
321 {                                                 321 {
322   out << "Polarized model for positron annihil    322   out << "Polarized model for positron annihilation into 2 photons.\n";
323 }                                                 323 }
324                                                   324