74{
77
78 struct {
79 double tcal = 0.001;
80 double pcal = 0.001;
81 double rcal = 0.001;
82 double acal = 0.001;
83 double ccal = 0.001;
84 int scal = 0xFFFFFFFF;
85 } precision;
86
87 string detectorFile_a;
88 string detectorFile_b;
91
92 try {
93
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
110
111 zap(argc, argv);
112 }
113 catch(const exception &error) {
114 FATAL(error.what() << endl);
115 }
116
117
120
121 try {
122 load(detectorFile_a, detector_a);
123 }
126 }
127
128 try {
129 load(detectorFile_b, detector_b);
130 }
133 }
134
135 size_t numberOfPMTs = 0;
136
137 bool is_equal = true;
138
141
143
144
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
157
159
160 DEBUG(
" * different UTM position: "
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
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 }
197
198
199
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
217
218 if (module_a->getLocation() == module_b->
getLocation()) {
219
221
222 } else {
223
225 DEBUG(
" * different location: "
227 <<
getLabel(*module_b) <<
" (B)" << endl);
228
229 is_equal = false;
230 }
231
232
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
245
246 if (module_a->getQuaternion() !=
JQuaternion3D(0.0, 0.0, 0.0, 0.0) &&
249
250 DEBUG(
" * different quaternion calibration: "
251 << module_a->getQuaternion() << " (A), "
254
255 is_equal = false;
256 }
257
258
259
261
262 DEBUG(
" * different position: "
263 << module_a->getPosition() << " (A), "
266
267 is_equal = false;
268 }
269
270
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
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: "
292 << endl);
293
294 is_equal = false;
295 }
296 }
297
298
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), "
305 << endl);
306
307 is_equal = false;
308 }
309
310
311
312
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);
318
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
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);
337
338 if (fabs(pmt_a.
getT0() - pmt_b.
getT0()) > precision.tcal) {
339
340 DEBUG(
" * different T0 PMT " << setw(2) << pmt <<
": "
343 <<
", B - A " << pmt_b.
getT0() - pmt_a.
getT0()
344 << endl);
345
346 is_equal = false;
347 }
348 }
349
350
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);
356
357
358
360
361 DEBUG(
" * different PMT position: "
365
366 is_equal = false;
367 }
368 }
369
370
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);
376
377
378
380
381 if ((1.0 - dot) > precision.rcal) {
382
383 DEBUG(
" * different PMT direction: "
387
388 is_equal = false;
389 }
390 }
391
392
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);
398
400
401 DEBUG(
" * different status PMT " << setw(2) << pmt <<
": "
404 << endl);
405
406 is_equal = false;
407 }
408 }
409 }
411
412
414
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
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
436 }
437
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 &&
488
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
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
521
522 return 0;
523}
void setFormat(const JFormat_t &format)
Set format for given type.
#define DEBUG(A)
Message macros.
#define ASSERT(A,...)
Assert macro.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
#define MAKE_CSTRING(A)
Make C-string.
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.
const JLocation & getLocation() const
Get location.
int getFloor() const
Get floor number.
Router for direct addressing of module data in detector data structure.
Data structure for a composite optical module.
const JPMT & getPMT(const int index) const
Get PMT.
Data structure for PMT geometry, calibration and status.
Utility class to parse parameter values.
const JDirection3D & getDirection() const
Get direction.
double getDot(const JAngle3D &angle) const
Get dot product.
Data structure for position in three dimensions.
const JPosition3D & getPosition() const
Get position.
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.
double getDistance(const JVector3D &pos) const
Get distance to point.
double getZ() const
Get z position.
double getDot(const JVector3D &vector) const
Get dot product.
double getX() const
Get x position.
int getID() const
Get identifier.
Utility class to parse command line options.
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys.
JPosition3D getPosition() const
Get position.
std::string getLabel(const JLocation &location)
Get module label for monitoring and other applications.
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)
Auxiliary data structure for sequence of same character.
Auxiliary data structure for floating point format specification.
int getStatus() const
Get status.
Auxiliary data structure for decomposition of quaternion in twist and swing quaternions.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...