TrackBasis.cxx
Go to the documentation of this file.
1 ///
2 /// \file TrackBasis.cxx
3 /// \brief Construct a track-based (x,y,z) coordinate system
4 /// \author messier@indiana.edu
5 /// \version $Id:$
6 ///
9 #include "TVector3.h"
10 #include "Utilities/func/MathUtil.h"
11 using namespace bpfit;
12 
13 #define SQR(x) ((x)*(x))
14 
15 void TrackBasis::FindAxis(const std::vector<double>& x,
16  const std::vector<double>& y,
17  const std::vector<double>& z,
18  const std::vector<double>& dx,
19  const std::vector<double>& dy,
20  const std::vector<double>& dz,
21  TVector3& z1)
22 {
23  unsigned int i;
24  std::vector<double> w(z.size()); // Hit weights in x/z then y/z views
25 
26  //
27  // Space point vectors must all have same number of entries
28  //
29  if(x.size() != z.size() || y.size() != z.size() ||
30  dx.size() != x.size() || dy.size() != y.size() ||
31  dz.size() != z.size()) {
32  throw BPException(__FILE__, __LINE__, kINDEX);
33  }
34 
35  //
36  // Make line fit in z-x projection
37  //
38  double mzx; // Slope in x-z view
39  double bzx; // Intercept in x-z view
40  for (i=0; i<w.size(); ++i) {
41  w[i] = dx[i]*dx[i] + dz[i]+dz[i];
42  }
43  util::LinFitTS(z, x, w, mzx, bzx, 2);
44 
45  //
46  // Make line fit in z-y projection
47  //
48  double mzy; // Slope in y-z view
49  double bzy; // Intercept in y-z view
50  for (i=0; i<w.size(); ++i) {
51  w[i] = dy[i]*dy[i] + dz[i]+dz[i];
52  }
53  util::LinFitTS(z, y, w, mzy, bzy, 2);
54 
55  z1.SetXYZ( mzx, mzy, 1.0);
56  z1.SetMag(1.0);
57 }
58 
59 //......................................................................
60 
61 void TrackBasis::MakeBasis(const std::vector<double>& x,
62  const std::vector<double>& y,
63  const std::vector<double>& z,
64  const std::vector<double>& dx,
65  const std::vector<double>& dy,
66  const std::vector<double>& dz,
67  double vx, double vy, double vz)
68 {
69  this->SetOrigin(vx,vy,vz);
70 
71  //
72  // The detector based basis vectors
73  //
74  static const TVector3 x0(1,0,0);
75  static const TVector3 y0(0,1,0);
76  static const TVector3 z0(0,0,1);
77 
78  //
79  // Set the z-axis to follow the track
80  //
81  TrackBasis::FindAxis(x, y, z, dx, dy, dz, fZ1);
82 
83  //
84  // Compute some cross products to find the detector axis that is
85  // most complementary
86  //
87  TVector3 ZxX(x0.Cross(fZ1));
88  TVector3 YxZ(y0.Cross(fZ1));
89  TVector3 ZxZ(z0.Cross(fZ1));
90 
91  double ZxXmag = ZxX.Mag2();
92  double YxZmag = YxZ.Mag2();
93  double ZxZmag = ZxZ.Mag2();
94 
95  //
96  // Choose the other two axes to maximize complementarity with the
97  // detector coordinates
98  //
99  if ((YxZmag>=ZxXmag) && (YxZmag>=ZxZmag)) {
100  fX1 = YxZ.Unit();
101  fY1 = fZ1.Cross(fX1);
102  }
103  else if ((ZxXmag>=YxZmag) && (ZxXmag>=ZxZmag)) {
104  fY1 = ZxX.Unit();
105  fX1 = fY1.Cross(fZ1);
106  }
107  else {
108  fX1 = ZxZ.Unit();
109  fY1 = fZ1.Cross(fX1);
110  }
111 
112  //
113  // Use the unit vectors to find the rotation matrix
114  //
115  fU[0][0] = fUT[0][0] = x0.Dot(fX1);
116  fU[0][1] = fUT[1][0] = y0.Dot(fX1);
117  fU[0][2] = fUT[2][0] = z0.Dot(fX1);
118 
119  fU[1][0] = fUT[0][1] = x0.Dot(fY1);
120  fU[1][1] = fUT[1][1] = y0.Dot(fY1);
121  fU[1][2] = fUT[2][1] = z0.Dot(fY1);
122 
123  fU[2][0] = fUT[0][2] = x0.Dot(fZ1);
124  fU[2][1] = fUT[1][2] = y0.Dot(fZ1);
125  fU[2][2] = fUT[2][2] = z0.Dot(fZ1);
126 
127  //
128  // Check that that the basis vectors point in the intended direction
129  // of the track (based on the vertex location) and adjust the rotation
130  // matrix if needed.
131  //
132  this->AdjustBasisDir(x,y,z,dx,dy,dz);
133 }
134 
135 //......................................................................
136 
137 void TrackBasis::AdjustBasisDir(const std::vector<double>& x,
138  const std::vector<double>& y,
139  const std::vector<double>& z,
140  const std::vector<double>& dx,
141  const std::vector<double>& dy,
142  const std::vector<double>& dz)
143 {
144  //
145  // Compute the average Z coordinate in the new basis (weighted by dz.)
146  //
147  unsigned int i;
148  double z1ave = 0.0;
149  double sumdz = 0.0;
150  double tol = 0.1;
151  double w;
152  double x1, y1, z1;
153  double dx1, dy1, dz1;
154 
155  for(i=0; i<x.size(); ++i) {
156  this->DetToTrk(x[i], y[i], z[i], dx[i], dy[i], dz[i],
157  &x1, &y1, &z1, &dx1, &dy1, &dz1);
158  if (dz1<tol) dz1 = tol;
159  w = 1.0/(dz1*dz1);
160  z1ave += w*z1;
161  sumdz += w;
162  }
163  z1ave /= sumdz;
164 
165  //
166  // If z1ave is behind the vertex (<0) then the new Z axis points away
167  // from the intended track direction. So reverse the Z axis in the
168  // rotation matrix, and reverse the Y axis as well to maintain a
169  // right handed coordinate system.
170  //
171  if(z1ave < 0.0) {
172  fY1 = -fY1;
173  fZ1 = -fZ1;
174 
175  fU[1][0] = fUT[0][1] = -fU[1][0];
176  fU[1][1] = fUT[1][1] = -fU[1][1];
177  fU[1][2] = fUT[2][1] = -fU[1][2];
178 
179  fU[2][0] = fUT[0][2] = -fU[2][0];
180  fU[2][1] = fUT[1][2] = -fU[2][1];
181  fU[2][2] = fUT[2][2] = -fU[2][2];
182  }
183 }
184 
185 //......................................................................
186 
187 void TrackBasis::SetOrigin(double x, double y, double z)
188 {
189  fO.SetXYZ(x, y, z);
190 }
191 
192 //......................................................................
193 
194 void TrackBasis::DetToTrk(double x, double y, double z,
195  double dx, double dy, double dz,
196  double* x1, double* y1, double* z1,
197  double* dx1, double* dy1, double* dz1)
198 {
199  x -= fO.X();
200  y -= fO.Y();
201  z -= fO.Z();
202 
203  *x1 = fU[0][0]*x + fU[0][1]*y + fU[0][2]*z;
204  *y1 = fU[1][0]*x + fU[1][1]*y + fU[1][2]*z;
205  *z1 = fU[2][0]*x + fU[2][1]*y + fU[2][2]*z;
206 
207  *dx1 = SQR(fU[0][0]*dx) + SQR(fU[0][1]*dy) + SQR(fU[0][2]*dz);
208  *dy1 = SQR(fU[1][0]*dx) + SQR(fU[1][1]*dy) + SQR(fU[1][2]*dz);
209  *dz1 = SQR(fU[2][0]*dx) + SQR(fU[2][1]*dy) + SQR(fU[2][2]*dz);
210 
211  *dx1 = sqrt(*dx1);
212  *dy1 = sqrt(*dy1);
213  *dz1 = sqrt(*dz1);
214 }
215 
216 //......................................................................
217 
218 void TrackBasis::TrkToDet(double x, double y, double z,
219  double dx, double dy, double dz,
220  double* x1, double* y1, double* z1,
221  double* dx1, double* dy1, double* dz1)
222 {
223  *x1 = fUT[0][0]*x + fUT[0][1]*y + fUT[0][2]*z;
224  *y1 = fUT[1][0]*x + fUT[1][1]*y + fUT[1][2]*z;
225  *z1 = fUT[2][0]*x + fUT[2][1]*y + fUT[2][2]*z;
226 
227  *x1 += fO.X();
228  *y1 += fO.Y();
229  *z1 += fO.Z();
230 
231  *dx1 = SQR(fUT[0][0]*dx) + SQR(fUT[0][1]*dy) + SQR(fUT[0][2]*dz);
232  *dy1 = SQR(fUT[1][0]*dx) + SQR(fUT[1][1]*dy) + SQR(fUT[1][2]*dz);
233  *dz1 = SQR(fUT[2][0]*dx) + SQR(fUT[2][1]*dy) + SQR(fUT[2][2]*dz);
234 
235  *dx1 = sqrt(*dx1);
236  *dy1 = sqrt(*dy1);
237  *dz1 = sqrt(*dz1);
238 }
239 
240 ////////////////////////////////////////////////////////////////////////
241 
TVector3 fO
System origin in detector coordinates.
Definition: TrackBasis.h:137
void MakeBasis(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &z, const std::vector< double > &dx, const std::vector< double > &dy, const std::vector< double > &dz, double vx, double vy, double vz)
Definition: TrackBasis.cxx:61
Float_t y1[n_points_granero]
Definition: compare.C:5
Float_t x1[n_points_granero]
Definition: compare.C:5
T sqrt(T number)
Definition: d0nt_math.hpp:156
double fU[3][3]
Transformation matrix.
Definition: TrackBasis.h:141
TVector3 fX1
x basis vector in detector coordinates
Definition: TrackBasis.h:138
double fUT[3][3]
Transpose of transformation matrix.
Definition: TrackBasis.h:142
double dy[NP][NC]
double dx[NP][NC]
double dz[NP][NC]
void DetToTrk(double x, double y, double z, double dx, double dy, double dz, double *x1, double *y1, double *z1, double *dx1, double *dy1, double *dz1)
Definition: TrackBasis.cxx:194
z
Definition: test.py:28
static void FindAxis(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &z, const std::vector< double > &dx, const std::vector< double > &dy, const std::vector< double > &dz, TVector3 &z1)
Definition: TrackBasis.cxx:15
const unsigned int kINDEX
index problem
Definition: BPException.h:21
#define SQR(x)
Definition: TrackBasis.cxx:13
void AdjustBasisDir(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &z, const std::vector< double > &dx, const std::vector< double > &dy, const std::vector< double > &dz)
Definition: TrackBasis.cxx:137
Construct a track-based (x,y,z) coordinate system.
void TrkToDet(double x, double y, double z, double dx, double dy, double dz, double *x1, double *y1, double *z1, double *dx1, double *dy1, double *dz1)
Definition: TrackBasis.cxx:218
void SetOrigin(double x, double y, double z)
Definition: TrackBasis.cxx:187
TVector3 fZ1
z basis vector in detector coordinates
Definition: TrackBasis.h:140
Float_t w
Definition: plot.C:20
TVector3 fY1
y basis vector in detector coordinates
Definition: TrackBasis.h:139
void LinFitTS(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &wt, double &m, double &b, int wmode)
Best fit line to points using Theil-Sens median method.
Definition: MathUtil.cxx:104