CMS 3D CMS Logo
Transform3DPJ.cc
Go to the documentation of this file.
1 // @(#)root/mathcore:$Name: V02-02-29 $:$Id: Transform3DPJ.cc,v 1.2 2008/01/22 20:41:27 muzaffar Exp $
2 // Authors: W. Brown, M. Fischler, L. Moneta 2005
3 
4 /**********************************************************************
5  * *
6  * Copyright (c) 2005 , LCG ROOT MathLib Team *
7  * *
8  * *
9  **********************************************************************/
10 
11 // implementation file for class Transform3D
12 //
13 // Created by: Lorenzo Moneta October 27 2005
14 //
15 //
16 
17 #include "FastSimulation/CaloGeometryTools/interface/Transform3DPJ.h"
18 
19 #include <cmath>
20 #include <algorithm>
21 
22 namespace ROOT {
23 
24  namespace Math {
25 
26  typedef Transform3DPJ::Point XYZPoint;
27  typedef Transform3DPJ::Vector XYZVector;
28 
29  // ========== Constructors and Assignment =====================
30 
31  // construct from two ref frames
32  Transform3DPJ::Transform3DPJ(const XYZPoint &fr0,
33  const XYZPoint &fr1,
34  const XYZPoint &fr2,
35  const XYZPoint &to0,
36  const XYZPoint &to1,
37  const XYZPoint &to2) {
38  // takes impl. from CLHEP ( E.Chernyaev). To be checked
39 
40  XYZVector x1, y1, z1, x2, y2, z2;
41  x1 = (fr1 - fr0).Unit();
42  y1 = (fr2 - fr0).Unit();
43  x2 = (to1 - to0).Unit();
44  y2 = (to2 - to0).Unit();
45 
46  // C H E C K A N G L E S
47 
48  double cos1, cos2;
49  cos1 = x1.Dot(y1);
50  cos2 = x2.Dot(y2);
51 
52  if (std::fabs(1.0 - cos1) <= 0.000001 || std::fabs(1.0 - cos2) <= 0.000001) {
53  std::cerr << "Transform3DPJ: Error : zero angle between axes" << std::endl;
54  SetIdentity();
55  } else {
56  if (std::fabs(cos1 - cos2) > 0.000001) {
57  std::cerr << "Transform3DPJ: Warning: angles between axes are not equal" << std::endl;
58  }
59 
60  // F I N D R O T A T I O N M A T R I X
61 
62  z1 = (x1.Cross(y1)).Unit();
63  y1 = z1.Cross(x1);
64 
65  z2 = (x2.Cross(y2)).Unit();
66  y2 = z2.Cross(x2);
67 
68  double x1x = x1.X(), x1y = x1.Y(), x1z = x1.Z();
69  double y1x = y1.X(), y1y = y1.Y(), y1z = y1.Z();
70  double z1x = z1.X(), z1y = z1.Y(), z1z = z1.Z();
71 
72  double detxx = (y1y * z1z - z1y * y1z);
73  double detxy = -(y1x * z1z - z1x * y1z);
74  double detxz = (y1x * z1y - z1x * y1y);
75  double detyx = -(x1y * z1z - z1y * x1z);
76  double detyy = (x1x * z1z - z1x * x1z);
77  double detyz = -(x1x * z1y - z1x * x1y);
78  double detzx = (x1y * y1z - y1y * x1z);
79  double detzy = -(x1x * y1z - y1x * x1z);
80  double detzz = (x1x * y1y - y1x * x1y);
81 
82  double x2x = x2.X(), x2y = x2.Y(), x2z = x2.Z();
83  double y2x = y2.X(), y2y = y2.Y(), y2z = y2.Z();
84  double z2x = z2.X(), z2y = z2.Y(), z2z = z2.Z();
85 
86  double txx = x2x * detxx + y2x * detyx + z2x * detzx;
87  double txy = x2x * detxy + y2x * detyy + z2x * detzy;
88  double txz = x2x * detxz + y2x * detyz + z2x * detzz;
89  double tyx = x2y * detxx + y2y * detyx + z2y * detzx;
90  double tyy = x2y * detxy + y2y * detyy + z2y * detzy;
91  double tyz = x2y * detxz + y2y * detyz + z2y * detzz;
92  double tzx = x2z * detxx + y2z * detyx + z2z * detzx;
93  double tzy = x2z * detxy + y2z * detyy + z2z * detzy;
94  double tzz = x2z * detxz + y2z * detyz + z2z * detzz;
95 
96  // S E T T R A N S F O R M A T I O N
97 
98  double dx1 = fr0.X(), dy1 = fr0.Y(), dz1 = fr0.Z();
99  double dx2 = to0.X(), dy2 = to0.Y(), dz2 = to0.Z();
100 
101  SetComponents(txx,
102  txy,
103  txz,
104  dx2 - txx * dx1 - txy * dy1 - txz * dz1,
105  tyx,
106  tyy,
107  tyz,
108  dy2 - tyx * dx1 - tyy * dy1 - tyz * dz1,
109  tzx,
110  tzy,
111  tzz,
112  dz2 - tzx * dx1 - tzy * dy1 - tzz * dz1);
113  }
114  }
115 
116  // inversion (from CLHEP)
117  void Transform3DPJ::Invert() {
118  //
119  // Name: Transform3DPJ::inverse Date: 24.09.96
120  // Author: E.Chernyaev (IHEP/Protvino) Revised:
121  //
122  // Function: Find inverse affine transformation.
123 
124  double detxx = fM[kYY] * fM[kZZ] - fM[kYZ] * fM[kZY];
125  double detxy = fM[kYX] * fM[kZZ] - fM[kYZ] * fM[kZX];
126  double detxz = fM[kYX] * fM[kZY] - fM[kYY] * fM[kZX];
127  double det = fM[kXX] * detxx - fM[kXY] * detxy + fM[kXZ] * detxz;
128  if (det == 0) {
129  std::cerr << "Transform3DPJ::inverse error: zero determinant" << std::endl;
130  return;
131  }
132  det = 1. / det;
133  detxx *= det;
134  detxy *= det;
135  detxz *= det;
136  double detyx = (fM[kXY] * fM[kZZ] - fM[kXZ] * fM[kZY]) * det;
137  double detyy = (fM[kXX] * fM[kZZ] - fM[kXZ] * fM[kZX]) * det;
138  double detyz = (fM[kXX] * fM[kZY] - fM[kXY] * fM[kZX]) * det;
139  double detzx = (fM[kXY] * fM[kYZ] - fM[kXZ] * fM[kYY]) * det;
140  double detzy = (fM[kXX] * fM[kYZ] - fM[kXZ] * fM[kYX]) * det;
141  double detzz = (fM[kXX] * fM[kYY] - fM[kXY] * fM[kYX]) * det;
142  SetComponents(detxx,
143  -detyx,
144  detzx,
145  -detxx * fM[kDX] + detyx * fM[kDY] - detzx * fM[kDZ],
146  -detxy,
147  detyy,
148  -detzy,
149  detxy * fM[kDX] - detyy * fM[kDY] + detzy * fM[kDZ],
150  detxz,
151  -detyz,
152  detzz,
153  -detxz * fM[kDX] + detyz * fM[kDY] - detzz * fM[kDZ]);
154  }
155 
156  // get rotations and translations
157  void Transform3DPJ::GetDecomposition(Rotation3D &r, XYZVector &v) const {
158  // decompose a trasfomation in a 3D rotation and in a 3D vector (cartesian coordinates)
159  r.SetComponents(fM[kXX], fM[kXY], fM[kXZ], fM[kYX], fM[kYY], fM[kYZ], fM[kZX], fM[kZY], fM[kZZ]);
160 
161  v.SetCoordinates(fM[kDX], fM[kDY], fM[kDZ]);
162  }
163 
164  // transformation on Position Vector (rotation + translations)
165  XYZPoint Transform3DPJ::operator()(const XYZPoint &p) const {
166  // pass through rotation class (could be implemented directly to be faster)
167 
168  Rotation3D r;
169  XYZVector t;
170  GetDecomposition(r, t);
171  XYZPoint pnew = r(p);
172  pnew += t;
173  return pnew;
174  }
175 
176  // transformation on Displacement Vector (only rotation)
177  XYZVector Transform3DPJ::operator()(const XYZVector &v) const {
178  // pass through rotation class ( could be implemented directly to be faster)
179 
180  Rotation3D r;
181  XYZVector t;
182  GetDecomposition(r, t);
183  // only rotation
184  return r(v);
185  }
186 
187  Transform3DPJ &Transform3DPJ::operator*=(const Transform3DPJ &t) {
188  // combination of transformations
189 
190  SetComponents(fM[kXX] * t.fM[kXX] + fM[kXY] * t.fM[kYX] + fM[kXZ] * t.fM[kZX],
191  fM[kXX] * t.fM[kXY] + fM[kXY] * t.fM[kYY] + fM[kXZ] * t.fM[kZY],
192  fM[kXX] * t.fM[kXZ] + fM[kXY] * t.fM[kYZ] + fM[kXZ] * t.fM[kZZ],
193  fM[kXX] * t.fM[kDX] + fM[kXY] * t.fM[kDY] + fM[kXZ] * t.fM[kDZ] + fM[kDX],
194 
195  fM[kYX] * t.fM[kXX] + fM[kYY] * t.fM[kYX] + fM[kYZ] * t.fM[kZX],
196  fM[kYX] * t.fM[kXY] + fM[kYY] * t.fM[kYY] + fM[kYZ] * t.fM[kZY],
197  fM[kYX] * t.fM[kXZ] + fM[kYY] * t.fM[kYZ] + fM[kYZ] * t.fM[kZZ],
198  fM[kYX] * t.fM[kDX] + fM[kYY] * t.fM[kDY] + fM[kYZ] * t.fM[kDZ] + fM[kDY],
199 
200  fM[kZX] * t.fM[kXX] + fM[kZY] * t.fM[kYX] + fM[kZZ] * t.fM[kZX],
201  fM[kZX] * t.fM[kXY] + fM[kZY] * t.fM[kYY] + fM[kZZ] * t.fM[kZY],
202  fM[kZX] * t.fM[kXZ] + fM[kZY] * t.fM[kYZ] + fM[kZZ] * t.fM[kZZ],
203  fM[kZX] * t.fM[kDX] + fM[kZY] * t.fM[kDY] + fM[kZZ] * t.fM[kDZ] + fM[kDZ]);
204 
205  return *this;
206  }
207 
208  void Transform3DPJ::SetIdentity() {
209  //set identity ( identity rotation and zero translation)
210  fM[kXX] = 1.0;
211  fM[kXY] = 0.0;
212  fM[kXZ] = 0.0;
213  fM[kDX] = 0.0;
214  fM[kYX] = 0.0;
215  fM[kYY] = 1.0;
216  fM[kYZ] = 0.0;
217  fM[kDY] = 0.0;
218  fM[kZX] = 0.0;
219  fM[kZY] = 0.0;
220  fM[kZZ] = 1.0;
221  fM[kDZ] = 0.0;
222  }
223 
224  void Transform3DPJ::AssignFrom(const Rotation3D &r, const XYZVector &v) {
225  // assignment from rotation + translation
226 
227  double rotData[9];
228  r.GetComponents(rotData, rotData + 9);
229  // first raw
230  for (int i = 0; i < 3; ++i)
231  fM[i] = rotData[i];
232  // second raw
233  for (int i = 0; i < 3; ++i)
234  fM[kYX + i] = rotData[3 + i];
235  // third raw
236  for (int i = 0; i < 3; ++i)
237  fM[kZX + i] = rotData[6 + i];
238 
239  // translation data
240  double vecData[3];
241  v.GetCoordinates(vecData, vecData + 3);
242  fM[kDX] = vecData[0];
243  fM[kDY] = vecData[1];
244  fM[kDZ] = vecData[2];
245  }
246 
247  void Transform3DPJ::AssignFrom(const Rotation3D &r) {
248  // assign from only a rotation (null translation)
249  double rotData[9];
250  r.GetComponents(rotData, rotData + 9);
251  for (int i = 0; i < 3; ++i) {
252  for (int j = 0; j < 3; ++j)
253  fM[4 * i + j] = rotData[3 * i + j];
254  // empty vector data
255  fM[4 * i + 3] = 0;
256  }
257  }
258 
259  void Transform3DPJ::AssignFrom(const XYZVector &v) {
260  // assign from a translation only (identity rotations)
261  fM[kXX] = 1.0;
262  fM[kXY] = 0.0;
263  fM[kXZ] = 0.0;
264  fM[kDX] = v.X();
265  fM[kYX] = 0.0;
266  fM[kYY] = 1.0;
267  fM[kYZ] = 0.0;
268  fM[kDY] = v.Y();
269  fM[kZX] = 0.0;
270  fM[kZY] = 0.0;
271  fM[kZZ] = 1.0;
272  fM[kDZ] = v.Z();
273  }
274 
275  Plane3D Transform3DPJ::operator()(const Plane3D &plane) const {
276  // transformations on a 3D plane
277  XYZVector n = plane.Normal();
278  // take a point on the plane. Use origin projection on the plane
279  // ( -ad, -bd, -cd) if (a**2 + b**2 + c**2 ) = 1
280  double d = plane.HesseDistance();
281  XYZPoint p(-d * n.X(), -d * n.Y(), -d * n.Z());
282  return Plane3D(operator()(n), operator()(p));
283  }
284 
285  std::ostream &operator<<(std::ostream &os, const Transform3DPJ &t) {
286  // TODO - this will need changing for machine-readable issues
287  // and even the human readable form needs formatiing improvements
288 
289  double m[12];
290  t.GetComponents(m, m + 12);
291  os << "\n" << m[0] << " " << m[1] << " " << m[2] << " " << m[3];
292  os << "\n" << m[4] << " " << m[5] << " " << m[6] << " " << m[7];
293  os << "\n" << m[8] << " " << m[9] << " " << m[10] << " " << m[11] << "\n";
294  return os;
295  }
296 
297  } // end namespace Math
298 } // end namespace ROOT
ROOT::Math::Transform3DPJ::fM
double fM[12]
Definition: Transform3DPJ.h:522
ROOT::Math::Transform3DPJ::Vector
DisplacementVector3D< Cartesian3D< double >, DefaultCoordinateSystemTag > Vector
Definition: Transform3DPJ.h:64
mps_fire.i
i
Definition: mps_fire.py:428
dqmiodumpmetadata.n
n
Definition: dqmiodumpmetadata.py:28
ROOT::Math::Transform3DPJ::AssignFrom
void AssignFrom(const Rotation3D &r, const Vector &v)
Definition: Transform3DPJ.cc:231
ROOT::Math::XYZPoint
Transform3DPJ::Point XYZPoint
Definition: Transform3DPJ.cc:33
ROOT::Math::Transform3DPJ::kXX
Definition: Transform3DPJ.h:68
testProducerWithPsetDescEmpty_cfi.x2
x2
Definition: testProducerWithPsetDescEmpty_cfi.py:28
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
ROOT::Math::Transform3DPJ::operator*=
Transform3DPJ & operator*=(const Transform3DPJ &t)
Definition: Transform3DPJ.cc:194
XYZVector
math::XYZVector XYZVector
Definition: RawParticle.h:26
ROOT::Math::Transform3DPJ::GetDecomposition
void GetDecomposition(Rotation3D &r, Vector &v) const
Definition: Transform3DPJ.cc:164
ROOT::Math::Transform3DPJ::kYZ
Definition: Transform3DPJ.h:74
findQualityFiles.v
v
Definition: findQualityFiles.py:179
ROOT::Math::Transform3DPJ::kXZ
Definition: Transform3DPJ.h:70
Plane3D
ROOT::Math::Plane3D Plane3D
Definition: PreshowerHitMaker.cc:13
testProducerWithPsetDescEmpty_cfi.z2
z2
Definition: testProducerWithPsetDescEmpty_cfi.py:41
ROOT::Math::Transform3DPJ::kZZ
Definition: Transform3DPJ.h:78
testProducerWithPsetDescEmpty_cfi.x1
x1
Definition: testProducerWithPsetDescEmpty_cfi.py:33
testProducerWithPsetDescEmpty_cfi.y1
y1
Definition: testProducerWithPsetDescEmpty_cfi.py:29
visualization-live-secondInstance_cfg.m
m
Definition: visualization-live-secondInstance_cfg.py:72
ROOT::Math::Transform3DPJ::Transform3DPJ
Transform3DPJ()
Definition: Transform3DPJ.h:85
ROOT::Math::Transform3DPJ::Point
PositionVector3D< Cartesian3D< double >, DefaultCoordinateSystemTag > Point
Definition: Transform3DPJ.h:65
testProducerWithPsetDescEmpty_cfi.y2
y2
Definition: testProducerWithPsetDescEmpty_cfi.py:30
ROOT::Math::Transform3DPJ::SetIdentity
void SetIdentity()
Definition: Transform3DPJ.cc:215
ROOT::Math::operator<<
std::ostream & operator<<(std::ostream &os, const Transform3DPJ &t)
Definition: Transform3DPJ.cc:292
ROOT::Math::Transform3DPJ::SetComponents
void SetComponents(IT begin, IT end)
Definition: Transform3DPJ.h:227
ROOT::Math::Transform3DPJ::operator()
Point operator()(const Point &p) const
ROOT::Math::Transform3DPJ::kYY
Definition: Transform3DPJ.h:73
ROOT::Math::Transform3DPJ
Definition: Transform3DPJ.h:62
XYZPoint
math::XYZVector XYZPoint
Definition: CalorimetryManager.cc:69
alignCSCRings.r
r
Definition: alignCSCRings.py:93
ROOT::Math::Transform3DPJ::kZY
Definition: Transform3DPJ.h:77
Transform3DPJ.h
ROOT::Math::Transform3DPJ::kYX
Definition: Transform3DPJ.h:72
ROOT::Math::Transform3DPJ::kXY
Definition: Transform3DPJ.h:69
ROOT::Math::Transform3DPJ::Invert
void Invert()
Definition: Transform3DPJ.cc:124
ROOT::Math::Transform3DPJ::kDY
Definition: Transform3DPJ.h:75
ROOT::Math::XYZVector
Transform3DPJ::Vector XYZVector
Definition: Transform3DPJ.cc:34
ztail.d
d
Definition: ztail.py:151
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
submitPVValidationJobs.t
string t
Definition: submitPVValidationJobs.py:644
EcnaPython_AdcPeg12_S1_10_R170298_1_0_150_Dee0.cerr
cerr
Definition: EcnaPython_AdcPeg12_S1_10_R170298_1_0_150_Dee0.py:8
ROOT::Math::Transform3DPJ::kDX
Definition: Transform3DPJ.h:71
ROOT
Definition: Transform3DPJ.h:35
ROOT::Math::Transform3DPJ::kZX
Definition: Transform3DPJ.h:76
ROOT::Math::Transform3DPJ::kDZ
Definition: Transform3DPJ.h:79
Generated for CMSSW Reference Manual by doxygen 1.8.16