Jpp
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 
10 #include "JGeometry3D/JQuaternion3D.hh"
11 #include "JTools/JRange.hh"
12 #include "JTools/JConstants.hh"
13 #include "JTools/JQuantile.hh"
14 #include "JDetector/JDetector.hh"
15 #include "JDetector/JDetectorToolkit.hh"
16 #include "JDetector/JModuleRouter.hh"
17 
18 #include "Jeep/JPrint.hh"
19 #include "Jeep/JParser.hh"
20 #include "Jeep/JMessage.hh"
21 
22 
23 namespace {
24 
25  using JDETECTOR::JModule;
26  using namespace JPP ;
27  /**
28  * Compare the t0s of two optical modules.
29  *
30  * \param module_a module
31  * \param module_b module
32  * \return q
33  */
34  inline JQuantile compareT0(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 /**
50  * \file
51  *
52  * Auxiliary program to find differences between two detector files.\n
53  * A ROOT outout file with histograms is optionally produced.
54  * \author mjongen
55  */
56 int main(int argc, char **argv)
57 {
58  using namespace std;
59  using namespace JPP;
60 
61  string detectorFile_a;
62  string detectorFile_b;
63  string outputFile;
64  double precision;
65  int debug;
66 
67  try {
68 
69  JParser<> zap("Auxiliary program to find differences between two detector files.");
70 
71  zap['a'] = make_field(detectorFile_a);
72  zap['b'] = make_field(detectorFile_b);
73  zap['o'] = make_field(outputFile) = "";
74  zap['p'] = make_field(precision) = 0.001; // [ns,m]
75  zap['d'] = make_field(debug) = 3;
76 
77  zap(argc, argv);
78  }
79  catch(const exception &error) {
80  FATAL(error.what() << endl);
81  }
82 
83 
84  JDetector detector_a;
85  JDetector detector_b;
86 
87  try {
88  load(detectorFile_a, detector_a);
89  }
90  catch(const JException& error) {
91  FATAL(error);
92  }
93 
94  try {
95  load(detectorFile_b, detector_b);
96  }
97  catch(const JException& error) {
98  FATAL(error);
99  }
100 
101 
102  bool is_equal = true;
103 
104  const JModuleRouter module_router_a(detector_a);
105  const JModuleRouter module_router_b(detector_b);
106 
107 
108  // compare detector IDs
109 
110  if (detector_a.getID() != detector_b.getID()) {
111 
112  DEBUG("* Different detector identifiers "
113  << setw(5) << detector_a.getID() << " (A) and " << endl
114  << setw(5) << detector_b.getID() << " (B)." << endl
115  << endl);
116 
117  is_equal = false;
118  }
119 
120 
121  // check whether the same modules are present in the file
122 
123  for (JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
124 
125  if (!module_router_b.hasModule(module->getID())) {
126 
127  DEBUG("* Module " << setw(10) << module->getID() << " is in A " << getModuleLabel(*module) << " but not in B" << endl);
128 
129  is_equal = false;
130  }
131  }
132 
133  for (JDetector::iterator module = detector_b.begin(); module != detector_b.end(); ++module) {
134 
135  if (!module_router_a.hasModule(module->getID())) {
136 
137  DEBUG("* Module " << setw(10) << module->getID() << " is in B " << getModuleLabel(*module) << " but not in A" << endl);
138 
139  is_equal = false;
140  }
141  }
142  DEBUG(endl);
143 
144 
145  // compare the modules that the files have in common
146 
147  DEBUG("Comparing module by module." << endl);
148 
149  for (JDetector::iterator module_a = detector_a.begin(); module_a != detector_a.end(); ++module_a) {
150 
151  if (!module_router_b.hasModule(module_a->getID())) {
152 
153  continue;
154 
155  is_equal = false;
156  }
157 
158  const JModule* module_b = &module_router_b.getModule(module_a->getID());
159 
160  DEBUG(" Module " << setw(10) << module_a->getID());
161 
162  // string and floor numbers
163 
164  if (module_a->getLocation() == module_b->getLocation()) {
165 
166  DEBUG(" " << getModuleLabel(*module_a) << endl);
167 
168  } else {
169 
170  DEBUG(endl);
171  DEBUG(" * different location: "
172  << getModuleLabel(*module_a) << " (A), "
173  << getModuleLabel(*module_b) << " (B)" << endl);
174 
175  is_equal = false;
176  }
177 
178  // average DOM positions
179 
180  if (module_a->getPosition().getDistance(module_b->getPosition()) > precision) {
181 
182  DEBUG(" * different position: "
183  << module_a->getPosition() << " (A), "
184  << module_b->getPosition() << " (B)"
185  << ", B - A " << JPosition3D(module_b->getPosition()-module_a->getPosition()) << endl);
186 
187  is_equal = false;
188  }
189 
190  // number of PMTs
191 
192  if (module_a->size() != module_b->size()) {
193 
194  DEBUG(" * different number of PMTs: "
195  << module_a->size() << " (A), "
196  << module_b->size() << " (B)" << endl);
197 
198  is_equal = false;
199  }
200 
201  // average t0
202 
203  const JQuantile q = compareT0(*module_a, *module_b);
204 
205  if (fabs(q.getMean()) > precision) {
206 
207  DEBUG(" * different average t0: "
208  << ", B - A " << q.getMean()
209  << endl);
210 
211  is_equal = false;
212  }
213 
214  // comparing by PMT
215 
216  // t0
217 
218  for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
219 
220  const JPMT& pmt_a = module_a->getPMT(pmt);
221  const JPMT& pmt_b = module_b->getPMT(pmt);
222 
223  if (fabs(pmt_a.getT0() - pmt_b.getT0()) > precision) {
224 
225  DEBUG(" * different T0 PMT " << pmt << ": "
226  << pmt_a.getT0() << " (A" << FILL(2,'0') << pmt << "), "
227  << pmt_b.getT0() << " (B" << FILL(2,'0') << pmt << ")"
228  << ", B - A " << pmt_b.getT0() - pmt_a.getT0()
229  << endl);
230 
231  is_equal = false;
232  }
233  }
234 
235  // positions
236 
237  for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
238 
239  const JPMT& pmt_a = module_a->getPMT(pmt);
240  const JPMT& pmt_b = module_b->getPMT(pmt);
241 
242  // PMT positions
243 
244  if (pmt_a.getPosition().getDistance(pmt_b.getPosition()) > precision) {
245 
246  DEBUG(" * different PMT position: "
247  << pmt_a.getPosition() << " (A" << FILL(2,'0') << pmt << "), "
248  << pmt_b.getPosition() << " (B" << FILL(2,'0') << pmt << ")"
249  << ", B - A " << JPosition3D(pmt_b.getPosition() - pmt_a.getPosition()) << endl);
250 
251  is_equal = false;
252  }
253  }
254 
255  // status
256 
257  for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
258 
259  const JPMT& pmt_a = module_a->getPMT(pmt);
260  const JPMT& pmt_b = module_b->getPMT(pmt);
261 
262  if (pmt_a.getStatus() != pmt_b.getStatus()) {
263 
264  DEBUG(" * different status PMT " << pmt << ": "
265  << pmt_a.getStatus() << " (A" << FILL(2,'0') << pmt << "), "
266  << pmt_b.getStatus() << " (B" << FILL(2,'0') << pmt << ")"
267  << endl);
268 
269  is_equal = false;
270  }
271  }
272  }
273  DEBUG(endl);
274 
275 
276  if (outputFile != "") {
277 
278  typedef JRange<int> JRange_t;
279 
280  const JRange_t string = combine(JRange_t(detector_a.begin(), detector_a.end(), &JModule::getString),
281  JRange_t(detector_b.begin(), detector_b.end(), &JModule::getString));
282  const JRange_t floor = combine(JRange_t(detector_a.begin(), detector_a.end(), &JModule::getFloor),
283  JRange_t(detector_b.begin(), detector_b.end(), &JModule::getFloor));
284 
285  TFile out(outputFile.c_str(), "recreate");
286 
287  TH2D M2("M2", NULL,
288  string.getLength() + 1,
289  string.getLowerLimit() - 0.5,
290  string.getUpperLimit() + 0.5,
291  floor.getLength() + 1,
292  floor.getLowerLimit() - 0.5,
293  floor.getUpperLimit() + 0.5);
294 
295  TH2D* X2 = (TH2D*) M2.Clone("X2" );
296  TH2D* Y2 = (TH2D*) M2.Clone("Y2" );
297  TH2D* Z2 = (TH2D*) M2.Clone("Z2" );
298  TH2D* T2 = (TH2D*) M2.Clone("T2" );
299  TH2D* RM2 = (TH2D*) M2.Clone("RM2");
300  TH2D* R2 = (TH2D*) M2.Clone("R2" );
301  TH2D* QA = (TH2D*) M2.Clone("QA" );
302  TH2D* QB = (TH2D*) M2.Clone("QB" );
303  TH2D* QC = (TH2D*) M2.Clone("QC" );
304  TH2D* QD = (TH2D*) M2.Clone("QD" );
305 
306  for( JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
307  if( !module_router_b.hasModule(module->getID()) ) {
308  M2.Fill(module->getString(), module->getFloor(), -1.0);
309  }
310  }
311 
312  for( JDetector::iterator module = detector_b.begin(); module != detector_b.end(); ++module) {
313  if( !module_router_a.hasModule(module->getID()) ) {
314  M2.Fill(module->getString(), module->getFloor(), +1.0);
315  }
316  }
317 
318  for( JDetector::iterator module_a = detector_a.begin(); module_a != detector_a.end(); ++module_a) {
319 
320  if (!module_router_b.hasModule(module_a->getID())) {
321  continue;
322  }
323 
324  const JModule* module_b = &module_router_b.getModule(module_a->getID());
325 
326  double phi = 0.0;
327 
328  for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
329 
330  double x = module_a->getPMT(pmt).getPhi() - module_b->getPMT(pmt).getPhi();
331 
332  while (x > +PI) { x -= PI; }
333  while (x < -PI) { x += PI; }
334 
335  phi += x;
336  }
337 
338  phi /= min(module_a->size(),
339  module_b->size());
340 
341  JQuaternion3D Q = getRotation(*module_a, *module_b);
342 
343  JQuantile q = compareT0(*module_a, *module_b);
344 
345  X2 ->Fill(module_a->getString(), module_a->getFloor(), module_a->getX() - module_b->getX());
346  Y2 ->Fill(module_a->getString(), module_a->getFloor(), module_a->getY() - module_b->getY());
347  Z2 ->Fill(module_a->getString(), module_a->getFloor(), module_a->getZ() - module_b->getZ());
348  T2 ->Fill(module_a->getString(), module_a->getFloor(), q.getMean());
349  RM2->Fill(module_a->getString(), module_a->getFloor(), q.getRMS());
350  R2 ->Fill(module_a->getString(), module_a->getFloor(), phi * 180.0/PI);
351  QA ->Fill(module_a->getString(), module_a->getFloor(), Q.getA());
352  QB ->Fill(module_a->getString(), module_a->getFloor(), Q.getB());
353  QC ->Fill(module_a->getString(), module_a->getFloor(), Q.getC());
354  QD ->Fill(module_a->getString(), module_a->getFloor(), Q.getD());
355  }
356 
357  out.Write();
358  out.Close();
359  }
360 
361  ASSERT(is_equal);
362 
363  return 0;
364 }
JDETECTOR::getModuleLabel
std::string getModuleLabel(const JModuleLocation &location)
Get module label for monitoring and other applications.
Definition: JDetectorToolkit.hh:688
JMessage.hh
ASSERT
#define ASSERT(A,...)
Assert macro.
Definition: JMessage.hh:90
JPrint.hh
JDETECTOR::JModuleRouter::getModule
const JModule & getModule(const JObjectID &id) const
Get module parameters.
Definition: JModuleRouter.hh:89
JGEOMETRY3D::JVector3D::getZ
double getZ() const
Get z position.
Definition: JVector3D.hh:114
JDETECTOR::load
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:456
JGEOMETRY3D::JVector3D::getDistance
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition: JVector3D.hh:268
JDETECTOR::JCalibration::getT0
double getT0() const
Get time offset.
Definition: JDetector/JCalibration.hh:83
JTOOLS::JRange< int >
JEEP::JStatus::getStatus
int getStatus() const
Get status.
Definition: JStatus.hh:56
JQuaternion3D.hh
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1493
JQuantile.hh
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:37
JRange.hh
JDETECTOR::JModule::getPMT
const JPMT & getPMT(const int index) const
Get PMT.
Definition: JModule.hh:173
JTOOLS::JQuantile::put
void put(const double x, const double w=1.0)
Put value.
Definition: JQuantile.hh:123
debug
int debug
debug level
Definition: JSirene.cc:59
JTOOLS::JQuantile::getRMS
double getRMS() const
Get RMS.
Definition: JQuantile.hh:231
FILL
Auxiliary data structure for sequence of same character.
Definition: JPrint.hh:361
JConstants.hh
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition: JPosition3D.hh:35
JLANG::JObjectID::getID
int getID() const
Get identifier.
Definition: JObjectID.hh:55
JModuleRouter.hh
JDETECTOR::JModule
Data structure for a composite optical module.
Definition: JModule.hh:49
JGEOMETRY3D::JQuaternion3D
Data structure for quaternion in three dimensions.
Definition: JQuaternion3D.hh:240
JDETECTOR::JPMT
Data structure for PMT geometry and calibration.
Definition: JPMT.hh:53
JTOOLS::JQuantile::getMean
double getMean() const
Get mean value.
Definition: JQuantile.hh:217
JTOOLS::combine
JRange< T, JComparator_t > combine(const JRange< T, JComparator_t > &first, const JRange< T, JComparator_t > &second)
Combine ranges.
Definition: JRange.hh:688
JParser.hh
JTOOLS::JQuantile
Quantile calculator.
Definition: JQuantile.hh:34
JDetectorToolkit.hh
JDETECTOR::JModuleLocation::getLocation
const JModuleLocation & getLocation() const
Get location.
Definition: JModuleLocation.hh:68
JDETECTOR::JModuleRouter
Router for direct addressing of module data in detector data structure.
Definition: JModuleRouter.hh:34
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1954
JDETECTOR::JDetector
Detector data structure.
Definition: JDetector.hh:80
JGEOMETRY3D::JPosition3D::getPosition
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:129
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
std
Definition: jaanetDictionary.h:36
JDETECTOR::JModuleRouter::hasModule
bool hasModule(const JObjectID &id) const
Has module.
Definition: JModuleRouter.hh:101
JGEOMETRY3D::JVector3D::getY
double getY() const
Get y position.
Definition: JVector3D.hh:103
JTOOLS::PI
static const double PI
Constants.
Definition: JConstants.hh:20
main
int main(int argc, char **argv)
Definition: JCompareDetector.cc:56
JDetector.hh
JGEOMETRY3D::JVector3D::getX
double getX() const
Get x position.
Definition: JVector3D.hh:93
JDETECTOR::JModuleLocation::getFloor
int getFloor() const
Get floor number.
Definition: JModuleLocation.hh:144
JDETECTOR::getRotation
static JRotation getRotation
Function object to get rotation matrix to go from first to second module.
Definition: JDetectorToolkit.hh:823
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JLANG::JException
General exception.
Definition: JException.hh:40
JDETECTOR::JModuleLocation::getString
int getString() const
Get string number.
Definition: JModuleLocation.hh:133
JGEOMETRY3D::JVersor3D::getPhi
double getPhi() const
Get phi angle.
Definition: JVersor3D.hh:140
Generated by   doxygen 1.8.16