Jpp
 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 "JDetector/JDetector.hh"
17 
18 #include "Jeep/JPrint.hh"
19 #include "Jeep/JParser.hh"
20 #include "Jeep/JMessage.hh"
21 
22 
23 namespace {
24 
25  using namespace JPP;
26 
27  /**
28  * Compare time calibration of PMTs between two optical modules.
29  *
30  * \param module_a module
31  * \param module_b module
32  * \return quantile
33  */
34  inline JQuantile getQuantile(const JModule& module_a, const JModule& module_b)
35  {
36  JQuantile Q;
37 
38  JModule::const_iterator pmt_a = module_a.begin();
39  JModule::const_iterator pmt_b = module_b.begin();
40 
41  for ( ; pmt_a != module_a.end() && pmt_b != module_b.end(); ++pmt_a, ++pmt_b) {
42  Q.put(pmt_a->getT0() - pmt_b->getT0());
43  }
44  return Q;
45  }
46 
47 
48  /**
49  * Get ROOT histogram abscissa value.
50  *
51  * \param buffer buffer
52  * \param value value
53  * \return value
54  */
55  inline Double_t getBin(const std::set<int>& buffer, const int value)
56  {
57  return std::distance(buffer.begin(), buffer.find(value));
58  }
59 }
60 
61 
62 /**
63  * \file
64  *
65  * Auxiliary program to find differences between two detector files.\n
66  * A ROOT output file with histograms is optionally produced.
67  * \author mjongen
68  */
69 int main(int argc, char **argv)
70 {
71  using namespace std;
72  using namespace JPP;
73 
74  string detectorFile_a;
75  string detectorFile_b;
76  string outputFile;
77  double precision;
78  int debug;
79 
80  try {
81 
82  JParser<> zap("Auxiliary program to find differences between two detector files.");
83 
84  zap['a'] = make_field(detectorFile_a);
85  zap['b'] = make_field(detectorFile_b);
86  zap['o'] = make_field(outputFile) = "";
87  zap['p'] = make_field(precision) = 0.001; // [ns,m]
88  zap['d'] = make_field(debug) = 3;
89 
90  zap(argc, argv);
91  }
92  catch(const exception &error) {
93  FATAL(error.what() << endl);
94  }
95 
96 
97  JDetector detector_a;
98  JDetector detector_b;
99 
100  try {
101  load(detectorFile_a, detector_a);
102  }
103  catch(const JException& error) {
104  FATAL(error);
105  }
106 
107  try {
108  load(detectorFile_b, detector_b);
109  }
110  catch(const JException& error) {
111  FATAL(error);
112  }
113 
114 
115  bool is_equal = true;
116 
117  const JModuleRouter module_router_a(detector_a);
118  const JModuleRouter module_router_b(detector_b);
119 
120  setFormat<JPosition3D> (JFormat_t(15, 9, std::ios::fixed | std::ios::showpos));
121 
122  // compare detector identifiers
123 
124  if (detector_a.getID() != detector_b.getID()) {
125 
126  DEBUG("* Different detector identifiers "
127  << setw(5) << detector_a.getID() << " (A) and " << endl
128  << setw(5) << detector_b.getID() << " (B)." << endl
129  << endl);
130 
131  is_equal = false;
132  }
133 
134 
135  // check whether the same modules are present in the file
136 
137  for (JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
138 
139  if (!module_router_b.hasModule(module->getID())) {
140 
141  DEBUG("* Module " << setw(10) << module->getID() << " is in A " << getLabel(*module) << " but not in B" << endl);
142 
143  is_equal = false;
144  }
145  }
146 
147  for (JDetector::iterator module = detector_b.begin(); module != detector_b.end(); ++module) {
148 
149  if (!module_router_a.hasModule(module->getID())) {
150 
151  DEBUG("* Module " << setw(10) << module->getID() << " is in B " << getLabel(*module) << " but not in A" << endl);
152 
153  is_equal = false;
154  }
155  }
156  DEBUG(endl);
157 
158 
159  // compare the modules that the files have in common
160 
161  DEBUG("Comparing module by module." << endl);
162 
163  for (JDetector::iterator module_a = detector_a.begin(); module_a != detector_a.end(); ++module_a) {
164 
165  if (!module_router_b.hasModule(module_a->getID())) {
166 
167  continue;
168 
169  is_equal = false;
170  }
171 
172  const JModule* module_b = &module_router_b.getModule(module_a->getID());
173 
174  DEBUG(" Module " << setw(10) << module_a->getID());
175 
176  // string and floor numbers
177 
178  if (module_a->getLocation() == module_b->getLocation()) {
179 
180  DEBUG(" " << getLabel(*module_a) << endl);
181 
182  } else {
183 
184  DEBUG(endl);
185  DEBUG(" * different location: "
186  << getLabel(*module_a) << " (A), "
187  << getLabel(*module_b) << " (B)" << endl);
188 
189  is_equal = false;
190  }
191 
192  // time offset
193 
194  if (fabs(module_a->getT0() - module_b->getT0()) > precision) {
195 
196  DEBUG(" * different T0: "
197  << module_a->getT0() << " (A), "
198  << module_b->getT0() << " (B) "
199  << ", B - A " << module_b->getT0() - module_a->getT0() << endl);
200 
201  is_equal = false;
202  }
203 
204  // quaternion calibration
205 
206  if (getAngle(module_a->getQuaternion(), module_b->getQuaternion()) > precision) {
207 
208  DEBUG(" * different quaternion calibration: "
209  << module_a->getQuaternion() << " (A), "
210  << module_b->getQuaternion() << " (B) "
211  << ", B - A " << getAngle(module_b->getQuaternion(), module_a->getQuaternion()) << endl);
212 
213  is_equal = false;
214  }
215 
216  // average DOM positions
217 
218  if (getDistance(module_a->getPosition(), module_b->getPosition()) > precision) {
219 
220  DEBUG(" * different position: "
221  << module_a->getPosition() << " (A), "
222  << module_b->getPosition() << " (B)"
223  << ", B - A " << JPosition3D(module_b->getPosition() - module_a->getPosition()) << endl);
224 
225  is_equal = false;
226  }
227 
228  // number of PMTs
229 
230  if (module_a->size() != module_b->size()) {
231 
232  DEBUG(" * different number of PMTs: "
233  << module_a->size() << " (A), "
234  << module_b->size() << " (B)" << endl);
235 
236  is_equal = false;
237  }
238 
239  // average t0
240 
241  if (!module_a->empty() &&
242  !module_b->empty()) {
243 
244  const JQuantile q = getQuantile(*module_a, *module_b);
245 
246  if (fabs(q.getMean()) > precision) {
247 
248  DEBUG(" * different average t0: "
249  << ", B - A " << q.getMean()
250  << endl);
251 
252  is_equal = false;
253  }
254  }
255 
256  // global t0
257 
258  if (fabs(module_a->getT0() - module_b->getT0()) > precision) {
259 
260  DEBUG(" * different global t0: "
261  << module_a->getT0() << " (A), "
262  << module_b->getT0() << " (B)"
263  << ", B - A " << module_b->getT0() - module_a->getT0()
264  << endl);
265 
266  is_equal = false;
267  }
268 
269  // comparing by PMT
270 
271  // t0
272 
273  for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
274 
275  const JPMT& pmt_a = module_a->getPMT(pmt);
276  const JPMT& pmt_b = module_b->getPMT(pmt);
277 
278  if (fabs(pmt_a.getT0() - pmt_b.getT0()) > precision) {
279 
280  DEBUG(" * different T0 PMT " << pmt << ": "
281  << pmt_a.getT0() << " (A" << FILL(2,'0') << pmt << "), "
282  << pmt_b.getT0() << " (B" << FILL(2,'0') << pmt << ")"
283  << ", B - A " << pmt_b.getT0() - pmt_a.getT0()
284  << endl);
285 
286  is_equal = false;
287  }
288  }
289 
290  // positions
291 
292  for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
293 
294  const JPMT& pmt_a = module_a->getPMT(pmt);
295  const JPMT& pmt_b = module_b->getPMT(pmt);
296 
297  // PMT positions
298 
299  if (pmt_a.getPosition().getDistance(pmt_b.getPosition()) > precision) {
300 
301  DEBUG(" * different PMT position: "
302  << pmt_a.getPosition() << " (A" << FILL(2,'0') << pmt << "), "
303  << pmt_b.getPosition() << " (B" << FILL(2,'0') << pmt << ")"
304  << ", B - A " << JPosition3D(pmt_b.getPosition() - pmt_a.getPosition()) << endl);
305 
306  is_equal = false;
307  }
308  }
309 
310  // status
311 
312  for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
313 
314  const JPMT& pmt_a = module_a->getPMT(pmt);
315  const JPMT& pmt_b = module_b->getPMT(pmt);
316 
317  if (pmt_a.getStatus() != pmt_b.getStatus()) {
318 
319  DEBUG(" * different status PMT " << pmt << ": "
320  << pmt_a.getStatus() << " (A" << FILL(2,'0') << pmt << "), "
321  << pmt_b.getStatus() << " (B" << FILL(2,'0') << pmt << ")"
322  << endl);
323 
324  is_equal = false;
325  }
326  }
327  }
328  DEBUG(endl);
329 
330 
331  if (outputFile != "") {
332 
333  TFile out(outputFile.c_str(), "recreate");
334 
335  set<int> string;
336  set<int> floor;
337 
338  for (JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
339  string.insert(module->getString());
340  floor .insert(module->getFloor ());
341  }
342 
343  for (JDetector::iterator module = detector_b.begin(); module != detector_b.end(); ++module) {
344  string.insert(module->getString());
345  floor .insert(module->getFloor ());
346  }
347 
348 
349  TH2D M2("M2", NULL,
350  string.size(), -0.5, string.size() - 0.5,
351  floor .size(), -0.5, floor .size() - 0.5);
352 
353  for (set<int>::const_iterator i = string.begin(); i != string.end(); ++i) {
354  M2.GetXaxis()->SetBinLabel(distance(string.begin(), i) + 1, MAKE_CSTRING(*i));
355  }
356 
357  for (set<int>::const_iterator i = floor.begin(); i != floor.end(); ++i) {
358  M2.GetYaxis()->SetBinLabel(distance(floor.begin(), i) + 1, MAKE_CSTRING(*i));
359  }
360 
361  TH2D* X2 = (TH2D*) M2.Clone("X2" );
362  TH2D* Y2 = (TH2D*) M2.Clone("Y2" );
363  TH2D* Z2 = (TH2D*) M2.Clone("Z2" );
364  TH2D* T2 = (TH2D*) M2.Clone("T2" );
365  TH2D* RM2 = (TH2D*) M2.Clone("RM2");
366  TH2D* R2 = (TH2D*) M2.Clone("R2" );
367  TH2D* QA = (TH2D*) M2.Clone("QA" );
368  TH2D* QB = (TH2D*) M2.Clone("QB" );
369  TH2D* QC = (TH2D*) M2.Clone("QC" );
370  TH2D* QD = (TH2D*) M2.Clone("QD" );
371 
372 
373  for (JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
374  if (!module_router_b.hasModule(module->getID()) ) {
375  M2.Fill(module->getString(), module->getFloor(), -1.0);
376  }
377  }
378 
379  for (JDetector::iterator module = detector_b.begin(); module != detector_b.end(); ++module) {
380  if (!module_router_a.hasModule(module->getID()) ) {
381  M2.Fill(module->getString(), module->getFloor(), +1.0);
382  }
383  }
384 
385  for (JDetector::iterator module_a = detector_a.begin(); module_a != detector_a.end(); ++module_a) {
386 
387  if (!module_router_b.hasModule(module_a->getID())) {
388  continue;
389  }
390 
391  const JModule* module_b = &module_router_b.getModule(module_a->getID());
392 
393  X2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), module_a->getX() - module_b->getX());
394  Y2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), module_a->getY() - module_b->getY());
395  Z2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), module_a->getZ() - module_b->getZ());
396 
397  if (module_a->getFloor() != 0 &&
398  module_b->getFloor() != 0) {
399 
400  double phi = 0.0;
401 
402  for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
403 
404  double x = module_a->getPMT(pmt).getPhi() - module_b->getPMT(pmt).getPhi();
405 
406  while (x > +PI) { x -= PI; }
407  while (x < -PI) { x += PI; }
408 
409  phi += x;
410  }
411 
412  phi /= min(module_a->size(),
413  module_b->size());
414 
415  const JQuaternion3D Q = getRotation(*module_a, *module_b);
416  const JQuantile q = getQuantile(*module_a, *module_b);
417 
418  R2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), phi * 180.0/PI);
419  QA ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), Q.getA());
420  QB ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), Q.getB());
421  QC ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), Q.getC());
422  QD ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), Q.getD());
423  T2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), q.getMean());
424  RM2->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), q.getRMS());
425  }
426  }
427 
428  out.Write();
429  out.Close();
430  }
431 
432  ASSERT(is_equal);
433 
434  return 0;
435 }
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
const JStatus & getStatus() const
Get status.
Definition: JStatus.hh:68
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:57
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Quantile calculator.
Definition: JQuantile.hh:88
std::string getLabel(const JLocation &location)
Get module label for monitoring and other applications.
Definition: JLocation.hh:246
Detector data structure.
Definition: JDetector.hh:80
Router for direct addressing of module data in detector data structure.
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
double getRMS() const
Get RMS.
Definition: JQuantile.hh:296
double getPhi() const
Get phi angle.
Definition: JVersor3D.hh:144
string outputFile
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition: JVector3D.hh:270
Data structure for detector geometry and calibration.
double getMean() const
Get mean value.
Definition: JQuantile.hh:282
double getDistance(const JFirst_t &first, const JSecond_t &second)
Get distance between objects.
const JQuaternion3D & getQuaternion() const
Get quaternion.
#define ASSERT(A,...)
Assert macro.
Definition: JMessage.hh:90
I/O formatting auxiliaries.
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 and calibration.
Definition: JPMT.hh:47
void put(const double x, const double w=1.0)
Put value.
Definition: JQuantile.hh:173
static const double PI
Mathematical constants.
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:211
General purpose messaging.
Auxiliary data structure for sequence of same character.
Definition: JManip.hh:327
#define FATAL(A)
Definition: JMessage.hh:67
Data structure for unit quaternion in three dimensions.
Direct access to module in detector data structure.
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.
bool hasModule(const JObjectID &id) const
Has module.
double getX() const
Get x position.
Definition: JVector3D.hh:94
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
Data structure for format specifications.
Definition: JManip.hh:521
double getZ() const
Get z position.
Definition: JVector3D.hh:115
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
double getT0() const
Get time offset.
int main(int argc, char *argv[])
Definition: Main.cpp:15