Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4JAEAPolarizedElasticScatteringModel.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/G4JAEAPolarizedElasticScatteringModel.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4JAEAPolarizedElasticScatteringModel.cc (Version 11.2.2)


  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:                                         27   Authors:
 28   M. Omer and R. Hajima  on 15 November 2019       28   M. Omer and R. Hajima  on 15 November 2019
 29   contact:                                         29   contact:
 30   omer.mohamed@jaea.go.jp and hajima.ryoichi@q     30   omer.mohamed@jaea.go.jp and hajima.ryoichi@qst.go.jp
 31   Publication Information:                         31   Publication Information:
 32   1- M. Omer, R. Hajima, Validating polarizati     32   1- M. Omer, R. Hajima, Validating polarization effects in gamma-rays elastic scattering by Monte
 33   Carlo simulation, New J. Phys., vol. 21, 201     33   Carlo simulation, New J. Phys., vol. 21, 2019, pp. 113006 (1-10),
 34   https://doi.org/10.1088/1367-2630/ab4d8a         34   https://doi.org/10.1088/1367-2630/ab4d8a
 35 */                                                 35 */
 36                                                    36 
 37 #include "G4JAEAPolarizedElasticScatteringMode     37 #include "G4JAEAPolarizedElasticScatteringModel.hh"
 38 #include "G4SystemOfUnits.hh"                      38 #include "G4SystemOfUnits.hh"
 39 #include "G4AutoLock.hh"                           39 #include "G4AutoLock.hh"
 40 namespace { G4Mutex G4JAEAPolarizedElasticScat     40 namespace { G4Mutex G4JAEAPolarizedElasticScatteringModelMutex = G4MUTEX_INITIALIZER; }
 41 using namespace std;                               41 using namespace std;
 42                                                    42 
 43 //....oooOO0OOooo........oooOO0OOooo........oo     43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 44 //....oooOO0OOooo........oooOO0OOooo........oo     44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 45                                                    45 
 46 G4PhysicsFreeVector* G4JAEAPolarizedElasticSca     46 G4PhysicsFreeVector* G4JAEAPolarizedElasticScatteringModel::dataCS[] = {nullptr};
 47 G4DataVector* G4JAEAPolarizedElasticScattering     47 G4DataVector* G4JAEAPolarizedElasticScatteringModel::Polarized_ES_Data[] = {nullptr};
 48                                                    48 
 49 G4JAEAPolarizedElasticScatteringModel::G4JAEAP     49 G4JAEAPolarizedElasticScatteringModel::G4JAEAPolarizedElasticScatteringModel()
 50   :G4VEmModel("G4JAEAPolarizedElasticScatterin     50   :G4VEmModel("G4JAEAPolarizedElasticScatteringModel"),isInitialised(false)
 51 {                                                  51 {
 52   fParticleChange = 0;                             52   fParticleChange = 0;
 53   lowEnergyLimit  = 100 * keV;   //low energy      53   lowEnergyLimit  = 100 * keV;   //low energy limit for JAEAElasticScattering cross section data
 54   fLinearPolarizationSensitvity1=1;                54   fLinearPolarizationSensitvity1=1;
 55   fLinearPolarizationSensitvity2=1;                55   fLinearPolarizationSensitvity2=1;
 56   fCircularPolarizationSensitvity=1;               56   fCircularPolarizationSensitvity=1;
 57                                                    57 
 58   verboseLevel= 0;                                 58   verboseLevel= 0;
 59   // Verbosity scale for debugging purposes:       59   // Verbosity scale for debugging purposes:
 60   // 0 = nothing                                   60   // 0 = nothing
 61   // 1 = calculation of cross sections, file o     61   // 1 = calculation of cross sections, file openings...
 62   // 2 = entering in methods                       62   // 2 = entering in methods
 63                                                    63 
 64   if(verboseLevel > 0)                             64   if(verboseLevel > 0)
 65     {                                              65     {
 66       G4cout << "G4JAEAPolarizedElasticScatter     66       G4cout << "G4JAEAPolarizedElasticScatteringModel is constructed " << G4endl;
 67     }                                              67     }
 68                                                    68 
 69 }                                                  69 }
 70                                                    70 
 71 //....oooOO0OOooo........oooOO0OOooo........oo     71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 72                                                    72 
 73 G4JAEAPolarizedElasticScatteringModel::~G4JAEA     73 G4JAEAPolarizedElasticScatteringModel::~G4JAEAPolarizedElasticScatteringModel()
 74 {                                                  74 {
 75   if(IsMaster()) {                                 75   if(IsMaster()) {
 76     for(G4int i=0; i<=maxZ; ++i) {                 76     for(G4int i=0; i<=maxZ; ++i) {
 77       if(dataCS[i]) {                              77       if(dataCS[i]) {
 78   delete dataCS[i];                                78   delete dataCS[i];
 79   dataCS[i] = nullptr;                             79   dataCS[i] = nullptr;
 80       }                                            80       }
 81       if (Polarized_ES_Data[i]){                   81       if (Polarized_ES_Data[i]){
 82   delete Polarized_ES_Data[i];                     82   delete Polarized_ES_Data[i];
 83   Polarized_ES_Data[i] = nullptr;                  83   Polarized_ES_Data[i] = nullptr;
 84       }                                            84       }
 85     }                                              85     }
 86   }                                                86   }
 87 }                                                  87 }
 88                                                    88 
 89 //....oooOO0OOooo........oooOO0OOooo........oo     89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 90                                                    90 
 91 void G4JAEAPolarizedElasticScatteringModel::In     91 void G4JAEAPolarizedElasticScatteringModel::Initialise(const G4ParticleDefinition* particle,
 92                    const G4DataVector& cuts)       92                    const G4DataVector& cuts)
 93 {                                                  93 {
 94   if (verboseLevel > 1)                            94   if (verboseLevel > 1)
 95     {                                              95     {
 96       G4cout << "Calling Initialise() of G4JAE     96       G4cout << "Calling Initialise() of G4JAEAPolarizedElasticScatteringModel." << G4endl
 97        << "Energy range: "                         97        << "Energy range: "
 98        << LowEnergyLimit() / eV << " eV - "        98        << LowEnergyLimit() / eV << " eV - "
 99        << HighEnergyLimit() / GeV << " GeV"        99        << HighEnergyLimit() / GeV << " GeV"
100        << G4endl;                                 100        << G4endl;
101     }                                             101     }
102                                                   102 
103   if(IsMaster()) {                                103   if(IsMaster()) {
104                                                   104 
105     // Initialise element selector                105     // Initialise element selector
106     InitialiseElementSelectors(particle, cuts)    106     InitialiseElementSelectors(particle, cuts);
107                                                   107 
108     // Access to elements                         108     // Access to elements
109     const char* path = G4FindDataDir("G4LEDATA    109     const char* path = G4FindDataDir("G4LEDATA");
110     G4ProductionCutsTable* theCoupleTable =       110     G4ProductionCutsTable* theCoupleTable =
111       G4ProductionCutsTable::GetProductionCuts    111       G4ProductionCutsTable::GetProductionCutsTable();
112     G4int numOfCouples = (G4int)theCoupleTable    112     G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
113                                                   113 
114     for(G4int i=0; i<numOfCouples; ++i)           114     for(G4int i=0; i<numOfCouples; ++i)
115       {                                           115       {
116   const G4MaterialCutsCouple* couple =            116   const G4MaterialCutsCouple* couple =
117     theCoupleTable->GetMaterialCutsCouple(i);     117     theCoupleTable->GetMaterialCutsCouple(i);
118   const G4Material* material = couple->GetMate    118   const G4Material* material = couple->GetMaterial();
119   const G4ElementVector* theElementVector = ma    119   const G4ElementVector* theElementVector = material->GetElementVector();
120   std::size_t nelm = material->GetNumberOfElem    120   std::size_t nelm = material->GetNumberOfElements();
121                                                   121 
122   for (std::size_t j=0; j<nelm; ++j)              122   for (std::size_t j=0; j<nelm; ++j)
123     {                                             123     {
124       G4int Z = G4lrint((*theElementVector)[j]    124       G4int Z = G4lrint((*theElementVector)[j]->GetZ());
125       if(Z < 1)          { Z = 1; }               125       if(Z < 1)          { Z = 1; }
126       else if(Z > maxZ)  { Z = maxZ; }            126       else if(Z > maxZ)  { Z = maxZ; }
127       if( (!dataCS[Z]) ) { ReadData(Z, path);     127       if( (!dataCS[Z]) ) { ReadData(Z, path); }
128     }                                             128     }
129       }                                           129       }
130   }                                               130   }
131                                                   131 
132   if(isInitialised) { return; }                   132   if(isInitialised) { return; }
133   fParticleChange = GetParticleChangeForGamma(    133   fParticleChange = GetParticleChangeForGamma();
134   isInitialised = true;                           134   isInitialised = true;
135                                                   135 
136 }                                                 136 }
137                                                   137 
138 //....oooOO0OOooo........oooOO0OOooo........oo    138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
139                                                   139 
140 void G4JAEAPolarizedElasticScatteringModel::In    140 void G4JAEAPolarizedElasticScatteringModel::InitialiseLocal(const G4ParticleDefinition*,
141                   G4VEmModel* masterModel)        141                   G4VEmModel* masterModel)
142 {                                                 142 {
143   SetElementSelectors(masterModel->GetElementS    143   SetElementSelectors(masterModel->GetElementSelectors());
144 }                                                 144 }
145                                                   145 
146 //....oooOO0OOooo........oooOO0OOooo........oo    146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
147                                                   147 
148 void G4JAEAPolarizedElasticScatteringModel::Re    148 void G4JAEAPolarizedElasticScatteringModel::ReadData(std::size_t Z, const char* path)
149 {                                                 149 {
150   if (verboseLevel > 1)                           150   if (verboseLevel > 1)
151     {                                             151     {
152       G4cout << "Calling ReadData() of G4JAEAP    152       G4cout << "Calling ReadData() of G4JAEAPolarizedElasticScatteringModel"
153        << G4endl;                                 153        << G4endl;
154     }                                             154     }
155                                                   155 
156   if(dataCS[Z]) { return; }                       156   if(dataCS[Z]) { return; }
157                                                   157 
158   const char* datadir = path;                     158   const char* datadir = path;
159   if(!datadir)                                    159   if(!datadir)
160     {                                             160     {
161       datadir = G4FindDataDir("G4LEDATA");        161       datadir = G4FindDataDir("G4LEDATA");
162       if(!datadir)                                162       if(!datadir)
163   {                                               163   {
164     G4Exception("G4JAEAPolarizedElasticScatter    164     G4Exception("G4JAEAPolarizedElasticScatteringModel::ReadData()","em0006",
165           FatalException,                         165           FatalException,
166           "Environment variable G4LEDATA not d    166           "Environment variable G4LEDATA not defined");
167     return;                                       167     return;
168   }                                               168   }
169     }                                             169     }
170                                                   170   
171   std::ostringstream ostCS;                       171   std::ostringstream ostCS;
172   ostCS << datadir << "/JAEAESData/amp_Z_" <<     172   ostCS << datadir << "/JAEAESData/amp_Z_" << Z ;
173   std::ifstream ES_Data_Buffer(ostCS.str().c_s    173   std::ifstream ES_Data_Buffer(ostCS.str().c_str(),ios::binary);
174   if( !ES_Data_Buffer.is_open() )                 174   if( !ES_Data_Buffer.is_open() )
175     {                                             175     {
176       G4ExceptionDescription ed;                  176       G4ExceptionDescription ed;
177       ed << "G4JAEAPolarizedElasticScattering     177       ed << "G4JAEAPolarizedElasticScattering Model data file <" << ostCS.str().c_str()
178    << "> is not opened!" << G4endl;               178    << "> is not opened!" << G4endl;
179       G4Exception("G4JAEAPolarizedElasticScatt    179       G4Exception("G4JAEAPolarizedElasticScatteringModel::ReadData()","em0003",FatalException,
180       ed,"G4LEDATA version should be G4EMLOW7.    180       ed,"G4LEDATA version should be G4EMLOW7.11 or later. Polarized Elastic Scattering Data are not loaded");
181       return;                                     181       return;
182     }                                             182     }
183   else                                            183   else
184     {                                             184     {
185       if(verboseLevel > 3) {                      185       if(verboseLevel > 3) {
186   G4cout << "File " << ostCS.str()                186   G4cout << "File " << ostCS.str()
187          << " is opened by G4JAEAPolarizedElas    187          << " is opened by G4JAEAPolarizedElasticScatteringModel" << G4endl;
188       }                                           188       }
189     }                                             189     }
190                                                   190   
191                                                   191   
192   if (!Polarized_ES_Data[Z])                      192   if (!Polarized_ES_Data[Z])
193     Polarized_ES_Data[Z] = new G4DataVector();    193     Polarized_ES_Data[Z] = new G4DataVector();
194                                                   194   
195   G4float buffer_var;                             195   G4float buffer_var;
196   while (ES_Data_Buffer.read(reinterpret_cast<    196   while (ES_Data_Buffer.read(reinterpret_cast<char*>(&buffer_var),sizeof(float)))
197     {                                             197     {
198       Polarized_ES_Data[Z]->push_back(buffer_v    198       Polarized_ES_Data[Z]->push_back(buffer_var);
199     }                                             199     }
200                                                   200   
201   dataCS[Z] = new G4PhysicsFreeVector(300,0.01    201   dataCS[Z] = new G4PhysicsFreeVector(300,0.01,3.,/*spline=*/true);
202                                                   202   
203   for (G4int i=0;i<300;++i)                       203   for (G4int i=0;i<300;++i)
204     dataCS[Z]->PutValues(i,10.*i*1e-3,Polarize    204     dataCS[Z]->PutValues(i,10.*i*1e-3,Polarized_ES_Data[Z]->at(i)*1e-22);
205                                                   205  
206   // Activation of spline interpolation           206   // Activation of spline interpolation
207   dataCS[Z] ->FillSecondDerivatives();            207   dataCS[Z] ->FillSecondDerivatives();
208 }                                                 208 }
209                                                   209 
210 //....oooOO0OOooo........oooOO0OOooo........oo    210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
211                                                   211 
212 G4double G4JAEAPolarizedElasticScatteringModel    212 G4double G4JAEAPolarizedElasticScatteringModel::ComputeCrossSectionPerAtom(
213                      const G4ParticleDefinitio    213                      const G4ParticleDefinition*,
214                      G4double GammaEnergy,        214                      G4double GammaEnergy,
215                      G4double Z, G4double,        215                      G4double Z, G4double,
216                      G4double, G4double)          216                      G4double, G4double)
217 {                                                 217 {
218   //Select the energy-grid point closest to th    218   //Select the energy-grid point closest to the photon energy
219   //    G4double *whichenergy = lower_bound(ES    219   //    G4double *whichenergy = lower_bound(ESdata[0],ESdata[0]+300,GammaEnergy);
220   //    int energyindex = max(0,(int)(whichene    220   //    int energyindex = max(0,(int)(whichenergy-ESdata[0]-1));
221                                                   221   
222   if (verboseLevel > 1)                           222   if (verboseLevel > 1)
223     {                                             223     {
224       G4cout << "G4JAEAPolarizedElasticScatter    224       G4cout << "G4JAEAPolarizedElasticScatteringModel::ComputeCrossSectionPerAtom()"
225        << G4endl;                                 225        << G4endl;
226     }                                             226     }
227                                                   227   
228   if(GammaEnergy < lowEnergyLimit) { return 0.    228   if(GammaEnergy < lowEnergyLimit) { return 0.0; }
229                                                   229   
230   G4double xs = 0.0;                              230   G4double xs = 0.0;
231                                                   231   
232   G4int intZ = G4lrint(Z);                        232   G4int intZ = G4lrint(Z);
233                                                   233   
234   if(intZ < 1 || intZ > maxZ) { return xs; }      234   if(intZ < 1 || intZ > maxZ) { return xs; }
235                                                   235 
236   G4PhysicsFreeVector* pv = dataCS[intZ];         236   G4PhysicsFreeVector* pv = dataCS[intZ];
237                                                   237 
238   // if element was not initialised               238   // if element was not initialised
239   // do initialisation safely for MT mode         239   // do initialisation safely for MT mode
240   if(!pv) {                                       240   if(!pv) {
241     InitialiseForElement(0, intZ);                241     InitialiseForElement(0, intZ);
242     pv = dataCS[intZ];                            242     pv = dataCS[intZ];
243     if(!pv) { return xs; }                        243     if(!pv) { return xs; }
244   }                                               244   }
245                                                   245 
246   std::size_t n = pv->GetVectorLength() - 1;      246   std::size_t n = pv->GetVectorLength() - 1;
247                                                   247 
248   G4double e = GammaEnergy;                       248   G4double e = GammaEnergy;
249   if(e >= pv->Energy(n)) {                        249   if(e >= pv->Energy(n)) {
250     xs = (*pv)[n];                                250     xs = (*pv)[n];
251   } else if(e >= pv->Energy(0)) {                 251   } else if(e >= pv->Energy(0)) {
252     xs = pv->Value(e);                            252     xs = pv->Value(e);
253   }                                               253   }
254                                                   254 
255   if(verboseLevel > 0)                            255   if(verboseLevel > 0)
256     {                                             256     {
257       G4cout  <<  "****** DEBUG: tcs value for    257       G4cout  <<  "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)="
258         << e << G4endl;                           258         << e << G4endl;
259       G4cout  <<  "  cs (Geant4 internal unit)    259       G4cout  <<  "  cs (Geant4 internal unit)=" << xs << G4endl;
260       G4cout  <<  "    -> first E*E*cs value i    260       G4cout  <<  "    -> first E*E*cs value in CS data file (iu) =" << (*pv)[0]
261         << G4endl;                                261         << G4endl;
262       G4cout  <<  "    -> last  E*E*cs value i    262       G4cout  <<  "    -> last  E*E*cs value in CS data file (iu) =" << (*pv)[n]
263         << G4endl;                                263         << G4endl;
264       G4cout  <<  "***************************    264       G4cout  <<  "*********************************************************"
265         << G4endl;                                265         << G4endl;
266     }                                             266     }
267                                                   267 
268   return (xs);                                    268   return (xs);
269 }                                                 269 }
270                                                   270 
271 //....oooOO0OOooo........oooOO0OOooo........oo    271 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
272                                                   272 
273 void G4JAEAPolarizedElasticScatteringModel::Sa    273 void G4JAEAPolarizedElasticScatteringModel::SampleSecondaries(
274                     std::vector<G4DynamicParti    274                     std::vector<G4DynamicParticle*>*,
275                     const G4MaterialCutsCouple    275                     const G4MaterialCutsCouple* couple,
276                     const G4DynamicParticle* a    276                     const G4DynamicParticle* aDynamicGamma,
277                     G4double, G4double)           277                     G4double, G4double)
278 {                                                 278 {
279   if (verboseLevel > 1) {                         279   if (verboseLevel > 1) {
280                                                   280 
281     G4cout << "Calling SampleSecondaries() of     281     G4cout << "Calling SampleSecondaries() of G4JAEAPolarizedElasticScatteringModel."
282      << G4endl;                                   282      << G4endl;
283   }                                               283   }
284   G4double photonEnergy0 = aDynamicGamma->GetK    284   G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
285                                                   285 
286   // absorption of low-energy gamma               286   // absorption of low-energy gamma
287   if (photonEnergy0 <= lowEnergyLimit)            287   if (photonEnergy0 <= lowEnergyLimit)
288     {                                             288     {
289       fParticleChange->ProposeTrackStatus(fSto    289       fParticleChange->ProposeTrackStatus(fStopAndKill);
290       fParticleChange->SetProposedKineticEnerg    290       fParticleChange->SetProposedKineticEnergy(0.);
291       fParticleChange->ProposeLocalEnergyDepos    291       fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
292       return ;                                    292       return ;
293     }                                             293     }
294                                                   294 
295   const G4ParticleDefinition* particle =  aDyn    295   const G4ParticleDefinition* particle =  aDynamicGamma->GetDefinition();
296   const G4Element* elm = SelectRandomAtom(coup    296   const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
297   G4int Z = G4lrint(elm->GetZ());                 297   G4int Z = G4lrint(elm->GetZ());
298                                                   298 
299   //Getting the corresponding distrbution         299   //Getting the corresponding distrbution
300   G4int energyindex=round(100*photonEnergy0)-1    300   G4int energyindex=round(100*photonEnergy0)-1;
301   G4double a1=0, a2=0, a3=0,a4=0;                 301   G4double a1=0, a2=0, a3=0,a4=0;
302   for (G4int i=0;i<=180;++i)                      302   for (G4int i=0;i<=180;++i)
303     {                                             303     {
304       a1=Polarized_ES_Data[Z]->at(4*i+300+181*    304       a1=Polarized_ES_Data[Z]->at(4*i+300+181*4*(energyindex));
305       a2=Polarized_ES_Data[Z]->at(4*i+1+300+18    305       a2=Polarized_ES_Data[Z]->at(4*i+1+300+181*4*(energyindex));
306       a3=Polarized_ES_Data[Z]->at(4*i+2+300+18    306       a3=Polarized_ES_Data[Z]->at(4*i+2+300+181*4*(energyindex));
307       a4=Polarized_ES_Data[Z]->at(4*i+3+300+18    307       a4=Polarized_ES_Data[Z]->at(4*i+3+300+181*4*(energyindex));
308       distribution[i]=a1*a1+a2*a2+a3*a3+a4*a4;    308       distribution[i]=a1*a1+a2*a2+a3*a3+a4*a4;
309     }                                             309     }
310                                                   310 
311   CLHEP::RandGeneral GenThetaDist(distribution    311   CLHEP::RandGeneral GenThetaDist(distribution,180);
312   //Intial sampling of the scattering angle. T    312   //Intial sampling of the scattering angle. To be updated for the circular polarization
313   G4double theta = CLHEP::pi*GenThetaDist.shoo    313   G4double theta = CLHEP::pi*GenThetaDist.shoot();
314   //G4double theta =45.*CLHEP::pi/180.;           314   //G4double theta =45.*CLHEP::pi/180.;
315   //Theta is in degree to call scattering ampl    315   //Theta is in degree to call scattering amplitudes
316   G4int theta_in_degree =round(theta*180./CLHE    316   G4int theta_in_degree =round(theta*180./CLHEP::pi);
317                                                   317 
318   //theta_in_degree=45;                           318   //theta_in_degree=45;
319                                                   319 
320   G4double am1=0,am2=0,am3=0,am4=0,aparaSquare    320   G4double am1=0,am2=0,am3=0,am4=0,aparaSquare=0,aperpSquare=0,
321     apara_aper_Asterisk=0,img_apara_aper_Aster    321     apara_aper_Asterisk=0,img_apara_aper_Asterisk=0;
322   am1=Polarized_ES_Data[Z]->at(4*theta_in_degr    322   am1=Polarized_ES_Data[Z]->at(4*theta_in_degree+300+181*4*(energyindex));
323   am2=Polarized_ES_Data[Z]->at(4*theta_in_degr    323   am2=Polarized_ES_Data[Z]->at(4*theta_in_degree+1+300+181*4*(energyindex));
324   am3=Polarized_ES_Data[Z]->at(4*theta_in_degr    324   am3=Polarized_ES_Data[Z]->at(4*theta_in_degree+2+300+181*4*(energyindex));
325   am4=Polarized_ES_Data[Z]->at(4*theta_in_degr    325   am4=Polarized_ES_Data[Z]->at(4*theta_in_degree+3+300+181*4*(energyindex));
326   aparaSquare=am1*am1+am2*am2;                    326   aparaSquare=am1*am1+am2*am2;
327   aperpSquare=am3*am3+am4*am4;                    327   aperpSquare=am3*am3+am4*am4;
328   apara_aper_Asterisk=2*a1*a3+2*a2*a4;            328   apara_aper_Asterisk=2*a1*a3+2*a2*a4;
329   img_apara_aper_Asterisk=2*a1*a4-2*a2*a3;        329   img_apara_aper_Asterisk=2*a1*a4-2*a2*a3;
330                                                   330 
331   G4ThreeVector Direction_Unpolarized(0.,0.,0.    331   G4ThreeVector Direction_Unpolarized(0.,0.,0.);
332   G4ThreeVector Direction_Linear1(0.,0.,0.);      332   G4ThreeVector Direction_Linear1(0.,0.,0.);
333   G4ThreeVector Direction_Linear2(0.,0.,0.);      333   G4ThreeVector Direction_Linear2(0.,0.,0.);
334   G4ThreeVector Direction_Circular(0.,0.,0.);     334   G4ThreeVector Direction_Circular(0.,0.,0.);
335   G4ThreeVector Polarization_Unpolarized(0.,0.    335   G4ThreeVector Polarization_Unpolarized(0.,0.,0.);
336   G4ThreeVector Polarization_Linear1(0.,0.,0.)    336   G4ThreeVector Polarization_Linear1(0.,0.,0.);
337   G4ThreeVector Polarization_Linear2(0.,0.,0.)    337   G4ThreeVector Polarization_Linear2(0.,0.,0.);
338   G4ThreeVector Polarization_Circular(0.,0.,0.    338   G4ThreeVector Polarization_Circular(0.,0.,0.);
339                                                   339 
340   //Stokes parameters for the incoming and out    340   //Stokes parameters for the incoming and outgoing photon
341   G4double Xi1=0, Xi2=0, Xi3=0, Xi1_Prime=0,Xi    341   G4double Xi1=0, Xi2=0, Xi3=0, Xi1_Prime=0,Xi2_Prime=0,Xi3_Prime=0;
342                                                   342     
343   //Getting the Stokes parameters for the inco    343   //Getting the Stokes parameters for the incoming photon
344   G4ThreeVector gammaPolarization0 = aDynamicG    344   G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
345   Xi1=gammaPolarization0.x();                     345   Xi1=gammaPolarization0.x();
346   Xi2=gammaPolarization0.y();                     346   Xi2=gammaPolarization0.y();
347   Xi3=gammaPolarization0.z();                     347   Xi3=gammaPolarization0.z();
348                                                   348     
349   //Polarization vector must be unit vector (5    349   //Polarization vector must be unit vector (5% tolerance)    
350   if ((gammaPolarization0.mag())>1.05 || (Xi1*    350   if ((gammaPolarization0.mag())>1.05 || (Xi1*Xi1>1.05) || (Xi2*Xi2>1.05) || (Xi3*Xi3>1.05))
351     {                                             351     {
352       G4Exception("G4JAEAPolarizedElasticScatt    352       G4Exception("G4JAEAPolarizedElasticScatteringModel::SampleSecondaries()","em1006",
353       JustWarning,                                353       JustWarning,
354       "WARNING: G4JAEAPolarizedElasticScatteri    354       "WARNING: G4JAEAPolarizedElasticScatteringModel is only compatible with a unit polarization vector.");
355       return;                                     355       return;
356     }                                             356     }
357   //Unpolarized gamma rays                        357   //Unpolarized gamma rays
358   if (Xi1==0 && Xi2==0 && Xi3==0)                 358   if (Xi1==0 && Xi2==0 && Xi3==0)
359     {                                             359     {
360       G4double Phi_Unpolarized=0;                 360       G4double Phi_Unpolarized=0;
361       if (fLinearPolarizationSensitvity1)         361       if (fLinearPolarizationSensitvity1)
362   Phi_Unpolarized=GeneratePolarizedPhi(aparaSq    362   Phi_Unpolarized=GeneratePolarizedPhi(aparaSquare,aperpSquare,0.);
363       else                                        363       else
364   Phi_Unpolarized=CLHEP::twopi*G4UniformRand()    364   Phi_Unpolarized=CLHEP::twopi*G4UniformRand();
365       Direction_Unpolarized.setX(sin(theta)*co    365       Direction_Unpolarized.setX(sin(theta)*cos(Phi_Unpolarized));
366       Direction_Unpolarized.setY(sin(theta)*si    366       Direction_Unpolarized.setY(sin(theta)*sin(Phi_Unpolarized));
367       Direction_Unpolarized.setZ(cos(theta));     367       Direction_Unpolarized.setZ(cos(theta));
368       Direction_Unpolarized.rotateUz(aDynamicG    368       Direction_Unpolarized.rotateUz(aDynamicGamma->GetMomentumDirection());
369       Xi1_Prime=(aparaSquare-aperpSquare)/(apa    369       Xi1_Prime=(aparaSquare-aperpSquare)/(aparaSquare+aperpSquare);
370       Polarization_Unpolarized.setX(Xi1_Prime)    370       Polarization_Unpolarized.setX(Xi1_Prime);
371       Polarization_Unpolarized.setY(0.);          371       Polarization_Unpolarized.setY(0.);
372       Polarization_Unpolarized.setZ(0.);          372       Polarization_Unpolarized.setZ(0.);
373       fParticleChange->ProposeMomentumDirectio    373       fParticleChange->ProposeMomentumDirection(Direction_Unpolarized);
374       fParticleChange->ProposePolarization(Pol    374       fParticleChange->ProposePolarization(Polarization_Unpolarized);
375       return;                                     375       return;
376     }                                             376     }
377                                                   377 
378   //Linear polarization defined by first Stoke    378   //Linear polarization defined by first Stokes parameter
379   G4double InitialAzimuth=aDynamicGamma->GetMo    379   G4double InitialAzimuth=aDynamicGamma->GetMomentumDirection().phi();
380   if(InitialAzimuth<0) InitialAzimuth=InitialA    380   if(InitialAzimuth<0) InitialAzimuth=InitialAzimuth+CLHEP::twopi;
381                                                   381     
382   G4double Phi_Linear1=0.;                        382   G4double Phi_Linear1=0.;
383                                                   383 
384   Phi_Linear1 = GeneratePolarizedPhi(aparaSqua    384   Phi_Linear1 = GeneratePolarizedPhi(aparaSquare+aperpSquare+Xi1*(aparaSquare-aperpSquare),
385              aparaSquare+aperpSquare-Xi1*(apar    385              aparaSquare+aperpSquare-Xi1*(aparaSquare-aperpSquare),InitialAzimuth);
386                                                   386 
387   Xi1_Prime=((aparaSquare-aperpSquare)+Xi1*(ap    387   Xi1_Prime=((aparaSquare-aperpSquare)+Xi1*(aparaSquare+aperpSquare)*cos(2*Phi_Linear1))/
388     ((aparaSquare+aperpSquare)+Xi1*(aparaSquar    388     ((aparaSquare+aperpSquare)+Xi1*(aparaSquare-aperpSquare)*cos(2*Phi_Linear1));
389   Xi2_Prime=(-Xi1*apara_aper_Asterisk*sin(2*Ph    389   Xi2_Prime=(-Xi1*apara_aper_Asterisk*sin(2*Phi_Linear1))/
390     ((aparaSquare+aperpSquare)+Xi1*(aparaSquar    390     ((aparaSquare+aperpSquare)+Xi1*(aparaSquare-aperpSquare)*cos(2*Phi_Linear1));
391   Xi3_Prime=(-Xi1*img_apara_aper_Asterisk*sin(    391   Xi3_Prime=(-Xi1*img_apara_aper_Asterisk*sin(2*Phi_Linear1))/
392     ((aparaSquare+aperpSquare)+Xi1*(aparaSquar    392     ((aparaSquare+aperpSquare)+Xi1*(aparaSquare-aperpSquare)*cos(2*Phi_Linear1));
393   //Store momentum direction and po;arization     393   //Store momentum direction and po;arization
394   Direction_Linear1.setX(sin(theta)*cos(Phi_Li    394   Direction_Linear1.setX(sin(theta)*cos(Phi_Linear1));
395   Direction_Linear1.setY(sin(theta)*sin(Phi_Li    395   Direction_Linear1.setY(sin(theta)*sin(Phi_Linear1));
396   Direction_Linear1.setZ(cos(theta));             396   Direction_Linear1.setZ(cos(theta));
397   Polarization_Linear1.setX(Xi1_Prime);           397   Polarization_Linear1.setX(Xi1_Prime);
398   Polarization_Linear1.setY(Xi2_Prime);           398   Polarization_Linear1.setY(Xi2_Prime);
399   Polarization_Linear1.setZ(Xi3_Prime);           399   Polarization_Linear1.setZ(Xi3_Prime);
400                                                   400     
401   //Set scattered photon polarization sensitiv    401   //Set scattered photon polarization sensitivity
402   Xi1_Prime=Xi1_Prime*fLinearPolarizationSensi    402   Xi1_Prime=Xi1_Prime*fLinearPolarizationSensitvity1;
403   Xi2_Prime=Xi2_Prime*fLinearPolarizationSensi    403   Xi2_Prime=Xi2_Prime*fLinearPolarizationSensitvity2;
404   Xi3_Prime=Xi3_Prime*fCircularPolarizationSen    404   Xi3_Prime=Xi3_Prime*fCircularPolarizationSensitvity;
405                                                   405     
406   G4double dsigmaL1=0.0;                          406   G4double dsigmaL1=0.0;
407   if(abs(Xi1)>0.0) dsigmaL1=0.25*((aparaSquare    407   if(abs(Xi1)>0.0) dsigmaL1=0.25*((aparaSquare+aperpSquare)*
408           (1+Xi1*Xi1_Prime*cos(2*Phi_Linear1))    408           (1+Xi1*Xi1_Prime*cos(2*Phi_Linear1))+
409           (aparaSquare-aperpSquare)*(Xi1*cos(2    409           (aparaSquare-aperpSquare)*(Xi1*cos(2*Phi_Linear1)+Xi1_Prime)
410           -Xi1*Xi2_Prime*apara_aper_Asterisk*s    410           -Xi1*Xi2_Prime*apara_aper_Asterisk*sin(2*Phi_Linear1)-
411           Xi1*Xi3_Prime*img_apara_aper_Asteris    411           Xi1*Xi3_Prime*img_apara_aper_Asterisk*sin(2*Phi_Linear1));
412                                                   412 
413   //Linear polarization defined by second Stok    413   //Linear polarization defined by second Stokes parameter
414   //G4double IntialAzimuth=aDynamicGamma->GetM    414   //G4double IntialAzimuth=aDynamicGamma->GetMomentumDirection().phi();
415   G4double Phi_Linear2=0.;                        415   G4double Phi_Linear2=0.;
416                                                   416 
417   InitialAzimuth=InitialAzimuth-CLHEP::pi/4.;     417   InitialAzimuth=InitialAzimuth-CLHEP::pi/4.;
418   if(InitialAzimuth<0) InitialAzimuth=InitialA    418   if(InitialAzimuth<0) InitialAzimuth=InitialAzimuth+CLHEP::twopi;
419                                                   419 
420   Phi_Linear2 = GeneratePolarizedPhi(aparaSqua    420   Phi_Linear2 = GeneratePolarizedPhi(aparaSquare+aperpSquare+Xi1*(aparaSquare-aperpSquare)
421              ,aparaSquare+aperpSquare-Xi1*(apa    421              ,aparaSquare+aperpSquare-Xi1*(aparaSquare-aperpSquare),InitialAzimuth);
422                                                   422     
423   Xi1_Prime=((aparaSquare-aperpSquare)+Xi2*(ap    423   Xi1_Prime=((aparaSquare-aperpSquare)+Xi2*(aparaSquare+aperpSquare)*sin(2*Phi_Linear2))/
424     ((aparaSquare+aperpSquare)+Xi2*(aparaSquar    424     ((aparaSquare+aperpSquare)+Xi2*(aparaSquare-aperpSquare)*sin(2*Phi_Linear2));
425   Xi2_Prime=(Xi2*apara_aper_Asterisk*cos(2*Phi    425   Xi2_Prime=(Xi2*apara_aper_Asterisk*cos(2*Phi_Linear2))/
426     ((aparaSquare+aperpSquare)+Xi2*(aparaSquar    426     ((aparaSquare+aperpSquare)+Xi2*(aparaSquare-aperpSquare)*sin(2*Phi_Linear2));
427   Xi3_Prime=(Xi2*img_apara_aper_Asterisk*cos(2    427   Xi3_Prime=(Xi2*img_apara_aper_Asterisk*cos(2*Phi_Linear2))/
428     ((aparaSquare+aperpSquare)+Xi2*(aparaSquar    428     ((aparaSquare+aperpSquare)+Xi2*(aparaSquare-aperpSquare)*sin(2*Phi_Linear2));
429   //Store momentum direction and polarization     429   //Store momentum direction and polarization
430   Direction_Linear2.setX(sin(theta)*cos(Phi_Li    430   Direction_Linear2.setX(sin(theta)*cos(Phi_Linear2));
431   Direction_Linear2.setY(sin(theta)*sin(Phi_Li    431   Direction_Linear2.setY(sin(theta)*sin(Phi_Linear2));
432   Direction_Linear2.setZ(cos(theta));             432   Direction_Linear2.setZ(cos(theta));
433   Polarization_Linear2.setX(Xi1_Prime);           433   Polarization_Linear2.setX(Xi1_Prime);
434   Polarization_Linear2.setY(Xi2_Prime);           434   Polarization_Linear2.setY(Xi2_Prime);
435   Polarization_Linear2.setZ(Xi3_Prime);           435   Polarization_Linear2.setZ(Xi3_Prime);
436                                                   436 
437   //Set scattered photon polarization sensitiv    437   //Set scattered photon polarization sensitivity
438   Xi1_Prime=Xi1_Prime*fLinearPolarizationSensi    438   Xi1_Prime=Xi1_Prime*fLinearPolarizationSensitvity1;
439   Xi2_Prime=Xi2_Prime*fLinearPolarizationSensi    439   Xi2_Prime=Xi2_Prime*fLinearPolarizationSensitvity2;
440   Xi3_Prime=Xi3_Prime*fCircularPolarizationSen    440   Xi3_Prime=Xi3_Prime*fCircularPolarizationSensitvity;
441                                                   441 
442   G4double dsigmaL2=0.0;                          442   G4double dsigmaL2=0.0;
443   if(abs(Xi2)>0.0)                                443   if(abs(Xi2)>0.0)
444     dsigmaL2=0.25*((aparaSquare+aperpSquare)*(    444     dsigmaL2=0.25*((aparaSquare+aperpSquare)*(1+Xi2*Xi1_Prime*sin(2*Phi_Linear2))+
445        (aparaSquare-aperpSquare)*(Xi2*sin(2*Ph    445        (aparaSquare-aperpSquare)*(Xi2*sin(2*Phi_Linear2)+Xi1_Prime)
446        +Xi2*Xi2_Prime*apara_aper_Asterisk*cos(    446        +Xi2*Xi2_Prime*apara_aper_Asterisk*cos(2*Phi_Linear2)-
447        Xi2*Xi3_Prime*img_apara_aper_Asterisk*c    447        Xi2*Xi3_Prime*img_apara_aper_Asterisk*cos(2*Phi_Linear2));
448                                                   448 
449   //Circular polarization                         449   //Circular polarization
450   G4double Phi_Circular = CLHEP::twopi*G4Unifo    450   G4double Phi_Circular = CLHEP::twopi*G4UniformRand();
451   G4double Theta_Circular = 0.;                   451   G4double Theta_Circular = 0.;
452                                                   452     
453   Xi1_Prime=(aparaSquare-aperpSquare)/(aparaSq    453   Xi1_Prime=(aparaSquare-aperpSquare)/(aparaSquare+aperpSquare);
454   Xi2_Prime=(-Xi3*img_apara_aper_Asterisk)/(ap    454   Xi2_Prime=(-Xi3*img_apara_aper_Asterisk)/(aparaSquare+aperpSquare);
455   Xi3_Prime=(Xi3*apara_aper_Asterisk)/(aparaSq    455   Xi3_Prime=(Xi3*apara_aper_Asterisk)/(aparaSquare+aperpSquare);
456                                                   456     
457   Polarization_Circular.setX(Xi1_Prime);          457   Polarization_Circular.setX(Xi1_Prime);
458   Polarization_Circular.setY(Xi2_Prime);          458   Polarization_Circular.setY(Xi2_Prime);
459   Polarization_Circular.setZ(Xi3_Prime);          459   Polarization_Circular.setZ(Xi3_Prime);
460                                                   460     
461   //Set scattered photon polarization sensitiv    461   //Set scattered photon polarization sensitivity
462   Xi1_Prime=Xi1_Prime*fLinearPolarizationSensi    462   Xi1_Prime=Xi1_Prime*fLinearPolarizationSensitvity1;
463   Xi2_Prime=Xi2_Prime*fLinearPolarizationSensi    463   Xi2_Prime=Xi2_Prime*fLinearPolarizationSensitvity2;
464   Xi3_Prime=Xi3_Prime*fCircularPolarizationSen    464   Xi3_Prime=Xi3_Prime*fCircularPolarizationSensitvity;
465                                                   465     
466   G4double dsigmaC=0.0;                           466   G4double dsigmaC=0.0;
467   if(abs(Xi3)>0.0)                                467   if(abs(Xi3)>0.0)
468     dsigmaC=0.25*(aparaSquare+aperpSquare+Xi1_    468     dsigmaC=0.25*(aparaSquare+aperpSquare+Xi1_Prime*(aparaSquare-aperpSquare)-
469       Xi3*Xi2_Prime*img_apara_aper_Asterisk       469       Xi3*Xi2_Prime*img_apara_aper_Asterisk
470       +Xi3*Xi3_Prime*apara_aper_Asterisk);        470       +Xi3*Xi3_Prime*apara_aper_Asterisk);
471                                                   471 
472   if (abs(Xi3)==0.0 && abs(Xi1_Prime)==0.0)       472   if (abs(Xi3)==0.0 && abs(Xi1_Prime)==0.0)
473     {                                             473     {
474       Direction_Circular.setX(sin(theta)*cos(P    474       Direction_Circular.setX(sin(theta)*cos(Phi_Circular));
475       Direction_Circular.setY(sin(theta)*sin(P    475       Direction_Circular.setY(sin(theta)*sin(Phi_Circular));
476       Direction_Circular.setZ(cos(theta));        476       Direction_Circular.setZ(cos(theta));
477     }                                             477     }
478   else                                            478   else
479     {                                             479     {
480       G4double c1=0, c2=0, c3=0,c4=0;             480       G4double c1=0, c2=0, c3=0,c4=0;
481       for (G4int i=0;i<=180;++i)                  481       for (G4int i=0;i<=180;++i)
482   {                                               482   {
483     c1=Polarized_ES_Data[Z]->at(4*i+300+181*4*    483     c1=Polarized_ES_Data[Z]->at(4*i+300+181*4*(energyindex));
484     c2=Polarized_ES_Data[Z]->at(4*i+1+300+181*    484     c2=Polarized_ES_Data[Z]->at(4*i+1+300+181*4*(energyindex));
485     c3=Polarized_ES_Data[Z]->at(4*i+2+300+181*    485     c3=Polarized_ES_Data[Z]->at(4*i+2+300+181*4*(energyindex));
486     c4=Polarized_ES_Data[Z]->at(4*i+3+300+181*    486     c4=Polarized_ES_Data[Z]->at(4*i+3+300+181*4*(energyindex));
487     cdistribution[i]=0.25*((c1*c1+c2*c2+c3*c3+    487     cdistribution[i]=0.25*((c1*c1+c2*c2+c3*c3+c4*c4)+
488          Xi1_Prime*(c1*c1+c2*c2-c3*c3-c4*c4)-     488          Xi1_Prime*(c1*c1+c2*c2-c3*c3-c4*c4)-
489          Xi3*Xi2_Prime*(2*c1*c4-2*c2*c3)          489          Xi3*Xi2_Prime*(2*c1*c4-2*c2*c3)
490          +Xi3*Xi3_Prime*(2*c1*c4-2*c2*c3));       490          +Xi3*Xi3_Prime*(2*c1*c4-2*c2*c3));
491   }                                               491   }
492       CLHEP::RandGeneral GenTheta_Circ_Dist(cd    492       CLHEP::RandGeneral GenTheta_Circ_Dist(cdistribution,180);
493       Theta_Circular=CLHEP::pi*GenTheta_Circ_D    493       Theta_Circular=CLHEP::pi*GenTheta_Circ_Dist.shoot();
494       Direction_Circular.setX(sin(Theta_Circul    494       Direction_Circular.setX(sin(Theta_Circular)*cos(Phi_Circular));
495       Direction_Circular.setY(sin(Theta_Circul    495       Direction_Circular.setY(sin(Theta_Circular)*sin(Phi_Circular));
496       Direction_Circular.setZ(cos(Theta_Circul    496       Direction_Circular.setZ(cos(Theta_Circular));
497     }                                             497     }
498                                                   498 
499   // Sampling scattered photon direction based    499   // Sampling scattered photon direction based on asymmetry arising from polarization mixing
500   G4double totalSigma= dsigmaL1+dsigmaL2+dsigm    500   G4double totalSigma= dsigmaL1+dsigmaL2+dsigmaC;
501   G4double prob1=dsigmaL1/totalSigma;             501   G4double prob1=dsigmaL1/totalSigma;
502   G4double prob2=dsigmaL2/totalSigma;             502   G4double prob2=dsigmaL2/totalSigma;
503   G4double probc=1-(prob1+prob2);                 503   G4double probc=1-(prob1+prob2);
504                                                   504     
505   //Check the Probability of polarization mixi    505   //Check the Probability of polarization mixing
506   if (abs(probc - dsigmaC/totalSigma)>=0.0001)    506   if (abs(probc - dsigmaC/totalSigma)>=0.0001)
507     {                                             507     {
508       G4Exception("G4JAEAPolarizedElasticScatt    508       G4Exception("G4JAEAPolarizedElasticScatteringModel::SampleSecondaries()","em1007",
509       JustWarning,                                509       JustWarning,
510       "WARNING: Polarization mixing might be i    510       "WARNING: Polarization mixing might be incorrect.");
511     }                                             511     }
512                                                   512     
513   // Generate outgoing photon direction           513   // Generate outgoing photon direction
514   G4ThreeVector finaldirection(0.0,0.0,0.0);      514   G4ThreeVector finaldirection(0.0,0.0,0.0);
515   G4ThreeVector outcomingPhotonPolarization(0.    515   G4ThreeVector outcomingPhotonPolarization(0.0,0.0,0.0);
516                                                   516     
517   //Polarization mixing                           517   //Polarization mixing
518   G4double polmix=G4UniformRand();                518   G4double polmix=G4UniformRand();
519   if (polmix<=prob1)                              519   if (polmix<=prob1)
520     {                                             520     {
521       finaldirection.setX(Direction_Linear1.x(    521       finaldirection.setX(Direction_Linear1.x());
522       finaldirection.setY(Direction_Linear1.y(    522       finaldirection.setY(Direction_Linear1.y());
523       finaldirection.setZ(Direction_Linear1.z(    523       finaldirection.setZ(Direction_Linear1.z());
524       outcomingPhotonPolarization.setX(Polariz    524       outcomingPhotonPolarization.setX(Polarization_Linear1.x());
525       outcomingPhotonPolarization.setY(Polariz    525       outcomingPhotonPolarization.setY(Polarization_Linear1.y());
526       outcomingPhotonPolarization.setZ(Polariz    526       outcomingPhotonPolarization.setZ(Polarization_Linear1.z());
527     }                                             527     }
528   else if ((polmix>prob1) && (polmix<=prob1+pr    528   else if ((polmix>prob1) && (polmix<=prob1+prob2))
529     {                                             529     {
530       finaldirection.setX(Direction_Linear2.x(    530       finaldirection.setX(Direction_Linear2.x());
531       finaldirection.setY(Direction_Linear2.y(    531       finaldirection.setY(Direction_Linear2.y());
532       finaldirection.setZ(Direction_Linear2.z(    532       finaldirection.setZ(Direction_Linear2.z());
533       outcomingPhotonPolarization.setX(Polariz    533       outcomingPhotonPolarization.setX(Polarization_Linear2.x());
534       outcomingPhotonPolarization.setY(Polariz    534       outcomingPhotonPolarization.setY(Polarization_Linear2.y());
535       outcomingPhotonPolarization.setZ(Polariz    535       outcomingPhotonPolarization.setZ(Polarization_Linear2.z());
536     }                                             536     }
537   else if (polmix>prob1+prob2)                    537   else if (polmix>prob1+prob2)
538     {                                             538     {
539       finaldirection.setX(Direction_Circular.x    539       finaldirection.setX(Direction_Circular.x());
540       finaldirection.setY(Direction_Circular.y    540       finaldirection.setY(Direction_Circular.y());
541       finaldirection.setZ(Direction_Circular.z    541       finaldirection.setZ(Direction_Circular.z());
542       outcomingPhotonPolarization.setX(Polariz    542       outcomingPhotonPolarization.setX(Polarization_Circular.x());
543       outcomingPhotonPolarization.setY(Polariz    543       outcomingPhotonPolarization.setY(Polarization_Circular.y());
544       outcomingPhotonPolarization.setZ(Polariz    544       outcomingPhotonPolarization.setZ(Polarization_Circular.z());
545     }                                             545     }
546                                                   546 
547   //Sampling the Final State                      547   //Sampling the Final State
548   finaldirection.rotateUz(aDynamicGamma->GetMo    548   finaldirection.rotateUz(aDynamicGamma->GetMomentumDirection());
549   fParticleChange->ProposeMomentumDirection(fi    549   fParticleChange->ProposeMomentumDirection(finaldirection);
550   fParticleChange->SetProposedKineticEnergy(ph    550   fParticleChange->SetProposedKineticEnergy(photonEnergy0);
551   fParticleChange->ProposePolarization(outcomi    551   fParticleChange->ProposePolarization(outcomingPhotonPolarization);
552                                                   552     
553 }                                                 553 }
554                                                   554 
555 //....oooOO0OOooo........oooOO0OOooo........oo    555 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
556                                                   556 
557 G4double G4JAEAPolarizedElasticScatteringModel    557 G4double G4JAEAPolarizedElasticScatteringModel::GeneratePolarizedPhi(G4double Sigma_para,
558                      G4double Sigma_perp,         558                      G4double Sigma_perp, 
559                      G4double initial_Pol_Plan    559                      G4double initial_Pol_Plane)
560 {                                                 560 {
561   G4double phi;                                   561   G4double phi;
562   G4double phiProbability;                        562   G4double phiProbability;
563   G4double Probability=Sigma_perp/(Sigma_para+    563   G4double Probability=Sigma_perp/(Sigma_para+Sigma_perp);
564   if (Probability<=G4UniformRand())               564   if (Probability<=G4UniformRand())
565     {                                             565     {
566       do                                          566       do
567   {                                               567   {
568     phi = CLHEP::twopi * G4UniformRand();         568     phi = CLHEP::twopi * G4UniformRand();
569     phiProbability = cos(phi+initial_Pol_Plane    569     phiProbability = cos(phi+initial_Pol_Plane)*cos(phi+initial_Pol_Plane);
570   }                                               570   }
571       while (phiProbability < G4UniformRand())    571       while (phiProbability < G4UniformRand());
572                                                   572 
573     }                                             573     }
574   else                                            574   else
575     {                                             575     {
576       do                                          576       do
577   {                                               577   {
578     phi = CLHEP::twopi * G4UniformRand();         578     phi = CLHEP::twopi * G4UniformRand();
579     phiProbability = sin(phi+initial_Pol_Plane    579     phiProbability = sin(phi+initial_Pol_Plane)*sin(phi+initial_Pol_Plane);
580   }                                               580   }
581       while (phiProbability < G4UniformRand())    581       while (phiProbability < G4UniformRand());
582     }                                             582     }
583   return phi;                                     583   return phi;
584                                                   584 
585 }                                                 585 }
586 //....oooOO0OOooo........oooOO0OOooo........oo    586 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
587                                                   587 
588 void                                              588 void
589 G4JAEAPolarizedElasticScatteringModel::Initial    589 G4JAEAPolarizedElasticScatteringModel::InitialiseForElement(const G4ParticleDefinition*,
590                   G4int Z)                        590                   G4int Z)
591 {                                                 591 {
592   G4AutoLock l(&G4JAEAPolarizedElasticScatteri    592   G4AutoLock l(&G4JAEAPolarizedElasticScatteringModelMutex);
593   //  G4cout << "G4JAEAPolarizedElasticScatter    593   //  G4cout << "G4JAEAPolarizedElasticScatteringModel::InitialiseForElement Z= "
594   //   << Z << G4endl;                            594   //   << Z << G4endl;
595   if(!dataCS[Z]) { ReadData(Z); }                 595   if(!dataCS[Z]) { ReadData(Z); }
596   l.unlock();                                     596   l.unlock();
597 }                                                 597 }
598                                                   598 
599 //....oooOO0OOooo........oooOO0OOooo........oo    599 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
600                                                   600