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