Jpp 20.0.0-27-g39925593c-D
the software that should make you happy
Loading...
Searching...
No Matches
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 "JMath/JConstants.hh"
12#include "JTools/JRange.hh"
13#include "JTools/JStats.hh"
14#include "JROOT/JManager.hh"
15#include "JROOT/JRootToolkit.hh"
16#include "JDetector/JDetector.hh"
17#include "JDetector/JDetectorToolkit.hh"
18#include "JDetector/JDetectorCalibration.hh"
19#include "JDetector/JModuleRouter.hh"
20
21#include "Jeep/JProperties.hh"
22#include "Jeep/JParser.hh"
23#include "Jeep/JMessage.hh"
24
25
26namespace {
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 JStats getQuantile(const JModule& module_a, const JModule& module_b)
38 {
39 JStats 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_b->getT0() - pmt_a->getT0());
46 }
47
48 return Q;
49 }
50
51
52 /**
53 * Get ROOT histogram abscissa value.
54 *
55 * \param buffer buffer
56 * \param value value
57 * \return value
58 */
59 inline Double_t getBin(const std::set<int>& buffer, const int value)
60 {
61 return std::distance(buffer.begin(), buffer.find(value));
62 }
63}
64
65
66/**
67 * \file
68 *
69 * Auxiliary program to find differences between two detector files.\n
70 * A ROOT output file with histograms is optionally produced.
71 * \author mjongen
72 */
73int main(int argc, char **argv)
74{
75 using namespace std;
76 using namespace JPP;
77
78 struct {
79 double tcal = 0.001; // [ns]
80 double pcal = 0.001; // [m]
81 double rcal = 0.001; // [rad]
82 double acal = 0.001; // [ns]
83 double ccal = 0.001; // [deg]
84 int scal = 0xFFFFFFFF; // [status]
85 } precision;
86
87 string detectorFile_a;
88 string detectorFile_b;
89 string outputFile;
90 int debug;
91
92 try {
93
94 JProperties properties;
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
105 zap['a'] = make_field(detectorFile_a);
106 zap['b'] = make_field(detectorFile_b);
107 zap['o'] = make_field(outputFile) = "";
108 zap['p'] = make_field(properties) = JPARSER::initialised();
109 zap['d'] = make_field(debug) = 3;
110
111 zap(argc, argv);
112 }
113 catch(const exception &error) {
114 FATAL(error.what() << endl);
115 }
116
117
118 JDetector detector_a;
119 JDetector detector_b;
120
121 try {
122 load(detectorFile_a, detector_a);
123 }
124 catch(const JException& error) {
125 FATAL(error);
126 }
127
128 try {
129 load(detectorFile_b, detector_b);
130 }
131 catch(const JException& error) {
132 FATAL(error);
133 }
134
135 size_t numberOfPMTs = 0;
136
137 bool is_equal = true;
138
139 const JModuleRouter module_router_a(detector_a);
140 const JModuleRouter module_router_b(detector_b);
141
142 setFormat<JPosition3D> (JFormat_t(15, 9, std::ios::fixed | std::ios::showpos));
143
144 // compare detector identifiers
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 // UTM positions
157
158 if (getDistance(detector_a.getPosition(), detector_b.getPosition()) > precision.pcal) {
159
160 DEBUG(" * different UTM position: "
161 << detector_a.getPosition() << " (A), "
162 << detector_b.getPosition() << " (B)"
163 << ", B - A " << JPosition3D(detector_b.getPosition() - detector_a.getPosition()) << endl);
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 // check whether the same modules are present in the file
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 }
196 DEBUG(endl);
197
198
199 // compare the modules that the files have in common
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 // string and floor numbers
217
218 if (module_a->getLocation() == module_b->getLocation()) {
219
220 DEBUG(" " << getLabel(*module_a) << endl);
221
222 } else {
223
224 DEBUG(endl);
225 DEBUG(" * different location: "
226 << getLabel(*module_a) << " (A), "
227 << getLabel(*module_b) << " (B)" << endl);
228
229 is_equal = false;
230 }
231
232 // time offset
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 // quaternion calibration
245
246 if (module_a->getQuaternion() != JQuaternion3D(0.0, 0.0, 0.0, 0.0) &&
247 module_b->getQuaternion() != JQuaternion3D(0.0, 0.0, 0.0, 0.0) &&
248 getAngle(module_a->getQuaternion(), module_b->getQuaternion()) > precision.ccal) {
249
250 DEBUG(" * different quaternion calibration: "
251 << module_a->getQuaternion() << " (A), "
252 << module_b->getQuaternion() << " (B) "
253 << ", B - A " << getAngle(module_b->getQuaternion(), module_a->getQuaternion()) << endl);
254
255 is_equal = false;
256 }
257
258 // average DOM positions
259
260 if (getDistance(module_a->getPosition(), module_b->getPosition()) > precision.pcal) {
261
262 DEBUG(" * different position: "
263 << module_a->getPosition() << " (A), "
264 << module_b->getPosition() << " (B)"
265 << ", B - A " << JPosition3D(module_b->getPosition() - module_a->getPosition()) << endl);
266
267 is_equal = false;
268 }
269
270 // number of PMTs
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 // average t0
282
283 if (!module_a->empty() &&
284 !module_b->empty()) {
285
286 const JStats q = getQuantile(*module_a, *module_b);
287
288 if (fabs(q.getMean()) > precision.tcal) {
289
290 DEBUG(" * different average t0: "
291 << ", B - A " << q.getMean()
292 << endl);
293
294 is_equal = false;
295 }
296 }
297
298 // status
299
300 if (module_a->getAND(precision.scal) != module_b->getAND(precision.scal)) {
301
302 DEBUG(" * different status module " << module_a->getID() << ": "
303 << module_a->getStatus() << " (A), "
304 << module_b->getStatus() << " (B)"
305 << endl);
306
307 is_equal = false;
308 }
309
310 // comparing by PMT
311
312 // identifier
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);
317 const JPMT& pmt_b = module_b->getPMT(pmt);
318
319 if (pmt_a.getID() != pmt_b.getID()) {
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 // t0
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);
336 const JPMT& pmt_b = module_b->getPMT(pmt);
337
338 if (fabs(pmt_a.getT0() - pmt_b.getT0()) > precision.tcal) {
339
340 DEBUG(" * different T0 PMT " << setw(2) << pmt << ": "
341 << FIXED(12,3) << pmt_a.getT0() << " (A" << FILL(2,'0') << pmt << "), " << FILL()
342 << FIXED(12,3) << pmt_b.getT0() << " (B" << FILL(2,'0') << pmt << ")" << FILL()
343 << ", B - A " << pmt_b.getT0() - pmt_a.getT0()
344 << endl);
345
346 is_equal = false;
347 }
348 }
349
350 // positions
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);
355 const JPMT& pmt_b = module_b->getPMT(pmt);
356
357 // PMT positions
358
359 if (pmt_a.getPosition().getDistance(pmt_b.getPosition()) > precision.pcal) {
360
361 DEBUG(" * different PMT position: "
362 << pmt_a.getPosition() << " (A" << FILL(2,'0') << pmt << "), " << FILL()
363 << pmt_b.getPosition() << " (B" << FILL(2,'0') << pmt << ")" << FILL()
364 << ", B - A " << JPosition3D(pmt_b.getPosition() - pmt_a.getPosition()) << endl);
365
366 is_equal = false;
367 }
368 }
369
370 // orientations
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);
375 const JPMT& pmt_b = module_b->getPMT(pmt);
376
377 // PMT positions
378
379 const double dot = pmt_a.getDirection().getDot(pmt_b.getDirection());
380
381 if ((1.0 - dot) > precision.rcal) {
382
383 DEBUG(" * different PMT direction: "
384 << pmt_a.getDirection() << " (A" << FILL(2,'0') << pmt << "), " << FILL()
385 << pmt_b.getDirection() << " (B" << FILL(2,'0') << pmt << ")" << FILL()
386 << ", B - A " << JPosition3D(JPosition3D(pmt_b.getDirection()) - JPosition3D(pmt_a.getDirection())) << ' ' << dot << endl);
387
388 is_equal = false;
389 }
390 }
391
392 // status
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);
397 const JPMT& pmt_b = module_b->getPMT(pmt);
398
399 if (pmt_a.getAND(precision.scal) != pmt_b.getAND(precision.scal)) {
400
401 DEBUG(" * different status PMT " << setw(2) << pmt << ": "
402 << pmt_a.getStatus() << " (A" << FILL(2,'0') << pmt << "), " << FILL()
403 << pmt_b.getStatus() << " (B" << FILL(2,'0') << pmt << ")" << FILL()
404 << endl);
405
406 is_equal = false;
407 }
408 }
409 }
410 DEBUG(endl);
411
412
413 if (outputFile != "") {
414
415 set<int> string;
416 set<int> floor;
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
428 JManager<int, TH1D> H1(new TH1D("M[%]", NULL, numberOfPMTs, -0.5, numberOfPMTs - 0.5));
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
434 for (set<int>::const_iterator i = string.begin(); i != string.end(); ++i) {
435 M2.GetXaxis()->SetBinLabel(distance(string.begin(), i) + 1, MAKE_CSTRING(*i));
436 }
437
438 for (set<int>::const_iterator i = floor.begin(); i != floor.end(); ++i) {
439 M2.GetYaxis()->SetBinLabel(distance(floor.begin(), i) + 1, MAKE_CSTRING(*i));
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 &&
487 module_b->getFloor() != 0) {
488
489 const JQuaternion3D Q = getRotation(*module_a, *module_b);
490 const JStats q = getQuantile(*module_a, *module_b);
491
492 const JQuaternion3D::decomposition q1(Q, JVector3Z_t);
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
508 TFile out(outputFile.c_str(), "recreate");
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
520 ASSERT(is_equal);
521
522 return 0;
523}
main
int main(int argc, char **argv)
Definition JCompareDetector.cc:73
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
JDetector.hh
Data structure for detector geometry and calibration.
JManager.hh
Dynamic ROOT object management.
JConstants.hh
Mathematical constants.
JMessage.hh
General purpose messaging.
DEBUG
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
ASSERT
#define ASSERT(A,...)
Assert macro.
Definition JMessage.hh:90
FATAL
#define FATAL(A)
Definition JMessage.hh:67
debug
int debug
debug level
Definition JSirene.cc:72
JModuleRouter.hh
Direct access to module in detector data structure.
JParser.hh
Utility class to parse command line options.
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
MAKE_CSTRING
#define MAKE_CSTRING(A)
Make C-string.
Definition JPrint.hh:72
JProperties.hh
Utility class to parse parameter values.
JRange.hh
Auxiliary class to define a range between two values.
Y2
#define Y2(d1, c12, d2)
Definition LambertW.hh:82
distance
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Definition build/antares-dataformat/PhysicsEvent.hh:434
JDETECTOR::JDetector
Detector data structure.
Definition JDetector.hh:96
JDETECTOR::JModuleRouter
Router for direct addressing of module data in detector data structure.
Definition JModuleRouter.hh:38
JDETECTOR::JModule
Data structure for a composite optical module.
Definition JModule.hh:76
JDETECTOR::JPMT
Data structure for PMT geometry, calibration and status.
Definition JPMT.hh:49
JEEP::JProperties
Utility class to parse parameter values.
Definition JProperties.hh:503
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition JPosition3D.hh:38
JGEOMETRY3D::JQuaternion3D
Data structure for unit quaternion in three dimensions.
Definition JQuaternion3D.hh:309
JGEOMETRY3D::JQuaternion3D::getAngle
double getAngle() const
Get rotation angle.
Definition JQuaternion3D.hh:531
JGEOMETRY3D::JQuaternion3D::getB
double getB() const
Get b value.
Definition JQuaternion3D.hh:556
JGEOMETRY3D::JQuaternion3D::getD
double getD() const
Get d value.
Definition JQuaternion3D.hh:578
JGEOMETRY3D::JQuaternion3D::getC
double getC() const
Get c value.
Definition JQuaternion3D.hh:567
JGEOMETRY3D::JQuaternion3D::getA
double getA() const
Get a value.
Definition JQuaternion3D.hh:545
JLANG::JException
General exception.
Definition JException.hh:24
JOSCPROB::JOscProbInterpolator
Template definition of a multi-dimensional oscillation probability interpolation table.
Definition JOscProbInterpolator.hh:71
JTOOLS::JMultiFunction::insert
void insert(const JMultiFunction< __JFunction_t, __JMaplist_t, __JDistance_t > &input)
Insert multidimensional input.
Definition JMultiFunction.hh:121
JDETECTOR::getLabel
std::string getLabel(const JLocation &location)
Get module label for monitoring and other applications.
Definition JLocation.hh:247
JDETECTOR::getRotation
static JRotation getRotation
Function object to get rotation matrix to go from first to second module.
Definition JDetectorToolkit.hh:801
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition JDetectorToolkit.hh:520
JGEOMETRY3D::getAngle
double getAngle(const JQuaternion3D &first, const JQuaternion3D &second)
Get space angle between quanternions.
Definition JQuaternion3D.hh:1165
JGEOMETRY3D::JVector3Z_t
static const JVector3D JVector3Z_t(0, 0, 1)
unit z-vector
JMATH::getDistance
double getDistance(const JFirst_t &first, const JSecond_t &second)
Get distance between objects.
Definition JMathToolkit.hh:73
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:45
JSON::TCAL
static const std::string TCAL
PMT time offsets.
Definition JSon/JSupport.hh:32
JSON::PCAL
static const std::string PCAL
(optical|base) module positions
Definition JSon/JSupport.hh:33
JSON::SCAL
static const std::string SCAL
(module|PMT) status
Definition JSon/JSupport.hh:37
JSON::RCAL
static const std::string RCAL
optical module orientations
Definition JSon/JSupport.hh:34
JSON::ACAL
static const std::string ACAL
acoustic time offsets (piezo sensor or hydrophone)
Definition JSon/JSupport.hh:35
JSON::CCAL
static const std::string CCAL
compass alignment (a.k.a. quaternion calibration)
Definition JSon/JSupport.hh:36
FILL
Auxiliary data structure for sequence of same character.
Definition JManip.hh:330
FIXED
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
JFormat_t
Data structure for format specifications.
Definition JManip.hh:524
JGEOMETRY3D::JQuaternion3D::decomposition
Auxiliary data structure for decomposition of quaternion in twist and swing quaternions.
Definition JQuaternion3D.hh:1127
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
JTOOLS::JStats
Auxiliary data structure for running average and standard deviation.
Definition JStats.hh:44
JTOOLS::JStats::put
void put(const double x, const double w=1.0)
Put value.
Definition JStats.hh:119
JTOOLS::JStats::getMean
double getMean() const
Get mean value.
Definition JStats.hh:234
JTOOLS::JStats::getSTDev
double getSTDev() const
Get standard deviation.
Definition JStats.hh:263
Generated by doxygen 1.10.0