Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4LivermorePolarizedComptonModel.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/lowenergy/src/G4LivermorePolarizedComptonModel.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4LivermorePolarizedComptonModel.cc (Version 11.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 // Authors: G.Depaola & F.Longo                    27 // Authors: G.Depaola & F.Longo
 28 //                                                 28 //
 29 // History:                                        29 // History:
 30 // -------                                         30 // -------
 31 //                                                 31 //
 32 // 05 Apr 2021   J Allison added quantum entan     32 // 05 Apr 2021   J Allison added quantum entanglement of e+ annihilation.
 33 // If the photons have been "tagged" as "quant     33 // If the photons have been "tagged" as "quantum-entangled", for example by
 34 // G4eplusAnnihilation for annihilation into 2     34 // G4eplusAnnihilation for annihilation into 2 photons, they are "analysed"
 35 // here if - and only if - both photons suffer     35 // here if - and only if - both photons suffer Compton scattering. Theoretical
 36 // predictions from Pryce and Ward, Nature No      36 // predictions from Pryce and Ward, Nature No 4065 (1947) p.435, and Snyder et al,
 37 // Physical Review 73 (1948) p.440. Experiment     37 // Physical Review 73 (1948) p.440. Experimental validation in "Photon quantum 
 38 // entanglement in the MeV regime and its appl     38 // entanglement in the MeV regime and its application in PET imaging",
 39 // D. Watts, J. Allison et al., Nature Communi     39 // D. Watts, J. Allison et al., Nature Communications (2021)12:2646
 40 // https://doi.org/10.1038/s41467-021-22907-5.     40 // https://doi.org/10.1038/s41467-021-22907-5.
 41 //                                                 41 //
 42 // 02 May 2009   S Incerti as V. Ivanchenko pr     42 // 02 May 2009   S Incerti as V. Ivanchenko proposed in G4LivermoreComptonModel.cc
 43 //                                                 43 //
 44 // Cleanup initialisation and generation of se     44 // Cleanup initialisation and generation of secondaries:
 45 //                  - apply internal high-ener     45 //                  - apply internal high-energy limit only in constructor 
 46 //                  - do not apply low-energy      46 //                  - do not apply low-energy limit (default is 0)
 47 //                  - remove GetMeanFreePath m     47 //                  - remove GetMeanFreePath method and table
 48 //                  - added protection against     48 //                  - added protection against numerical problem in energy sampling 
 49 //                  - use G4ElementSelector        49 //                  - use G4ElementSelector
 50                                                    50 
 51 #include "G4LivermorePolarizedComptonModel.hh"     51 #include "G4LivermorePolarizedComptonModel.hh"
 52 #include "G4PhysicalConstants.hh"                  52 #include "G4PhysicalConstants.hh"
 53 #include "G4SystemOfUnits.hh"                      53 #include "G4SystemOfUnits.hh"
 54 #include "G4AutoLock.hh"                           54 #include "G4AutoLock.hh"
 55 #include "G4Electron.hh"                           55 #include "G4Electron.hh"
 56 #include "G4ParticleChangeForGamma.hh"             56 #include "G4ParticleChangeForGamma.hh"
 57 #include "G4LossTableManager.hh"                   57 #include "G4LossTableManager.hh"
 58 #include "G4VAtomDeexcitation.hh"                  58 #include "G4VAtomDeexcitation.hh"
 59 #include "G4AtomicShell.hh"                        59 #include "G4AtomicShell.hh"
 60 #include "G4Gamma.hh"                              60 #include "G4Gamma.hh"
 61 #include "G4ShellData.hh"                          61 #include "G4ShellData.hh"
 62 #include "G4DopplerProfile.hh"                     62 #include "G4DopplerProfile.hh"
 63 #include "G4Log.hh"                                63 #include "G4Log.hh"
 64 #include "G4Exp.hh"                                64 #include "G4Exp.hh"
 65 #include "G4Pow.hh"                                65 #include "G4Pow.hh"
 66 #include "G4LogLogInterpolation.hh"                66 #include "G4LogLogInterpolation.hh"
 67 #include "G4PhysicsModelCatalog.hh"                67 #include "G4PhysicsModelCatalog.hh"
 68 #include "G4EntanglementAuxInfo.hh"                68 #include "G4EntanglementAuxInfo.hh"
 69 #include "G4eplusAnnihilationEntanglementClipB     69 #include "G4eplusAnnihilationEntanglementClipBoard.hh"
 70                                                    70 
 71 //....oooOO0OOooo........oooOO0OOooo........oo     71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 72                                                    72 
 73 using namespace std;                               73 using namespace std;
 74 namespace { G4Mutex LivermorePolarizedComptonM     74 namespace { G4Mutex LivermorePolarizedComptonModelMutex = G4MUTEX_INITIALIZER; }
 75                                                    75 
 76                                                    76 
 77 G4PhysicsFreeVector* G4LivermorePolarizedCompt     77 G4PhysicsFreeVector* G4LivermorePolarizedComptonModel::data[] = {nullptr};
 78 G4ShellData*       G4LivermorePolarizedCompton     78 G4ShellData*       G4LivermorePolarizedComptonModel::shellData = nullptr;
 79 G4DopplerProfile*  G4LivermorePolarizedCompton     79 G4DopplerProfile*  G4LivermorePolarizedComptonModel::profileData = nullptr; 
 80 G4CompositeEMDataSet* G4LivermorePolarizedComp     80 G4CompositeEMDataSet* G4LivermorePolarizedComptonModel::scatterFunctionData = nullptr;
 81                                                    81 
 82 //....oooOO0OOooo........oooOO0OOooo........oo     82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 83                                                    83 
 84 G4LivermorePolarizedComptonModel::G4LivermoreP     84 G4LivermorePolarizedComptonModel::G4LivermorePolarizedComptonModel(const G4ParticleDefinition*, const G4String& nam)
 85   :G4VEmModel(nam),isInitialised(false)            85   :G4VEmModel(nam),isInitialised(false)
 86 {                                                  86 { 
 87   verboseLevel= 1;                                 87   verboseLevel= 1;
 88   // Verbosity scale:                              88   // Verbosity scale:
 89   // 0 = nothing                                   89   // 0 = nothing 
 90   // 1 = warning for energy non-conservation       90   // 1 = warning for energy non-conservation 
 91   // 2 = details of energy budget                  91   // 2 = details of energy budget
 92   // 3 = calculation of cross sections, file o     92   // 3 = calculation of cross sections, file openings, sampling of atoms
 93   // 4 = entering in methods                       93   // 4 = entering in methods
 94                                                    94 
 95   if( verboseLevel>1 )                             95   if( verboseLevel>1 )  
 96     G4cout << "Livermore Polarized Compton is      96     G4cout << "Livermore Polarized Compton is constructed " << G4endl;
 97                                                    97         
 98   //Mark this model as "applicable" for atomic     98   //Mark this model as "applicable" for atomic deexcitation
 99   SetDeexcitationFlag(true);                       99   SetDeexcitationFlag(true);
100                                                   100   
101   fParticleChange = nullptr;                      101   fParticleChange = nullptr;
102   fAtomDeexcitation = nullptr;                    102   fAtomDeexcitation = nullptr;
103   fEntanglementModelID = G4PhysicsModelCatalog    103   fEntanglementModelID = G4PhysicsModelCatalog::GetModelID("model_GammaGammaEntanglement");
104 }                                                 104 }
105                                                   105 
106 //....oooOO0OOooo........oooOO0OOooo........oo    106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107                                                   107 
108 G4LivermorePolarizedComptonModel::~G4Livermore    108 G4LivermorePolarizedComptonModel::~G4LivermorePolarizedComptonModel()
109 {                                                 109 {  
110   if(IsMaster()) {                                110   if(IsMaster()) {
111     delete shellData;                             111     delete shellData;
112     shellData = nullptr;                          112     shellData = nullptr;
113     delete profileData;                           113     delete profileData;
114     profileData = nullptr;                        114     profileData = nullptr;
115     delete scatterFunctionData;                   115     delete scatterFunctionData;
116     scatterFunctionData = nullptr;                116     scatterFunctionData = nullptr;
117     for(G4int i=0; i<maxZ; ++i) {                 117     for(G4int i=0; i<maxZ; ++i) {
118       if(data[i]) {                               118       if(data[i]) { 
119   delete data[i];                                 119   delete data[i];
120   data[i] = nullptr;                              120   data[i] = nullptr;
121       }                                           121       }
122     }                                             122     }
123   }                                               123   }
124 }                                                 124 }
125                                                   125 
126 //....oooOO0OOooo........oooOO0OOooo........oo    126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
127                                                   127 
128 void G4LivermorePolarizedComptonModel::Initial    128 void G4LivermorePolarizedComptonModel::Initialise(const G4ParticleDefinition* particle,
129                                        const G    129                                        const G4DataVector& cuts)
130 {                                                 130 {
131   if (verboseLevel > 1)                           131   if (verboseLevel > 1)
132     G4cout << "Calling G4LivermorePolarizedCom    132     G4cout << "Calling G4LivermorePolarizedComptonModel::Initialise()" << G4endl;
133                                                   133 
134   // Initialise element selector                  134   // Initialise element selector 
135   if(IsMaster()) {                                135   if(IsMaster()) {
136     // Access to elements                         136     // Access to elements 
137     const char* path = G4FindDataDir("G4LEDATA    137     const char* path = G4FindDataDir("G4LEDATA");
138                                                   138 
139     G4ProductionCutsTable* theCoupleTable =       139     G4ProductionCutsTable* theCoupleTable = 
140       G4ProductionCutsTable::GetProductionCuts    140       G4ProductionCutsTable::GetProductionCutsTable();
141                                                   141 
142     G4int numOfCouples = (G4int)theCoupleTable    142     G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
143                                                   143     
144     for(G4int i=0; i<numOfCouples; ++i) {         144     for(G4int i=0; i<numOfCouples; ++i) {
145       const G4Material* material =                145       const G4Material* material = 
146         theCoupleTable->GetMaterialCutsCouple(    146         theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
147       const G4ElementVector* theElementVector     147       const G4ElementVector* theElementVector = material->GetElementVector();
148       std::size_t nelm = material->GetNumberOf    148       std::size_t nelm = material->GetNumberOfElements();
149                                                   149       
150       for (std::size_t j=0; j<nelm; ++j) {        150       for (std::size_t j=0; j<nelm; ++j) {
151         G4int Z = G4lrint((*theElementVector)[    151         G4int Z = G4lrint((*theElementVector)[j]->GetZ());
152         if(Z < 1)        { Z = 1; }               152         if(Z < 1)        { Z = 1; }
153         else if(Z > maxZ){ Z = maxZ; }            153         else if(Z > maxZ){ Z = maxZ; }
154                                                   154 
155         if( (!data[Z]) ) { ReadData(Z, path);     155         if( (!data[Z]) ) { ReadData(Z, path); }
156       }                                           156       }
157     }                                             157     }
158                                                   158     
159     // For Doppler broadening                     159     // For Doppler broadening
160     if(!shellData) {                              160     if(!shellData) {
161       shellData = new G4ShellData();              161       shellData = new G4ShellData(); 
162       shellData->SetOccupancyData();              162       shellData->SetOccupancyData();
163       G4String file = "/doppler/shell-doppler"    163       G4String file = "/doppler/shell-doppler";
164       shellData->LoadData(file);                  164       shellData->LoadData(file);
165     }                                             165     }
166     if(!profileData) { profileData = new G4Dop    166     if(!profileData) { profileData = new G4DopplerProfile(); }
167                                                   167 
168     // Scattering Function                        168     // Scattering Function 
169     if(!scatterFunctionData)                      169     if(!scatterFunctionData)
170       {                                           170       {
171                                                   171   
172   G4VDataSetAlgorithm* scatterInterpolation =     172   G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation;
173   G4String scatterFile = "comp/ce-sf-";           173   G4String scatterFile = "comp/ce-sf-";
174   scatterFunctionData = new G4CompositeEMDataS    174   scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.);
175   scatterFunctionData->LoadData(scatterFile);     175   scatterFunctionData->LoadData(scatterFile);
176       }                                           176       }
177                                                   177     
178     InitialiseElementSelectors(particle, cuts)    178     InitialiseElementSelectors(particle, cuts);
179   }                                               179   }
180                                                   180  
181   if (verboseLevel > 2) {                         181   if (verboseLevel > 2) {
182     G4cout << "Loaded cross section files" <<     182     G4cout << "Loaded cross section files" << G4endl;
183   }                                               183   }
184                                                   184   
185   if( verboseLevel>1 ) {                          185   if( verboseLevel>1 ) { 
186     G4cout << "G4LivermoreComptonModel is init    186     G4cout << "G4LivermoreComptonModel is initialized " << G4endl
187      << "Energy range: "                          187      << "Energy range: "
188      << LowEnergyLimit() / eV << " eV - "         188      << LowEnergyLimit() / eV << " eV - "
189      << HighEnergyLimit() / GeV << " GeV"         189      << HighEnergyLimit() / GeV << " GeV"
190      << G4endl;                                   190      << G4endl;
191   }                                               191   }
192   //                                              192   //  
193   if(isInitialised) { return; }                   193   if(isInitialised) { return; }
194                                                   194   
195   fParticleChange = GetParticleChangeForGamma(    195   fParticleChange = GetParticleChangeForGamma();
196   fAtomDeexcitation  = G4LossTableManager::Ins    196   fAtomDeexcitation  = G4LossTableManager::Instance()->AtomDeexcitation();
197   isInitialised = true;                           197   isInitialised = true;
198 }                                                 198 }
199                                                   199 
200                                                   200 
201 void G4LivermorePolarizedComptonModel::Initial    201 void G4LivermorePolarizedComptonModel::InitialiseLocal(const G4ParticleDefinition*,
202                 G4VEmModel* masterModel)          202                 G4VEmModel* masterModel)
203 {                                                 203 {
204   SetElementSelectors(masterModel->GetElementS    204   SetElementSelectors(masterModel->GetElementSelectors());
205 }                                                 205 }
206                                                   206 
207 //....oooOO0OOooo........oooOO0OOooo........oo    207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
208                                                   208 
209 void G4LivermorePolarizedComptonModel::ReadDat    209 void G4LivermorePolarizedComptonModel::ReadData(std::size_t Z, const char* path)
210 {                                                 210 {
211   if (verboseLevel > 1)                           211   if (verboseLevel > 1) 
212     {                                             212     {
213       G4cout << "G4LivermorePolarizedComptonMo    213       G4cout << "G4LivermorePolarizedComptonModel::ReadData()" 
214        << G4endl;                                 214        << G4endl;
215     }                                             215     }
216   if(data[Z]) { return; }                         216   if(data[Z]) { return; }  
217   const char* datadir = path;                     217   const char* datadir = path;
218   if(!datadir)                                    218   if(!datadir) 
219     {                                             219     {
220       datadir = G4FindDataDir("G4LEDATA");        220       datadir = G4FindDataDir("G4LEDATA");
221       if(!datadir)                                221       if(!datadir) 
222   {                                               222   {
223     G4Exception("G4LivermorePolarizedComptonMo    223     G4Exception("G4LivermorePolarizedComptonModel::ReadData()",
224           "em0006",FatalException,                224           "em0006",FatalException,
225           "Environment variable G4LEDATA not d    225           "Environment variable G4LEDATA not defined");
226     return;                                       226     return;
227   }                                               227   }
228     }                                             228     }
229                                                   229   
230   data[Z] = new G4PhysicsFreeVector();            230   data[Z] = new G4PhysicsFreeVector();
231                                                   231     
232   std::ostringstream ost;                         232   std::ostringstream ost;
233   ost << datadir << "/livermore/comp/ce-cs-" <    233   ost << datadir << "/livermore/comp/ce-cs-" << Z <<".dat";
234   std::ifstream fin(ost.str().c_str());           234   std::ifstream fin(ost.str().c_str());
235                                                   235   
236   if( !fin.is_open())                             236   if( !fin.is_open()) 
237     {                                             237     {
238       G4ExceptionDescription ed;                  238       G4ExceptionDescription ed;
239       ed << "G4LivermorePolarizedComptonModel     239       ed << "G4LivermorePolarizedComptonModel data file <" << ost.str().c_str()
240    << "> is not opened!" << G4endl;               240    << "> is not opened!" << G4endl;
241       G4Exception("G4LivermoreComptonModel::Re    241       G4Exception("G4LivermoreComptonModel::ReadData()",
242       "em0003",FatalException,                    242       "em0003",FatalException,
243       ed,"G4LEDATA version should be G4EMLOW8.    243       ed,"G4LEDATA version should be G4EMLOW8.0 or later");
244       return;                                     244       return;
245     } else {                                      245     } else {
246     if(verboseLevel > 3) {                        246     if(verboseLevel > 3) {
247       G4cout << "File " << ost.str()              247       G4cout << "File " << ost.str() 
248        << " is opened by G4LivermorePolarizedC    248        << " is opened by G4LivermorePolarizedComptonModel" << G4endl;
249     }                                             249     }   
250     data[Z]->Retrieve(fin, true);                 250     data[Z]->Retrieve(fin, true);
251     data[Z]->ScaleVector(MeV, MeV*barn);          251     data[Z]->ScaleVector(MeV, MeV*barn);
252   }                                               252   }   
253   fin.close();                                    253   fin.close();
254 }                                                 254 }
255                                                   255 
256 //....oooOO0OOooo........oooOO0OOooo........oo    256 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
257                                                   257 
258 G4double G4LivermorePolarizedComptonModel::Com    258 G4double G4LivermorePolarizedComptonModel::ComputeCrossSectionPerAtom(
259                                        const G    259                                        const G4ParticleDefinition*,
260                                              G    260                                              G4double GammaEnergy,
261                                              G    261                                              G4double Z, G4double,
262                                              G    262                                              G4double, G4double)
263 {                                                 263 {
264   if (verboseLevel > 3)                           264   if (verboseLevel > 3)
265     G4cout << "Calling ComputeCrossSectionPerA    265     G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePolarizedComptonModel" << G4endl;
266                                                   266 
267   G4double cs = 0.0;                              267   G4double cs = 0.0; 
268                                                   268   
269   if (GammaEnergy < LowEnergyLimit())             269   if (GammaEnergy < LowEnergyLimit()) 
270     return 0.0;                                   270     return 0.0;
271                                                   271 
272   G4int intZ = G4lrint(Z);                        272   G4int intZ = G4lrint(Z);
273   if(intZ < 1 || intZ > maxZ) { return cs; }      273   if(intZ < 1 || intZ > maxZ) { return cs; } 
274                                                   274   
275   G4PhysicsFreeVector* pv = data[intZ];           275   G4PhysicsFreeVector* pv = data[intZ];
276                                                   276   
277   // if element was not initialised               277   // if element was not initialised
278   // do initialisation safely for MT mode         278   // do initialisation safely for MT mode
279   if(!pv)                                         279   if(!pv) 
280     {                                             280     {
281       InitialiseForElement(0, intZ);              281       InitialiseForElement(0, intZ);
282       pv = data[intZ];                            282       pv = data[intZ];
283       if(!pv) { return cs; }                      283       if(!pv) { return cs; }
284     }                                             284     }
285                                                   285   
286   G4int n = G4int(pv->GetVectorLength() - 1);     286   G4int n = G4int(pv->GetVectorLength() - 1);   
287   G4double e1 = pv->Energy(0);                    287   G4double e1 = pv->Energy(0);
288   G4double e2 = pv->Energy(n);                    288   G4double e2 = pv->Energy(n);
289                                                   289   
290   if(GammaEnergy <= e1)      { cs = GammaEnerg    290   if(GammaEnergy <= e1)      { cs = GammaEnergy/(e1*e1)*pv->Value(e1); }
291   else if(GammaEnergy <= e2) { cs = pv->Value(    291   else if(GammaEnergy <= e2) { cs = pv->Value(GammaEnergy)/GammaEnergy; }
292   else if(GammaEnergy > e2)  { cs = pv->Value(    292   else if(GammaEnergy > e2)  { cs = pv->Value(e2)/GammaEnergy; }
293                                                   293   
294   return cs;                                      294   return cs;
295                                                   295 
296 }                                                 296 }
297                                                   297 
298 //....oooOO0OOooo........oooOO0OOooo........oo    298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
299                                                   299 
300 void G4LivermorePolarizedComptonModel::SampleS    300 void G4LivermorePolarizedComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
301                 const G4MaterialCutsCouple* co    301                 const G4MaterialCutsCouple* couple,
302                 const G4DynamicParticle* aDyna    302                 const G4DynamicParticle* aDynamicGamma,
303                 G4double,                         303                 G4double,
304                 G4double)                         304                 G4double)
305 {                                                 305 {
306   // The scattered gamma energy is sampled acc    306   // The scattered gamma energy is sampled according to Klein - Nishina formula.
307   // The random number techniques of Butcher &    307   // The random number techniques of Butcher & Messel are used (Nuc Phys 20(1960),15).
308   // GEANT4 internal units                        308   // GEANT4 internal units
309   //                                              309   //
310   // Note : Effects due to binding of atomic e    310   // Note : Effects due to binding of atomic electrons are negliged.
311                                                   311 
312   if (verboseLevel > 3)                           312   if (verboseLevel > 3)
313     G4cout << "Calling SampleSecondaries() of     313     G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedComptonModel" << G4endl;
314                                                   314 
315   G4double gammaEnergy0 = aDynamicGamma->GetKi    315   G4double gammaEnergy0 = aDynamicGamma->GetKineticEnergy();
316                                                   316  
317   // do nothing below the threshold               317   // do nothing below the threshold
318   // should never get here because the XS is z    318   // should never get here because the XS is zero below the limit
319   if (gammaEnergy0 < LowEnergyLimit())            319   if (gammaEnergy0 < LowEnergyLimit())     
320     return ;                                      320     return ; 
321                                                   321 
322   G4ThreeVector gammaPolarization0 = aDynamicG    322   G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
323                                                   323 
324   // Protection: a polarisation parallel to th    324   // Protection: a polarisation parallel to the
325   // direction causes problems;                   325   // direction causes problems;
326   // in that case find a random polarization      326   // in that case find a random polarization
327   G4ThreeVector gammaDirection0 = aDynamicGamm    327   G4ThreeVector gammaDirection0 = aDynamicGamma->GetMomentumDirection();
328                                                   328 
329   // Make sure that the polarization vector is    329   // Make sure that the polarization vector is perpendicular to the
330   // gamma direction. If not                      330   // gamma direction. If not
331   if(!(gammaPolarization0.isOrthogonal(gammaDi    331   if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
332     { // only for testing now                     332     { // only for testing now
333       gammaPolarization0 = GetRandomPolarizati    333       gammaPolarization0 = GetRandomPolarization(gammaDirection0);
334     }                                             334     }
335   else                                            335   else
336     {                                             336     {
337       if ( gammaPolarization0.howOrthogonal(ga    337       if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
338   {                                               338   {
339     gammaPolarization0 = GetPerpendicularPolar    339     gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
340   }                                               340   }
341     }                                             341     }
342   // End of Protection                            342   // End of Protection
343                                                   343 
344   G4double E0_m = gammaEnergy0 / electron_mass    344   G4double E0_m = gammaEnergy0 / electron_mass_c2 ;
345                                                   345 
346   // Select randomly one element in the curren    346   // Select randomly one element in the current material
347   //G4int Z = crossSectionHandler->SelectRando    347   //G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0);
348   const G4ParticleDefinition* particle =  aDyn    348   const G4ParticleDefinition* particle =  aDynamicGamma->GetDefinition();
349   const G4Element* elm = SelectRandomAtom(coup    349   const G4Element* elm = SelectRandomAtom(couple,particle,gammaEnergy0);
350   G4int Z = (G4int)elm->GetZ();                   350   G4int Z = (G4int)elm->GetZ();
351                                                   351 
352   // Sample the energy and the polarization of    352   // Sample the energy and the polarization of the scattered photon
353   G4double epsilon, epsilonSq, onecost, sinThe    353   G4double epsilon, epsilonSq, onecost, sinThetaSqr, greject ;
354                                                   354 
355   G4double epsilon0Local = 1./(1. + 2*E0_m);      355   G4double epsilon0Local = 1./(1. + 2*E0_m);
356   G4double epsilon0Sq = epsilon0Local*epsilon0    356   G4double epsilon0Sq = epsilon0Local*epsilon0Local;
357   G4double alpha1   = - G4Log(epsilon0Local);     357   G4double alpha1   = - G4Log(epsilon0Local);
358   G4double alpha2 = 0.5*(1.- epsilon0Sq);         358   G4double alpha2 = 0.5*(1.- epsilon0Sq);
359                                                   359 
360   G4double wlGamma = h_Planck*c_light/gammaEne    360   G4double wlGamma = h_Planck*c_light/gammaEnergy0;
361   G4double gammaEnergy1;                          361   G4double gammaEnergy1;
362   G4ThreeVector gammaDirection1;                  362   G4ThreeVector gammaDirection1;
363                                                   363 
364   do {                                            364   do {
365     if ( alpha1/(alpha1+alpha2) > G4UniformRan    365     if ( alpha1/(alpha1+alpha2) > G4UniformRand() )
366       {                                           366       {
367   epsilon   = G4Exp(-alpha1*G4UniformRand());     367   epsilon   = G4Exp(-alpha1*G4UniformRand());  
368   epsilonSq = epsilon*epsilon;                    368   epsilonSq = epsilon*epsilon; 
369       }                                           369       }
370     else                                          370     else 
371       {                                           371       {
372   epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4    372   epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4UniformRand();
373   epsilon   = std::sqrt(epsilonSq);               373   epsilon   = std::sqrt(epsilonSq);
374       }                                           374       }
375                                                   375 
376     onecost = (1.- epsilon)/(epsilon*E0_m);       376     onecost = (1.- epsilon)/(epsilon*E0_m);
377     sinThetaSqr   = onecost*(2.-onecost);         377     sinThetaSqr   = onecost*(2.-onecost);
378                                                   378 
379     // Protection                                 379     // Protection
380     if (sinThetaSqr > 1.)                         380     if (sinThetaSqr > 1.)
381       {                                           381       {
382   G4cout                                          382   G4cout
383     << " -- Warning -- G4LivermorePolarizedCom    383     << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
384     << "sin(theta)**2 = "                         384     << "sin(theta)**2 = "
385     << sinThetaSqr                                385     << sinThetaSqr
386     << "; set to 1"                               386     << "; set to 1"
387     << G4endl;                                    387     << G4endl;
388   sinThetaSqr = 1.;                               388   sinThetaSqr = 1.;
389       }                                           389       }
390     if (sinThetaSqr < 0.)                         390     if (sinThetaSqr < 0.)
391       {                                           391       {
392   G4cout                                          392   G4cout
393     << " -- Warning -- G4LivermorePolarizedCom    393     << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
394     << "sin(theta)**2 = "                         394     << "sin(theta)**2 = "
395     << sinThetaSqr                                395     << sinThetaSqr
396     << "; set to 0"                               396     << "; set to 0"
397     << G4endl;                                    397     << G4endl;
398   sinThetaSqr = 0.;                               398   sinThetaSqr = 0.;
399       }                                           399       }
400     // End protection                             400     // End protection
401                                                   401 
402     G4double x =  std::sqrt(onecost/2.) / (wlG    402     G4double x =  std::sqrt(onecost/2.) / (wlGamma/cm);;
403     G4double scatteringFunction = scatterFunct    403     G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
404     greject = (1. - epsilon*sinThetaSqr/(1.+ e    404     greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction;
405                                                   405 
406   } while(greject < G4UniformRand()*Z);           406   } while(greject < G4UniformRand()*Z);
407                                                   407 
408                                                   408 
409   // *****************************************    409   // ****************************************************
410   //    Phi determination                         410   //    Phi determination
411   // *****************************************    411   // ****************************************************
412   G4double phi = SetPhi(epsilon,sinThetaSqr);     412   G4double phi = SetPhi(epsilon,sinThetaSqr);
413                                                   413 
414   //                                              414   //
415   // scattered gamma angles. ( Z - axis along     415   // scattered gamma angles. ( Z - axis along the parent gamma)
416   //                                              416   //
417   G4double cosTheta = 1. - onecost;               417   G4double cosTheta = 1. - onecost;
418                                                   418 
419   // Protection                                   419   // Protection
420   if (cosTheta > 1.)                              420   if (cosTheta > 1.)
421     {                                             421     {
422       G4cout                                      422       G4cout
423   << " -- Warning -- G4LivermorePolarizedCompt    423   << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
424   << "cosTheta = "                                424   << "cosTheta = "
425   << cosTheta                                     425   << cosTheta
426   << "; set to 1"                                 426   << "; set to 1"
427   << G4endl;                                      427   << G4endl;
428       cosTheta = 1.;                              428       cosTheta = 1.;
429     }                                             429     }
430   if (cosTheta < -1.)                             430   if (cosTheta < -1.)
431     {                                             431     {
432       G4cout                                      432       G4cout 
433   << " -- Warning -- G4LivermorePolarizedCompt    433   << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
434   << "cosTheta = "                                434   << "cosTheta = " 
435   << cosTheta                                     435   << cosTheta
436   << "; set to -1"                                436   << "; set to -1"
437   << G4endl;                                      437   << G4endl;
438       cosTheta = -1.;                             438       cosTheta = -1.;
439     }                                             439     }
440   // End protection                               440   // End protection      
441                                                   441  
442   G4double sinTheta = std::sqrt (sinThetaSqr);    442   G4double sinTheta = std::sqrt (sinThetaSqr);
443                                                   443   
444   // Protection                                   444   // Protection
445   if (sinTheta > 1.)                              445   if (sinTheta > 1.)
446     {                                             446     {
447       G4cout                                      447       G4cout 
448   << " -- Warning -- G4LivermorePolarizedCompt    448   << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
449   << "sinTheta = "                                449   << "sinTheta = " 
450   << sinTheta                                     450   << sinTheta
451   << "; set to 1"                                 451   << "; set to 1"
452   << G4endl;                                      452   << G4endl;
453       sinTheta = 1.;                              453       sinTheta = 1.;
454     }                                             454     }
455   if (sinTheta < -1.)                             455   if (sinTheta < -1.)
456     {                                             456     {
457       G4cout                                      457       G4cout 
458   << " -- Warning -- G4LivermorePolarizedCompt    458   << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
459   << "sinTheta = "                                459   << "sinTheta = " 
460   << sinTheta                                     460   << sinTheta
461   << "; set to -1"                                461   << "; set to -1" 
462   << G4endl;                                      462   << G4endl;
463       sinTheta = -1.;                             463       sinTheta = -1.;
464     }                                             464     }
465   // End protection                               465   // End protection
466                                                   466   
467   // Check for entanglement and re-sample phi     467   // Check for entanglement and re-sample phi if necessary
468                                                   468   
469   const auto* auxInfo                             469   const auto* auxInfo
470   = fParticleChange->GetCurrentTrack()->GetAux    470   = fParticleChange->GetCurrentTrack()->GetAuxiliaryTrackInformation(fEntanglementModelID);
471   if (auxInfo) {                                  471   if (auxInfo) {
472     const auto* entanglementAuxInfo = dynamic_    472     const auto* entanglementAuxInfo = dynamic_cast<const G4EntanglementAuxInfo*>(auxInfo);
473     if (entanglementAuxInfo) {                    473     if (entanglementAuxInfo) {
474       auto* clipBoard = dynamic_cast<G4eplusAn    474       auto* clipBoard = dynamic_cast<G4eplusAnnihilationEntanglementClipBoard*>
475       (entanglementAuxInfo->GetEntanglementCli    475       (entanglementAuxInfo->GetEntanglementClipBoard());
476       if (clipBoard) {                            476       if (clipBoard) {
477         // This is an entangled photon from ep    477         // This is an entangled photon from eplus annihilation at rest.
478         // If this is the first scatter of the    478         // If this is the first scatter of the first photon, place theta and
479         // phi on the clipboard.                  479         // phi on the clipboard.
480         // If this is the first scatter of the    480         // If this is the first scatter of the second photon, use theta and
481         // phi of the first scatter of the fir    481         // phi of the first scatter of the first photon, together with the
482         // theta of the second photon, to samp    482         // theta of the second photon, to sample phi.
483         if (clipBoard->IsTrack1Measurement())     483         if (clipBoard->IsTrack1Measurement()) {
484           // Check we have the relevant track.    484           // Check we have the relevant track. Not sure this is strictly
485           // necessary but I want to be sure t    485           // necessary but I want to be sure tracks from, say, more than one
486           // entangled system are properly pai    486           // entangled system are properly paired.
487           // Note: the tracking manager pops t    487           // Note: the tracking manager pops the tracks in the reverse order. We
488           // will rely on that.  (If not, the     488           // will rely on that.  (If not, the logic here would have to be a bit
489           // more complicated to ensure we mat    489           // more complicated to ensure we matched the right tracks.)
490           // So our track 1 is clipboard track    490           // So our track 1 is clipboard track B.
491           if (clipBoard->GetTrackB() == fParti    491           if (clipBoard->GetTrackB() == fParticleChange->GetCurrentTrack()) {
492             // This is the first scatter of th    492             // This is the first scatter of the first photon. Reset flag.
493             //            // Debug                493             //            // Debug
494             //            auto* track1 = fPart    494             //            auto* track1 = fParticleChange->GetCurrentTrack();
495             //            G4cout                  495             //            G4cout
496             //            << "This is the firs    496             //            << "This is the first scatter of the first photon. Reset flag."
497             //            << "\nTrack: " << tr    497             //            << "\nTrack: " << track1->GetTrackID()
498             //            << ", Parent: " << t    498             //            << ", Parent: " << track1->GetParentID()
499             //            << ", Name: " << cli    499             //            << ", Name: " << clipBoard->GetParentParticleDefinition()->GetParticleName()
500             //            << G4endl;              500             //            << G4endl;
501             //            // End debug            501             //            // End debug
502             clipBoard->ResetTrack1Measurement(    502             clipBoard->ResetTrack1Measurement();
503             // Store cos(theta),phi of first p    503             // Store cos(theta),phi of first photon.
504             clipBoard->SetComptonCosTheta1(cos    504             clipBoard->SetComptonCosTheta1(cosTheta);
505             clipBoard->SetComptonPhi1(phi);       505             clipBoard->SetComptonPhi1(phi);
506           }                                       506           }
507         } else if (clipBoard->IsTrack2Measurem    507         } else if (clipBoard->IsTrack2Measurement()) {
508           // Check we have the relevant track.    508           // Check we have the relevant track.
509           // Remember our track 2 is clipboard    509           // Remember our track 2 is clipboard track A.
510           if (clipBoard->GetTrackA() == fParti    510           if (clipBoard->GetTrackA() == fParticleChange->GetCurrentTrack()) {
511             // This is the first scatter of th    511             // This is the first scatter of the second photon. Reset flag.
512             //            // Debug                512             //            // Debug
513             //            auto* track2 = fPart    513             //            auto* track2 = fParticleChange->GetCurrentTrack();
514             //            G4cout                  514             //            G4cout
515             //            << "This is the firs    515             //            << "This is the first scatter of the second photon. Reset flag."
516             //            << "\nTrack: " << tr    516             //            << "\nTrack: " << track2->GetTrackID()
517             //            << ", Parent: " << t    517             //            << ", Parent: " << track2->GetParentID()
518             //            << ", Name: " << cli    518             //            << ", Name: " << clipBoard->GetParentParticleDefinition()->GetParticleName()
519             //            << G4endl;              519             //            << G4endl;
520             //            // End debug            520             //            // End debug
521             clipBoard->ResetTrack2Measurement(    521             clipBoard->ResetTrack2Measurement();
522                                                   522             
523             // Get cos(theta),phi of first pho    523             // Get cos(theta),phi of first photon.
524             const G4double& cosTheta1 = clipBo    524             const G4double& cosTheta1 = clipBoard->GetComptonCosTheta1();
525             const G4double& phi1 = clipBoard->    525             const G4double& phi1 = clipBoard->GetComptonPhi1();
526             // For clarity make aliases for th    526             // For clarity make aliases for the current cos(theta),phi.
527             const G4double& cosTheta2 = cosThe    527             const G4double& cosTheta2 = cosTheta;
528             G4double& phi2 = phi;                 528             G4double& phi2 = phi;
529             //            G4cout << "cosTheta1    529             //            G4cout << "cosTheta1,phi1: " << cosTheta1 << ',' << phi1 << G4endl;
530             //            G4cout << "cosTheta2    530             //            G4cout << "cosTheta2,phi2: " << cosTheta2 << ',' << phi2 << G4endl;
531                                                   531             
532             // Re-sample phi                      532             // Re-sample phi
533             // Draw the difference of azimutha    533             // Draw the difference of azimuthal angles, deltaPhi, from
534             // A + B * cos(2*deltaPhi), or rat    534             // A + B * cos(2*deltaPhi), or rather C + D * cos(2*deltaPhi), where
535             // C = A / (A + |B|) and D = B / (    535             // C = A / (A + |B|) and D = B / (A + |B|), so that maximum is 1.
536             const G4double sin2Theta1 = 1.-cos    536             const G4double sin2Theta1 = 1.-cosTheta1*cosTheta1;
537             const G4double sin2Theta2 = 1.-cos    537             const G4double sin2Theta2 = 1.-cosTheta2*cosTheta2;
538                                                   538             
539             // Pryce and Ward, Nature No 4065     539             // Pryce and Ward, Nature No 4065 (1947) p.435.
540             auto* g4Pow = G4Pow::GetInstance()    540             auto* g4Pow = G4Pow::GetInstance();
541             const G4double A =                    541             const G4double A =
542             ((g4Pow->powN(1.-cosTheta1,3))+2.)    542             ((g4Pow->powN(1.-cosTheta1,3))+2.)*(g4Pow->powN(1.-cosTheta2,3)+2.)/
543             ((g4Pow->powN(2.-cosTheta1,3)*g4Po    543             ((g4Pow->powN(2.-cosTheta1,3)*g4Pow->powN(2.-cosTheta2,3)));
544             const G4double B = -(sin2Theta1*si    544             const G4double B = -(sin2Theta1*sin2Theta2)/
545             ((g4Pow->powN(2.-cosTheta1,2)*g4Po    545             ((g4Pow->powN(2.-cosTheta1,2)*g4Pow->powN(2.-cosTheta2,2)));
546                                                   546 
547             //    // Snyder et al, Physical Re    547             //    // Snyder et al, Physical Review 73 (1948) p.440.
548             //    // (This is an alternative f    548             //    // (This is an alternative formulation but result is identical.)
549             //    const G4double& k0 = gammaEn    549             //    const G4double& k0 = gammaEnergy0;
550             //    const G4double k1 = k0/(2.-c    550             //    const G4double k1 = k0/(2.-cosTheta1);
551             //    const G4double k2 = k0/(2.-c    551             //    const G4double k2 = k0/(2.-cosTheta2);
552             //    const G4double gamma1 = k1/k    552             //    const G4double gamma1 = k1/k0+k0/k1;
553             //    const G4double gamma2 = k2/k    553             //    const G4double gamma2 = k2/k0+k0/k2;
554             //    const G4double A1 = gamma1*g    554             //    const G4double A1 = gamma1*gamma2-gamma1*sin2Theta2-gamma2*sin2Theta1;
555             //    const G4double B1 = 2.*sin2T    555             //    const G4double B1 = 2.*sin2Theta1*sin2Theta2;
556             //    // That's A1 + B1*sin2(delta    556             //    // That's A1 + B1*sin2(deltaPhi) = A1 + B1*(0.5*(1.-cos(2.*deltaPhi).
557             //    const G4double A = A1 + 0.5*    557             //    const G4double A = A1 + 0.5*B1;
558             //    const G4double B = -0.5*B1;     558             //    const G4double B = -0.5*B1;
559                                                   559             
560             const G4double maxValue = A + std:    560             const G4double maxValue = A + std::abs(B);
561             const G4double C = A / maxValue;      561             const G4double C = A / maxValue;
562             const G4double D = B / maxValue;      562             const G4double D = B / maxValue;
563             //    G4cout << "A,B,C,D: " << A <    563             //    G4cout << "A,B,C,D: " << A << ',' << B  << ',' << C << ',' << D << G4endl;
564                                                   564             
565             // Sample delta phi                   565             // Sample delta phi
566             G4double deltaPhi;                    566             G4double deltaPhi;
567             const G4int maxCount = 999999;        567             const G4int maxCount = 999999;
568             G4int iCount = 0;                     568             G4int iCount = 0;
569             for (; iCount < maxCount; ++iCount    569             for (; iCount < maxCount; ++iCount) {
570               deltaPhi = twopi * G4UniformRand    570               deltaPhi = twopi * G4UniformRand();
571               if (G4UniformRand() < C + D * co    571               if (G4UniformRand() < C + D * cos(2.*deltaPhi)) break;
572             }                                     572             }
573             if (iCount >= maxCount ) {            573             if (iCount >= maxCount ) {
574               G4cout << "G4LivermorePolarizedC    574               G4cout << "G4LivermorePolarizedComptonModel::SampleSecondaries: "
575               << "Re-sampled delta phi not fou    575               << "Re-sampled delta phi not found in " << maxCount
576               << " tries - carrying on anyway.    576               << " tries - carrying on anyway." << G4endl;
577             }                                     577             }
578                                                   578             
579             // Thus, the desired second photon    579             // Thus, the desired second photon azimuth
580             phi2 = deltaPhi - phi1 + halfpi;      580             phi2 = deltaPhi - phi1 + halfpi;
581             // The minus sign is in above stat    581             // The minus sign is in above statement because, since the two
582             // annihilation photons are in opp    582             // annihilation photons are in opposite directions, their phi's
583             // are measured in the opposite di    583             // are measured in the opposite direction.
584             // halfpi is added for the followi    584             // halfpi is added for the following reason:
585             // In this function phi is relativ    585             // In this function phi is relative to the polarisation - see
586             // SystemOfRefChange below. We kno    586             // SystemOfRefChange below. We know from G4eplusAnnihilation that
587             // the polarisations of the two an    587             // the polarisations of the two annihilation photons are perpendicular
588             // to each other, i.e., halfpi dif    588             // to each other, i.e., halfpi different.
589             // Furthermore, only sin(phi) and     589             // Furthermore, only sin(phi) and cos(phi) are used below so no
590             // need to place any range constra    590             // need to place any range constraints.
591             //            if (phi2 > pi) {        591             //            if (phi2 > pi) {
592             //              phi2 -= twopi;        592             //              phi2 -= twopi;
593             //            }                       593             //            }
594             //            if (phi2 < -pi) {       594             //            if (phi2 < -pi) {
595             //              phi2 += twopi;        595             //              phi2 += twopi;
596             //            }                       596             //            }
597           }                                       597           }
598         }                                         598         }
599       }                                           599       }
600     }                                             600     }
601   }                                               601   }
602                                                   602   
603   // End of entanglement                          603   // End of entanglement
604                                                   604       
605   G4double dirx = sinTheta*std::cos(phi);         605   G4double dirx = sinTheta*std::cos(phi);
606   G4double diry = sinTheta*std::sin(phi);         606   G4double diry = sinTheta*std::sin(phi);
607   G4double dirz = cosTheta ;                      607   G4double dirz = cosTheta ;
608                                                   608  
609   // oneCosT , eom                                609   // oneCosT , eom
610                                                   610 
611   // Doppler broadening -  Method based on:       611   // Doppler broadening -  Method based on:
612   // Y. Namito, S. Ban and H. Hirayama,           612   // Y. Namito, S. Ban and H. Hirayama, 
613   // "Implementation of the Doppler Broadening    613   // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code" 
614   // NIM A 349, pp. 489-494, 1994                 614   // NIM A 349, pp. 489-494, 1994
615                                                   615   
616   // Maximum number of sampling iterations        616   // Maximum number of sampling iterations
617   static G4int maxDopplerIterations = 1000;       617   static G4int maxDopplerIterations = 1000;
618   G4double bindingE = 0.;                         618   G4double bindingE = 0.;
619   G4double photonEoriginal = epsilon * gammaEn    619   G4double photonEoriginal = epsilon * gammaEnergy0;
620   G4double photonE = -1.;                         620   G4double photonE = -1.;
621   G4int iteration = 0;                            621   G4int iteration = 0;
622   G4double eMax = gammaEnergy0;                   622   G4double eMax = gammaEnergy0;
623                                                   623 
624   G4int shellIdx = 0;                             624   G4int shellIdx = 0;
625                                                   625 
626   if (verboseLevel > 3) {                         626   if (verboseLevel > 3) {
627     G4cout << "Started loop to sample broading    627     G4cout << "Started loop to sample broading" << G4endl;
628   }                                               628   }
629                                                   629 
630   do                                              630   do
631     {                                             631     {
632       iteration++;                                632       iteration++;
633       // Select shell based on shell occupancy    633       // Select shell based on shell occupancy
634       shellIdx = shellData->SelectRandomShell(    634       shellIdx = shellData->SelectRandomShell(Z);
635       bindingE = shellData->BindingEnergy(Z,sh    635       bindingE = shellData->BindingEnergy(Z,shellIdx);
636                                                   636       
637       if (verboseLevel > 3) {                     637       if (verboseLevel > 3) {
638   G4cout << "Shell ID= " << shellIdx              638   G4cout << "Shell ID= " << shellIdx 
639          << " Ebind(keV)= " << bindingE/keV <<    639          << " Ebind(keV)= " << bindingE/keV << G4endl;
640       }                                           640       }
641       eMax = gammaEnergy0 - bindingE;             641       eMax = gammaEnergy0 - bindingE;
642                                                   642       
643       // Randomly sample bound electron moment    643       // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
644       G4double pSample = profileData->RandomSe    644       G4double pSample = profileData->RandomSelectMomentum(Z,shellIdx);
645                                                   645 
646       if (verboseLevel > 3) {                     646       if (verboseLevel > 3) {
647        G4cout << "pSample= " << pSample << G4e    647        G4cout << "pSample= " << pSample << G4endl;
648      }                                            648      }
649       // Rescale from atomic units                649       // Rescale from atomic units
650       G4double pDoppler = pSample * fine_struc    650       G4double pDoppler = pSample * fine_structure_const;
651       G4double pDoppler2 = pDoppler * pDoppler    651       G4double pDoppler2 = pDoppler * pDoppler;
652       G4double var2 = 1. + onecost * E0_m;        652       G4double var2 = 1. + onecost * E0_m;
653       G4double var3 = var2*var2 - pDoppler2;      653       G4double var3 = var2*var2 - pDoppler2;
654       G4double var4 = var2 - pDoppler2 * cosTh    654       G4double var4 = var2 - pDoppler2 * cosTheta;
655       G4double var = var4*var4 - var3 + pDoppl    655       G4double var = var4*var4 - var3 + pDoppler2 * var3;
656       if (var > 0.)                               656       if (var > 0.)
657   {                                               657   {
658     G4double varSqrt = std::sqrt(var);            658     G4double varSqrt = std::sqrt(var);        
659     G4double scale = gammaEnergy0 / var3;         659     G4double scale = gammaEnergy0 / var3;  
660           // Random select either root            660           // Random select either root
661     if (G4UniformRand() < 0.5) photonE = (var4    661     if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;               
662     else photonE = (var4 + varSqrt) * scale;      662     else photonE = (var4 + varSqrt) * scale;
663   }                                               663   } 
664       else                                        664       else
665   {                                               665   {
666     photonE = -1.;                                666     photonE = -1.;
667   }                                               667   }
668    } while ( iteration <= maxDopplerIterations    668    } while ( iteration <= maxDopplerIterations && 
669        (photonE < 0. || photonE > eMax || phot    669        (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) );
670                                                   670 
671   // End of recalculation of photon energy wit    671   // End of recalculation of photon energy with Doppler broadening
672   // Revert to original if maximum number of i    672   // Revert to original if maximum number of iterations threshold has been reached
673   if (iteration >= maxDopplerIterations)          673   if (iteration >= maxDopplerIterations)
674     {                                             674     {
675       photonE = photonEoriginal;                  675       photonE = photonEoriginal;
676       bindingE = 0.;                              676       bindingE = 0.;
677     }                                             677     }
678                                                   678 
679   gammaEnergy1 = photonE;                         679   gammaEnergy1 = photonE;
680                                                   680  
681   //                                              681   //
682   // update G4VParticleChange for the scattere    682   // update G4VParticleChange for the scattered photon 
683   //                                              683   //
684   // New polarization                             684   // New polarization
685   G4ThreeVector gammaPolarization1 = SetNewPol    685   G4ThreeVector gammaPolarization1 = SetNewPolarization(epsilon,
686               sinThetaSqr,                        686               sinThetaSqr,
687               phi,                                687               phi,
688               cosTheta);                          688               cosTheta);
689                                                   689 
690   // Set new direction                            690   // Set new direction
691   G4ThreeVector tmpDirection1( dirx,diry,dirz     691   G4ThreeVector tmpDirection1( dirx,diry,dirz );
692   gammaDirection1 = tmpDirection1;                692   gammaDirection1 = tmpDirection1;
693                                                   693 
694   // Change reference frame.                      694   // Change reference frame.
695   SystemOfRefChange(gammaDirection0,gammaDirec    695   SystemOfRefChange(gammaDirection0,gammaDirection1,
696         gammaPolarization0,gammaPolarization1)    696         gammaPolarization0,gammaPolarization1);
697                                                   697 
698   if (gammaEnergy1 > 0.)                          698   if (gammaEnergy1 > 0.)
699     {                                             699     {
700       fParticleChange->SetProposedKineticEnerg    700       fParticleChange->SetProposedKineticEnergy( gammaEnergy1 ) ;
701       fParticleChange->ProposeMomentumDirectio    701       fParticleChange->ProposeMomentumDirection( gammaDirection1 );
702       fParticleChange->ProposePolarization( ga    702       fParticleChange->ProposePolarization( gammaPolarization1 );
703     }                                             703     }
704   else                                            704   else
705     {                                             705     {
706       gammaEnergy1 = 0.;                          706       gammaEnergy1 = 0.;
707       fParticleChange->SetProposedKineticEnerg    707       fParticleChange->SetProposedKineticEnergy(0.) ;
708       fParticleChange->ProposeTrackStatus(fSto    708       fParticleChange->ProposeTrackStatus(fStopAndKill);
709     }                                             709     }
710                                                   710 
711   //                                              711   //
712   // kinematic of the scattered electron          712   // kinematic of the scattered electron
713   //                                              713   //
714   G4double ElecKineEnergy = gammaEnergy0 - gam    714   G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;
715                                                   715 
716   // SI -protection against negative final ene    716   // SI -protection against negative final energy: no e- is created
717   // like in G4LivermoreComptonModel.cc           717   // like in G4LivermoreComptonModel.cc
718   if(ElecKineEnergy < 0.0) {                      718   if(ElecKineEnergy < 0.0) {
719     fParticleChange->ProposeLocalEnergyDeposit    719     fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0 - gammaEnergy1);
720     return;                                       720     return;
721   }                                               721   }
722                                                   722  
723   G4double ElecMomentum = std::sqrt(ElecKineEn    723   G4double ElecMomentum = std::sqrt(ElecKineEnergy*(ElecKineEnergy+2.*electron_mass_c2));
724                                                   724 
725   G4ThreeVector ElecDirection((gammaEnergy0 *     725   G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 -
726            gammaEnergy1 * gammaDirection1) * (    726            gammaEnergy1 * gammaDirection1) * (1./ElecMomentum));
727                                                   727 
728   G4DynamicParticle* dp =                         728   G4DynamicParticle* dp = 
729     new G4DynamicParticle (G4Electron::Electro    729     new G4DynamicParticle (G4Electron::Electron(),ElecDirection.unit(),ElecKineEnergy) ;
730   fvect->push_back(dp);                           730   fvect->push_back(dp);
731                                                   731 
732   // sample deexcitation                          732   // sample deexcitation
733   //                                              733   //  
734   if (verboseLevel > 3) {                         734   if (verboseLevel > 3) {
735     G4cout << "Started atomic de-excitation "     735     G4cout << "Started atomic de-excitation " << fAtomDeexcitation << G4endl;
736   }                                               736   }
737                                                   737   
738   if(fAtomDeexcitation && iteration < maxDoppl    738   if(fAtomDeexcitation && iteration < maxDopplerIterations) {
739     G4int index = couple->GetIndex();             739     G4int index = couple->GetIndex();
740     if(fAtomDeexcitation->CheckDeexcitationAct    740     if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
741       std::size_t nbefore = fvect->size();        741       std::size_t nbefore = fvect->size();
742       G4AtomicShellEnumerator as = G4AtomicShe    742       G4AtomicShellEnumerator as = G4AtomicShellEnumerator(shellIdx);
743       const G4AtomicShell* shell = fAtomDeexci    743       const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
744       fAtomDeexcitation->GenerateParticles(fve    744       fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
745       std::size_t nafter = fvect->size();         745       std::size_t nafter = fvect->size();
746       if(nafter > nbefore) {                      746       if(nafter > nbefore) {
747   for (std::size_t i=nbefore; i<nafter; ++i) {    747   for (std::size_t i=nbefore; i<nafter; ++i) {
748     //Check if there is enough residual energy    748     //Check if there is enough residual energy 
749     if (bindingE >= ((*fvect)[i])->GetKineticE    749     if (bindingE >= ((*fvect)[i])->GetKineticEnergy())
750             {                                     750             {
751               //Ok, this is a valid secondary:    751               //Ok, this is a valid secondary: keep it
752         bindingE -= ((*fvect)[i])->GetKineticE    752         bindingE -= ((*fvect)[i])->GetKineticEnergy();
753             }                                     753             }
754     else                                          754     else
755             {                                     755             {
756         //Invalid secondary: not enough energy    756         //Invalid secondary: not enough energy to create it!
757         //Keep its energy in the local deposit    757         //Keep its energy in the local deposit
758         delete (*fvect)[i];                       758         delete (*fvect)[i]; 
759               (*fvect)[i]=0;                      759               (*fvect)[i]=0;
760             }                                     760             }
761   }                                               761   } 
762       }                                           762       }
763     }                                             763     }
764   }                                               764   }
765   //This should never happen                      765   //This should never happen
766   if(bindingE < 0.0)                              766   if(bindingE < 0.0) 
767     G4Exception("G4LivermoreComptonModel::Samp    767     G4Exception("G4LivermoreComptonModel::SampleSecondaries()",
768     "em2050",FatalException,"Negative local en    768     "em2050",FatalException,"Negative local energy deposit");   
769                                                   769   
770   fParticleChange->ProposeLocalEnergyDeposit(b    770   fParticleChange->ProposeLocalEnergyDeposit(bindingE);
771                                                   771 
772 }                                                 772 }
773                                                   773 
774 //....oooOO0OOooo........oooOO0OOooo........oo    774 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
775                                                   775 
776 G4double G4LivermorePolarizedComptonModel::Set    776 G4double G4LivermorePolarizedComptonModel::SetPhi(G4double energyRate,
777                G4double sinSqrTh)                 777                G4double sinSqrTh)
778 {                                                 778 {
779   G4double rand1;                                 779   G4double rand1;
780   G4double rand2;                                 780   G4double rand2;
781   G4double phiProbability;                        781   G4double phiProbability;
782   G4double phi;                                   782   G4double phi;
783   G4double a, b;                                  783   G4double a, b;
784                                                   784 
785   do                                              785   do
786     {                                             786     {
787       rand1 = G4UniformRand();                    787       rand1 = G4UniformRand();
788       rand2 = G4UniformRand();                    788       rand2 = G4UniformRand();
789       phiProbability=0.;                          789       phiProbability=0.;
790       phi = twopi*rand1;                          790       phi = twopi*rand1;
791                                                   791       
792       a = 2*sinSqrTh;                             792       a = 2*sinSqrTh;
793       b = energyRate + 1/energyRate;              793       b = energyRate + 1/energyRate;
794                                                   794       
795       phiProbability = 1 - (a/b)*(std::cos(phi    795       phiProbability = 1 - (a/b)*(std::cos(phi)*std::cos(phi));
796     }                                             796     }
797   while ( rand2 > phiProbability );               797   while ( rand2 > phiProbability );
798   return phi;                                     798   return phi;
799 }                                                 799 }
800                                                   800 
801                                                   801 
802 //....oooOO0OOooo........oooOO0OOooo........oo    802 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
803                                                   803 
804 G4ThreeVector G4LivermorePolarizedComptonModel    804 G4ThreeVector G4LivermorePolarizedComptonModel::SetPerpendicularVector(G4ThreeVector& a)
805 {                                                 805 {
806   G4double dx = a.x();                            806   G4double dx = a.x();
807   G4double dy = a.y();                            807   G4double dy = a.y();
808   G4double dz = a.z();                            808   G4double dz = a.z();
809   G4double x = dx < 0.0 ? -dx : dx;               809   G4double x = dx < 0.0 ? -dx : dx;
810   G4double y = dy < 0.0 ? -dy : dy;               810   G4double y = dy < 0.0 ? -dy : dy;
811   G4double z = dz < 0.0 ? -dz : dz;               811   G4double z = dz < 0.0 ? -dz : dz;
812   if (x < y) {                                    812   if (x < y) {
813     return x < z ? G4ThreeVector(-dy,dx,0) : G    813     return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
814   }else{                                          814   }else{
815     return y < z ? G4ThreeVector(dz,0,-dx) : G    815     return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
816   }                                               816   }
817 }                                                 817 }
818                                                   818 
819 //....oooOO0OOooo........oooOO0OOooo........oo    819 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
820                                                   820 
821 G4ThreeVector G4LivermorePolarizedComptonModel    821 G4ThreeVector G4LivermorePolarizedComptonModel::GetRandomPolarization(G4ThreeVector& direction0)
822 {                                                 822 {
823   G4ThreeVector d0 = direction0.unit();           823   G4ThreeVector d0 = direction0.unit();
824   G4ThreeVector a1 = SetPerpendicularVector(d0    824   G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
825   G4ThreeVector a0 = a1.unit(); // unit vector    825   G4ThreeVector a0 = a1.unit(); // unit vector
826                                                   826 
827   G4double rand1 = G4UniformRand();               827   G4double rand1 = G4UniformRand();
828                                                   828   
829   G4double angle = twopi*rand1; // random pola    829   G4double angle = twopi*rand1; // random polar angle
830   G4ThreeVector b0 = d0.cross(a0); // cross pr    830   G4ThreeVector b0 = d0.cross(a0); // cross product
831                                                   831   
832   G4ThreeVector c;                                832   G4ThreeVector c;
833                                                   833   
834   c.setX(std::cos(angle)*(a0.x())+std::sin(ang    834   c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
835   c.setY(std::cos(angle)*(a0.y())+std::sin(ang    835   c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
836   c.setZ(std::cos(angle)*(a0.z())+std::sin(ang    836   c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
837                                                   837   
838   G4ThreeVector c0 = c.unit();                    838   G4ThreeVector c0 = c.unit();
839                                                   839 
840   return c0;                                      840   return c0;
841 }                                                 841 }
842                                                   842 
843 //....oooOO0OOooo........oooOO0OOooo........oo    843 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
844                                                   844 
845 G4ThreeVector G4LivermorePolarizedComptonModel    845 G4ThreeVector G4LivermorePolarizedComptonModel::GetPerpendicularPolarization
846 (const G4ThreeVector& gammaDirection, const G4    846 (const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
847 {                                                 847 {
848   //                                              848   // 
849   // The polarization of a photon is always pe    849   // The polarization of a photon is always perpendicular to its momentum direction.
850   // Therefore this function removes those vec    850   // Therefore this function removes those vector component of gammaPolarization, which
851   // points in direction of gammaDirection        851   // points in direction of gammaDirection
852   //                                              852   //
853   // Mathematically we search the projection o    853   // Mathematically we search the projection of the vector a on the plane E, where n is the
854   // plains normal vector.                        854   // plains normal vector.
855   // The basic equation can be found in each g    855   // The basic equation can be found in each geometry book (e.g. Bronstein):
856   // p = a - (a o n)/(n o n)*n                    856   // p = a - (a o n)/(n o n)*n
857                                                   857   
858   return gammaPolarization - gammaPolarization    858   return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
859 }                                                 859 }
860                                                   860 
861 //....oooOO0OOooo........oooOO0OOooo........oo    861 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
862                                                   862 
863 G4ThreeVector G4LivermorePolarizedComptonModel    863 G4ThreeVector G4LivermorePolarizedComptonModel::SetNewPolarization(G4double epsilon,
864                     G4double sinSqrTh,            864                     G4double sinSqrTh, 
865                     G4double phi,                 865                     G4double phi,
866                     G4double costheta)            866                     G4double costheta) 
867 {                                                 867 {
868   G4double rand1;                                 868   G4double rand1;
869   G4double rand2;                                 869   G4double rand2;
870   G4double cosPhi = std::cos(phi);                870   G4double cosPhi = std::cos(phi);
871   G4double sinPhi = std::sin(phi);                871   G4double sinPhi = std::sin(phi);
872   G4double sinTheta = std::sqrt(sinSqrTh);        872   G4double sinTheta = std::sqrt(sinSqrTh);
873   G4double cosSqrPhi = cosPhi*cosPhi;             873   G4double cosSqrPhi = cosPhi*cosPhi;
874   //  G4double cossqrth = 1.-sinSqrTh;            874   //  G4double cossqrth = 1.-sinSqrTh;
875   //  G4double sinsqrphi = sinPhi*sinPhi;         875   //  G4double sinsqrphi = sinPhi*sinPhi;
876   G4double normalisation = std::sqrt(1. - cosS    876   G4double normalisation = std::sqrt(1. - cosSqrPhi*sinSqrTh);
877                                                   877  
878   // Determination of Theta                       878   // Determination of Theta 
879   G4double theta;                                 879   G4double theta;
880                                                   880 
881   // Dan Xu method (IEEE TNS, 52, 1160 (2005))    881   // Dan Xu method (IEEE TNS, 52, 1160 (2005))
882   rand1 = G4UniformRand();                        882   rand1 = G4UniformRand();
883   rand2 = G4UniformRand();                        883   rand2 = G4UniformRand();
884                                                   884 
885   if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsi    885   if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsilon+1.0/epsilon)-4.0*sinSqrTh*cosSqrPhi))
886     {                                             886     {
887       if (rand2<0.5)                              887       if (rand2<0.5)
888   theta = pi/2.0;                                 888   theta = pi/2.0;
889       else                                        889       else
890   theta = 3.0*pi/2.0;                             890   theta = 3.0*pi/2.0;
891     }                                             891     }
892   else                                            892   else
893     {                                             893     {
894       if (rand2<0.5)                              894       if (rand2<0.5)
895   theta = 0;                                      895   theta = 0;
896       else                                        896       else
897   theta = pi;                                     897   theta = pi;
898     }                                             898     }
899   G4double cosBeta = std::cos(theta);             899   G4double cosBeta = std::cos(theta);
900   G4double sinBeta = std::sqrt(1-cosBeta*cosBe    900   G4double sinBeta = std::sqrt(1-cosBeta*cosBeta);
901                                                   901   
902   G4ThreeVector gammaPolarization1;               902   G4ThreeVector gammaPolarization1;
903                                                   903 
904   G4double xParallel = normalisation*cosBeta;     904   G4double xParallel = normalisation*cosBeta;
905   G4double yParallel = -(sinSqrTh*cosPhi*sinPh    905   G4double yParallel = -(sinSqrTh*cosPhi*sinPhi)*cosBeta/normalisation;
906   G4double zParallel = -(costheta*sinTheta*cos    906   G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation;
907   G4double xPerpendicular = 0.;                   907   G4double xPerpendicular = 0.;
908   G4double yPerpendicular = (costheta)*sinBeta    908   G4double yPerpendicular = (costheta)*sinBeta/normalisation;
909   G4double zPerpendicular = -(sinTheta*sinPhi)    909   G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation;
910                                                   910 
911   G4double xTotal = (xParallel + xPerpendicula    911   G4double xTotal = (xParallel + xPerpendicular);
912   G4double yTotal = (yParallel + yPerpendicula    912   G4double yTotal = (yParallel + yPerpendicular);
913   G4double zTotal = (zParallel + zPerpendicula    913   G4double zTotal = (zParallel + zPerpendicular);
914                                                   914   
915   gammaPolarization1.setX(xTotal);                915   gammaPolarization1.setX(xTotal);
916   gammaPolarization1.setY(yTotal);                916   gammaPolarization1.setY(yTotal);
917   gammaPolarization1.setZ(zTotal);                917   gammaPolarization1.setZ(zTotal);
918                                                   918   
919   return gammaPolarization1;                      919   return gammaPolarization1;
920 }                                                 920 }
921                                                   921 
922 //....oooOO0OOooo........oooOO0OOooo........oo    922 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
923                                                   923 
924 void G4LivermorePolarizedComptonModel::SystemO    924 void G4LivermorePolarizedComptonModel::SystemOfRefChange(G4ThreeVector& direction0,
925                 G4ThreeVector& direction1,        925                 G4ThreeVector& direction1,
926                 G4ThreeVector& polarization0,     926                 G4ThreeVector& polarization0,
927                 G4ThreeVector& polarization1)     927                 G4ThreeVector& polarization1)
928 {                                                 928 {
929   // direction0 is the original photon directi    929   // direction0 is the original photon direction ---> z
930   // polarization0 is the original photon pola    930   // polarization0 is the original photon polarization ---> x
931   // need to specify y axis in the real refere    931   // need to specify y axis in the real reference frame ---> y 
932   G4ThreeVector Axis_Z0 = direction0.unit();      932   G4ThreeVector Axis_Z0 = direction0.unit();
933   G4ThreeVector Axis_X0 = polarization0.unit()    933   G4ThreeVector Axis_X0 = polarization0.unit();
934   G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_    934   G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
935                                                   935 
936   G4double direction_x = direction1.getX();       936   G4double direction_x = direction1.getX();
937   G4double direction_y = direction1.getY();       937   G4double direction_y = direction1.getY();
938   G4double direction_z = direction1.getZ();       938   G4double direction_z = direction1.getZ();
939                                                   939   
940   direction1 = (direction_x*Axis_X0 + directio    940   direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
941   G4double polarization_x = polarization1.getX    941   G4double polarization_x = polarization1.getX();
942   G4double polarization_y = polarization1.getY    942   G4double polarization_y = polarization1.getY();
943   G4double polarization_z = polarization1.getZ    943   G4double polarization_z = polarization1.getZ();
944                                                   944 
945   polarization1 = (polarization_x*Axis_X0 + po    945   polarization1 = (polarization_x*Axis_X0 + polarization_y*Axis_Y0 + polarization_z*Axis_Z0).unit();
946                                                   946 
947 }                                                 947 }
948                                                   948 
949                                                   949 
950 //....oooOO0OOooo........oooOO0OOooo........oo    950 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
951 //....oooOO0OOooo........oooOO0OOooo........oo    951 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
952                                                   952 
953 void                                              953 void 
954 G4LivermorePolarizedComptonModel::InitialiseFo    954 G4LivermorePolarizedComptonModel::InitialiseForElement(const G4ParticleDefinition*, 
955                 G4int Z)                          955                 G4int Z)
956 {                                                 956 {
957   G4AutoLock l(&LivermorePolarizedComptonModel    957   G4AutoLock l(&LivermorePolarizedComptonModelMutex);
958   if(!data[Z]) { ReadData(Z); }                   958   if(!data[Z]) { ReadData(Z); }
959   l.unlock();                                     959   l.unlock();
960 }                                                 960 }
961                                                   961