File indexing completed on 2024-05-10 02:20:35
0001 #include <cmath>
0002 #include <iomanip>
0003 #include <iostream>
0004
0005 #include <CLHEP/Units/SystemOfUnits.h>
0006 #include <CLHEP/Vector/Rotation.h>
0007 #include <CLHEP/Vector/RotationInterfaces.h>
0008 #include <CLHEP/Vector/ThreeVector.h>
0009 #include "DetectorDescription/Core/interface/DDRotationMatrix.h"
0010 #include "DetectorDescription/Core/interface/DDTranslation.h"
0011 #include "Math/GenVector/Cartesian3D.h"
0012 #include "Math/GenVector/DisplacementVector3D.h"
0013 #include "Math/GenVector/Rotation3D.h"
0014
0015 typedef CLHEP::Hep3Vector H3V;
0016 typedef CLHEP::HepRotation HRM;
0017 using CLHEP::deg;
0018
0019 using namespace std;
0020
0021 void hrmOut(const HRM& r) {
0022 cout << "row 1 = " << std::setw(12) << std::fixed << std::setprecision(5) << r.xx();
0023 cout << ", " << std::setw(12) << std::fixed << std::setprecision(5) << r.xy();
0024 cout << ", " << std::setw(12) << std::fixed << std::setprecision(5) << r.xz() << endl;
0025 cout << "row 2 = " << std::setw(12) << std::fixed << std::setprecision(5) << r.yx();
0026 cout << ", " << std::setw(12) << std::fixed << std::setprecision(5) << r.yy();
0027 cout << ", " << std::setw(12) << std::fixed << std::setprecision(5) << r.yz() << endl;
0028 cout << "row 3 = " << std::setw(12) << std::fixed << std::setprecision(5) << r.zx();
0029 cout << ", " << std::setw(12) << std::fixed << std::setprecision(5) << r.zy();
0030 cout << ", " << std::setw(12) << std::fixed << std::setprecision(5) << r.zz() << endl;
0031 }
0032
0033 void hepSetAxes(
0034 H3V& x, H3V& y, H3V& z, double thetaX, double phiX, double thetaY, double phiY, double thetaZ, double phiZ) {
0035 x[0] = sin(thetaX) * cos(phiX);
0036 x[1] = sin(thetaX) * sin(phiX);
0037 x[2] = cos(thetaX);
0038 y[0] = sin(thetaY) * cos(phiY);
0039 y[1] = sin(thetaY) * sin(phiY);
0040 y[2] = cos(thetaY);
0041 z[0] = sin(thetaZ) * cos(phiZ);
0042 z[1] = sin(thetaZ) * sin(phiZ);
0043 z[2] = cos(thetaZ);
0044 }
0045
0046 void hepOutVecs(const H3V& x, const H3V& y, const H3V& z) {
0047 cout << "Vectors used in construction:" << endl;
0048 cout << "x vector = " << std::setw(12) << std::fixed << std::setprecision(5) << x[0] << ", " << std::setw(12)
0049 << std::fixed << std::setprecision(5) << y[0] << ", " << std::setw(12) << std::fixed << std::setprecision(5)
0050 << z[0] << endl;
0051 cout << "y vector = " << std::setw(12) << std::fixed << std::setprecision(5) << x[1] << ", " << std::setw(12)
0052 << std::fixed << std::setprecision(5) << y[1] << ", " << std::setw(12) << std::fixed << std::setprecision(5)
0053 << z[1] << endl;
0054 cout << "z vector = " << std::setw(12) << std::fixed << std::setprecision(5) << x[2] << ", " << std::setw(12)
0055 << std::fixed << std::setprecision(5) << y[2] << ", " << std::setw(12) << std::fixed << std::setprecision(5)
0056 << z[2] << endl;
0057 }
0058
0059 void ddRotOut(const DDRotationMatrix& r) {
0060 double xx, xy, xz, yx, yy, yz, zx, zy, zz;
0061 r.GetComponents(xx, xy, xz, yx, yy, yz, zx, zy, zz);
0062 cout << "row 1 = " << std::setw(12) << std::fixed << std::setprecision(5) << xx;
0063 cout << ", " << std::setw(12) << std::fixed << std::setprecision(5) << xy;
0064 cout << ", " << std::setw(12) << std::fixed << std::setprecision(5) << xz << endl;
0065 cout << "row 2 = " << std::setw(12) << std::fixed << std::setprecision(5) << yx;
0066 cout << ", " << std::setw(12) << std::fixed << std::setprecision(5) << yy;
0067 cout << ", " << std::setw(12) << std::fixed << std::setprecision(5) << yz << endl;
0068 cout << "row 3 = " << std::setw(12) << std::fixed << std::setprecision(5) << zx;
0069 cout << ", " << std::setw(12) << std::fixed << std::setprecision(5) << zy;
0070 cout << ", " << std::setw(12) << std::fixed << std::setprecision(5) << zz << endl;
0071 cout << "USING INTERNAL << OPERATOR" << endl;
0072 cout << r << endl;
0073 cout << "USING VECTOR COMPONENTS" << endl;
0074 DD3Vector x, y, z;
0075 r.GetComponents(x, y, z);
0076 cout << "col 1 = " << std::setw(12) << std::fixed << std::setprecision(5) << x.X();
0077 cout << ", " << std::setw(12) << std::fixed << std::setprecision(5) << x.Y();
0078 cout << ", " << std::setw(12) << std::fixed << std::setprecision(5) << x.Z() << endl;
0079 cout << "col 2 = " << std::setw(12) << std::fixed << std::setprecision(5) << y.X();
0080 cout << ", " << std::setw(12) << std::fixed << std::setprecision(5) << y.Y();
0081 cout << ", " << std::setw(12) << std::fixed << std::setprecision(5) << y.Z() << endl;
0082 cout << "col 3 = " << std::setw(12) << std::fixed << std::setprecision(5) << z.X();
0083 cout << ", " << std::setw(12) << std::fixed << std::setprecision(5) << z.Y();
0084 cout << ", " << std::setw(12) << std::fixed << std::setprecision(5) << z.Z() << endl;
0085 }
0086
0087 void ddSetAxes(DD3Vector& x,
0088 DD3Vector& y,
0089 DD3Vector& z,
0090 double thetaX,
0091 double phiX,
0092 double thetaY,
0093 double phiY,
0094 double thetaZ,
0095 double phiZ) {
0096 x.SetX(sin(thetaX) * cos(phiX));
0097 x.SetY(sin(thetaX) * sin(phiX));
0098 x.SetZ(cos(thetaX));
0099 y.SetX(sin(thetaY) * cos(phiY));
0100 y.SetY(sin(thetaY) * sin(phiY));
0101 y.SetZ(cos(thetaY));
0102 z.SetX(sin(thetaZ) * cos(phiZ));
0103 z.SetY(sin(thetaZ) * sin(phiZ));
0104 z.SetZ(cos(thetaZ));
0105 }
0106
0107 void ddOutVecs(const DD3Vector& x, const DD3Vector& y, const DD3Vector& z) {
0108 cout << "Vectors used in construction:" << endl;
0109 cout << "x vector = " << std::setw(12) << std::fixed << std::setprecision(5) << x.X() << ", " << std::setw(12)
0110 << std::fixed << std::setprecision(5) << x.Y() << ", " << std::setw(12) << std::fixed << std::setprecision(5)
0111 << x.Z() << endl;
0112 cout << "y vector = " << std::setw(12) << std::fixed << std::setprecision(5) << y.X() << ", " << std::setw(12)
0113 << std::fixed << std::setprecision(5) << y.Y() << ", " << std::setw(12) << std::fixed << std::setprecision(5)
0114 << y.Z() << endl;
0115 cout << "z vector = " << std::setw(12) << std::fixed << std::setprecision(5) << z.X() << ", " << std::setw(12)
0116 << std::fixed << std::setprecision(5) << z.Y() << ", " << std::setw(12) << std::fixed << std::setprecision(5)
0117 << z.Z() << endl;
0118 }
0119
0120 void checkNorm(double check) {
0121 double tol = 1.0e-3;
0122 if (1.0 - std::abs(check) > tol) {
0123 cout << "NOT orthonormal!" << endl;
0124 } else if (1.0 + check <= tol) {
0125 cout << "IS Left-handed (reflection)" << endl;
0126 } else {
0127 cout << "IS Right-handed (proper)" << endl;
0128 }
0129 }
0130
0131 int main(int , char** ) {
0132 cout << "====================== CLHEP: ========================" << endl;
0133
0134
0135 {
0136 H3V x, y, z;
0137 hepSetAxes(x, y, z, 90 * deg, 0 * deg, 90 * deg, 90 * deg, 180 * deg, 0 * deg);
0138 CLHEP::HepRotation R;
0139 R.rotateAxes(x, y, z);
0140 HRM ddr(R);
0141 cout << " *** REFLECTION *** " << endl;
0142 checkNorm((x.cross(y)) * z);
0143 hepOutVecs(x, y, z);
0144 cout << "Matrix output built up from vectors:" << endl;
0145 hrmOut(ddr);
0146 cout << "Matrix build from CLHEP::HepRep3x3 to preserve left-handedness:" << endl;
0147 CLHEP::HepRep3x3 temp(x.x(), y.x(), z.x(), x.y(), y.y(), z.y(), x.z(), y.z(), z.z());
0148 HRM ddr2(temp);
0149 hrmOut(ddr2);
0150 }
0151
0152 {
0153 H3V x, y, z;
0154 hepSetAxes(x, y, z, 90 * deg, -51.39999 * deg, 90 * deg, 38.60001 * deg, 0 * deg, 0 * deg);
0155 CLHEP::HepRotation R;
0156 R.rotateAxes(x, y, z);
0157 HRM ddr(R);
0158 cout << " *** ROTATION *** " << endl;
0159 checkNorm((x.cross(y)) * z);
0160 hepOutVecs(x, y, z);
0161 cout << "Matrix output built up from vectors:" << endl;
0162 hrmOut(ddr);
0163 cout << "Matrix build from CLHEP::HepRep3x3 to preserve left-handedness:" << endl;
0164 CLHEP::HepRep3x3 temp(x.x(), y.x(), z.x(), x.y(), y.y(), z.y(), x.z(), y.z(), z.z());
0165 HRM ddr2(temp);
0166 hrmOut(ddr2);
0167 }
0168
0169 cout << "====================== ROOT::Math ========================" << endl;
0170
0171 {
0172 DD3Vector x, y, z;
0173 ddSetAxes(x, y, z, 90 * deg, 0 * deg, 90 * deg, 90 * deg, 180 * deg, 0 * deg);
0174 DDRotationMatrix R(x, y, z);
0175 cout << " *** REFLECTION *** " << endl;
0176 checkNorm((x.Cross(y)).Dot(z));
0177 ddOutVecs(x, y, z);
0178 cout << "Matrix output built up from vectors:" << endl;
0179 ddRotOut(R);
0180 cout << "Matrix built to preserve left-handedness:" << endl;
0181 DDRotationMatrix temp(x.x(), y.x(), z.x(), x.y(), y.y(), z.y(), x.z(), y.z(), z.z());
0182 ddRotOut(temp);
0183 }
0184
0185 {
0186 DD3Vector x, y, z;
0187 ddSetAxes(x, y, z, 90 * deg, -51.39999 * deg, 90 * deg, 38.60001 * deg, 0 * deg, 0 * deg);
0188 DDRotationMatrix R(x, y, z);
0189 cout << " *** ROTATION *** " << endl;
0190 checkNorm((x.Cross(y)).Dot(z));
0191 ddOutVecs(x, y, z);
0192 cout << "Matrix output built up from vectors:" << endl;
0193 ddRotOut(R);
0194 cout << "Matrix built to preserve left-handedness:" << endl;
0195 DDRotationMatrix temp(x.x(), y.x(), z.x(), x.y(), y.y(), z.y(), x.z(), y.z(), z.z());
0196 ddRotOut(temp);
0197 }
0198 }