CoordinateTransformation.cxx
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////////
2 // \file CoordinateTransformation.cxx
3 // \brief
4 // \version $Id: CoordinateTransformation.cxx,v 1.3 2012-11-09 17:22:59 jpaley Exp $
5 // \author denis@fnal.gov
6 /////////////////////////////////////////////////////////////////////////
7 
9 #include "NovaDAQConventions/DAQConventions.h"
10 
11 #include <boost/numeric/ublas/io.hpp>
12 
14 
15 #include "boost/format.hpp"
16 
17 #include <cstdio>
18 #include <iostream>
19 #include <iomanip>
20 
21 namespace geo{
22 
23  //---------------------------------------------------------------------
25  _detectorID(0),
26  _beamType(0),
27  _rotationMatrix(3,3),
28  _rotationMatrixInverse(3,3),
29  _simpleTranslation(3)
30  {
31  }
32 
33  //------------------------------------------------------------------------------
35  _detectorID(detector_id),
36  _beamType(beam_type),
37  _rotationMatrix(3,3),
40  {
42  }
43 
44 
45  //------------------------------------------------------------------------------
46  bool CoordinateTransformation::setDetectorAndBeam(uint32_t detector_id, uint32_t beam_type){
47 
48  /// Trying to set the detector parameters
49  /// If unsuccessfull, set the detector and beam ID to zero
50  if(!setCoordinateParameters(detector_id, beam_type)){
51  _detectorID = 0;
52  _beamType = 0;
53  return false;
54  }// end of setting the detector parameters
55 
56 
57  /// Have to reset the detector ID and beam type
58  _detectorID = detector_id;
60 
61  return true;
62  }
63 
64 
65  //------------------------------------------------------------------------------
67  LOG_DEBUG("CoordinateTransformation")
68  << "~CoordinateTransformation\n";
69  }
70 
71 
72  //------------------------------------------------------------------------------
74 
75  double ca = cos(_alpha);
76  double sa = sin(_alpha);
77  double cb = cos(_beta);
78  double sb = sin(_beta);
79 
80  _rotationMatrix(0, 0) = cb;
81  _rotationMatrix(0, 1) = -sb*sa;
82  _rotationMatrix(0, 2) = sb*ca;
83 
84  _rotationMatrix(1, 0) = 0.0;
85  _rotationMatrix(1, 1) = ca;
86  _rotationMatrix(1, 2) = sa;
87 
88  _rotationMatrix(2, 0) = -sb;
89  _rotationMatrix(2, 1) = -cb*sa;
90  _rotationMatrix(2, 2) = cb*ca;
91 
92 
93  _rotationMatrixInverse(0, 0) = cb;
94  _rotationMatrixInverse(0, 1) = 0.0;
95  _rotationMatrixInverse(0, 2) = -sb;
96 
97  _rotationMatrixInverse(1, 0) = -sa*sb;
98  _rotationMatrixInverse(1, 1) = ca;
99  _rotationMatrixInverse(1, 2) = -sa*cb;
100 
101  _rotationMatrixInverse(2, 0) = ca*sb;
102  _rotationMatrixInverse(2, 1) = sa;
103  _rotationMatrixInverse(2, 2) = ca*cb;
104 
105  return true;
106  }
107 
108  //------------------------------------------------------------------------------
110  (double* input_beam_coords,
111  double* output_detector_coords) const {
112 
113 
114  /// Check if the beam type and detector id are set
115  if(!_beamType || !_detectorID) return false;
116 
117  /// Calculate the beam coordinates
118  /// x`` = M^{-1} * (x-x0)
119  for(int i = 0; i < 3; i++){
120  double sum = 0.0;
121  for(int j = 0; j < 3; j++)
122  sum += _rotationMatrixInverse(i,j)*(input_beam_coords[j] - _simpleTranslation[j]);
123 
124  output_detector_coords[i] = sum;
125  }// end of loop over coordinates
126 
127  return true;
128  }// end of transformBeamToDetectorCoordinates
129 
130  //------------------------------------------------------------------------------
132  (double* input_detector_coords,
133  double* output_beam_coords) const {
134 
135 
136  /// Check if the beam type and detector id are set
137  if(!_beamType || !_detectorID) return false;
138 
139 
140  /// Calculate the detector coordinates
141  /// x = x0 + M * x''
142  for(int i = 0; i < 3; i++) {
143  double sum = 0.0;
144  for(int j = 0; j < 3; j++)
145  sum += _rotationMatrix(i,j)*input_detector_coords[j];
146 
147  output_beam_coords[i] = _simpleTranslation[i] + sum;
148  }// end of loop over coordinates
149 
150  return true;
151  }// end of transformDetectorToBeamCoordinates
152 
153 
154  //------------------------------------------------------------------------------
156 
157  std::cout << "Coordinate system transformations for "
158  << "Beam Type = " << _beamType << " "
159  << "Detector = " << _detectorID << ":\n";
160 
161  std::cout << "Translation:\n";
162  printf(" |%0.10E|\n", _simpleTranslation[0]);
163  printf("S = |%0.10E|\n", _simpleTranslation[1]);
164  printf(" |%0.10E|\n", _simpleTranslation[2]);
165 
166 
167 
168  std::cout << "Beam->Detector\n";
169  printf("|x_det| | %0.10E %0.10E %0.10E | | x_beam - %0.10E |\n", _rotationMatrixInverse(0,0), _rotationMatrixInverse(0,1), _rotationMatrixInverse(0,2) , _simpleTranslation[0]);
170  printf("|y_det| = | %0.10E %0.10E %0.10E | * | y_beam - %0.10E |\n", _rotationMatrixInverse(1,0), _rotationMatrixInverse(1,1), _rotationMatrixInverse(1,2) , _simpleTranslation[1]);
171  printf("|z_det| | %0.10E %0.10E %0.10E | | z_beam - %0.10E |\n", _rotationMatrixInverse(2,0), _rotationMatrixInverse(2,1), _rotationMatrixInverse(2,2) , _simpleTranslation[2]);
172 
173 
174  std::cout << "\n";
175  std::cout << "Detector->Beam\n";
176  printf("|x_beam| | %0.10E | | %0.10E %0.10E %0.10E | | x_det |\n", _simpleTranslation[0], _rotationMatrix(0,0), _rotationMatrix(0,1), _rotationMatrix(0,2));
177  printf("|y_beam| = | %0.10E | + | %0.10E %0.10E %0.10E | * | y_det |\n", _simpleTranslation[1], _rotationMatrix(1,0), _rotationMatrix(1,1), _rotationMatrix(1,2));
178  printf("|z_beam| | %0.10E | | %0.10E %0.10E %0.10E | | z_det |\n", _simpleTranslation[2], _rotationMatrix(2,0), _rotationMatrix(2,1), _rotationMatrix(2,2));
179 
180  /// Now calculating the angle at which point S sees the front surface of the detector
181 
182  /// First figuring out the detector UZ in beam coordinates
183  double p0_det[3] = {0,0,0};
184  double p0_beam[3];
185  transformDetectorToBeamCoordinates(p0_det, p0_beam);
186 
187  double p1_det[3] = {0,0,1};
188  double p1_beam[3];
189  transformDetectorToBeamCoordinates(p1_det, p1_beam);
190 
191  /// Unit vector of detector's Uz in beam coordinates
192  double dr_detector_uz[3] = {
193  p1_beam[0] - p0_beam[0],
194  p1_beam[1] - p0_beam[1],
195  p1_beam[2] - p0_beam[2]
196  };
197 
198 
199  /// Distance from point S to detector zero
200  double dr_from_pointS_to_detector_zero[3] = {
203  _simpleTranslation[2] - pointS_Z
204  };
205 
206  /// Now need to normalize it
207  {
208  double drmod_from_pointS_to_detector_zero = sqrt(
209  dr_from_pointS_to_detector_zero[0]*dr_from_pointS_to_detector_zero[0]
210  + dr_from_pointS_to_detector_zero[1]*dr_from_pointS_to_detector_zero[1]
211  + dr_from_pointS_to_detector_zero[2]*dr_from_pointS_to_detector_zero[2]
212  );
213  for(int i=0; i<3; i++)
214  dr_from_pointS_to_detector_zero[i] /= drmod_from_pointS_to_detector_zero;
215  }// end of normalization
216 
217  /// cosine of the angle at which point-S sees the front surface of the detector
218  double cos_detector_front_surface = dr_detector_uz[0] * dr_from_pointS_to_detector_zero[0]
219  + dr_detector_uz[1] * dr_from_pointS_to_detector_zero[1]
220  + dr_detector_uz[2] * dr_from_pointS_to_detector_zero[2];
221 
222  LOG_DEBUG("CoordinateTransformation")
223  << "Cosine of the angle at which point S sees the front surface of the detector: "
224  << cos_detector_front_surface
225  << "\n";
226 
227  return true;
228  }
229 
230 
231  //------------------------------------------------------------------------------
233  uint32_t beam_type) {
234 
235 
236  /// See if we support this beam type
237  if(beam_type != BEAMTYPE_NUMI) {
238  LOG_DEBUG("CoordinateTransformation")
239  << "From CoordinateTransformation::setCoordinateParameters. "
240  << "Beam type " << beam_type << " is not supported.\n";
241  return false;
242  }// end of checking the beam configuration
243 
244 
245  /// See if we support this detector type
246  if (detector_id == novadaq::cnv::kNEARDET) {
247  /// From NOvA-doc-5485
248  _simpleTranslation[0] = 11.414 * CLHEP::m;
249  _simpleTranslation[1] = -3.456 * CLHEP::m;
250  _simpleTranslation[2] = 994.665 * CLHEP::m;
251 
254  }
255  else if (detector_id == novadaq::cnv::kFARDET) {
256  /// From NOvA-doc-5485
257  _simpleTranslation[0] = 11.037 * CLHEP::km;
258  _simpleTranslation[1] = -4.1626 * CLHEP::km;
259  _simpleTranslation[2] = 810.42 * CLHEP::km;
260 
263  }
264  else if (detector_id == novadaq::cnv::kNDOS) {
265  /// From NOvA-doc-5485
266  _simpleTranslation[0] = -0.29 * CLHEP::m;
267  _simpleTranslation[1] = 92.21 * CLHEP::m;
268  _simpleTranslation[2] = 841.76* CLHEP::m;
269 
271  _beta = 0.0;
272  }
273 
274  else {
275  LOG_DEBUG("CoordinateTransformation")
276  << "From geo::CoordinateTransformation::setCoordinateParameters. "
277  << "Detector ID " << detector_id << " is not supported.\n";
278  return false;
279  }// end of checking the detector configuration
280 
281  return init();
282  }
283 
284 }
285 ////////////////////////////////////////////////////////////////////////////
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
boost::numeric::ublas::matrix< double > _rotationMatrixInverse
Rotational matrix inverse.
beam_type
Definition: Constants.h:26
bool setDetectorAndBeam(uint32_t detector_id, uint32_t beam_type)
Try to set the detector and beam.
double _alpha
Angle between NuMI beamline and the earth&#39;s surface.
bool transformBeamToDetectorCoordinates(double *input_beam_coords, double *output_detector_coords) const
T sqrt(T number)
Definition: d0nt_math.hpp:156
bool init()
Initializing matrices.
bool transformDetectorToBeamCoordinates(double *input_detector_coords, double *output_beam_coords) const
bool setCoordinateParameters(uint32_t detector_id, uint32_t beam_type)
Far Detector at Ash River, MN.
Prototype Near Detector on the surface at FNAL.
printf("%d Experimental points found\n", nlines)
Near Detector in the NuMI cavern.
boost::numeric::ublas::matrix< double > _rotationMatrix
Rotational matrix.
const double j
Definition: BetheBloch.cxx:29
uint32_t _beamType
Current Beam Type.
static constexpr double m
OStream cout
Definition: OStream.cxx:6
T sin(T number)
Definition: d0nt_math.hpp:132
double _beta
Angle between NuMi beamline and the detector looking from top.
T cos(T number)
Definition: d0nt_math.hpp:78
boost::numeric::ublas::vector< double > _simpleTranslation
coordinates of the detector zero in the beam coordinates
Double_t sum
Definition: plot.C:31
Helper for AttenCurve.
Definition: Path.h:10
static constexpr double km
const double NOvAneardetNuMiAngle
Angle between NuMi beamline and the NOvA near detectors.
uint32_t _detectorID
Current Detector ID.
const double pointS_Z
const double NuMiBeamLineAngle