ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Rotation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Rotation.cc
1 // -*- C++ -*-
2 // ---------------------------------------------------------------------------
3 //
4 // This file is a part of the CLHEP - a Class Library for High Energy Physics.
5 //
6 // This is the implementation of the parts of the the HepRotation class which
7 // were present in the original CLHEP before the merge with ZOOM PhysicsVectors.
8 //
9 
10 #ifdef GNUPRAGMA
11 #pragma implementation
12 #endif
13 
14 #include "CLHEP/Vector/Rotation.h"
16 
17 #include <iostream>
18 #include <cmath>
19 
20 namespace CLHEP {
21 
22 static inline double safe_acos (double x) {
23  if (std::abs(x) <= 1.0) return std::acos(x);
24  return ( (x>0) ? 0 : CLHEP::pi );
25 }
26 
27 double HepRotation::operator() (int i, int j) const {
28  if (i == 0) {
29  if (j == 0) { return xx(); }
30  if (j == 1) { return xy(); }
31  if (j == 2) { return xz(); }
32  } else if (i == 1) {
33  if (j == 0) { return yx(); }
34  if (j == 1) { return yy(); }
35  if (j == 2) { return yz(); }
36  } else if (i == 2) {
37  if (j == 0) { return zx(); }
38  if (j == 1) { return zy(); }
39  if (j == 2) { return zz(); }
40  }
41  std::cerr << "HepRotation subscripting: bad indices "
42  << "(" << i << "," << j << ")" << std::endl;
43  return 0.0;
44 }
45 
46 HepRotation & HepRotation::rotate(double a, const Hep3Vector& aaxis) {
47  if (a != 0.0) {
48  double ll = aaxis.mag();
49  if (ll == 0.0) {
50  std::cerr << "HepRotation::rotate() - "
51  << "HepRotation: zero axis" << std::endl;
52  }else{
53  double sa = std::sin(a), ca = std::cos(a);
54  double dx = aaxis.x()/ll, dy = aaxis.y()/ll, dz = aaxis.z()/ll;
55  HepRotation m1(
56  ca+(1-ca)*dx*dx, (1-ca)*dx*dy-sa*dz, (1-ca)*dx*dz+sa*dy,
57  (1-ca)*dy*dx+sa*dz, ca+(1-ca)*dy*dy, (1-ca)*dy*dz-sa*dx,
58  (1-ca)*dz*dx-sa*dy, (1-ca)*dz*dy+sa*dx, ca+(1-ca)*dz*dz );
59  transform(m1);
60  }
61  }
62  return *this;
63 }
64 
66  double c1 = std::cos(a);
67  double s1 = std::sin(a);
68  double x1 = ryx, y1 = ryy, z1 = ryz;
69  ryx = c1*x1 - s1*rzx;
70  ryy = c1*y1 - s1*rzy;
71  ryz = c1*z1 - s1*rzz;
72  rzx = s1*x1 + c1*rzx;
73  rzy = s1*y1 + c1*rzy;
74  rzz = s1*z1 + c1*rzz;
75  return *this;
76 }
77 
79  double c1 = std::cos(a);
80  double s1 = std::sin(a);
81  double x1 = rzx, y1 = rzy, z1 = rzz;
82  rzx = c1*x1 - s1*rxx;
83  rzy = c1*y1 - s1*rxy;
84  rzz = c1*z1 - s1*rxz;
85  rxx = s1*x1 + c1*rxx;
86  rxy = s1*y1 + c1*rxy;
87  rxz = s1*z1 + c1*rxz;
88  return *this;
89 }
90 
92  double c1 = std::cos(a);
93  double s1 = std::sin(a);
94  double x1 = rxx, y1 = rxy, z1 = rxz;
95  rxx = c1*x1 - s1*ryx;
96  rxy = c1*y1 - s1*ryy;
97  rxz = c1*z1 - s1*ryz;
98  ryx = s1*x1 + c1*ryx;
99  ryy = s1*y1 + c1*ryy;
100  ryz = s1*z1 + c1*ryz;
101  return *this;
102 }
103 
105  const Hep3Vector &newY,
106  const Hep3Vector &newZ) {
107  double del = 0.001;
108  Hep3Vector w = newX.cross(newY);
109 
110  if (std::abs(newZ.x()-w.x()) > del ||
111  std::abs(newZ.y()-w.y()) > del ||
112  std::abs(newZ.z()-w.z()) > del ||
113  std::abs(newX.mag2()-1.) > del ||
114  std::abs(newY.mag2()-1.) > del ||
115  std::abs(newZ.mag2()-1.) > del ||
116  std::abs(newX.dot(newY)) > del ||
117  std::abs(newY.dot(newZ)) > del ||
118  std::abs(newZ.dot(newX)) > del) {
119  std::cerr << "HepRotation::rotateAxes: bad axis vectors" << std::endl;
120  return *this;
121  }else{
122  return transform(HepRotation(newX.x(), newY.x(), newZ.x(),
123  newX.y(), newY.y(), newZ.y(),
124  newX.z(), newY.z(), newZ.z()));
125  }
126 }
127 
128 double HepRotation::phiX() const {
129  return (yx() == 0.0 && xx() == 0.0) ? 0.0 : std::atan2(yx(),xx());
130 }
131 
132 double HepRotation::phiY() const {
133  return (yy() == 0.0 && xy() == 0.0) ? 0.0 : std::atan2(yy(),xy());
134 }
135 
136 double HepRotation::phiZ() const {
137  return (yz() == 0.0 && xz() == 0.0) ? 0.0 : std::atan2(yz(),xz());
138 }
139 
140 double HepRotation::thetaX() const {
141  return safe_acos(zx());
142 }
143 
144 double HepRotation::thetaY() const {
145  return safe_acos(zy());
146 }
147 
148 double HepRotation::thetaZ() const {
149  return safe_acos(zz());
150 }
151 
152 void HepRotation::getAngleAxis(double &angle, Hep3Vector &aaxis) const {
153  double cosa = 0.5*(xx()+yy()+zz()-1);
154  double cosa1 = 1-cosa;
155  if (cosa1 <= 0) {
156  angle = 0;
157  aaxis = Hep3Vector(0,0,1);
158  }else{
159  double x=0, y=0, z=0;
160  if (xx() > cosa) x = std::sqrt((xx()-cosa)/cosa1);
161  if (yy() > cosa) y = std::sqrt((yy()-cosa)/cosa1);
162  if (zz() > cosa) z = std::sqrt((zz()-cosa)/cosa1);
163  if (zy() < yz()) x = -x;
164  if (xz() < zx()) y = -y;
165  if (yx() < xy()) z = -z;
166  angle = (cosa < -1.) ? std::acos(-1.) : std::acos(cosa);
167  aaxis = Hep3Vector(x,y,z);
168  }
169 }
170 
172  return (rxx == 1.0 && rxy == 0.0 && rxz == 0.0 &&
173  ryx == 0.0 && ryy == 1.0 && ryz == 0.0 &&
174  rzx == 0.0 && rzy == 0.0 && rzz == 1.0) ? true : false;
175 }
176 
177 int HepRotation::compare ( const HepRotation & r ) const {
178  if (rzz<r.rzz) return -1; else if (rzz>r.rzz) return 1;
179  else if (rzy<r.rzy) return -1; else if (rzy>r.rzy) return 1;
180  else if (rzx<r.rzx) return -1; else if (rzx>r.rzx) return 1;
181  else if (ryz<r.ryz) return -1; else if (ryz>r.ryz) return 1;
182  else if (ryy<r.ryy) return -1; else if (ryy>r.ryy) return 1;
183  else if (ryx<r.ryx) return -1; else if (ryx>r.ryx) return 1;
184  else if (rxz<r.rxz) return -1; else if (rxz>r.rxz) return 1;
185  else if (rxy<r.rxy) return -1; else if (rxy>r.rxy) return 1;
186  else if (rxx<r.rxx) return -1; else if (rxx>r.rxx) return 1;
187  else return 0;
188 }
189 
190 
192 
193 } // namespace CLHEP
194 
195