Jpp 19.3.0-rc.1
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
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"
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 JQuantile getQuantile(const JModule& module_a, const JModule& module_b)
38 {
39 JQuantile 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_a->getT0() - pmt_b->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 JQuantile 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->getStatus(precision.scal) != module_b->getStatus(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.getStatus(precision.scal) != pmt_b.getStatus(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 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
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}
int main(int argc, char **argv)
string outputFile
Data structure for detector geometry and calibration.
Dynamic ROOT object management.
void setFormat(const JFormat_t &format)
Set format for given type.
Definition JManip.hh:714
Mathematical constants.
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define ASSERT(A,...)
Assert macro.
Definition JMessage.hh:90
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
Direct access to module in detector data structure.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
#define MAKE_CSTRING(A)
Make C-string.
Definition JPrint.hh:72
Utility class to parse parameter values.
Auxiliary class to define a range between two values.
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.
Detector data structure.
Definition JDetector.hh:96
const JLocation & getLocation() const
Get location.
Definition JLocation.hh:70
int getFloor() const
Get floor number.
Definition JLocation.hh:146
int getString() const
Get string number.
Definition JLocation.hh:135
Router for direct addressing of module data in detector data structure.
bool hasModule(const JObjectID &id) const
Has module.
const JModule & getModule(const JObjectID &id) const
Get module parameters.
Data structure for a composite optical module.
Definition JModule.hh:75
const JPMT & getPMT(const int index) const
Get PMT.
Definition JModule.hh:172
Data structure for PMT geometry, calibration and status.
Definition JPMT.hh:49
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.
Definition JVector3D.hh:104
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition JVector3D.hh:270
double getZ() const
Get z position.
Definition JVector3D.hh:115
double getDot(const JVector3D &vector) const
Get dot product.
Definition JVector3D.hh:282
double getX() const
Get x position.
Definition JVector3D.hh:94
General exception.
Definition JException.hh:24
int getID() const
Get identifier.
Definition JObjectID.hh:50
Utility class to parse command line options.
Definition JParser.hh:1698
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys.
Definition JManager.hh:47
void Write(TDirectory &out, const bool wm=false)
Write objects to file.
Definition JManager.hh:304
JPosition3D getPosition() const
Get position.
std::string getLabel(const JLocation &location)
Get module label for monitoring and other applications.
Definition JLocation.hh:247
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.
Definition JManip.hh:330
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
int getStatus() const
Get status.
Definition JStatus.hh:63
Data structure for format specifications.
Definition JManip.hh:524
Auxiliary data structure for decomposition of quaternion in twist and swing quaternions.
JQuaternion3D twist
rotation around parallel axis
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
Auxiliary data structure for running average, standard deviation and quantiles.
Definition JQuantile.hh:46
double getSTDev() const
Get standard deviation.
Definition JQuantile.hh:281
void put(const double x, const double w=1.0)
Put value.
Definition JQuantile.hh:133
double getMean() const
Get mean value.
Definition JQuantile.hh:252