00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044 #include "G4StokesVector.hh"
00045 #include "G4PolarizationHelper.hh"
00046 #include "G4PhysicalConstants.hh"
00047 #include "Randomize.hh"
00048
00049 const G4StokesVector G4StokesVector::ZERO=G4ThreeVector(0.,0.,0.);
00050 const G4StokesVector G4StokesVector::P1=G4ThreeVector(1.,0.,0.);
00051 const G4StokesVector G4StokesVector::P2=G4ThreeVector(0.,1.,0.);
00052 const G4StokesVector G4StokesVector::P3=G4ThreeVector(0.,0.,1.);
00053 const G4StokesVector G4StokesVector::M1=G4ThreeVector(-1.,0.,0.);
00054 const G4StokesVector G4StokesVector::M2=G4ThreeVector(0.,-1.,0.);
00055 const G4StokesVector G4StokesVector::M3=G4ThreeVector(0.,0.,-1.);
00056
00057 G4StokesVector::G4StokesVector()
00058 : G4ThreeVector(),isPhoton(false)
00059 {
00060 }
00061
00062 G4StokesVector::G4StokesVector(const G4ThreeVector & v)
00063 : G4ThreeVector(v),isPhoton(false)
00064 {
00065 }
00066
00067 G4StokesVector::~G4StokesVector()
00068 {
00069 }
00070
00071 void G4StokesVector::RotateAz(G4ThreeVector nInteractionFrame,
00072 G4ThreeVector particleDirection)
00073 {
00074 G4ThreeVector yParticleFrame =
00075 G4PolarizationHelper::GetParticleFrameY(particleDirection);
00076
00077
00078 G4double cosphi=yParticleFrame*nInteractionFrame;
00079 if (cosphi>(1.+1.e-8) || cosphi<(-1.-1.e-8)) {
00080 G4cout<<" warning G4StokesVector::RotateAz cosphi>1 or cosphi<-1\n"
00081 <<" cosphi="<<cosphi<<"\n"
00082 <<" zAxis="<<particleDirection<<" ("<<particleDirection.mag()<<")\n"
00083 <<" yAxis="<<yParticleFrame<<" ("<<yParticleFrame.mag()<<")\n"
00084 <<" nAxis="<<nInteractionFrame<<" ("
00085 <<nInteractionFrame.mag()<<")"<<G4endl;
00086 }
00087 if (cosphi>1.) cosphi=1.;
00088 else if (cosphi<-1.) cosphi=-1.;
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099 G4double hel=(yParticleFrame.cross(nInteractionFrame)*particleDirection)>0?1.:-1.;
00100
00101 G4double sinphi=hel*std::sqrt(1.-cosphi*cosphi);
00102
00103
00104 RotateAz(cosphi,sinphi);
00105 }
00106
00107
00108 void G4StokesVector::InvRotateAz(G4ThreeVector nInteractionFrame,
00109 G4ThreeVector particleDirection)
00110 {
00111
00112
00113
00114
00115
00116 G4ThreeVector yParticleFrame =
00117 G4PolarizationHelper::GetParticleFrameY(particleDirection);
00118 G4double cosphi=yParticleFrame*nInteractionFrame;
00119
00120 if (cosphi>1.+1.e-8 || cosphi<-1.-1.e-8) {
00121 G4cout<<" warning G4StokesVector::RotateAz cosphi>1 or cosphi<-1\n";
00122 }
00123 if (cosphi>1) cosphi=1.;
00124 else if (cosphi<-1)cosphi=-1.;
00125
00126
00127 G4double hel=(yParticleFrame.cross(nInteractionFrame)*particleDirection)>0?1.:-1.;
00128 G4double sinphi=hel*std::sqrt(std::fabs(1.-cosphi*cosphi));
00129 RotateAz(cosphi,-sinphi);
00130 }
00131
00132 void G4StokesVector::RotateAz(G4double cosphi, G4double sinphi)
00133 {
00134 if (!isPhoton) {
00135 G4double xsi1= cosphi*p1() + sinphi*p2();
00136 G4double xsi2= -sinphi*p1() + cosphi*p2();
00137 setX(xsi1);
00138 setY(xsi2);
00139 return;
00140 }
00141
00142 G4double sin2phi=2.*cosphi*sinphi;
00143 G4double cos2phi=cosphi*cosphi-sinphi*sinphi;
00144
00145 G4double xsi1= cos2phi*p1() + sin2phi*p2();
00146 G4double xsi2= -sin2phi*p1() + cos2phi*p2();
00147 setX(xsi1);
00148 setY(xsi2);
00149 }
00150
00151 G4double G4StokesVector::GetBeta()
00152 {
00153 G4double bet=getPhi();
00154 if (isPhoton) { bet *= 0.5; }
00155 return bet;
00156 }
00157
00158 void G4StokesVector::DiceUniform()
00159 {
00160 G4double costheta=2.*G4UniformRand()-1.;
00161 G4double sintheta=std::sqrt(1.-costheta*costheta);
00162 G4double aphi =2.*pi*G4UniformRand();
00163 setX(std::sin(aphi)*sintheta);
00164 setY(std::cos(aphi)*sintheta);
00165 setZ(costheta);
00166 }
00167
00168 void G4StokesVector::DiceP1()
00169 {
00170 if (G4UniformRand()>0.5) setX(1.);
00171 else setX(-1.);
00172 setY(0.);
00173 setZ(0.);
00174 }
00175
00176 void G4StokesVector::DiceP2()
00177 {
00178 setX(0.);
00179 if (G4UniformRand()>0.5) setY(1.);
00180 else setY(-1.);
00181 setZ(0.);
00182 }
00183
00184 void G4StokesVector::DiceP3()
00185 {
00186 setX(0.);
00187 setY(0.);
00188 if (G4UniformRand()>0.5) setZ(1.);
00189 else setZ(-1.);
00190 }
00191
00192 void G4StokesVector::FlipP3()
00193 {
00194 setZ(-z());
00195 }
00196
00197 G4ThreeVector G4StokesVector::PolError(const G4StokesVector & sum2, long n)
00198 {
00199
00200 G4StokesVector mean=(1./n)*(*this);
00201 return G4StokesVector((1./(n-1.)*((1./n)*sum2 - mean.PolSqr()))).PolSqrt();
00202 }
00203
00204 G4ThreeVector G4StokesVector::PolDiv(const G4StokesVector & b)
00205 {return G4ThreeVector(b.x()!=0. ? x()/b.x() : 11111.,
00206 b.y()!=0. ? y()/b.y() : 11111.,
00207 b.z()!=0. ? z()/b.z() : 11111.);}