Jpp  master_rocky-43-ge265d140c
the software that should make you happy
JAlignDetector.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <cmath>
5 #include <vector>
6 
7 #include "JMath/JMath.hh"
8 #include "JDetector/JDetector.hh"
11 #include "JDetector/JTripod.hh"
12 
13 #include "JFit/JSimplex.hh"
14 #include "JLang/JException.hh"
15 #include "JSupport/JMeta.hh"
16 #include "JTools/JRange.hh"
17 
18 #include "Jeep/JContainer.hh"
19 #include "Jeep/JParser.hh"
20 #include "Jeep/JMessage.hh"
21 
22 
23 namespace {
24 
25  using namespace JPP;
26 
27  typedef JRange<int> JRange_t;
28 
29 
30  /**
31  * Model for detector alignment.
32  *
33  * The model allows for a global offset of the detector position, a rotation around the z-axis and a tilt of the (x,y) plane.
34  */
35  struct JModel_t :
36  JMath<JModel_t>
37  {
38  /**
39  * Default constructor.
40  */
41  JModel_t() :
42  x (0.0),
43  y (0.0),
44  z (0.0),
45  phi(0.0),
46  tx (0.0),
47  ty (0.0)
48  {}
49 
50 
51  /**
52  * Constructor.
53  *
54  * \param x x
55  * \param y y
56  * \param z z
57  * \param phi phi
58  * \param tx Tx
59  * \param ty Ty
60  */
61  JModel_t(const double x,
62  const double y,
63  const double z,
64  const double phi,
65  const double tx,
66  const double ty) :
67  x (x),
68  y (y),
69  z (z),
70  phi(phi),
71  tx (tx),
72  ty (ty)
73  {}
74 
75 
76  /**
77  * Negate model.
78  *
79  * \return this model
80  */
81  JModel_t& negate()
82  {
83  x = -x;
84  y = -y;
85  z = -z;
86  phi = -phi;
87  tx = -tx;
88  ty = -ty;
89 
90  return *this;
91  }
92 
93 
94  /**
95  * Add model.
96  *
97  * \param model model
98  * \return this model
99  */
100  JModel_t& add(const JModel_t& model)
101  {
102  x += model.x;
103  y += model.y;
104  z += model.z;
105  phi += model.phi;
106  tx += model.tx;
107  ty += model.ty;
108 
109  return *this;
110  }
111 
112 
113  /**
114  * Subtract model.
115  *
116  * \param model model
117  * \return this model
118  */
119  JModel_t& sub(const JModel_t& model)
120  {
121  x -= model.x;
122  y -= model.y;
123  z -= model.z;
124  phi -= model.phi;
125  tx -= model.tx;
126  ty -= model.ty;
127 
128  return *this;
129  }
130 
131 
132  /**
133  * Scale model.
134  *
135  * \param factor multiplication factor
136  * \return this model
137  */
138  JModel_t& mul(const double factor)
139  {
140  x *= factor;
141  y *= factor;
142  z *= factor;
143  phi *= factor;
144  tx *= factor;
145  ty *= factor;
146 
147  return *this;
148  }
149 
150 
151  /**
152  * Scale model.
153  *
154  * \param factor division factor
155  * \return this model
156  */
157  JModel_t& div(const double factor)
158  {
159  x /= factor;
160  y /= factor;
161  z /= factor;
162  phi /= factor;
163  tx /= factor;
164  ty /= factor;
165 
166  return *this;
167  }
168 
169 
170  double x; //!< x-offset
171  double y; //!< y-offset
172  double z; //!< z-offset
173  double phi; //!< rotation around z-axis [rad]
174  double tx; //!< tilt angle [rad]
175  double ty; //!< tilt angle [rad]
176  };
177 
178 
179  /**
180  * Fit function.
181  */
182  class JFit_t :
183  public JModuleRouter
184  {
185  public:
186  /**
187  * Constructor.
188  *
189  * \param detector detector (reference)
190  * \param sigma_m resolution [m]
191  * \param option option to keep strings vertical
192  * \param range range of floors used in fit
193  */
194  JFit_t(const JDetector& detector,
195  const double sigma_m,
196  const bool option,
197  const JRange_t range) :
199  sigma (sigma_m),
200  option(option),
201  range (range)
202  {
203  using namespace JPP;
204 
205  center = getAverage(make_array(detector.begin(), detector.end(), &JModule::getPosition));
206 
207  center.setPosition(JPosition3D(center.getX(), center.getY(), 0.0));
208  }
209 
210 
211  /**
212  * Get chi2.
213  *
214  * \param model model
215  * \param module module
216  * \return chi2
217  */
218  double operator()(const JModel_t& model, const JModule& module) const
219  {
220  if (this->hasModule(module.getID()) && range(module.getFloor())) {
221 
222  const JPosition3D p1 = this->getModule(module.getID()).getPosition();
223  const JPosition3D p2 = this->getPosition(model, module.getPosition());
224 
225  return p2.getDistanceSquared(p1) / (sigma*sigma);
226  }
227 
228  return 0.0;
229  }
230 
231 
232  /**
233  * Get position according model.
234  *
235  * \param model model
236  * \param pos position
237  * \return position
238  */
239  JPosition3D getPosition(const JModel_t& model, const JPosition3D& pos) const
240  {
241  using namespace std;
242  using namespace JPP;
243 
244  JPosition3D p1 = pos;
245 
246  p1 -= center;
247 
248  p1.rotate(JRotation3Z(model.phi));
249 
250  {
251  const double tx = model.tx;
252  const double tz = sqrt((1.0 + tx) * (1.0 - tx));
253 
254  if (!option)
255  p1 = JPosition3D(tz * p1.getX() + tx * p1.getZ(), p1.getY(), tz * p1.getZ() - tx * p1.getX());
256  else
257  p1 = JPosition3D(tz * p1.getX(), p1.getY(), p1.getZ() - tx * p1.getX());
258  }
259  {
260  const double ty = model.ty;
261  const double tz = sqrt((1.0 + ty) * (1.0 - ty));
262 
263  if (!option)
264  p1 = JPosition3D(p1.getX(), tz * p1.getY() + ty * p1.getZ(), tz * p1.getZ() - ty * p1.getY());
265  else
266  p1 = JPosition3D(p1.getX(), tz * p1.getY(), p1.getZ() - ty * p1.getY());
267  }
268 
269  p1 += JPosition3D(model.x, model.y, model.z);
270  p1 += center;
271 
272  return p1;
273  }
274 
275  double sigma;
276  bool option;
277  JRange_t range;
278  JPosition3D center;
279  };
280 }
281 
282 
283 /**
284  * \file
285  *
286  * Auxiliary program to align two detectors.\n
287  * \author mjongen
288  */
289 int main(int argc, char **argv)
290 {
291  using namespace std;
292  using namespace JPP;
293 
294  string detectorFile_a;
295  string detectorFile_b;
296  bool overwriteDetector;
297  string tripodFile;
298  double sigma_m;
299  bool option;
300  JRange_t range;
301  int debug;
302 
303  try {
304 
305  JParser<> zap("Auxiliary program to align two detectors.");
306 
307  zap['a'] = make_field(detectorFile_a, "detector - subject to alignment (option -A)");
308  zap['b'] = make_field(detectorFile_b, "detector - reference for alignment");
309  zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with modified positions.");
310  zap['T'] = make_field(tripodFile, "tripods") = "";
311  zap['s'] = make_field(sigma_m) = 0.2;
312  zap['O'] = make_field(option, "keep strings vertical");
313  zap['r'] = make_field(range, "range of floors used in fit") = JRange_t();
314  zap['d'] = make_field(debug) = 2;
315 
316  zap(argc, argv);
317  }
318  catch(const exception &error) {
319  FATAL(error.what() << endl);
320  }
321 
322 
323  if (overwriteDetector) {
324  if (tripodFile == "") {
325  FATAL("No tripod file.");
326  }
327  }
328 
329 
330  JDetector detector_a;
331  JDetector detector_b;
332 
333  try {
334  load(detectorFile_a, detector_a);
335  }
336  catch(const JException& error) {
337  FATAL(error);
338  }
339 
340  try {
341  load(detectorFile_b, detector_b);
342  }
343  catch(const JException& error) {
344  FATAL(error);
345  }
346 
347 
348  const JFit_t fit(detector_b, sigma_m, option, range);
349 
351 
352  for (JDetector::const_iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
353  if (fit.hasModule(module->getID())) {
354  data.push_back(*module);
355  }
356  }
357 
358 
359  JSimplex<JModel_t> simplex;
360 
363 
364  simplex.value = JModel_t(0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
365 
366  simplex.step.resize(6);
367 
368  simplex.step[0] = JModel_t(0.01, 0.00, 0.00, 0.0, 0.0, 0.0);
369  simplex.step[1] = JModel_t(0.00, 0.01, 0.00, 0.0, 0.0, 0.0);
370  simplex.step[2] = JModel_t(0.00, 0.00, 0.01, 0.0, 0.0, 0.0);
371  simplex.step[3] = JModel_t(0.00, 0.00, 0.00, 5.0e-4, 0.0, 0.0);
372  simplex.step[4] = JModel_t(0.00, 0.00, 0.00, 0.0, 1.0e-4, 0.0);
373  simplex.step[5] = JModel_t(0.00, 0.00, 0.00, 0.0, 0.0, 1.0e-4);
374 
375  const double chi2 = simplex(fit, data.begin(), data.end());
376 
377  cout << "Number of iterations " << simplex.numberOfIterations << endl;
378  cout << "chi2/NDF " << FIXED(7,3) << chi2 << '/' << (detector_a.size() - simplex.step.size()) << endl;
379 
380  cout << "model:" << endl;
381 
382  cout << "x " << FIXED(7,3) << simplex.value.x << endl;
383  cout << "y " << FIXED(7,3) << simplex.value.y << endl;
384  cout << "z " << FIXED(7,3) << simplex.value.z << endl;
385  cout << "phi " << FIXED(9,5) << simplex.value.phi << endl;
386  cout << "Tx " << FIXED(9,5) << simplex.value.tx << endl;
387  cout << "Ty " << FIXED(9,5) << simplex.value.ty << endl;
388 
389 
390  if (overwriteDetector) {
391 
392  NOTICE("Store alignment data on files " << detectorFile_a << " and " << tripodFile << endl);
393 
394  detector_a.comment.add(JMeta(argc,argv));
395 
396  for (JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
397 
398  const JPosition3D p1 = module->getPosition();
399  const JPosition3D p2 = fit.getPosition(simplex.value, p1);
400 
401  module->add(p2 - p1);
402  }
403 
404  try {
405  store(detectorFile_a, detector_a);
406  }
407  catch(const JException& error) {
408  FATAL(error);
409  }
410 
411  JContainer< vector<JTripod> > tripods;
412 
413  tripods.load(tripodFile.c_str());
414 
415  tripods.comment.add(JMeta(argc,argv));
416 
417  for (vector<JTripod>::iterator tripod = tripods.begin(); tripod != tripods.end(); ++tripod) {
418 
419  const JPosition3D p1 = tripod->getUTMPosition() - detector_a.getUTMPosition();
420  const JPosition3D p2 = fit.getPosition(simplex.value, p1);
421 
422  tripod->add(p2 - p1);
423  }
424 
425  tripods.store(tripodFile.c_str());
426  }
427 }
int main(int argc, char **argv)
const JModule & getModule(const JDetector &detector, const JModuleLocation &location)
find module with a given string and floor number
Container I/O.
Data structure for detector geometry and calibration.
TPaveText * p1
Exceptions.
Base class for data structures with artithmetic capabilities.
General purpose messaging.
#define NOTICE(A)
Definition: JMessage.hh:64
#define FATAL(A)
Definition: JMessage.hh:67
int debug
debug level
Definition: JSirene.cc:69
ROOT I/O of application specific meta data.
Direct access to module in detector data structure.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2142
Auxiliary class to define a range between two values.
Data structure for tripod.
Detector data structure.
Definition: JDetector.hh:96
int getFloor() const
Get floor number.
Definition: JLocation.hh:146
Router for direct addressing of module data in detector data structure.
Data structure for a composite optical module.
Definition: JModule.hh:75
JModule & add(const JVector3D &pos)
Add position.
Definition: JModule.hh:419
Simple fit method based on Powell's algorithm, see reference: Numerical Recipes in C++,...
Definition: JSimplex.hh:44
JModel_t value
Definition: JSimplex.hh:240
std::vector< JModel_t > step
Definition: JSimplex.hh:241
int numberOfIterations
Definition: JSimplex.hh:242
Data structure for position in three dimensions.
Definition: JPosition3D.hh:38
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
Rotation around Z-axis.
Definition: JRotation3D.hh:87
double getDistanceSquared(const JVector3D &pos) const
Get squared of distance to point.
Definition: JVector3D.hh:258
General exception.
Definition: JException.hh:24
int getID() const
Get identifier.
Definition: JObjectID.hh:50
Utility class to parse command line options.
Definition: JParser.hh:1698
const JUTMPosition & getUTMPosition() const
Get UTM position.
Definition: JUTMPosition.hh:84
const double sigma[]
Definition: JQuadrature.cc:74
JPosition3D getPosition(const Vec &pos)
Get position.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
const array_type< JValue_t > & make_array(const JValue_t(&array)[N])
Method to create array of values.
Definition: JVectorize.hh:54
std::iterator_traits< T >::value_type getAverage(T __begin, T __end)
Get average.
Definition: JMath.hh:494
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JSTDTypes.hh:14
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:448
Type definition of range.
Definition: JHead.hh:43
Detector file.
Definition: JHead.hh:227
JComment & add(const std::string &comment)
Add comment.
Definition: JComment.hh:100
Auxiliary wrapper for I/O of container with optional comment (see JComment).
Definition: JContainer.hh:42
JComment comment
Definition: JContainer.hh:90
void store(const char *file_name) const
Store to output file.
void load(const char *file_name)
Load from input file.
Auxiliary base class for aritmetic operations of derived class types.
Definition: JMath.hh:347
Fit model.
Definition: JMath/JModel.hh:32
JModel_t & sub(const JModel_t &model)
Subtract model.
JModel_t & add(const JModel_t &model)
Add model.
JModel_t & div(const double factor)
Scale model.
JModel_t & negate()
Negate model.
Definition: JMath/JModel.hh:94
JModel_t & mul(const double factor)
Scale model.
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:72