Back to home page

Project CMSSW displayed by LXR

 
 

    


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 /*argc*/, char** /*argv[]*/) {
0132   cout << "====================== CLHEP: ========================" << endl;
0133   // Examples from DD XML
0134   // <ReflectionRotation name="180R" thetaX="90*deg" phiX="0*deg" thetaY="90*deg" phiY="90*deg" thetaZ="180*deg" phiZ="0*deg" />
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());  //matrix representation
0148     HRM ddr2(temp);
0149     hrmOut(ddr2);
0150   }
0151   // <Rotation name="RM1509" thetaX="90*deg" phiX="-51.39999*deg" thetaY="90*deg" phiY="38.60001*deg" thetaZ="0*deg" phiZ="0*deg" />
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());  //matrix representation
0165     HRM ddr2(temp);
0166     hrmOut(ddr2);
0167   }
0168 
0169   cout << "====================== ROOT::Math ========================" << endl;
0170   // <ReflectionRotation name="180R" thetaX="90*deg" phiX="0*deg" thetaY="90*deg" phiY="90*deg" thetaZ="180*deg" phiZ="0*deg" />
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());  //matrix representation
0182     ddRotOut(temp);
0183   }
0184   // <Rotation name="RM1509" thetaX="90*deg" phiX="-51.39999*deg" thetaY="90*deg" phiY="38.60001*deg" thetaZ="0*deg" phiZ="0*deg" />
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());  //matrix representation
0196     ddRotOut(temp);
0197   }
0198 }