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)...