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 #include "G3toG4.hh"
00030 #include "G3RotTable.hh"
00031 #include "G3toG4RotationMatrix.hh"
00032 #include "G4PhysicalConstants.hh"
00033 #include "G4ThreeVector.hh"
00034
00035 void PG4gsrotm(G4String *tokens)
00036 {
00037
00038 G3fillParams(tokens,PTgsrotm);
00039
00040
00041 G4int irot = Ipar[0];
00042
00043
00044 G4double theta1 = Rpar[0];
00045 G4double phi1 = Rpar[1];
00046 G4double theta2 = Rpar[2];
00047 G4double phi2 = Rpar[3];
00048 G4double theta3 = Rpar[4];
00049 G4double phi3 = Rpar[5];
00050
00051 G4gsrotm(irot, theta1,phi1, theta2,phi2, theta3,phi3);
00052 }
00053
00054 void G4gsrotm(G4int irot, G4double theta1, G4double phi1,
00055 G4double theta2, G4double phi2, G4double theta3, G4double phi3)
00056 {
00057 G4double degrad = pi/180;
00058
00059 G4double th1r = theta1*degrad;
00060 G4double th2r = theta2*degrad;
00061 G4double th3r = theta3*degrad;
00062
00063 G4double phi1r = phi1*degrad;
00064 G4double phi2r = phi2*degrad;
00065 G4double phi3r = phi3*degrad;
00066
00067
00068
00069 G4ThreeVector x(std::sin(th1r)*std::cos(phi1r), std::sin(th1r)*std::sin(phi1r), std::cos(th1r));
00070 G4ThreeVector y(std::sin(th2r)*std::cos(phi2r), std::sin(th2r)*std::sin(phi2r), std::cos(th2r));
00071 G4ThreeVector z(std::sin(th3r)*std::cos(phi3r), std::sin(th3r)*std::sin(phi3r), std::cos(th3r));
00072
00073
00074
00075 G4double check = (x.cross(y))*z;
00076 G4double tol = 1.0e-3;
00077
00078 if (1-std::abs(check)>tol) {
00079 G4cerr << "Coordinate axes forming rotation matrix "
00080 << irot << " are not orthonormal.(" << 1-std::abs(check) << ")"
00081 << G4endl;
00082 G4cerr << " theta1=" << theta1;
00083 G4cerr << " phi1=" << phi1;
00084 G4cerr << " theta2=" << theta2;
00085 G4cerr << " phi2=" << phi2;
00086 G4cerr << " theta3=" << theta3;
00087 G4cerr << " phi3=" << phi3;
00088 G4cerr << G4endl;
00089 G4Exception("G4gsrotm()", "G3toG40023", FatalException,
00090 "Non orthogonal axes!");
00091 return;
00092 }
00093
00094
00095
00096
00097
00098 G3toG4RotationMatrix* rotp = new G3toG4RotationMatrix;
00099
00100 rotp->SetRotationMatrixByRow(x, y, z);
00101
00102
00103
00104 G3Rot.Put(irot, rotp);
00105 }