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