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