Jpp  master_rocky
the software that should make you happy
JCompareDetector.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 
5 #include "TROOT.h"
6 #include "TFile.h"
7 #include "TH1D.h"
8 #include "TH2D.h"
9 
11 #include "JMath/JConstants.hh"
12 #include "JTools/JRange.hh"
13 #include "JTools/JQuantile.hh"
14 #include "JROOT/JManager.hh"
15 #include "JROOT/JRootToolkit.hh"
16 #include "JDetector/JDetector.hh"
20 
21 #include "Jeep/JProperties.hh"
22 #include "Jeep/JParser.hh"
23 #include "Jeep/JMessage.hh"
24 
25 
26 namespace {
27 
28  using namespace JPP;
29 
30  /**
31  * Compare time calibration of PMTs between two optical modules.
32  *
33  * \param module_a module
34  * \param module_b module
35  * \return quantile
36  */
37  inline JQuantile getQuantile(const JModule& module_a, const JModule& module_b)
38  {
39  JQuantile Q;
40 
41  JModule::const_iterator pmt_a = module_a.begin();
42  JModule::const_iterator pmt_b = module_b.begin();
43 
44  for ( ; pmt_a != module_a.end() && pmt_b != module_b.end(); ++pmt_a, ++pmt_b) {
45  Q.put(pmt_a->getT0() - pmt_b->getT0());
46  }
47 
48  return Q;
49  }
50 
51 
52  /**
53  * Get ROOT histogram abscissa value.
54  *
55  * \param buffer buffer
56  * \param value value
57  * \return value
58  */
59  inline Double_t getBin(const std::set<int>& buffer, const int value)
60  {
61  return std::distance(buffer.begin(), buffer.find(value));
62  }
63 }
64 
65 
66 /**
67  * \file
68  *
69  * Auxiliary program to find differences between two detector files.\n
70  * A ROOT output file with histograms is optionally produced.
71  * \author mjongen
72  */
73 int main(int argc, char **argv)
74 {
75  using namespace std;
76  using namespace JPP;
77 
78  struct {
79  double tcal = 0.001; // [ns]
80  double pcal = 0.001; // [m]
81  double rcal = 0.001; // [rad]
82  double acal = 0.001; // [ns]
83  double ccal = 0.001; // [deg]
84  int scal = 0xFFFFFFFF; // [status]
85  } precision;
86 
87  string detectorFile_a;
88  string detectorFile_b;
89  string outputFile;
90  int debug;
91 
92  try {
93 
94  JProperties properties;
95 
96  properties[TCAL] = precision.tcal;
97  properties[PCAL] = precision.pcal;
98  properties[RCAL] = precision.rcal;
99  properties[ACAL] = precision.acal;
100  properties[CCAL] = precision.ccal;
101  properties[SCAL] = precision.scal;
102 
103  JParser<> zap("Auxiliary program to find differences between two detector files.");
104 
105  zap['a'] = make_field(detectorFile_a);
106  zap['b'] = make_field(detectorFile_b);
107  zap['o'] = make_field(outputFile) = "";
108  zap['p'] = make_field(properties) = JPARSER::initialised();
109  zap['d'] = make_field(debug) = 3;
110 
111  zap(argc, argv);
112  }
113  catch(const exception &error) {
114  FATAL(error.what() << endl);
115  }
116 
117 
118  JDetector detector_a;
119  JDetector detector_b;
120 
121  try {
122  load(detectorFile_a, detector_a);
123  }
124  catch(const JException& error) {
125  FATAL(error);
126  }
127 
128  try {
129  load(detectorFile_b, detector_b);
130  }
131  catch(const JException& error) {
132  FATAL(error);
133  }
134 
135  size_t numberOfPMTs = 0;
136 
137  bool is_equal = true;
138 
139  const JModuleRouter module_router_a(detector_a);
140  const JModuleRouter module_router_b(detector_b);
141 
142  setFormat<JPosition3D> (JFormat_t(15, 9, std::ios::fixed | std::ios::showpos));
143 
144  // compare detector identifiers
145 
146  if (detector_a.getID() != detector_b.getID()) {
147 
148  DEBUG("* Different detector identifiers "
149  << setw(5) << detector_a.getID() << " (A) and " << endl
150  << setw(5) << detector_b.getID() << " (B)." << endl
151  << endl);
152 
153  is_equal = false;
154  }
155 
156  // UTM positions
157 
158  if (getDistance(detector_a.getPosition(), detector_b.getPosition()) > precision.pcal) {
159 
160  DEBUG(" * different UTM position: "
161  << detector_a.getPosition() << " (A), "
162  << detector_b.getPosition() << " (B)"
163  << ", B - A " << JPosition3D(detector_b.getPosition() - detector_a.getPosition()) << endl);
164 
165  is_equal = false;
166  }
167 
168 
169  for (JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
170  if (module->size() > numberOfPMTs) {
171  numberOfPMTs = module->size();
172  }
173  }
174 
175  // check whether the same modules are present in the file
176 
177  for (JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
178 
179  if (!module_router_b.hasModule(module->getID())) {
180 
181  DEBUG("* Module " << setw(10) << module->getID() << " is in A " << getLabel(*module) << " but not in B" << endl);
182 
183  is_equal = false;
184  }
185  }
186 
187  for (JDetector::iterator module = detector_b.begin(); module != detector_b.end(); ++module) {
188 
189  if (!module_router_a.hasModule(module->getID())) {
190 
191  DEBUG("* Module " << setw(10) << module->getID() << " is in B " << getLabel(*module) << " but not in A" << endl);
192 
193  is_equal = false;
194  }
195  }
196  DEBUG(endl);
197 
198 
199  // compare the modules that the files have in common
200 
201  DEBUG("Comparing module by module." << endl);
202 
203  for (JDetector::iterator module_a = detector_a.begin(); module_a != detector_a.end(); ++module_a) {
204 
205  if (!module_router_b.hasModule(module_a->getID())) {
206 
207  continue;
208 
209  is_equal = false;
210  }
211 
212  const JModule* module_b = &module_router_b.getModule(module_a->getID());
213 
214  DEBUG(" Module " << setw(10) << module_a->getID());
215 
216  // string and floor numbers
217 
218  if (module_a->getLocation() == module_b->getLocation()) {
219 
220  DEBUG(" " << getLabel(*module_a) << endl);
221 
222  } else {
223 
224  DEBUG(endl);
225  DEBUG(" * different location: "
226  << getLabel(*module_a) << " (A), "
227  << getLabel(*module_b) << " (B)" << endl);
228 
229  is_equal = false;
230  }
231 
232  // time offset
233 
234  if (fabs(module_a->getT0() - module_b->getT0()) > precision.acal) {
235 
236  DEBUG(" * different T0: "
237  << FIXED(12,3) << module_a->getT0() << " (A), "
238  << FIXED(12,3) << module_b->getT0() << " (B) "
239  << ", B - A " << module_b->getT0() - module_a->getT0() << endl);
240 
241  is_equal = false;
242  }
243 
244  // quaternion calibration
245 
246  if (module_a->getQuaternion() != JQuaternion3D(0.0, 0.0, 0.0, 0.0) &&
247  module_b->getQuaternion() != JQuaternion3D(0.0, 0.0, 0.0, 0.0) &&
248  getAngle(module_a->getQuaternion(), module_b->getQuaternion()) > precision.ccal) {
249 
250  DEBUG(" * different quaternion calibration: "
251  << module_a->getQuaternion() << " (A), "
252  << module_b->getQuaternion() << " (B) "
253  << ", B - A " << getAngle(module_b->getQuaternion(), module_a->getQuaternion()) << endl);
254 
255  is_equal = false;
256  }
257 
258  // average DOM positions
259 
260  if (getDistance(module_a->getPosition(), module_b->getPosition()) > precision.pcal) {
261 
262  DEBUG(" * different position: "
263  << module_a->getPosition() << " (A), "
264  << module_b->getPosition() << " (B)"
265  << ", B - A " << JPosition3D(module_b->getPosition() - module_a->getPosition()) << endl);
266 
267  is_equal = false;
268  }
269 
270  // number of PMTs
271 
272  if (module_a->size() != module_b->size()) {
273 
274  DEBUG(" * different number of PMTs: "
275  << module_a->size() << " (A), "
276  << module_b->size() << " (B)" << endl);
277 
278  is_equal = false;
279  }
280 
281  // average t0
282 
283  if (!module_a->empty() &&
284  !module_b->empty()) {
285 
286  const JQuantile q = getQuantile(*module_a, *module_b);
287 
288  if (fabs(q.getMean()) > precision.tcal) {
289 
290  DEBUG(" * different average t0: "
291  << ", B - A " << q.getMean()
292  << endl);
293 
294  is_equal = false;
295  }
296  }
297 
298  // status
299 
300  if (module_a->getStatus(precision.scal) != module_b->getStatus(precision.scal)) {
301 
302  DEBUG(" * different status module " << module_a->getID() << ": "
303  << module_a->getStatus() << " (A), "
304  << module_b->getStatus() << " (B)"
305  << endl);
306 
307  is_equal = false;
308  }
309 
310  // comparing by PMT
311 
312  // identifier
313 
314  for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
315 
316  const JPMT& pmt_a = module_a->getPMT(pmt);
317  const JPMT& pmt_b = module_b->getPMT(pmt);
318 
319  if (pmt_a.getID() != pmt_b.getID()) {
320 
321  DEBUG(" * different identifier PMT " << setw(2) << pmt << ": "
322  << setw(8) << pmt_a.getID() << " (A" << FILL(2,'0') << pmt << "), " << FILL()
323  << setw(8) << pmt_b.getID() << " (B" << FILL(2,'0') << pmt << ")" << FILL()
324  << ", B - A " << pmt_b.getID() - pmt_a.getID()
325  << endl);
326 
327  is_equal = false;
328  }
329  }
330 
331  // t0
332 
333  for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
334 
335  const JPMT& pmt_a = module_a->getPMT(pmt);
336  const JPMT& pmt_b = module_b->getPMT(pmt);
337 
338  if (fabs(pmt_a.getT0() - pmt_b.getT0()) > precision.tcal) {
339 
340  DEBUG(" * different T0 PMT " << setw(2) << pmt << ": "
341  << FIXED(12,3) << pmt_a.getT0() << " (A" << FILL(2,'0') << pmt << "), " << FILL()
342  << FIXED(12,3) << pmt_b.getT0() << " (B" << FILL(2,'0') << pmt << ")" << FILL()
343  << ", B - A " << pmt_b.getT0() - pmt_a.getT0()
344  << endl);
345 
346  is_equal = false;
347  }
348  }
349 
350  // positions
351 
352  for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
353 
354  const JPMT& pmt_a = module_a->getPMT(pmt);
355  const JPMT& pmt_b = module_b->getPMT(pmt);
356 
357  // PMT positions
358 
359  if (pmt_a.getPosition().getDistance(pmt_b.getPosition()) > precision.pcal) {
360 
361  DEBUG(" * different PMT position: "
362  << pmt_a.getPosition() << " (A" << FILL(2,'0') << pmt << "), " << FILL()
363  << pmt_b.getPosition() << " (B" << FILL(2,'0') << pmt << ")" << FILL()
364  << ", B - A " << JPosition3D(pmt_b.getPosition() - pmt_a.getPosition()) << endl);
365 
366  is_equal = false;
367  }
368  }
369 
370  // orientations
371 
372  for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
373 
374  const JPMT& pmt_a = module_a->getPMT(pmt);
375  const JPMT& pmt_b = module_b->getPMT(pmt);
376 
377  // PMT positions
378 
379  const double dot = pmt_a.getDirection().getDot(pmt_b.getDirection());
380 
381  if ((1.0 - dot) > precision.rcal) {
382 
383  DEBUG(" * different PMT direction: "
384  << pmt_a.getDirection() << " (A" << FILL(2,'0') << pmt << "), " << FILL()
385  << pmt_b.getDirection() << " (B" << FILL(2,'0') << pmt << ")" << FILL()
386  << ", B - A " << JPosition3D(JPosition3D(pmt_b.getDirection()) - JPosition3D(pmt_a.getDirection())) << ' ' << dot << endl);
387 
388  is_equal = false;
389  }
390  }
391 
392  // status
393 
394  for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
395 
396  const JPMT& pmt_a = module_a->getPMT(pmt);
397  const JPMT& pmt_b = module_b->getPMT(pmt);
398 
399  if (pmt_a.getStatus(precision.scal) != pmt_b.getStatus(precision.scal)) {
400 
401  DEBUG(" * different status PMT " << setw(2) << pmt << ": "
402  << pmt_a.getStatus() << " (A" << FILL(2,'0') << pmt << "), " << FILL()
403  << pmt_b.getStatus() << " (B" << FILL(2,'0') << pmt << ")" << FILL()
404  << endl);
405 
406  is_equal = false;
407  }
408  }
409  }
410  DEBUG(endl);
411 
412 
413  if (outputFile != "") {
414 
415  set<int> string;
416  set<int> floor;
417 
418  for (JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
419  string.insert(module->getString());
420  floor .insert(module->getFloor ());
421  }
422 
423  for (JDetector::iterator module = detector_b.begin(); module != detector_b.end(); ++module) {
424  string.insert(module->getString());
425  floor .insert(module->getFloor ());
426  }
427 
428  JManager<int, TH1D> H1(new TH1D("M[%]", NULL, numberOfPMTs, -0.5, numberOfPMTs - 0.5));
429 
430  TH2D M2("M2", NULL,
431  string.size(), -0.5, string.size() - 0.5,
432  floor .size(), -0.5, floor .size() - 0.5);
433 
434  for (set<int>::const_iterator i = string.begin(); i != string.end(); ++i) {
435  M2.GetXaxis()->SetBinLabel(distance(string.begin(), i) + 1, MAKE_CSTRING(*i));
436  }
437 
438  for (set<int>::const_iterator i = floor.begin(); i != floor.end(); ++i) {
439  M2.GetYaxis()->SetBinLabel(distance(floor.begin(), i) + 1, MAKE_CSTRING(*i));
440  }
441 
442  TH2D* X2 = (TH2D*) M2.Clone("X2" );
443  TH2D* Y2 = (TH2D*) M2.Clone("Y2" );
444  TH2D* Z2 = (TH2D*) M2.Clone("Z2" );
445  TH2D* T2 = (TH2D*) M2.Clone("T2" );
446  TH2D* RMS = (TH2D*) M2.Clone("RMS");
447  TH2D* R2 = (TH2D*) M2.Clone("R2" );
448  TH2D* QA = (TH2D*) M2.Clone("QA" );
449  TH2D* QB = (TH2D*) M2.Clone("QB" );
450  TH2D* QC = (TH2D*) M2.Clone("QC" );
451  TH2D* QD = (TH2D*) M2.Clone("QD" );
452  TH2D* Q2 = (TH2D*) M2.Clone("Q2" );
453 
454  for (JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
455  if (!module_router_b.hasModule(module->getID()) ) {
456  M2.Fill(module->getString(), module->getFloor(), -1.0);
457  }
458  }
459 
460  for (JDetector::iterator module = detector_b.begin(); module != detector_b.end(); ++module) {
461  if (!module_router_a.hasModule(module->getID()) ) {
462  M2.Fill(module->getString(), module->getFloor(), +1.0);
463  }
464  }
465 
466 
467  for (JDetector::iterator module_a = detector_a.begin(); module_a != detector_a.end(); ++module_a) {
468 
469  if (!module_router_b.hasModule(module_a->getID())) {
470  continue;
471  }
472 
473  const JModule* module_b = &module_router_b.getModule(module_a->getID());
474 
475  for (size_t i = 0; i != numberOfPMTs; ++i) {
476 
477  if (module_a->getFloor() != 0) {
478  H1[module_a->getID()]->SetBinContent(i + 1, module_a->getPMT(i).getT0() - module_b->getPMT(i).getT0());
479  }
480  }
481 
482  X2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), module_a->getX() - module_b->getX() + numeric_limits<double>::min());
483  Y2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), module_a->getY() - module_b->getY() + numeric_limits<double>::min());
484  Z2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), module_a->getZ() - module_b->getZ() + numeric_limits<double>::min());
485 
486  if (module_a->getFloor() != 0 &&
487  module_b->getFloor() != 0) {
488 
489  const JQuaternion3D Q = getRotation(*module_a, *module_b);
490  const JQuantile q = getQuantile(*module_a, *module_b);
491 
493 
494  const double phi = (JVector3Z_t.getDot(q1.twist) >= 0.0 ? +1.0 : -1.0) * q1.twist.getAngle();
495 
496  R2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), phi);
497  QA ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), Q.getA());
498  QB ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), Q.getB());
499  QC ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), Q.getC());
500  QD ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), Q.getD());
501  Q2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), Q.getAngle());
502  T2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), q.getMean());
503  RMS->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), q.getSTDev());
504  }
505  }
506 
507 
508  TFile out(outputFile.c_str(), "recreate");
509 
510  for (TH2D* h2 : { &M2, X2, Y2, Z2, T2, RMS, R2, QA, QB, QC, QD, Q2 }) {
511  out << *h2;
512  }
513 
514  out << H1;
515 
516  out.Write();
517  out.Close();
518  }
519 
520  ASSERT(is_equal);
521 
522  return 0;
523 }
int main(int argc, char **argv)
string outputFile
Data structure for detector geometry and calibration.
Dynamic ROOT object management.
Mathematical constants.
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
#define ASSERT(A,...)
Assert macro.
Definition: JMessage.hh:90
#define FATAL(A)
Definition: JMessage.hh:67
int debug
debug level
Definition: JSirene.cc:69
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
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:72
Utility class to parse parameter values.
Auxiliary class to define a range between two values.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
double getT0() const
Get time offset.
Detector data structure.
Definition: JDetector.hh:96
int getFloor() const
Get floor number.
Definition: JLocation.hh:146
int getString() const
Get string number.
Definition: JLocation.hh:135
const JLocation & getLocation() const
Get location.
Definition: JLocation.hh:70
Router for direct addressing of module data in detector data structure.
const JModule & getModule(const JObjectID &id) const
Get module parameters.
bool hasModule(const JObjectID &id) const
Has module.
Data structure for a composite optical module.
Definition: JModule.hh:75
const JPMT & getPMT(const int index) const
Get PMT.
Definition: JModule.hh:172
Data structure for PMT geometry, calibration and status.
Definition: JPMT.hh:49
Utility class to parse parameter values.
Definition: JProperties.hh:501
const JDirection3D & getDirection() const
Get direction.
double getDot(const JAngle3D &angle) const
Get dot product.
Data structure for position in three dimensions.
Definition: JPosition3D.hh:38
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
Data structure for unit quaternion in three dimensions.
const JQuaternion3D & getQuaternion() const
Get quaternion.
double getAngle() const
Get rotation angle.
double getB() const
Get b value.
double getD() const
Get d value.
double getC() const
Get c value.
double getA() const
Get a value.
double getY() const
Get y position.
Definition: JVector3D.hh:104
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition: JVector3D.hh:270
double getZ() const
Get z position.
Definition: JVector3D.hh:115
double getDot(const JVector3D &vector) const
Get dot product.
Definition: JVector3D.hh:282
double getX() const
Get x position.
Definition: JVector3D.hh:94
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
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys.
Definition: JManager.hh:47
void Write(TDirectory &out, const bool wm=false)
Write objects to file.
Definition: JManager.hh:304
JPosition3D getPosition() const
Get position.
std::string getLabel(const JLocation &location)
Get module label for monitoring and other applications.
Definition: JLocation.hh:247
static JRotation getRotation
Function object to get rotation matrix to go from first to second module.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
double getAngle(const JQuaternion3D &first, const JQuaternion3D &second)
Get space angle between quanternions.
static const JVector3D JVector3Z_t(0, 0, 1)
unit z-vector
double getDistance(const JFirst_t &first, const JSecond_t &second)
Get distance between objects.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
static const std::string TCAL
PMT time offsets.
static const std::string PCAL
(optical|base) module positions
static const std::string SCAL
(module|PMT) status
static const std::string RCAL
optical module orientations
static const std::string ACAL
acoustic time offsets (piezo sensor or hydrophone)
static const std::string CCAL
compass alignment (a.k.a. quaternion calibration)
Definition: JSTDTypes.hh:14
Auxiliary data structure for sequence of same character.
Definition: JManip.hh:330
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:448
int getStatus() const
Get status.
Definition: JStatus.hh:63
Data structure for format specifications.
Definition: JManip.hh:524
Auxiliary data structure for decomposition of quaternion in twist and swing quaternions.
JQuaternion3D twist
rotation around parallel axis
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:68
Auxiliary data structure for running average, standard deviation and quantiles.
Definition: JQuantile.hh:46
double getSTDev() const
Get standard deviation.
Definition: JQuantile.hh:281
void put(const double x, const double w=1.0)
Put value.
Definition: JQuantile.hh:133
double getMean() const
Get mean value.
Definition: JQuantile.hh:252