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 #include "G4PolarizationHelper.hh"
00044 #include "G4PhysicalConstants.hh"
00045 #include "G4StokesVector.hh"
00046
00047
00048 G4ThreeVector G4PolarizationHelper::GetFrame(const G4ThreeVector & mom1, const G4ThreeVector & mom2)
00049 {
00050 G4ThreeVector normal = (mom1.cross(mom2)).unit();
00051 return normal;
00052
00053 }
00054
00055 G4ThreeVector G4PolarizationHelper::GetParticleFrameY(const G4ThreeVector &uZ)
00056 {
00057
00058
00059 if (uZ.x()==0. && uZ.y()==0.) {
00060 return G4ThreeVector(0.,1.,0.);
00061 }
00062
00063 G4double invPerp = 1./std::sqrt(sqr(uZ.x())+sqr(uZ.y()));
00064 return G4ThreeVector(-uZ.y()*invPerp,uZ.x()*invPerp,0);
00065 }
00066
00067 G4ThreeVector G4PolarizationHelper::GetParticleFrameX(const G4ThreeVector &uZ)
00068 {
00069
00070
00071 if (uZ.x()==0. && uZ.y()==0.) {
00072 if (uZ.z()>=0.) return G4ThreeVector(1.,0.,0.);
00073 return G4ThreeVector(-1.,0.,0.);
00074 }
00075
00076 G4double perp = std::sqrt(sqr(uZ.x())+sqr(uZ.y()));
00077 G4double invPerp = uZ.z()/perp;
00078 return G4ThreeVector(uZ.x()*invPerp,uZ.y()*invPerp,-perp);
00079 }
00080
00081 G4ThreeVector G4PolarizationHelper::GetRandomFrame(const G4ThreeVector & mom1)
00082 {
00083 G4double phi =2.*pi*G4UniformRand();
00084 G4ThreeVector normal = std::cos(phi)*GetParticleFrameX(mom1)
00085 + std::sin(phi)*G4PolarizationHelper::GetParticleFrameY(mom1);
00086 return normal;
00087 }
00088
00089
00090 G4ThreeVector G4PolarizationHelper::GetSpinInPRF(const G4ThreeVector &uZ, const G4ThreeVector & spin)
00091 {
00092
00093
00094 if (uZ.x()==0. && uZ.y()==0.) {
00095 if (uZ.z()>=0.) return spin;
00096 return G4ThreeVector(-spin.x(),spin.y(),-spin.z());
00097 }
00098
00099 G4double perp = std::sqrt(sqr(uZ.x())+sqr(uZ.y()));
00100 G4double invPerp = 1./perp;
00101
00102 G4ThreeVector uX(uZ.x()*uZ.z()*invPerp,uZ.y()*uZ.z()*invPerp,-perp);
00103 G4ThreeVector uY(-uZ.y()*invPerp,uZ.x()*invPerp,0);
00104
00105 return G4ThreeVector(spin*uX,spin*uY,spin*uZ);
00106 }
00107
00108 void G4PolarizationHelper::TestPolarizationTransformations()
00109 {
00110 G4double theta=0.;
00111 G4cout<<"========================================\n\n";
00112 for (G4int i=0; i<=10; ++i) {
00113 theta=pi*i/10.;
00114 G4ThreeVector zAxis = G4ThreeVector(std::sin(theta),0.,std::cos(theta));
00115 if (i==5) zAxis = G4ThreeVector(1.,0.,0.);
00116 if (i==10) zAxis = G4ThreeVector(0.,0.,-1.);
00117 G4ThreeVector yAxis = GetParticleFrameY(zAxis);
00118
00119 G4cout<<zAxis<<" "<<zAxis.mag()<<"\n";
00120 G4cout<<yAxis<<" "<<yAxis.mag()<<"\n";
00121 G4ThreeVector xAxis = yAxis.cross(zAxis);
00122 G4cout<<xAxis<<" "<<xAxis.mag()<<"\n\n";
00123 }
00124
00125 G4cout<<"========================================\n\n";
00126
00127 for (G4int i=0; i<=10; ++i) {
00128 theta=pi*i/10.;
00129 G4ThreeVector zAxis = G4ThreeVector(0.,std::sin(theta),std::cos(theta));
00130 if (i==5) zAxis = G4ThreeVector(0.,1.,0.);
00131 if (i==10) zAxis = G4ThreeVector(0.,0.,-1.);
00132 G4ThreeVector yAxis = GetParticleFrameY(zAxis);
00133
00134 G4cout<<zAxis<<" "<<zAxis.mag()<<"\n";
00135 G4cout<<yAxis<<" "<<yAxis.mag()<<"\n";
00136 G4ThreeVector xAxis = yAxis.cross(zAxis);
00137 G4cout<<xAxis<<" "<<xAxis.mag()<<"\n\n";
00138
00139 G4cout<<"spat : "<<xAxis*yAxis.cross(zAxis)<<"\n\n";
00140 }
00141 G4cout<<"========================================\n\n";
00142 }
00143
00144 void G4PolarizationHelper::TestInteractionFrame()
00145 {
00146
00147
00148
00149
00150
00151
00152 G4cout<<"========================================\n\n";
00153
00154 G4double theta=0.;
00155
00156 G4ThreeVector dir0=G4ThreeVector(0.,0.,1.);
00157 G4ThreeVector dir2=G4ThreeVector(std::sin(theta),0.,std::cos(theta));
00158
00159 G4StokesVector pol0=G4ThreeVector(0.,0.,1.);
00160 G4StokesVector pol1=G4ThreeVector(0.,0.,1.);
00161
00162 pol1.rotateUz(dir0);
00163
00164 G4cout<<"========================================\n\n";
00165
00166
00167 }