Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/eRosita/physics/src/G4LowEnergyPolarizedCompton.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 /examples/advanced/eRosita/physics/src/G4LowEnergyPolarizedCompton.cc (Version 11.3.0) and /examples/advanced/eRosita/physics/src/G4LowEnergyPolarizedCompton.cc (Version 10.1.p1)


  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 // $Id$
                                                   >>  28 // GEANT4 tag $Name:  $
 27 //                                                 29 //
 28 // -------------------------------------------     30 // ------------------------------------------------------------
 29 //      GEANT 4 class implementation file          31 //      GEANT 4 class implementation file
 30 //      CERN Geneva Switzerland                    32 //      CERN Geneva Switzerland
 31 //                                                 33 //
 32                                                    34 
 33 // --------- G4LowEnergyPolarizedCompton class     35 // --------- G4LowEnergyPolarizedCompton class -----
 34 //                                                 36 //
 35 //           by G.Depaola & F.Longo (21 may 20     37 //           by G.Depaola & F.Longo (21 may 2001)
 36 //                                                 38 //
 37 // 21 May 2001 - MGP      Modified to inherit      39 // 21 May 2001 - MGP      Modified to inherit from G4VDiscreteProcess
 38 //                        Applies same algorit     40 //                        Applies same algorithm as LowEnergyCompton
 39 //                        if the incoming phot     41 //                        if the incoming photon is not polarised
 40 //                        Temporary protection     42 //                        Temporary protection to avoid crash in the case 
 41 //                        of polarisation || i     43 //                        of polarisation || incident photon direction
 42 //                                                 44 //
 43 // 17 October 2001 - F.Longo - Revised accordi     45 // 17 October 2001 - F.Longo - Revised according to a design iteration
 44 //                                                 46 //
 45 // 21 February 2002 - F.Longo Revisions with A     47 // 21 February 2002 - F.Longo Revisions with A.Zoglauer and G.Depaola
 46 //                            - better descrip     48 //                            - better description of parallelism
 47 //                            - system of ref      49 //                            - system of ref change method improved
 48 // 22 January  2003 - V.Ivanchenko Cut per reg     50 // 22 January  2003 - V.Ivanchenko Cut per region
 49 // 24 April    2003 - V.Ivanchenko Cut per reg     51 // 24 April    2003 - V.Ivanchenko Cut per region mfpt
 50 //                                                 52 //
 51 //                                                 53 //
 52 // *******************************************     54 // ************************************************************
 53 //                                                 55 //
 54 // Corrections by Rui Curado da Silva (2000)       56 // Corrections by Rui Curado da Silva (2000)
 55 // New Implementation by G.Depaola & F.Longo       57 // New Implementation by G.Depaola & F.Longo
 56 //                                                 58 //
 57 // - sampling of phi                               59 // - sampling of phi
 58 // - polarization of scattered photon              60 // - polarization of scattered photon
 59 //                                                 61 //
 60 // -------------------------------------------     62 // --------------------------------------------------------------
 61                                                    63 
 62 #include "G4LowEnergyPolarizedCompton.hh"          64 #include "G4LowEnergyPolarizedCompton.hh"
 63 #include "Randomize.hh"                            65 #include "Randomize.hh"
 64 #include "G4PhysicalConstants.hh"                  66 #include "G4PhysicalConstants.hh"
 65 #include "G4SystemOfUnits.hh"                      67 #include "G4SystemOfUnits.hh"
 66 #include "G4ParticleDefinition.hh"                 68 #include "G4ParticleDefinition.hh"
 67 #include "G4Track.hh"                              69 #include "G4Track.hh"
 68 #include "G4Step.hh"                               70 #include "G4Step.hh"
 69 #include "G4ForceCondition.hh"                     71 #include "G4ForceCondition.hh"
 70 #include "G4Gamma.hh"                              72 #include "G4Gamma.hh"
 71 #include "G4Electron.hh"                           73 #include "G4Electron.hh"
 72 #include "G4DynamicParticle.hh"                    74 #include "G4DynamicParticle.hh"
 73 #include "G4VParticleChange.hh"                    75 #include "G4VParticleChange.hh"
 74 #include "G4ThreeVector.hh"                        76 #include "G4ThreeVector.hh"
 75 #include "G4RDVCrossSectionHandler.hh"             77 #include "G4RDVCrossSectionHandler.hh"
 76 #include "G4RDCrossSectionHandler.hh"              78 #include "G4RDCrossSectionHandler.hh"
 77 #include "G4RDVEMDataSet.hh"                       79 #include "G4RDVEMDataSet.hh"
 78 #include "G4RDCompositeEMDataSet.hh"               80 #include "G4RDCompositeEMDataSet.hh"
 79 #include "G4RDVDataSetAlgorithm.hh"                81 #include "G4RDVDataSetAlgorithm.hh"
 80 #include "G4RDLogLogInterpolation.hh"              82 #include "G4RDLogLogInterpolation.hh"
 81 #include "G4RDVRangeTest.hh"                       83 #include "G4RDVRangeTest.hh"
 82 #include "G4RDRangeTest.hh"                        84 #include "G4RDRangeTest.hh"
 83 #include "G4MaterialCutsCouple.hh"                 85 #include "G4MaterialCutsCouple.hh"
 84                                                    86 
 85 // constructor                                     87 // constructor
 86                                                    88 
 87 G4LowEnergyPolarizedCompton::G4LowEnergyPolari     89 G4LowEnergyPolarizedCompton::G4LowEnergyPolarizedCompton(const G4String& processName)
 88   : G4VDiscreteProcess(processName),               90   : G4VDiscreteProcess(processName),
 89     lowEnergyLimit (250*eV),              // i     91     lowEnergyLimit (250*eV),              // initialization
 90     highEnergyLimit(100*GeV),                      92     highEnergyLimit(100*GeV),
 91     intrinsicLowEnergyLimit(10*eV),                93     intrinsicLowEnergyLimit(10*eV),
 92     intrinsicHighEnergyLimit(100*GeV)              94     intrinsicHighEnergyLimit(100*GeV)
 93 {                                                  95 {
 94   if (lowEnergyLimit < intrinsicLowEnergyLimit     96   if (lowEnergyLimit < intrinsicLowEnergyLimit ||
 95       highEnergyLimit > intrinsicHighEnergyLim     97       highEnergyLimit > intrinsicHighEnergyLimit)
 96     {                                              98     {
 97       G4Exception("G4LowEnergyPolarizedCompton     99       G4Exception("G4LowEnergyPolarizedCompton::G4LowEnergyPolarizedCompton()",
 98                   "OutOfRange", FatalException    100                   "OutOfRange", FatalException,
 99                   "Energy outside intrinsic pr    101                   "Energy outside intrinsic process validity range!");
100     }                                             102     }
101                                                   103 
102   crossSectionHandler = new G4RDCrossSectionHa    104   crossSectionHandler = new G4RDCrossSectionHandler;
103                                                   105 
104                                                   106 
105   G4RDVDataSetAlgorithm* scatterInterpolation     107   G4RDVDataSetAlgorithm* scatterInterpolation = new G4RDLogLogInterpolation;
106   G4String scatterFile = "comp/ce-sf-";           108   G4String scatterFile = "comp/ce-sf-";
107   scatterFunctionData = new G4RDCompositeEMDat    109   scatterFunctionData = new G4RDCompositeEMDataSet(scatterInterpolation,1.,1.);
108   scatterFunctionData->LoadData(scatterFile);     110   scatterFunctionData->LoadData(scatterFile);
109                                                   111 
110   meanFreePathTable = 0;                          112   meanFreePathTable = 0;
111                                                   113 
112   rangeTest = new G4RDRangeTest;                  114   rangeTest = new G4RDRangeTest;
113                                                   115 
114   // For Doppler broadening                       116   // For Doppler broadening
115   shellData.SetOccupancyData();                   117   shellData.SetOccupancyData();
116                                                   118 
117                                                   119 
118    if (verboseLevel > 0)                          120    if (verboseLevel > 0)
119      {                                            121      {
120        G4cout << GetProcessName() << " is crea    122        G4cout << GetProcessName() << " is created " << G4endl
121         << "Energy range: "                       123         << "Energy range: "
122         << lowEnergyLimit / keV << " keV - "      124         << lowEnergyLimit / keV << " keV - "
123         << highEnergyLimit / GeV << " GeV"        125         << highEnergyLimit / GeV << " GeV"
124         << G4endl;                                126         << G4endl;
125      }                                            127      }
126 }                                                 128 }
127                                                   129 
128 // destructor                                     130 // destructor
129                                                   131 
130 G4LowEnergyPolarizedCompton::~G4LowEnergyPolar    132 G4LowEnergyPolarizedCompton::~G4LowEnergyPolarizedCompton()
131 {                                                 133 {
132   delete meanFreePathTable;                       134   delete meanFreePathTable;
133   delete crossSectionHandler;                     135   delete crossSectionHandler;
134   delete scatterFunctionData;                     136   delete scatterFunctionData;
135   delete rangeTest;                               137   delete rangeTest;
136 }                                                 138 }
137                                                   139 
138                                                   140 
139 void G4LowEnergyPolarizedCompton::BuildPhysics    141 void G4LowEnergyPolarizedCompton::BuildPhysicsTable(const G4ParticleDefinition& )
140 {                                                 142 {
141                                                   143 
142   crossSectionHandler->Clear();                   144   crossSectionHandler->Clear();
143   G4String crossSectionFile = "comp/ce-cs-";      145   G4String crossSectionFile = "comp/ce-cs-";
144   crossSectionHandler->LoadData(crossSectionFi    146   crossSectionHandler->LoadData(crossSectionFile);
145   delete meanFreePathTable;                       147   delete meanFreePathTable;
146   meanFreePathTable = crossSectionHandler->Bui    148   meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
147                                                   149 
148   // For Doppler broadening                       150   // For Doppler broadening
149   G4String file = "/doppler/shell-doppler";       151   G4String file = "/doppler/shell-doppler";
150   shellData.LoadData(file);                       152   shellData.LoadData(file);
151                                                   153 
152 }                                                 154 }
153                                                   155 
154 G4VParticleChange* G4LowEnergyPolarizedCompton    156 G4VParticleChange* G4LowEnergyPolarizedCompton::PostStepDoIt(const G4Track& aTrack,
155                    const G4Step&  aStep)          157                    const G4Step&  aStep)
156 {                                                 158 {
157   // The scattered gamma energy is sampled acc    159   // The scattered gamma energy is sampled according to Klein - Nishina formula.
158   // The random number techniques of Butcher &    160   // The random number techniques of Butcher & Messel are used (Nuc Phys 20(1960),15).
159   // GEANT4 internal units                        161   // GEANT4 internal units
160   //                                              162   //
161   // Note : Effects due to binding of atomic e    163   // Note : Effects due to binding of atomic electrons are negliged.
162                                                   164 
163   aParticleChange.Initialize(aTrack);             165   aParticleChange.Initialize(aTrack);
164                                                   166 
165   // Dynamic particle quantities                  167   // Dynamic particle quantities
166   const G4DynamicParticle* incidentPhoton = aT    168   const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
167   G4double gammaEnergy0 = incidentPhoton->GetK    169   G4double gammaEnergy0 = incidentPhoton->GetKineticEnergy();
168   G4ThreeVector gammaPolarization0 = incidentP    170   G4ThreeVector gammaPolarization0 = incidentPhoton->GetPolarization();
169                                                   171 
170   //  gammaPolarization0 = gammaPolarization0.    172   //  gammaPolarization0 = gammaPolarization0.unit(); //
171                                                   173 
172   // Protection: a polarisation parallel to th    174   // Protection: a polarisation parallel to the
173   // direction causes problems;                   175   // direction causes problems;
174   // in that case find a random polarization      176   // in that case find a random polarization
175                                                   177 
176   G4ThreeVector gammaDirection0 = incidentPhot    178   G4ThreeVector gammaDirection0 = incidentPhoton->GetMomentumDirection();
177   // ---- MGP ---- Next two lines commented ou    179   // ---- MGP ---- Next two lines commented out to remove compilation warnings
178   // G4double scalarproduct = gammaPolarizatio    180   // G4double scalarproduct = gammaPolarization0.dot(gammaDirection0);
179   // G4double angle = gammaPolarization0.angle    181   // G4double angle = gammaPolarization0.angle(gammaDirection0);
180                                                   182 
181   // Make sure that the polarization vector is    183   // Make sure that the polarization vector is perpendicular to the
182   // gamma direction. If not                      184   // gamma direction. If not
183                                                   185 
184   if(!(gammaPolarization0.isOrthogonal(gammaDi    186   if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
185     { // only for testing now                     187     { // only for testing now
186       gammaPolarization0 = GetRandomPolarizati    188       gammaPolarization0 = GetRandomPolarization(gammaDirection0);
187     }                                             189     }
188   else                                            190   else
189     {                                             191     {
190       if ( gammaPolarization0.howOrthogonal(ga    192       if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
191   {                                               193   {
192     gammaPolarization0 = GetPerpendicularPolar    194     gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
193   }                                               195   }
194     }                                             196     }
195                                                   197 
196   // End of Protection                            198   // End of Protection
197                                                   199 
198   // Within energy limit?                         200   // Within energy limit?
199                                                   201 
200   if(gammaEnergy0 <= lowEnergyLimit)              202   if(gammaEnergy0 <= lowEnergyLimit)
201     {                                             203     {
202       aParticleChange.ProposeTrackStatus(fStop    204       aParticleChange.ProposeTrackStatus(fStopAndKill);
203       aParticleChange.ProposeEnergy(0.);          205       aParticleChange.ProposeEnergy(0.);
204       aParticleChange.ProposeLocalEnergyDeposi    206       aParticleChange.ProposeLocalEnergyDeposit(gammaEnergy0);
205       return G4VDiscreteProcess::PostStepDoIt(    207       return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
206     }                                             208     }
207                                                   209 
208   G4double E0_m = gammaEnergy0 / electron_mass    210   G4double E0_m = gammaEnergy0 / electron_mass_c2 ;
209                                                   211 
210   // Select randomly one element in the curren    212   // Select randomly one element in the current material
211                                                   213 
212   const G4MaterialCutsCouple* couple = aTrack.    214   const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
213   G4int Z = crossSectionHandler->SelectRandomA    215   G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0);
214                                                   216 
215   // Sample the energy and the polarization of    217   // Sample the energy and the polarization of the scattered photon
216                                                   218 
217   G4double epsilon, epsilonSq, onecost, sinThe    219   G4double epsilon, epsilonSq, onecost, sinThetaSqr, greject ;
218                                                   220 
219   G4double epsilon0 = 1./(1. + 2*E0_m);           221   G4double epsilon0 = 1./(1. + 2*E0_m);
220   G4double epsilon0Sq = epsilon0*epsilon0;        222   G4double epsilon0Sq = epsilon0*epsilon0;
221   G4double alpha1   = - std::log(epsilon0);       223   G4double alpha1   = - std::log(epsilon0);
222   G4double alpha2 = 0.5*(1.- epsilon0Sq);         224   G4double alpha2 = 0.5*(1.- epsilon0Sq);
223                                                   225 
224   G4double wlGamma = h_Planck*c_light/gammaEne    226   G4double wlGamma = h_Planck*c_light/gammaEnergy0;
225   G4double gammaEnergy1;                          227   G4double gammaEnergy1;
226   G4ThreeVector gammaDirection1;                  228   G4ThreeVector gammaDirection1;
227                                                   229 
228   do {                                            230   do {
229     if ( alpha1/(alpha1+alpha2) > G4UniformRan    231     if ( alpha1/(alpha1+alpha2) > G4UniformRand() )
230       {                                           232       {
231   epsilon   = std::exp(-alpha1*G4UniformRand()    233   epsilon   = std::exp(-alpha1*G4UniformRand());  
232   epsilonSq = epsilon*epsilon;                    234   epsilonSq = epsilon*epsilon; 
233       }                                           235       }
234     else                                          236     else 
235       {                                           237       {
236   epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4    238   epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4UniformRand();
237   epsilon   = std::sqrt(epsilonSq);               239   epsilon   = std::sqrt(epsilonSq);
238       }                                           240       }
239                                                   241 
240     onecost = (1.- epsilon)/(epsilon*E0_m);       242     onecost = (1.- epsilon)/(epsilon*E0_m);
241     sinThetaSqr   = onecost*(2.-onecost);         243     sinThetaSqr   = onecost*(2.-onecost);
242                                                   244 
243     // Protection                                 245     // Protection
244     if (sinThetaSqr > 1.)                         246     if (sinThetaSqr > 1.)
245       {                                           247       {
246   if (verboseLevel>0) G4cout                      248   if (verboseLevel>0) G4cout
247     << " -- Warning -- G4LowEnergyPolarizedCom    249     << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
248     << "sin(theta)**2 = "                         250     << "sin(theta)**2 = "
249     << sinThetaSqr                                251     << sinThetaSqr
250     << "; set to 1"                               252     << "; set to 1"
251     << G4endl;                                    253     << G4endl;
252   sinThetaSqr = 1.;                               254   sinThetaSqr = 1.;
253       }                                           255       }
254     if (sinThetaSqr < 0.)                         256     if (sinThetaSqr < 0.)
255       {                                           257       {
256   if (verboseLevel>0) G4cout                      258   if (verboseLevel>0) G4cout
257     << " -- Warning -- G4LowEnergyPolarizedCom    259     << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
258     << "sin(theta)**2 = "                         260     << "sin(theta)**2 = "
259     << sinThetaSqr                                261     << sinThetaSqr
260     << "; set to 0"                               262     << "; set to 0"
261     << G4endl;                                    263     << G4endl;
262   sinThetaSqr = 0.;                               264   sinThetaSqr = 0.;
263       }                                           265       }
264     // End protection                             266     // End protection
265                                                   267 
266     G4double x =  std::sqrt(onecost/2.) / (wlG    268     G4double x =  std::sqrt(onecost/2.) / (wlGamma/cm);;
267     G4double scatteringFunction = scatterFunct    269     G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
268     greject = (1. - epsilon*sinThetaSqr/(1.+ e    270     greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction;
269                                                   271 
270   } while(greject < G4UniformRand()*Z);           272   } while(greject < G4UniformRand()*Z);
271                                                   273 
272                                                   274 
273   // *****************************************    275   // ****************************************************
274   //    Phi determination                         276   //    Phi determination
275   // *****************************************    277   // ****************************************************
276                                                   278 
277   G4double phi = SetPhi(epsilon,sinThetaSqr);     279   G4double phi = SetPhi(epsilon,sinThetaSqr);
278                                                   280 
279   //                                              281   //
280   // scattered gamma angles. ( Z - axis along     282   // scattered gamma angles. ( Z - axis along the parent gamma)
281   //                                              283   //
282                                                   284 
283   G4double cosTheta = 1. - onecost;               285   G4double cosTheta = 1. - onecost;
284                                                   286 
285   // Protection                                   287   // Protection
286                                                   288 
287   if (cosTheta > 1.)                              289   if (cosTheta > 1.)
288     {                                             290     {
289       if (verboseLevel>0) G4cout                  291       if (verboseLevel>0) G4cout
290   << " -- Warning -- G4LowEnergyPolarizedCompt    292   << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
291   << "cosTheta = "                                293   << "cosTheta = "
292   << cosTheta                                     294   << cosTheta
293   << "; set to 1"                                 295   << "; set to 1"
294   << G4endl;                                      296   << G4endl;
295       cosTheta = 1.;                              297       cosTheta = 1.;
296     }                                             298     }
297   if (cosTheta < -1.)                             299   if (cosTheta < -1.)
298     {                                             300     {
299       if (verboseLevel>0) G4cout                  301       if (verboseLevel>0) G4cout 
300   << " -- Warning -- G4LowEnergyPolarizedCompt    302   << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
301   << "cosTheta = "                                303   << "cosTheta = " 
302   << cosTheta                                     304   << cosTheta
303   << "; set to -1"                                305   << "; set to -1"
304   << G4endl;                                      306   << G4endl;
305       cosTheta = -1.;                             307       cosTheta = -1.;
306     }                                             308     }
307   // End protection                               309   // End protection      
308                                                   310   
309                                                   311   
310   G4double sinTheta = std::sqrt (sinThetaSqr);    312   G4double sinTheta = std::sqrt (sinThetaSqr);
311                                                   313   
312   // Protection                                   314   // Protection
313   if (sinTheta > 1.)                              315   if (sinTheta > 1.)
314     {                                             316     {
315       if (verboseLevel>0) G4cout                  317       if (verboseLevel>0) G4cout 
316   << " -- Warning -- G4LowEnergyPolarizedCompt    318   << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
317   << "sinTheta = "                                319   << "sinTheta = " 
318   << sinTheta                                     320   << sinTheta
319   << "; set to 1"                                 321   << "; set to 1"
320   << G4endl;                                      322   << G4endl;
321       sinTheta = 1.;                              323       sinTheta = 1.;
322     }                                             324     }
323   if (sinTheta < -1.)                             325   if (sinTheta < -1.)
324     {                                             326     {
325       if (verboseLevel>0) G4cout                  327       if (verboseLevel>0) G4cout 
326   << " -- Warning -- G4LowEnergyPolarizedCompt    328   << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
327   << "sinTheta = "                                329   << "sinTheta = " 
328   << sinTheta                                     330   << sinTheta
329   << "; set to -1"                                331   << "; set to -1" 
330   << G4endl;                                      332   << G4endl;
331       sinTheta = -1.;                             333       sinTheta = -1.;
332     }                                             334     }
333   // End protection                               335   // End protection
334                                                   336   
335                                                   337       
336   G4double dirx = sinTheta*std::cos(phi);         338   G4double dirx = sinTheta*std::cos(phi);
337   G4double diry = sinTheta*std::sin(phi);         339   G4double diry = sinTheta*std::sin(phi);
338   G4double dirz = cosTheta ;                      340   G4double dirz = cosTheta ;
339                                                   341   
340                                                   342 
341   // oneCosT , eom                                343   // oneCosT , eom
342                                                   344 
343                                                   345 
344                                                   346 
345   // Doppler broadening -  Method based on:       347   // Doppler broadening -  Method based on:
346   // Y. Namito, S. Ban and H. Hirayama,           348   // Y. Namito, S. Ban and H. Hirayama, 
347   // "Implementation of the Doppler Broadening    349   // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code" 
348   // NIM A 349, pp. 489-494, 1994                 350   // NIM A 349, pp. 489-494, 1994
349                                                   351   
350   // Maximum number of sampling iterations        352   // Maximum number of sampling iterations
351                                                   353 
352   G4int maxDopplerIterations = 1000;              354   G4int maxDopplerIterations = 1000;
353   G4double bindingE = 0.;                         355   G4double bindingE = 0.;
354   G4double photonEoriginal = epsilon * gammaEn    356   G4double photonEoriginal = epsilon * gammaEnergy0;
355   G4double photonE = -1.;                         357   G4double photonE = -1.;
356   G4int iteration = 0;                            358   G4int iteration = 0;
357   G4double eMax = gammaEnergy0;                   359   G4double eMax = gammaEnergy0;
358                                                   360 
359   do                                              361   do
360     {                                             362     {
361       iteration++;                                363       iteration++;
362       // Select shell based on shell occupancy    364       // Select shell based on shell occupancy
363       G4int shell = shellData.SelectRandomShel    365       G4int shell = shellData.SelectRandomShell(Z);
364       bindingE = shellData.BindingEnergy(Z,she    366       bindingE = shellData.BindingEnergy(Z,shell);
365                                                   367       
366       eMax = gammaEnergy0 - bindingE;             368       eMax = gammaEnergy0 - bindingE;
367                                                   369      
368       // Randomly sample bound electron moment    370       // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
369       G4double pSample = profileData.RandomSel    371       G4double pSample = profileData.RandomSelectMomentum(Z,shell);
370       // Rescale from atomic units                372       // Rescale from atomic units
371       G4double pDoppler = pSample * fine_struc    373       G4double pDoppler = pSample * fine_structure_const;
372       G4double pDoppler2 = pDoppler * pDoppler    374       G4double pDoppler2 = pDoppler * pDoppler;
373       G4double var2 = 1. + onecost * E0_m;        375       G4double var2 = 1. + onecost * E0_m;
374       G4double var3 = var2*var2 - pDoppler2;      376       G4double var3 = var2*var2 - pDoppler2;
375       G4double var4 = var2 - pDoppler2 * cosTh    377       G4double var4 = var2 - pDoppler2 * cosTheta;
376       G4double var = var4*var4 - var3 + pDoppl    378       G4double var = var4*var4 - var3 + pDoppler2 * var3;
377       if (var > 0.)                               379       if (var > 0.)
378   {                                               380   {
379     G4double varSqrt = std::sqrt(var);            381     G4double varSqrt = std::sqrt(var);        
380     G4double scale = gammaEnergy0 / var3;         382     G4double scale = gammaEnergy0 / var3;  
381           // Random select either root            383           // Random select either root
382     if (G4UniformRand() < 0.5) photonE = (var4    384     if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;               
383     else photonE = (var4 + varSqrt) * scale;      385     else photonE = (var4 + varSqrt) * scale;
384   }                                               386   } 
385       else                                        387       else
386   {                                               388   {
387     photonE = -1.;                                389     photonE = -1.;
388   }                                               390   }
389    } while ( iteration <= maxDopplerIterations    391    } while ( iteration <= maxDopplerIterations && 
390        (photonE < 0. || photonE > eMax || phot    392        (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) );
391                                                   393  
392   // End of recalculation of photon energy wit    394   // End of recalculation of photon energy with Doppler broadening
393   // Revert to original if maximum number of i    395   // Revert to original if maximum number of iterations threshold has been reached
394   if (iteration >= maxDopplerIterations)          396   if (iteration >= maxDopplerIterations)
395     {                                             397     {
396       photonE = photonEoriginal;                  398       photonE = photonEoriginal;
397       bindingE = 0.;                              399       bindingE = 0.;
398     }                                             400     }
399                                                   401 
400   gammaEnergy1 = photonE;                         402   gammaEnergy1 = photonE;
401                                                   403  
402   // G4cout << "--> PHOTONENERGY1 = " << photo    404   // G4cout << "--> PHOTONENERGY1 = " << photonE/keV << G4endl;
403                                                   405 
404                                                   406 
405   /// Doppler Broadeing                           407   /// Doppler Broadeing 
406                                                   408 
407                                                   409 
408                                                   410 
409                                                   411 
410   //                                              412   //
411   // update G4VParticleChange for the scattere    413   // update G4VParticleChange for the scattered photon 
412   //                                              414   //
413                                                   415 
414   //  gammaEnergy1 = epsilon*gammaEnergy0;        416   //  gammaEnergy1 = epsilon*gammaEnergy0;
415                                                   417 
416                                                   418 
417   // New polarization                             419   // New polarization
418                                                   420 
419   G4ThreeVector gammaPolarization1 = SetNewPol    421   G4ThreeVector gammaPolarization1 = SetNewPolarization(epsilon,
420               sinThetaSqr,                        422               sinThetaSqr,
421               phi,                                423               phi,
422               cosTheta);                          424               cosTheta);
423                                                   425 
424   // Set new direction                            426   // Set new direction
425   G4ThreeVector tmpDirection1( dirx,diry,dirz     427   G4ThreeVector tmpDirection1( dirx,diry,dirz );
426   gammaDirection1 = tmpDirection1;                428   gammaDirection1 = tmpDirection1;
427                                                   429 
428   // Change reference frame.                      430   // Change reference frame.
429                                                   431 
430   SystemOfRefChange(gammaDirection0,gammaDirec    432   SystemOfRefChange(gammaDirection0,gammaDirection1,
431         gammaPolarization0,gammaPolarization1)    433         gammaPolarization0,gammaPolarization1);
432                                                   434 
433   if (gammaEnergy1 > 0.)                          435   if (gammaEnergy1 > 0.)
434     {                                             436     {
435       aParticleChange.ProposeEnergy( gammaEner    437       aParticleChange.ProposeEnergy( gammaEnergy1 ) ;
436       aParticleChange.ProposeMomentumDirection    438       aParticleChange.ProposeMomentumDirection( gammaDirection1 );
437       aParticleChange.ProposePolarization( gam    439       aParticleChange.ProposePolarization( gammaPolarization1 );
438     }                                             440     }
439   else                                            441   else
440     {                                             442     {
441       aParticleChange.ProposeEnergy(0.) ;         443       aParticleChange.ProposeEnergy(0.) ;
442       aParticleChange.ProposeTrackStatus(fStop    444       aParticleChange.ProposeTrackStatus(fStopAndKill);
443     }                                             445     }
444                                                   446 
445   //                                              447   //
446   // kinematic of the scattered electron          448   // kinematic of the scattered electron
447   //                                              449   //
448                                                   450 
449   G4double ElecKineEnergy = gammaEnergy0 - gam    451   G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;
450                                                   452 
451                                                   453 
452   // Generate the electron only if with large     454   // Generate the electron only if with large enough range w.r.t. cuts and safety
453                                                   455 
454   G4double safety = aStep.GetPostStepPoint()->    456   G4double safety = aStep.GetPostStepPoint()->GetSafety();
455                                                   457 
456                                                   458 
457   if (rangeTest->Escape(G4Electron::Electron()    459   if (rangeTest->Escape(G4Electron::Electron(),couple,ElecKineEnergy,safety))
458     {                                             460     {
459       G4double ElecMomentum = std::sqrt(ElecKi    461       G4double ElecMomentum = std::sqrt(ElecKineEnergy*(ElecKineEnergy+2.*electron_mass_c2));
460       G4ThreeVector ElecDirection((gammaEnergy    462       G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 -
461            gammaEnergy1 * gammaDirection1) * (    463            gammaEnergy1 * gammaDirection1) * (1./ElecMomentum));
462       G4DynamicParticle* electron = new G4Dyna    464       G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),ElecDirection.unit(),ElecKineEnergy) ;
463       aParticleChange.SetNumberOfSecondaries(1    465       aParticleChange.SetNumberOfSecondaries(1);
464       aParticleChange.AddSecondary(electron);     466       aParticleChange.AddSecondary(electron);
465       //      aParticleChange.ProposeLocalEner    467       //      aParticleChange.ProposeLocalEnergyDeposit(0.); 
466       aParticleChange.ProposeLocalEnergyDeposi    468       aParticleChange.ProposeLocalEnergyDeposit(bindingE);
467     }                                             469     }
468   else                                            470   else
469     {                                             471     {
470       aParticleChange.SetNumberOfSecondaries(0    472       aParticleChange.SetNumberOfSecondaries(0);
471       aParticleChange.ProposeLocalEnergyDeposi    473       aParticleChange.ProposeLocalEnergyDeposit(ElecKineEnergy+bindingE);
472     }                                             474     }
473                                                   475   
474   return G4VDiscreteProcess::PostStepDoIt( aTr    476   return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep);
475                                                   477   
476 }                                                 478 }
477                                                   479 
478                                                   480 
479 G4double G4LowEnergyPolarizedCompton::SetPhi(G    481 G4double G4LowEnergyPolarizedCompton::SetPhi(G4double energyRate,
480                G4double sinSqrTh)                 482                G4double sinSqrTh)
481 {                                                 483 {
482   G4double rand1;                                 484   G4double rand1;
483   G4double rand2;                                 485   G4double rand2;
484   G4double phiProbability;                        486   G4double phiProbability;
485   G4double phi;                                   487   G4double phi;
486   G4double a, b;                                  488   G4double a, b;
487                                                   489 
488   do                                              490   do
489     {                                             491     {
490       rand1 = G4UniformRand();                    492       rand1 = G4UniformRand();
491       rand2 = G4UniformRand();                    493       rand2 = G4UniformRand();
492       phiProbability=0.;                          494       phiProbability=0.;
493       phi = twopi*rand1;                          495       phi = twopi*rand1;
494                                                   496       
495       a = 2*sinSqrTh;                             497       a = 2*sinSqrTh;
496       b = energyRate + 1/energyRate;              498       b = energyRate + 1/energyRate;
497                                                   499       
498       phiProbability = 1 - (a/b)*(std::cos(phi    500       phiProbability = 1 - (a/b)*(std::cos(phi)*std::cos(phi));
499                                                   501 
500                                                   502       
501                                                   503  
502     }                                             504     }
503   while ( rand2 > phiProbability );               505   while ( rand2 > phiProbability );
504   return phi;                                     506   return phi;
505 }                                                 507 }
506                                                   508 
507                                                   509 
508 G4ThreeVector G4LowEnergyPolarizedCompton::Set    510 G4ThreeVector G4LowEnergyPolarizedCompton::SetPerpendicularVector(G4ThreeVector& a)
509 {                                                 511 {
510   G4double dx = a.x();                            512   G4double dx = a.x();
511   G4double dy = a.y();                            513   G4double dy = a.y();
512   G4double dz = a.z();                            514   G4double dz = a.z();
513   G4double x = dx < 0.0 ? -dx : dx;               515   G4double x = dx < 0.0 ? -dx : dx;
514   G4double y = dy < 0.0 ? -dy : dy;               516   G4double y = dy < 0.0 ? -dy : dy;
515   G4double z = dz < 0.0 ? -dz : dz;               517   G4double z = dz < 0.0 ? -dz : dz;
516   if (x < y) {                                    518   if (x < y) {
517     return x < z ? G4ThreeVector(-dy,dx,0) : G    519     return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
518   }else{                                          520   }else{
519     return y < z ? G4ThreeVector(dz,0,-dx) : G    521     return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
520   }                                               522   }
521 }                                                 523 }
522                                                   524 
523 G4ThreeVector G4LowEnergyPolarizedCompton::Get    525 G4ThreeVector G4LowEnergyPolarizedCompton::GetRandomPolarization(G4ThreeVector& direction0)
524 {                                                 526 {
525   G4ThreeVector d0 = direction0.unit();           527   G4ThreeVector d0 = direction0.unit();
526   G4ThreeVector a1 = SetPerpendicularVector(d0    528   G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
527   G4ThreeVector a0 = a1.unit(); // unit vector    529   G4ThreeVector a0 = a1.unit(); // unit vector
528                                                   530 
529   G4double rand1 = G4UniformRand();               531   G4double rand1 = G4UniformRand();
530                                                   532   
531   G4double angle = twopi*rand1; // random pola    533   G4double angle = twopi*rand1; // random polar angle
532   G4ThreeVector b0 = d0.cross(a0); // cross pr    534   G4ThreeVector b0 = d0.cross(a0); // cross product
533                                                   535   
534   G4ThreeVector c;                                536   G4ThreeVector c;
535                                                   537   
536   c.setX(std::cos(angle)*(a0.x())+std::sin(ang    538   c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
537   c.setY(std::cos(angle)*(a0.y())+std::sin(ang    539   c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
538   c.setZ(std::cos(angle)*(a0.z())+std::sin(ang    540   c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
539                                                   541   
540   G4ThreeVector c0 = c.unit();                    542   G4ThreeVector c0 = c.unit();
541                                                   543 
542   return c0;                                      544   return c0;
543                                                   545   
544 }                                                 546 }
545                                                   547 
546                                                   548 
547 G4ThreeVector G4LowEnergyPolarizedCompton::Get    549 G4ThreeVector G4LowEnergyPolarizedCompton::GetPerpendicularPolarization
548 (const G4ThreeVector& gammaDirection, const G4    550 (const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
549 {                                                 551 {
550                                                   552 
551   //                                              553   // 
552   // The polarization of a photon is always pe    554   // The polarization of a photon is always perpendicular to its momentum direction.
553   // Therefore this function removes those vec    555   // Therefore this function removes those vector component of gammaPolarization, which
554   // points in direction of gammaDirection        556   // points in direction of gammaDirection
555   //                                              557   //
556   // Mathematically we search the projection o    558   // Mathematically we search the projection of the vector a on the plane E, where n is the
557   // plains normal vector.                        559   // plains normal vector.
558   // The basic equation can be found in each g    560   // The basic equation can be found in each geometry book (e.g. Bronstein):
559   // p = a - (a o n)/(n o n)*n                    561   // p = a - (a o n)/(n o n)*n
560                                                   562   
561   return gammaPolarization - gammaPolarization    563   return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
562 }                                                 564 }
563                                                   565 
564                                                   566 
565 G4ThreeVector G4LowEnergyPolarizedCompton::Set    567 G4ThreeVector G4LowEnergyPolarizedCompton::SetNewPolarization(G4double epsilon,
566                     G4double sinSqrTh,            568                     G4double sinSqrTh, 
567                     G4double phi,                 569                     G4double phi,
568                     G4double costheta)            570                     G4double costheta) 
569 {                                                 571 {
570   G4double rand1;                                 572   G4double rand1;
571   G4double rand2;                                 573   G4double rand2;
572   G4double cosPhi = std::cos(phi);                574   G4double cosPhi = std::cos(phi);
573   G4double sinPhi = std::sin(phi);                575   G4double sinPhi = std::sin(phi);
574   G4double sinTheta = std::sqrt(sinSqrTh);        576   G4double sinTheta = std::sqrt(sinSqrTh);
575   G4double cosSqrPhi = cosPhi*cosPhi;             577   G4double cosSqrPhi = cosPhi*cosPhi;
576   //  G4double cossqrth = 1.-sinSqrTh;            578   //  G4double cossqrth = 1.-sinSqrTh;
577   //  G4double sinsqrphi = sinPhi*sinPhi;         579   //  G4double sinsqrphi = sinPhi*sinPhi;
578   G4double normalisation = std::sqrt(1. - cosS    580   G4double normalisation = std::sqrt(1. - cosSqrPhi*sinSqrTh);
579                                                   581  
580                                                   582 
581   // Determination of Theta                       583   // Determination of Theta 
582                                                   584   
583   // ---- MGP ---- Commented out the following    585   // ---- MGP ---- Commented out the following 3 lines to avoid compilation 
584   // warnings (unused variables)                  586   // warnings (unused variables)
585   // G4double thetaProbability;                   587   // G4double thetaProbability;
586   G4double theta;                                 588   G4double theta;
587   // G4double a, b;                               589   // G4double a, b;
588   // G4double cosTheta;                           590   // G4double cosTheta;
589                                                   591 
590   /*                                              592   /*
591                                                   593 
592   depaola method                                  594   depaola method
593                                                   595   
594   do                                              596   do
595   {                                               597   {
596       rand1 = G4UniformRand();                    598       rand1 = G4UniformRand();
597       rand2 = G4UniformRand();                    599       rand2 = G4UniformRand();
598       thetaProbability=0.;                        600       thetaProbability=0.;
599       theta = twopi*rand1;                        601       theta = twopi*rand1;
600       a = 4*normalisation*normalisation;          602       a = 4*normalisation*normalisation;
601       b = (epsilon + 1/epsilon) - 2;              603       b = (epsilon + 1/epsilon) - 2;
602       thetaProbability = (b + a*std::cos(theta    604       thetaProbability = (b + a*std::cos(theta)*std::cos(theta))/(a+b);
603       cosTheta = std::cos(theta);                 605       cosTheta = std::cos(theta);
604     }                                             606     }
605   while ( rand2 > thetaProbability );             607   while ( rand2 > thetaProbability );
606                                                   608   
607   G4double cosBeta = cosTheta;                    609   G4double cosBeta = cosTheta;
608                                                   610 
609   */                                              611   */
610                                                   612 
611                                                   613 
612   // Dan Xu method (IEEE TNS, 52, 1160 (2005))    614   // Dan Xu method (IEEE TNS, 52, 1160 (2005))
613                                                   615 
614   rand1 = G4UniformRand();                        616   rand1 = G4UniformRand();
615   rand2 = G4UniformRand();                        617   rand2 = G4UniformRand();
616                                                   618 
617   if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsi    619   if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsilon+1.0/epsilon)-4.0*sinSqrTh*cosSqrPhi))
618     {                                             620     {
619       if (rand2<0.5)                              621       if (rand2<0.5)
620   theta = pi/2.0;                                 622   theta = pi/2.0;
621       else                                        623       else
622   theta = 3.0*pi/2.0;                             624   theta = 3.0*pi/2.0;
623     }                                             625     }
624   else                                            626   else
625     {                                             627     {
626       if (rand2<0.5)                              628       if (rand2<0.5)
627   theta = 0;                                      629   theta = 0;
628       else                                        630       else
629   theta = pi;                                     631   theta = pi;
630     }                                             632     }
631   G4double cosBeta = std::cos(theta);             633   G4double cosBeta = std::cos(theta);
632   G4double sinBeta = std::sqrt(1-cosBeta*cosBe    634   G4double sinBeta = std::sqrt(1-cosBeta*cosBeta);
633                                                   635   
634   G4ThreeVector gammaPolarization1;               636   G4ThreeVector gammaPolarization1;
635                                                   637 
636   G4double xParallel = normalisation*cosBeta;     638   G4double xParallel = normalisation*cosBeta;
637   G4double yParallel = -(sinSqrTh*cosPhi*sinPh    639   G4double yParallel = -(sinSqrTh*cosPhi*sinPhi)*cosBeta/normalisation;
638   G4double zParallel = -(costheta*sinTheta*cos    640   G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation;
639   G4double xPerpendicular = 0.;                   641   G4double xPerpendicular = 0.;
640   G4double yPerpendicular = (costheta)*sinBeta    642   G4double yPerpendicular = (costheta)*sinBeta/normalisation;
641   G4double zPerpendicular = -(sinTheta*sinPhi)    643   G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation;
642                                                   644 
643   G4double xTotal = (xParallel + xPerpendicula    645   G4double xTotal = (xParallel + xPerpendicular);
644   G4double yTotal = (yParallel + yPerpendicula    646   G4double yTotal = (yParallel + yPerpendicular);
645   G4double zTotal = (zParallel + zPerpendicula    647   G4double zTotal = (zParallel + zPerpendicular);
646                                                   648   
647   gammaPolarization1.setX(xTotal);                649   gammaPolarization1.setX(xTotal);
648   gammaPolarization1.setY(yTotal);                650   gammaPolarization1.setY(yTotal);
649   gammaPolarization1.setZ(zTotal);                651   gammaPolarization1.setZ(zTotal);
650                                                   652   
651   return gammaPolarization1;                      653   return gammaPolarization1;
652                                                   654 
653 }                                                 655 }
654                                                   656 
655                                                   657 
656 void G4LowEnergyPolarizedCompton::SystemOfRefC    658 void G4LowEnergyPolarizedCompton::SystemOfRefChange(G4ThreeVector& direction0,
657                 G4ThreeVector& direction1,        659                 G4ThreeVector& direction1,
658                 G4ThreeVector& polarization0,     660                 G4ThreeVector& polarization0,
659                 G4ThreeVector& polarization1)     661                 G4ThreeVector& polarization1)
660 {                                                 662 {
661   // direction0 is the original photon directi    663   // direction0 is the original photon direction ---> z
662   // polarization0 is the original photon pola    664   // polarization0 is the original photon polarization ---> x
663   // need to specify y axis in the real refere    665   // need to specify y axis in the real reference frame ---> y 
664   G4ThreeVector Axis_Z0 = direction0.unit();      666   G4ThreeVector Axis_Z0 = direction0.unit();
665   G4ThreeVector Axis_X0 = polarization0.unit()    667   G4ThreeVector Axis_X0 = polarization0.unit();
666   G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_    668   G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
667                                                   669 
668   G4double direction_x = direction1.getX();       670   G4double direction_x = direction1.getX();
669   G4double direction_y = direction1.getY();       671   G4double direction_y = direction1.getY();
670   G4double direction_z = direction1.getZ();       672   G4double direction_z = direction1.getZ();
671                                                   673   
672   direction1 = (direction_x*Axis_X0 + directio    674   direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
673   G4double polarization_x = polarization1.getX    675   G4double polarization_x = polarization1.getX();
674   G4double polarization_y = polarization1.getY    676   G4double polarization_y = polarization1.getY();
675   G4double polarization_z = polarization1.getZ    677   G4double polarization_z = polarization1.getZ();
676                                                   678 
677   polarization1 = (polarization_x*Axis_X0 + po    679   polarization1 = (polarization_x*Axis_X0 + polarization_y*Axis_Y0 + polarization_z*Axis_Z0).unit();
678                                                   680 
679 }                                                 681 }
680                                                   682 
681                                                   683 
682 G4bool G4LowEnergyPolarizedCompton::IsApplicab    684 G4bool G4LowEnergyPolarizedCompton::IsApplicable(const G4ParticleDefinition& particle)
683 {                                                 685 {
684   return ( &particle == G4Gamma::Gamma() );       686   return ( &particle == G4Gamma::Gamma() ); 
685 }                                                 687 }
686                                                   688 
687                                                   689 
688 G4double G4LowEnergyPolarizedCompton::GetMeanF    690 G4double G4LowEnergyPolarizedCompton::GetMeanFreePath(const G4Track& track,
689                   G4double,                       691                   G4double,
690                   G4ForceCondition*)              692                   G4ForceCondition*)
691 {                                                 693 {
692   const G4DynamicParticle* photon = track.GetD    694   const G4DynamicParticle* photon = track.GetDynamicParticle();
693   G4double energy = photon->GetKineticEnergy()    695   G4double energy = photon->GetKineticEnergy();
694   const G4MaterialCutsCouple* couple = track.G    696   const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
695   size_t materialIndex = couple->GetIndex();      697   size_t materialIndex = couple->GetIndex();
696   G4double meanFreePath;                          698   G4double meanFreePath;
697   if (energy > highEnergyLimit) meanFreePath =    699   if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
698   else if (energy < lowEnergyLimit) meanFreePa    700   else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
699   else meanFreePath = meanFreePathTable->FindV    701   else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
700   return meanFreePath;                            702   return meanFreePath;
701 }                                                 703 }
702                                                   704 
703                                                   705 
704                                                   706 
705                                                   707 
706                                                   708 
707                                                   709 
708                                                   710 
709                                                   711 
710                                                   712 
711                                                   713 
712                                                   714 
713                                                   715 
714                                                   716 
715                                                   717 
716                                                   718