Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/polarisation/src/G4PolarizedCompton.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/G4PolarizedCompton.cc (Version 11.3.0) and /processes/electromagnetic/polarisation/src/G4PolarizedCompton.cc (Version 10.2.p3)


  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 //
                                                   >>  27 // $Id: G4PolarizedCompton.cc 93113 2015-10-07 07:49:04Z gcosmo $
                                                   >>  28 // 
                                                   >>  29 //
 26 // File name:     G4PolarizedCompton               30 // File name:     G4PolarizedCompton
 27 //                                                 31 //
 28 // Author:        Andreas Schaelicke               32 // Author:        Andreas Schaelicke
 29 //                based on code by Michel Mair     33 //                based on code by Michel Maire / Vladimir IVANTCHENKO
 30 //                                             << 
 31 // Class description                               34 // Class description
 32 //   modified version respecting media and bea <<  35 //
 33 //   using the stokes formalism                <<  36 // modified version respecting media and beam polarization
                                                   >>  37 //     using the stokes formalism
                                                   >>  38 //
                                                   >>  39 // Creation date: 01.05.2005
                                                   >>  40 //
                                                   >>  41 // Modifications:
                                                   >>  42 //
                                                   >>  43 // 01-01-05, include polarization description (A.Stahl)
                                                   >>  44 // 01-01-05, create asymmetry table and determine interactionlength (A.Stahl)
                                                   >>  45 // 01-05-05, update handling of media polarization (A.Schalicke)
                                                   >>  46 // 01-05-05, update polarized differential cross section (A.Schalicke)
                                                   >>  47 // 20-05-05, added polarization transfer (A.Schalicke)
                                                   >>  48 // 10-06-05, transformation between different reference frames (A.Schalicke)
                                                   >>  49 // 17-10-05, correct reference frame dependence in GetMeanFreePath (A.Schalicke)
                                                   >>  50 // 26-07-06, cross section recalculated (P.Starovoitov)
                                                   >>  51 // 09-08-06, make it work under current geant4 release (A.Schalicke)
                                                   >>  52 // 11-06-07, add PostStepGetPhysicalInteractionLength (A.Schalicke)
                                                   >>  53 // -----------------------------------------------------------------------------
 34                                                    54 
 35 #include "G4PolarizedCompton.hh"               << 
 36                                                    55 
                                                   >>  56 #include "G4PolarizedCompton.hh"
                                                   >>  57 #include "G4SystemOfUnits.hh"
 37 #include "G4Electron.hh"                           58 #include "G4Electron.hh"
 38 #include "G4EmParameters.hh"                   <<  59 
 39 #include "G4KleinNishinaCompton.hh"            <<  60 #include "G4StokesVector.hh"
 40 #include "G4PhysicsTableHelper.hh"             << 
 41 #include "G4PolarizationManager.hh"                61 #include "G4PolarizationManager.hh"
 42 #include "G4PolarizedComptonModel.hh"              62 #include "G4PolarizedComptonModel.hh"
 43 #include "G4ProductionCutsTable.hh"                63 #include "G4ProductionCutsTable.hh"
 44 #include "G4StokesVector.hh"                   <<  64 #include "G4PhysicsTableHelper.hh"
 45 #include "G4SystemOfUnits.hh"                  <<  65 #include "G4KleinNishinaCompton.hh"
                                                   >>  66 #include "G4PolarizedComptonModel.hh"
                                                   >>  67 #include "G4EmParameters.hh"
 46                                                    68 
 47 //....oooOO0OOooo........oooOO0OOooo........oo     69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >>  70 
 48 G4PhysicsTable* G4PolarizedCompton::theAsymmet     71 G4PhysicsTable* G4PolarizedCompton::theAsymmetryTable = nullptr;
 49                                                    72 
 50 G4PolarizedCompton::G4PolarizedCompton(const G     73 G4PolarizedCompton::G4PolarizedCompton(const G4String& processName,
 51                                        G4Proce <<  74   G4ProcessType type):
 52   : G4VEmProcess(processName, type)            <<  75   G4VEmProcess (processName, type),
 53   , fType(10)                                  <<  76   buildAsymmetryTable(true),
 54   , fBuildAsymmetryTable(true)                 <<  77   useAsymmetryTable(true),
 55   , fUseAsymmetryTable(true)                   <<  78   isInitialised(false),
 56   , fIsInitialised(false)                      <<  79   mType(10),
                                                   >>  80   targetPolarization(0.0,0.0,0.0)
 57 {                                                  81 {
 58   SetStartFromNullFlag(true);                      82   SetStartFromNullFlag(true);
 59   SetBuildTableFlag(true);                         83   SetBuildTableFlag(true);
 60   SetSecondaryParticle(G4Electron::Electron())     84   SetSecondaryParticle(G4Electron::Electron());
 61   SetProcessSubType(fComptonScattering);           85   SetProcessSubType(fComptonScattering);
 62   SetMinKinEnergyPrim(1. * MeV);               <<  86   SetMinKinEnergyPrim(1*MeV);
 63   SetSplineFlag(true);                             87   SetSplineFlag(true);
 64   fEmModel = nullptr;                          <<  88   emModel = nullptr;
 65 }                                                  89 }
 66                                                    90 
 67 //....oooOO0OOooo........oooOO0OOooo........oo     91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 68 G4PolarizedCompton::~G4PolarizedCompton() { Cl <<  92  
 69                                                <<  93 G4PolarizedCompton::~G4PolarizedCompton()
 70 //....oooOO0OOooo........oooOO0OOooo........oo << 
 71 void G4PolarizedCompton::ProcessDescription(st << 
 72 {                                                  94 {
 73   out << "Polarized model for Compton scatteri <<  95   CleanTable();
 74                                                << 
 75   G4VEmProcess::ProcessDescription(out);       << 
 76 }                                                  96 }
 77                                                    97 
 78 //....oooOO0OOooo........oooOO0OOooo........oo     98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >>  99  
 79 void G4PolarizedCompton::CleanTable()             100 void G4PolarizedCompton::CleanTable()
 80 {                                                 101 {
 81   if(theAsymmetryTable)                        << 102   if( theAsymmetryTable) {
 82   {                                            << 
 83     theAsymmetryTable->clearAndDestroy();         103     theAsymmetryTable->clearAndDestroy();
 84     delete theAsymmetryTable;                     104     delete theAsymmetryTable;
 85     theAsymmetryTable = nullptr;                  105     theAsymmetryTable = nullptr;
 86   }                                               106   }
 87 }                                                 107 }
 88                                                   108 
 89 //....oooOO0OOooo........oooOO0OOooo........oo    109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 110 
 90 G4bool G4PolarizedCompton::IsApplicable(const     111 G4bool G4PolarizedCompton::IsApplicable(const G4ParticleDefinition& p)
 91 {                                                 112 {
 92   return (&p == G4Gamma::Gamma());                113   return (&p == G4Gamma::Gamma());
 93 }                                                 114 }
 94                                                   115 
 95 //....oooOO0OOooo........oooOO0OOooo........oo    116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 117 
 96 void G4PolarizedCompton::InitialiseProcess(con    118 void G4PolarizedCompton::InitialiseProcess(const G4ParticleDefinition*)
 97 {                                                 119 {
 98   if(!fIsInitialised)                          << 120   if(!isInitialised) {
 99   {                                            << 121     isInitialised = true;
100     fIsInitialised = true;                     << 122     if(0 == mType) {
101     if(0 == fType)                             << 123       if(!EmModel(1)) { SetEmModel(new G4KleinNishinaCompton(), 1); }
102     {                                          << 124     } else {
103       if(nullptr == EmModel(0))                << 125       emModel = new G4PolarizedComptonModel();
104       {                                        << 126       SetEmModel(emModel, 1); 
105         SetEmModel(new G4KleinNishinaCompton() << 
106       }                                        << 
107     }                                          << 
108     else                                       << 
109     {                                          << 
110       fEmModel = new G4PolarizedComptonModel() << 
111       SetEmModel(fEmModel);                    << 
112     }                                             127     }
113     G4EmParameters* param = G4EmParameters::In    128     G4EmParameters* param = G4EmParameters::Instance();
114     EmModel(0)->SetLowEnergyLimit(param->MinKi << 129     EmModel(1)->SetLowEnergyLimit(param->MinKinEnergy());
115     EmModel(0)->SetHighEnergyLimit(param->MaxK << 130     EmModel(1)->SetHighEnergyLimit(param->MaxKinEnergy());
116     AddEmModel(1, EmModel(0));                 << 131     AddEmModel(1, EmModel(1));
117   }                                            << 132   } 
118 }                                                 133 }
119                                                   134 
120 //....oooOO0OOooo........oooOO0OOooo........oo    135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 136 
                                                   >> 137 void G4PolarizedCompton::PrintInfo()
                                                   >> 138 {
                                                   >> 139   G4cout << " Total cross sections has a good parametrisation"
                                                   >> 140          << " from 10 KeV to (100/Z) GeV" 
                                                   >> 141          << "\n      Sampling according " <<  EmModel(1)->GetName() << " model" 
                                                   >> 142    << G4endl;
                                                   >> 143 }         
                                                   >> 144 
                                                   >> 145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 146 
121 void G4PolarizedCompton::SetModel(const G4Stri    147 void G4PolarizedCompton::SetModel(const G4String& ss)
122 {                                                 148 {
123   if(ss == "Klein-Nishina")                    << 149   if(ss == "Klein-Nishina")     { mType = 0; }
124   {                                            << 150   if(ss == "Polarized-Compton") { mType = 10; }
125     fType = 0;                                 << 
126   }                                            << 
127   if(ss == "Polarized-Compton")                << 
128   {                                            << 
129     fType = 10;                                << 
130   }                                            << 
131 }                                                 151 }
132                                                   152 
133 //....oooOO0OOooo........oooOO0OOooo........oo    153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 154 
134 G4double G4PolarizedCompton::GetMeanFreePath(c    155 G4double G4PolarizedCompton::GetMeanFreePath(const G4Track& aTrack,
135                                              G << 156                G4double   previousStepSize,
136                                              G << 157                G4ForceCondition* condition)
137 {                                                 158 {
138   // *** get unploarised mean free path from l    159   // *** get unploarised mean free path from lambda table ***
139   G4double mfp =                               << 160   G4double mfp = G4VEmProcess::GetMeanFreePath(aTrack, previousStepSize, condition);
140     G4VEmProcess::GetMeanFreePath(aTrack, prev << 
141                                                   161 
142   if(theAsymmetryTable && fUseAsymmetryTable & << 162   if (theAsymmetryTable && useAsymmetryTable && mfp < DBL_MAX) {
143   {                                            << 
144     mfp *= ComputeSaturationFactor(aTrack);       163     mfp *= ComputeSaturationFactor(aTrack);
145   }                                               164   }
146   if(verboseLevel >= 2)                        << 165   if (verboseLevel>=2) {
147   {                                            << 166     G4cout << "G4PolarizedCompton::MeanFreePath:  " << mfp / mm << " mm " << G4endl;
148     G4cout << "G4PolarizedCompton::MeanFreePat << 
149            << G4endl;                          << 
150   }                                               167   }
151   return mfp;                                     168   return mfp;
152 }                                                 169 }
153                                                   170 
154 //....oooOO0OOooo........oooOO0OOooo........oo    171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 172 
155 G4double G4PolarizedCompton::PostStepGetPhysic    173 G4double G4PolarizedCompton::PostStepGetPhysicalInteractionLength(
156   const G4Track& aTrack, G4double previousStep << 174            const G4Track& aTrack,
                                                   >> 175            G4double   previousStepSize,
                                                   >> 176            G4ForceCondition* condition)
157 {                                                 177 {
158   // save previous values                      << 178   // save previous value
159   G4double nLength = theNumberOfInteractionLen    179   G4double nLength = theNumberOfInteractionLengthLeft;
160   G4double iLength = currentInteractionLength; << 
161                                                   180 
162   // *** compute unpolarized step limit ***    << 181   // *** compute uppolarized step limit ***
163   // this changes theNumberOfInteractionLength << 182   G4double x = G4VEmProcess::PostStepGetPhysicalInteractionLength(aTrack, 
164   G4double x = G4VEmProcess::PostStepGetPhysic << 183                   previousStepSize, 
165     aTrack, previousStepSize, condition);      << 184                   condition);
166   G4double x0      = x;                        << 
167   G4double satFact = 1.0;                      << 
168                                                   185 
169   // *** add corrections on polarisation ***      186   // *** add corrections on polarisation ***
170   if(theAsymmetryTable && fUseAsymmetryTable & << 187   if (theAsymmetryTable && useAsymmetryTable && x < DBL_MAX) {
171   {                                            << 188     G4double curLength = currentInteractionLength*ComputeSaturationFactor(aTrack);
172     satFact            = ComputeSaturationFact << 189     if(nLength > 0.0) {
173     G4double curLength = currentInteractionLen << 190       theNumberOfInteractionLengthLeft = 
174     G4double prvLength = iLength * satFact;    << 191   std::max(nLength - previousStepSize/curLength, 0.0);
175     if(nLength > 0.0)                          << 
176     {                                          << 
177       theNumberOfInteractionLengthLeft =       << 
178         std::max(nLength - previousStepSize /  << 
179     }                                             192     }
180     x = theNumberOfInteractionLengthLeft * cur    193     x = theNumberOfInteractionLengthLeft * curLength;
181   }                                               194   }
182   if(verboseLevel >= 2)                        << 195   if (verboseLevel>=2) {
183   {                                            << 196     G4cout << "G4PolarizedCompton::PostStepGetPhysicalInteractionLength:  " 
184     G4cout << "G4PolarizedCompton::PostStepGPI << 197      << x/mm << " mm " << G4endl;
185            << x / mm << " mm;" << G4endl       << 
186            << "               unpolarized valu << 
187            << x0 / mm << " mm." << G4endl;     << 
188   }                                               198   }
189   return x;                                       199   return x;
190 }                                                 200 }
191                                                   201 
192 //....oooOO0OOooo........oooOO0OOooo........oo    202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 203 
193 G4double G4PolarizedCompton::ComputeSaturation    204 G4double G4PolarizedCompton::ComputeSaturationFactor(const G4Track& aTrack)
194 {                                                 205 {
195   G4double factor = 1.0;                          206   G4double factor = 1.0;
196                                                   207 
197   // *** get asymmetry, if target is polarized    208   // *** get asymmetry, if target is polarized ***
198   const G4DynamicParticle* aDynamicGamma = aTr    209   const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
199   const G4double GammaEnergy             = aDy << 210   const G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
200   const G4StokesVector GammaPolarization =     << 211   const G4StokesVector GammaPolarization = aTrack.GetPolarization();
201     G4StokesVector(aTrack.GetPolarization());  << 212   const G4ParticleMomentum GammaDirection0 = aDynamicGamma->GetMomentumDirection();
202   const G4ParticleMomentum GammaDirection0 =   << 213 
203     aDynamicGamma->GetMomentumDirection();     << 214   G4Material*         aMaterial = aTrack.GetMaterial();
204                                                << 215   G4VPhysicalVolume*  aPVolume  = aTrack.GetVolume();
205   const G4Material* aMaterial = aTrack.GetMate << 216   G4LogicalVolume*    aLVolume  = aPVolume->GetLogicalVolume();
206   G4VPhysicalVolume* aPVolume = aTrack.GetVolu << 217 
207   G4LogicalVolume* aLVolume   = aPVolume->GetL << 218   //   G4Material* bMaterial = aLVolume->GetMaterial();
208                                                << 219   G4PolarizationManager * polarizationManger = G4PolarizationManager::GetInstance();
209   G4PolarizationManager* polarizationManager = << 220 
210     G4PolarizationManager::GetInstance();      << 221   const G4bool VolumeIsPolarized = polarizationManger->IsPolarized(aLVolume);
211                                                << 222   G4StokesVector ElectronPolarization = polarizationManger->GetVolumePolarization(aLVolume);
212   const G4bool VolumeIsPolarized = polarizatio << 223 
213   G4StokesVector ElectronPolarization =        << 224   if (VolumeIsPolarized) {
214     polarizationManager->GetVolumePolarization << 225      
215                                                << 226     if (verboseLevel>=2) {
216   if(VolumeIsPolarized)                        << 
217   {                                            << 
218     if(verboseLevel >= 2)                      << 
219     {                                          << 
220       G4cout << "G4PolarizedCompton::ComputeSa    227       G4cout << "G4PolarizedCompton::ComputeSaturationFactor: " << G4endl;
221       G4cout << " Mom " << GammaDirection0 <<  << 228       G4cout << " Mom " << GammaDirection0  << G4endl;
222       G4cout << " Polarization " << GammaPolar << 229       G4cout << " Polarization " << GammaPolarization  << G4endl;
223       G4cout << " MaterialPol. " << ElectronPo << 230       G4cout << " MaterialPol. " << ElectronPolarization  << G4endl;
224       G4cout << " Phys. Volume " << aPVolume->    231       G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl;
225       G4cout << " Log. Volume  " << aLVolume->    232       G4cout << " Log. Volume  " << aLVolume->GetName() << G4endl;
226       G4cout << " Material     " << aMaterial  << 233       G4cout << " Material     " << aMaterial          << G4endl;
227     }                                             234     }
228                                                   235 
229     std::size_t midx               = CurrentMa << 236     size_t midx = CurrentMaterialCutsCoupleIndex();
230     const G4PhysicsVector* aVector = nullptr;     237     const G4PhysicsVector* aVector = nullptr;
231     if(midx < theAsymmetryTable->size())       << 238     if(midx < theAsymmetryTable->size()) { 
232     {                                          << 
233       aVector = (*theAsymmetryTable)(midx);       239       aVector = (*theAsymmetryTable)(midx);
234     }                                             240     }
235     if(aVector)                                << 241     if (aVector) {
236     {                                          << 
237       G4double asymmetry = aVector->Value(Gamm    242       G4double asymmetry = aVector->Value(GammaEnergy);
238                                                   243 
239       //  we have to determine angle between p << 244       //  we have to determine angle between particle motion 
240       //  and target polarisation here         << 245       //  and target polarisation here  
241       //      circ pol * Vec(ElectronPol)*Vec(    246       //      circ pol * Vec(ElectronPol)*Vec(PhotonMomentum)
242       //  both vectors in global reference fra    247       //  both vectors in global reference frame
243                                                << 248      
244       G4double pol        = ElectronPolarizati << 249       G4double pol = ElectronPolarization*GammaDirection0;     
245       G4double polProduct = GammaPolarization.    250       G4double polProduct = GammaPolarization.p3() * pol;
246       factor /= (1. + polProduct * asymmetry);    251       factor /= (1. + polProduct * asymmetry);
247       if(verboseLevel >= 2)                    << 252       if (verboseLevel>=2) {
248       {                                        << 253   G4cout << " Asymmetry:     " << asymmetry      << G4endl;
249         G4cout << " Asymmetry:     " << asymme << 254   G4cout << " PolProduct:    " << polProduct     << G4endl;
250         G4cout << " PolProduct:    " << polPro << 255   G4cout << " Factor:        " << factor         << G4endl;
251         G4cout << " Factor:        " << factor << 256       }   
252       }                                        << 257     } else {
253     }                                          << 
254     else                                       << 
255     {                                          << 
256       G4ExceptionDescription ed;                  258       G4ExceptionDescription ed;
257       ed << "Problem with asymmetry table: mat << 259       ed << "Problem with asymmetry table: material index " << midx 
258          << " is out of range or the table is  << 260    << " is out of range or the table is not filled";
259       G4Exception("G4PolarizedComptonModel::Co << 261       G4Exception("G4PolarizedComptonModel::ComputeSaturationFactor","em0048",
260                   JustWarning, ed, "");        << 262       JustWarning, ed, "");
261     }                                             263     }
262   }                                               264   }
263   return factor;                                  265   return factor;
264 }                                                 266 }
265                                                   267 
266 //....oooOO0OOooo........oooOO0OOooo........oo    268 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 269 
267 void G4PolarizedCompton::BuildPhysicsTable(con    270 void G4PolarizedCompton::BuildPhysicsTable(const G4ParticleDefinition& part)
268 {                                                 271 {
269   // *** build (unpolarized) cross section tab    272   // *** build (unpolarized) cross section tables (Lambda)
270   G4VEmProcess::BuildPhysicsTable(part);          273   G4VEmProcess::BuildPhysicsTable(part);
271   if(fBuildAsymmetryTable && fEmModel)         << 274   if(buildAsymmetryTable && emModel) { 
272   {                                            << 
273     G4bool isMaster = true;                       275     G4bool isMaster = true;
274     const G4PolarizedCompton* masterProcess =  << 276     const G4PolarizedCompton* masterProcess = 
275       static_cast<const G4PolarizedCompton*>(G    277       static_cast<const G4PolarizedCompton*>(GetMasterProcess());
276     if(masterProcess && masterProcess != this) << 278     if(masterProcess && masterProcess != this) { isMaster = false; }
277     {                                          << 279     if(isMaster) { BuildAsymmetryTable(part); }
278       isMaster = false;                        << 
279     }                                          << 
280     if(isMaster)                               << 
281     {                                          << 
282       BuildAsymmetryTable(part);               << 
283     }                                          << 
284   }                                               280   }
285 }                                                 281 }
286                                                   282 
287 //....oooOO0OOooo........oooOO0OOooo........oo    283 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 284 
288 void G4PolarizedCompton::BuildAsymmetryTable(c    285 void G4PolarizedCompton::BuildAsymmetryTable(const G4ParticleDefinition& part)
289 {                                                 286 {
290   // cleanup old, initialise new table            287   // cleanup old, initialise new table
291   CleanTable();                                   288   CleanTable();
292   theAsymmetryTable =                          << 289   theAsymmetryTable = 
293     G4PhysicsTableHelper::PreparePhysicsTable(    290     G4PhysicsTableHelper::PreparePhysicsTable(theAsymmetryTable);
294                                                   291 
295   // Access to materials                          292   // Access to materials
296   const G4ProductionCutsTable* theCoupleTable  << 293   const G4ProductionCutsTable* theCoupleTable=
297     G4ProductionCutsTable::GetProductionCutsTa << 294         G4ProductionCutsTable::GetProductionCutsTable();
298   G4int numOfCouples = (G4int)theCoupleTable-> << 295   size_t numOfCouples = theCoupleTable->GetTableSize();
299   if(!theAsymmetryTable)                       << 296   if(!theAsymmetryTable) { return; }
300   {                                            << 297   G4int nbins = LambdaBinning();
301     return;                                    << 298   G4double emin = MinKinEnergy();
302   }                                            << 299   G4double emax = MaxKinEnergy();
303   G4int nbins                 = LambdaBinning( << 300   G4PhysicsLogVector* aVector = 0;
304   G4double emin               = MinKinEnergy() << 301   G4PhysicsLogVector* bVector = 0;
305   G4double emax               = MaxKinEnergy() << 302 
306   G4PhysicsLogVector* aVector = nullptr;       << 303   for(size_t i=0; i<numOfCouples; ++i) {
307   G4PhysicsLogVector* bVector = nullptr;       << 304     if (theAsymmetryTable->GetFlag(i)) {
308                                                << 305 
309   for(G4int i = 0; i < numOfCouples; ++i)      << 
310   {                                            << 
311     if(theAsymmetryTable->GetFlag(i))          << 
312     {                                          << 
313       // create physics vector and fill it        306       // create physics vector and fill it
314       const G4MaterialCutsCouple* couple =     << 307       const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
315         theCoupleTable->GetMaterialCutsCouple( << 
316       // use same parameters as for lambda        308       // use same parameters as for lambda
317       if(!aVector)                             << 309       if(!aVector) { 
318       {                                        << 310   aVector = new G4PhysicsLogVector(emin, emax, nbins); 
319         aVector = new G4PhysicsLogVector(emin, << 311         aVector->SetSpline(true);
320         bVector = aVector;                        312         bVector = aVector;
321       }                                        << 313       } else {
322       else                                     << 314   bVector = new G4PhysicsLogVector(*aVector);
323       {                                        << 
324         bVector = new G4PhysicsLogVector(*aVec << 
325       }                                           315       }
326                                                   316 
327       for(G4int j = 0; j <= nbins; ++j)        << 317       for (G4int j = 0; j <= nbins; ++j ) {
328       {                                        << 318   G4double energy = bVector->Energy(j);
329         G4double energy = bVector->Energy(j);  << 319   G4double tasm=0.;
330         G4double tasm   = 0.;                  << 320   G4double asym = ComputeAsymmetry(energy, couple, part, 0., tasm);
331         G4double asym   = ComputeAsymmetry(ene << 321   bVector->PutValue(j,asym);
332         bVector->PutValue(j, asym);            << 
333       }                                           322       }
334       bVector->FillSecondDerivatives();        << 
335       G4PhysicsTableHelper::SetPhysicsVector(t    323       G4PhysicsTableHelper::SetPhysicsVector(theAsymmetryTable, i, bVector);
336     }                                             324     }
337   }                                               325   }
338 }                                                 326 }
339                                                   327 
340 //....oooOO0OOooo........oooOO0OOooo........oo    328 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
341 G4double G4PolarizedCompton::ComputeAsymmetry( << 329 
342   G4double energy, const G4MaterialCutsCouple* << 330 G4double G4PolarizedCompton::ComputeAsymmetry(G4double energy,
343   const G4ParticleDefinition& aParticle, G4dou << 331                 const G4MaterialCutsCouple* couple,
                                                   >> 332                 const G4ParticleDefinition& aParticle,
                                                   >> 333                 G4double cut,
                                                   >> 334                 G4double & tAsymmetry)
344 {                                                 335 {
345   G4double lAsymmetry = 0.0;                      336   G4double lAsymmetry = 0.0;
346   tAsymmetry          = 0;                     << 337   tAsymmetry=0;
347                                                   338 
                                                   >> 339   //
348   // calculate polarized cross section            340   // calculate polarized cross section
349   G4ThreeVector thePolarization = G4ThreeVecto << 341   //
350   fEmModel->SetTargetPolarization(thePolarizat << 342   G4ThreeVector thePolarization=G4ThreeVector(0.,0.,1.);
351   fEmModel->SetBeamPolarization(thePolarizatio << 343   emModel->SetTargetPolarization(thePolarization);
352   G4double sigma2 =                            << 344   emModel->SetBeamPolarization(thePolarization);
353     fEmModel->CrossSection(couple, &aParticle, << 345   G4double sigma2=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
354                                                   346 
                                                   >> 347   //
355   // calculate unpolarized cross section          348   // calculate unpolarized cross section
356   thePolarization = G4ThreeVector();           << 349   //
357   fEmModel->SetTargetPolarization(thePolarizat << 350   thePolarization=G4ThreeVector();
358   fEmModel->SetBeamPolarization(thePolarizatio << 351   emModel->SetTargetPolarization(thePolarization);
359   G4double sigma0 =                            << 352   emModel->SetBeamPolarization(thePolarization);
360     fEmModel->CrossSection(couple, &aParticle, << 353   G4double sigma0=emModel->CrossSection(couple,&aParticle,energy,cut,energy);
361                                                << 354 
362   // determine asymmetries                     << 355   // determine assymmetries
363   if(sigma0 > 0.)                              << 356   if (sigma0 > 0.) {
364   {                                            << 357     lAsymmetry = sigma2/sigma0-1.;
365     lAsymmetry = sigma2 / sigma0 - 1.;         << 
366   }                                               358   }
367   return lAsymmetry;                              359   return lAsymmetry;
368 }                                                 360 }
                                                   >> 361 
                                                   >> 362 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
369                                                   363