74 string detectorFile_a;
75 string detectorFile_b;
82 JParser<> zap(
"Auxiliary program to find differences between two detector files.");
92 catch(
const exception &error) {
93 FATAL(error.what() << endl);
101 load(detectorFile_a, detector_a);
108 load(detectorFile_b, detector_b);
115 bool is_equal =
true;
120 setFormat<JPosition3D> (
JFormat_t(15, 9, std::ios::fixed | std::ios::showpos));
124 if (detector_a.getID() != detector_b.getID()) {
126 DEBUG(
"* Different detector identifiers "
127 << setw(5) << detector_a.getID() <<
" (A) and " << endl
128 << setw(5) << detector_b.getID() <<
" (B)." << endl
137 for (JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
139 if (!module_router_b.hasModule(module->getID())) {
141 DEBUG(
"* Module " << setw(10) << module->getID() <<
" is in A " <<
getLabel(*module) <<
" but not in B" << endl);
147 for (JDetector::iterator module = detector_b.begin(); module != detector_b.end(); ++module) {
149 if (!module_router_a.hasModule(module->getID())) {
151 DEBUG(
"* Module " << setw(10) << module->getID() <<
" is in B " <<
getLabel(*module) <<
" but not in A" << endl);
161 DEBUG(
"Comparing module by module." << endl);
163 for (JDetector::iterator module_a = detector_a.begin(); module_a != detector_a.end(); ++module_a) {
165 if (!module_router_b.hasModule(module_a->getID())) {
172 const JModule* module_b = &module_router_b.getModule(module_a->getID());
174 DEBUG(
" Module " << setw(10) << module_a->getID());
178 if (module_a->getLocation() == module_b->
getLocation()) {
185 DEBUG(
" * different location: "
187 <<
getLabel(*module_b) <<
" (B)" << endl);
194 if (fabs(module_a->getT0() - module_b->
getT0()) > precision) {
196 DEBUG(
" * different T0: "
197 << module_a->getT0() <<
" (A), "
198 << module_b->
getT0() <<
" (B) "
199 <<
", B - A " << module_b->
getT0() - module_a->getT0() << endl);
208 DEBUG(
" * different quaternion calibration: "
209 << module_a->getQuaternion() <<
" (A), "
220 DEBUG(
" * different position: "
221 << module_a->getPosition() <<
" (A), "
230 if (module_a->size() != module_b->size()) {
232 DEBUG(
" * different number of PMTs: "
233 << module_a->size() <<
" (A), "
234 << module_b->size() <<
" (B)" << endl);
241 if (!module_a->empty() &&
242 !module_b->empty()) {
244 const JQuantile q = getQuantile(*module_a, *module_b);
246 if (fabs(q.
getMean()) > precision) {
248 DEBUG(
" * different average t0: "
258 if (fabs(module_a->getT0() - module_b->
getT0()) > precision) {
260 DEBUG(
" * different global t0: "
261 << module_a->getT0() <<
" (A), "
262 << module_b->
getT0() <<
" (B)"
263 <<
", B - A " << module_b->
getT0() - module_a->getT0()
273 for (
unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
275 const JPMT& pmt_a = module_a->getPMT(pmt);
278 if (fabs(pmt_a.
getT0() - pmt_b.
getT0()) > precision) {
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()
292 for (
unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
294 const JPMT& pmt_a = module_a->getPMT(pmt);
301 DEBUG(
" * different PMT position: "
312 for (
unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
314 const JPMT& pmt_a = module_a->getPMT(pmt);
319 DEBUG(
" * different status PMT " << pmt <<
": "
338 for (JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
339 string.insert(module->getString());
340 floor .insert(module->getFloor ());
343 for (JDetector::iterator module = detector_b.begin(); module != detector_b.end(); ++module) {
344 string.insert(module->getString());
345 floor .insert(module->getFloor ());
350 string.size(), -0.5,
string.size() - 0.5,
351 floor .size(), -0.5, floor .size() - 0.5);
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" );
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);
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);
385 for (JDetector::iterator module_a = detector_a.begin(); module_a != detector_a.end(); ++module_a) {
387 if (!module_router_b.hasModule(module_a->getID())) {
391 const JModule* module_b = &module_router_b.getModule(module_a->getID());
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());
397 if (module_a->getFloor() != 0 &&
402 for (
unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
406 while (x > +
PI) { x -=
PI; }
407 while (x < -
PI) { x +=
PI; }
412 phi /= min(module_a->size(),
416 const JQuantile q = getQuantile(*module_a, *module_b);
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());
Utility class to parse command line options.
double getAngle(const JQuaternion3D &first, const JQuaternion3D &second)
Get space angle between quanternions.
const JStatus & getStatus() const
Get status.
int getFloor() const
Get floor number.
Data structure for a composite optical module.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
std::string getLabel(const JLocation &location)
Get module label for monitoring and other applications.
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.
double getPhi() const
Get phi angle.
double getDistance(const JVector3D &pos) const
Get distance to point.
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.
const JLocation & getLocation() const
Get location.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Data structure for PMT geometry and calibration.
static const double PI
Mathematical constants.
double getY() const
Get y position.
const JPosition3D & getPosition() const
Get position.
const JPMT & getPMT(const int index) const
Get PMT.
Auxiliary data structure for sequence of same character.
Data structure for unit quaternion in three dimensions.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
double getX() const
Get x position.
Data structure for position in three dimensions.
double getZ() const
Get z position.
#define DEBUG(A)
Message macros.
double getT0() const
Get time offset.