Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/polarisation/src/G4StokesVector.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/polarisation/src/G4StokesVector.cc (Version 11.3.0) and /processes/electromagnetic/polarisation/src/G4StokesVector.cc (Version 11.0.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 // Geant4 Class file                               26 // Geant4 Class file
 27 //                                                 27 //
 28 // File name:     G4StokesVector                   28 // File name:     G4StokesVector
 29 //                                                 29 //
 30 // Author:        Andreas Schaelicke               30 // Author:        Andreas Schaelicke
 31 //                                                 31 //
 32 // Class Description:                              32 // Class Description:
 33 //   Provides Stokes vector representation emp     33 //   Provides Stokes vector representation employed in polarized processes.
 34                                                    34 
 35 #include "G4StokesVector.hh"                       35 #include "G4StokesVector.hh"
 36                                                    36 
 37 #include "G4PolarizationHelper.hh"                 37 #include "G4PolarizationHelper.hh"
 38 #include "Randomize.hh"                            38 #include "Randomize.hh"
 39                                                    39 
 40 const G4StokesVector G4StokesVector::ZERO =        40 const G4StokesVector G4StokesVector::ZERO =
 41   G4StokesVector(G4ThreeVector(0., 0., 0.));       41   G4StokesVector(G4ThreeVector(0., 0., 0.));
 42 const G4StokesVector G4StokesVector::P1 =          42 const G4StokesVector G4StokesVector::P1 =
 43   G4StokesVector(G4ThreeVector(1., 0., 0.));       43   G4StokesVector(G4ThreeVector(1., 0., 0.));
 44 const G4StokesVector G4StokesVector::P2 =          44 const G4StokesVector G4StokesVector::P2 =
 45   G4StokesVector(G4ThreeVector(0., 1., 0.));       45   G4StokesVector(G4ThreeVector(0., 1., 0.));
 46 const G4StokesVector G4StokesVector::P3 =          46 const G4StokesVector G4StokesVector::P3 =
 47   G4StokesVector(G4ThreeVector(0., 0., 1.));       47   G4StokesVector(G4ThreeVector(0., 0., 1.));
 48 const G4StokesVector G4StokesVector::M1 =          48 const G4StokesVector G4StokesVector::M1 =
 49   G4StokesVector(G4ThreeVector(-1., 0., 0.));      49   G4StokesVector(G4ThreeVector(-1., 0., 0.));
 50 const G4StokesVector G4StokesVector::M2 =          50 const G4StokesVector G4StokesVector::M2 =
 51   G4StokesVector(G4ThreeVector(0., -1., 0.));      51   G4StokesVector(G4ThreeVector(0., -1., 0.));
 52 const G4StokesVector G4StokesVector::M3 =          52 const G4StokesVector G4StokesVector::M3 =
 53   G4StokesVector(G4ThreeVector(0., 0., -1.));      53   G4StokesVector(G4ThreeVector(0., 0., -1.));
 54                                                    54 
 55 G4StokesVector::G4StokesVector()                   55 G4StokesVector::G4StokesVector()
 56   : G4ThreeVector()                                56   : G4ThreeVector()
 57   , fIsPhoton(false)                               57   , fIsPhoton(false)
 58 {}                                                 58 {}
 59                                                    59 
 60 G4StokesVector::G4StokesVector(const G4ThreeVe     60 G4StokesVector::G4StokesVector(const G4ThreeVector& v)
 61   : G4ThreeVector(v)                               61   : G4ThreeVector(v)
 62   , fIsPhoton(false)                               62   , fIsPhoton(false)
 63 {}                                                 63 {}
 64                                                    64 
 65 G4bool G4StokesVector::IsZero() const { return     65 G4bool G4StokesVector::IsZero() const { return *this == ZERO; }
 66                                                    66 
 67 void G4StokesVector::RotateAz(G4ThreeVector nI     67 void G4StokesVector::RotateAz(G4ThreeVector nInteractionFrame,
 68                               G4ThreeVector pa     68                               G4ThreeVector particleDirection)
 69 {                                                  69 {
 70   G4ThreeVector yParticleFrame =                   70   G4ThreeVector yParticleFrame =
 71     G4PolarizationHelper::GetParticleFrameY(pa     71     G4PolarizationHelper::GetParticleFrameY(particleDirection);
 72                                                    72 
 73   G4double cosphi = yParticleFrame * nInteract     73   G4double cosphi = yParticleFrame * nInteractionFrame;
 74   if(cosphi > (1. + 1.e-8) || cosphi < (-1. -      74   if(cosphi > (1. + 1.e-8) || cosphi < (-1. - 1.e-8))
 75   {                                                75   {
 76     G4ExceptionDescription ed;                     76     G4ExceptionDescription ed;
 77     ed << " warning G4StokesVector::RotateAz       77     ed << " warning G4StokesVector::RotateAz  cosphi>1 or cosphi<-1\n"
 78        << " cosphi=" << cosphi << "\n"             78        << " cosphi=" << cosphi << "\n"
 79        << " zAxis=" << particleDirection << "      79        << " zAxis=" << particleDirection << " (" << particleDirection.mag()
 80        << ")\n"                                    80        << ")\n"
 81        << " yAxis=" << yParticleFrame << " ("      81        << " yAxis=" << yParticleFrame << " (" << yParticleFrame.mag() << ")\n"
 82        << " nAxis=" << nInteractionFrame << "      82        << " nAxis=" << nInteractionFrame << " (" << nInteractionFrame.mag()
 83        << ")\n";                                   83        << ")\n";
 84     G4Exception("G4StokesVector::RotateAz", "p     84     G4Exception("G4StokesVector::RotateAz", "pol030", JustWarning, ed);
 85   }                                                85   }
 86   if(cosphi > 1.)                                  86   if(cosphi > 1.)
 87     cosphi = 1.;                                   87     cosphi = 1.;
 88   else if(cosphi < -1.)                            88   else if(cosphi < -1.)
 89     cosphi = -1.;                                  89     cosphi = -1.;
 90                                                    90 
 91   G4double hel =                                   91   G4double hel =
 92     (yParticleFrame.cross(nInteractionFrame) *     92     (yParticleFrame.cross(nInteractionFrame) * particleDirection) > 0. ? 1.
 93                                                    93                                                                        : -1.;
 94                                                    94 
 95   G4double sinphi = hel * std::sqrt(1. - cosph     95   G4double sinphi = hel * std::sqrt(1. - cosphi * cosphi);
 96                                                    96 
 97   RotateAz(cosphi, sinphi);                        97   RotateAz(cosphi, sinphi);
 98 }                                                  98 }
 99                                                    99 
100 void G4StokesVector::InvRotateAz(G4ThreeVector    100 void G4StokesVector::InvRotateAz(G4ThreeVector nInteractionFrame,
101                                  G4ThreeVector    101                                  G4ThreeVector particleDirection)
102 {                                                 102 {
103   // note if incoming particle is on z-axis,      103   // note if incoming particle is on z-axis,
104   // we might encounter some nummerical proble    104   // we might encounter some nummerical problems, since
105   // nInteratonFrame and yParticleFrame are ac    105   // nInteratonFrame and yParticleFrame are actually (almost) the same momentum
106   // and the normalization is only good to 10^    106   // and the normalization is only good to 10^-12 !
107                                                   107 
108   G4ThreeVector yParticleFrame =                  108   G4ThreeVector yParticleFrame =
109     G4PolarizationHelper::GetParticleFrameY(pa    109     G4PolarizationHelper::GetParticleFrameY(particleDirection);
110   G4double cosphi = yParticleFrame * nInteract    110   G4double cosphi = yParticleFrame * nInteractionFrame;
111                                                   111 
112   if(cosphi > 1. + 1.e-8 || cosphi < -1. - 1.e    112   if(cosphi > 1. + 1.e-8 || cosphi < -1. - 1.e-8)
113   {                                               113   {
114     G4ExceptionDescription ed;                    114     G4ExceptionDescription ed;
115     ed << " warning G4StokesVector::RotateAz      115     ed << " warning G4StokesVector::RotateAz  cosphi>1 or cosphi<-1\n";
116     G4Exception("G4StokesVector::InvRotateAz",    116     G4Exception("G4StokesVector::InvRotateAz", "pol030", JustWarning, ed);
117   }                                               117   }
118   if(cosphi > 1.)                                 118   if(cosphi > 1.)
119     cosphi = 1.;                                  119     cosphi = 1.;
120   else if(cosphi < -1.)                           120   else if(cosphi < -1.)
121     cosphi = -1.;                                 121     cosphi = -1.;
122                                                   122 
123   // check sign once more!                        123   // check sign once more!
124   G4double hel =                                  124   G4double hel =
125     (yParticleFrame.cross(nInteractionFrame) *    125     (yParticleFrame.cross(nInteractionFrame) * particleDirection) > 0. ? 1.
126                                                   126                                                                        : -1.;
127   G4double sinphi = hel * std::sqrt(std::fabs(    127   G4double sinphi = hel * std::sqrt(std::fabs(1. - cosphi * cosphi));
128   RotateAz(cosphi, -sinphi);                      128   RotateAz(cosphi, -sinphi);
129 }                                                 129 }
130                                                   130 
131 void G4StokesVector::RotateAz(G4double cosphi,    131 void G4StokesVector::RotateAz(G4double cosphi, G4double sinphi)
132 {                                                 132 {
133   if(!fIsPhoton)                                  133   if(!fIsPhoton)
134   {                                               134   {
135     G4double xsi1 = cosphi * p1() + sinphi * p    135     G4double xsi1 = cosphi * p1() + sinphi * p2();
136     G4double xsi2 = -sinphi * p1() + cosphi *     136     G4double xsi2 = -sinphi * p1() + cosphi * p2();
137     setX(xsi1);                                   137     setX(xsi1);
138     setY(xsi2);                                   138     setY(xsi2);
139     return;                                       139     return;
140   }                                               140   }
141                                                   141 
142   G4double sin2phi = 2. * cosphi * sinphi;        142   G4double sin2phi = 2. * cosphi * sinphi;
143   G4double cos2phi = cosphi * cosphi - sinphi     143   G4double cos2phi = cosphi * cosphi - sinphi * sinphi;
144                                                   144 
145   G4double xsi1 = cos2phi * p1() + sin2phi * p    145   G4double xsi1 = cos2phi * p1() + sin2phi * p2();
146   G4double xsi2 = -sin2phi * p1() + cos2phi *     146   G4double xsi2 = -sin2phi * p1() + cos2phi * p2();
147   setX(xsi1);                                     147   setX(xsi1);
148   setY(xsi2);                                     148   setY(xsi2);
149 }                                                 149 }
150                                                   150 
151 G4double G4StokesVector::GetBeta()                151 G4double G4StokesVector::GetBeta()
152 {                                                 152 {
153   G4double bet = getPhi();                        153   G4double bet = getPhi();
154   if(fIsPhoton)                                   154   if(fIsPhoton)
155   {                                               155   {
156     bet *= 0.5;                                   156     bet *= 0.5;
157   }                                               157   }
158   return bet;                                     158   return bet;
159 }                                                 159 }
160                                                   160 
161 void G4StokesVector::DiceUniform()                161 void G4StokesVector::DiceUniform()
162 {                                                 162 {
163   G4double costheta = 2. * G4UniformRand() - 1    163   G4double costheta = 2. * G4UniformRand() - 1.;
164   G4double sintheta = std::sqrt(1. - costheta     164   G4double sintheta = std::sqrt(1. - costheta * costheta);
165   G4double aphi     = 2. * CLHEP::pi * G4Unifo    165   G4double aphi     = 2. * CLHEP::pi * G4UniformRand();
166   setX(std::sin(aphi) * sintheta);                166   setX(std::sin(aphi) * sintheta);
167   setY(std::cos(aphi) * sintheta);                167   setY(std::cos(aphi) * sintheta);
168   setZ(costheta);                                 168   setZ(costheta);
169 }                                                 169 }
170                                                   170 
171 void G4StokesVector::DiceP1()                     171 void G4StokesVector::DiceP1()
172 {                                                 172 {
173   if(G4UniformRand() > 0.5)                       173   if(G4UniformRand() > 0.5)
174     setX(1.);                                     174     setX(1.);
175   else                                            175   else
176     setX(-1.);                                    176     setX(-1.);
177   setY(0.);                                       177   setY(0.);
178   setZ(0.);                                       178   setZ(0.);
179 }                                                 179 }
180                                                   180 
181 void G4StokesVector::DiceP2()                     181 void G4StokesVector::DiceP2()
182 {                                                 182 {
183   setX(0.);                                       183   setX(0.);
184   if(G4UniformRand() > 0.5)                       184   if(G4UniformRand() > 0.5)
185     setY(1.);                                     185     setY(1.);
186   else                                            186   else
187     setY(-1.);                                    187     setY(-1.);
188   setZ(0.);                                       188   setZ(0.);
189 }                                                 189 }
190                                                   190 
191 void G4StokesVector::DiceP3()                     191 void G4StokesVector::DiceP3()
192 {                                                 192 {
193   setX(0.);                                       193   setX(0.);
194   setY(0.);                                       194   setY(0.);
195   if(G4UniformRand() > 0.5)                       195   if(G4UniformRand() > 0.5)
196     setZ(1.);                                     196     setZ(1.);
197   else                                            197   else
198     setZ(-1.);                                    198     setZ(-1.);
199 }                                                 199 }
200                                                   200 
201 void G4StokesVector::FlipP3() { setZ(-z()); }     201 void G4StokesVector::FlipP3() { setZ(-z()); }
202                                                   202 
203 G4ThreeVector G4StokesVector::PolError(const G    203 G4ThreeVector G4StokesVector::PolError(const G4StokesVector& sum2, long n)
204 {                                                 204 {
205   // delta x = sqrt[ ( <x^2> - <x>^2 )/(n-1) ]    205   // delta x = sqrt[ ( <x^2> - <x>^2 )/(n-1) ]
206   G4ThreeVector mean   = (1. / n) * G4ThreeVec    206   G4ThreeVector mean   = (1. / n) * G4ThreeVector(*this);
207   G4ThreeVector polsqr = G4StokesVector(mean).    207   G4ThreeVector polsqr = G4StokesVector(mean).PolSqr();
208   G4ThreeVector result =                          208   G4ThreeVector result =
209     G4StokesVector((1. / (n - 1.) * ((1. / n)     209     G4StokesVector((1. / (n - 1.) * ((1. / n) * sum2 - polsqr))).PolSqrt();
210   return result;                                  210   return result;
211 }                                                 211 }
212                                                   212 
213 G4ThreeVector G4StokesVector::PolDiv(const G4S    213 G4ThreeVector G4StokesVector::PolDiv(const G4StokesVector& b)
214 {                                                 214 {
215   return G4ThreeVector(b.x() != 0. ? x() / b.x    215   return G4ThreeVector(b.x() != 0. ? x() / b.x() : 11111.,
216                        b.y() != 0. ? y() / b.y    216                        b.y() != 0. ? y() / b.y() : 11111.,
217                        b.z() != 0. ? z() / b.z    217                        b.z() != 0. ? z() / b.z() : 11111.);
218 }                                                 218 }
219                                                   219