JCompareDetector.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
 
#include "JGeometry3D/JQuaternion3D.hh"
#include "JMath/JConstants.hh"
#include "JTools/JRange.hh"
#include "JTools/JQuantile.hh"
#include "JROOT/JManager.hh"
#include "JROOT/JRootToolkit.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JDetectorCalibration.hh"
#include "JDetector/JModuleRouter.hh"
 
#include "Jeep/JProperties.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
namespace {
 
  using namespace JPP;
 
  /**
   * Compare time calibration of PMTs between two optical modules.
   *
   * \param  module_a        module
   * \param  module_b        module
   * \return                 quantile
   */
  inline JQuantile getQuantile(const JModule& module_a, const JModule& module_b)
  {
    JQuantile Q;
 
    JModule::const_iterator pmt_a = module_a.begin();
    JModule::const_iterator pmt_b = module_b.begin();
 
    for ( ; pmt_a != module_a.end() && pmt_b != module_b.end(); ++pmt_a, ++pmt_b) {
      Q.put(pmt_a->getT0() - pmt_b->getT0());
    }      
 
    return Q;
  }
 
 
  /**
   * Get ROOT histogram abscissa value.
   *
   * \param  buffer          buffer
   * \param  value           value
   * \return                 value
   */
  inline Double_t getBin(const std::set<int>& buffer, const int value)
  {
    return std::distance(buffer.begin(), buffer.find(value));
  }
}
 
 
/**
 * \file
 *
 * Auxiliary program to find differences between two detector files.\n
 * A ROOT output file with histograms is optionally produced.
 * \author mjongen
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  struct {
    double tcal  =  0.001;      // [ns]
    double pcal  =  0.001;      // [m]
    double rcal  =  0.001;      // [rad]
    double acal  =  0.001;      // [ns]
    double ccal  =  0.001;      // [deg]
    int    scal  =  0xFFFFFFFF; // [status]
  } precision;
 
  string  detectorFile_a;
  string  detectorFile_b;
  string  outputFile;
  int     debug;
 
  try {
 
    JProperties properties;
 
    properties[TCAL] = precision.tcal;
    properties[PCAL] = precision.pcal;
    properties[RCAL] = precision.rcal;
    properties[ACAL] = precision.acal;
    properties[CCAL] = precision.ccal;
    properties[SCAL] = precision.scal;
 
    JParser<> zap("Auxiliary program to find differences between two detector files.");
 
    zap['a'] = make_field(detectorFile_a);
    zap['b'] = make_field(detectorFile_b);
    zap['o'] = make_field(outputFile)            = "";
    zap['p'] = make_field(properties)            = JPARSER::initialised();
    zap['d'] = make_field(debug)                 = 3;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
  
  JDetector detector_a;
  JDetector detector_b;
 
  try {
    load(detectorFile_a, detector_a);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  try {
    load(detectorFile_b, detector_b);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  size_t numberOfPMTs = 0;
 
  bool is_equal = true;
 
  const JModuleRouter module_router_a(detector_a);
  const JModuleRouter module_router_b(detector_b);
 
  setFormat<JPosition3D>  (JFormat_t(15, 9, std::ios::fixed | std::ios::showpos));
 
  // compare detector identifiers
 
  if (detector_a.getID() != detector_b.getID()) {
 
    DEBUG("* Different detector identifiers " 
          << setw(5) << detector_a.getID() << " (A) and " << endl 
          << setw(5) << detector_b.getID() << " (B)."     << endl
          << endl);
 
    is_equal = false;
  }
 
  // UTM positions
 
  if (getDistance(detector_a.getPosition(), detector_b.getPosition()) > precision.pcal) {
 
    DEBUG("  * different UTM position: "
          << detector_a.getPosition() << " (A), " 
          << detector_b.getPosition() << " (B)"
          << ", B - A " << JPosition3D(detector_b.getPosition() - detector_a.getPosition()) << endl);
    
    is_equal = false;
  }
 
 
  for (JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
    if (module->size() > numberOfPMTs) {
       numberOfPMTs = module->size();
    }
  }
 
  // check whether the same modules are present in the file
 
  for (JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
 
    if (!module_router_b.hasModule(module->getID())) {
 
      DEBUG("* Module " << setw(10) << module->getID() << " is in A " << getLabel(*module) << " but not in B" << endl);
 
      is_equal = false;
    }
  }
 
  for (JDetector::iterator module = detector_b.begin(); module != detector_b.end(); ++module) {
 
    if (!module_router_a.hasModule(module->getID())) {
 
      DEBUG("* Module " << setw(10) << module->getID() << " is in B " << getLabel(*module) << " but not in A" << endl);
 
      is_equal = false;
    }
  }
  DEBUG(endl);
 
 
  // compare the modules that the files have in common
 
  DEBUG("Comparing module by module." << endl);
 
  for (JDetector::iterator module_a = detector_a.begin(); module_a != detector_a.end(); ++module_a) {
 
    if (!module_router_b.hasModule(module_a->getID())) {
 
      continue;
 
      is_equal = false;
    }
    
    const JModule* module_b = &module_router_b.getModule(module_a->getID());
    
    DEBUG("  Module " << setw(10) << module_a->getID());
           
    // string and floor numbers
 
    if (module_a->getLocation() == module_b->getLocation()) {
 
      DEBUG(" " << getLabel(*module_a) << endl);
 
    } else {
 
      DEBUG(endl);
      DEBUG("  * different location: "
            << getLabel(*module_a) << " (A), "
            << getLabel(*module_b) << " (B)" << endl);
 
      is_equal = false;
    }
 
    // time offset
 
    if (fabs(module_a->getT0() - module_b->getT0()) > precision.acal) {
 
      DEBUG("  * different T0: " 
            << FIXED(12,3) << module_a->getT0() << " (A), "
            << FIXED(12,3) << module_b->getT0() << " (B) "
            << ", B - A " << module_b->getT0() - module_a->getT0() << endl);
        
      is_equal = false;
    }
 
    // quaternion calibration
 
    if (module_a->getQuaternion() != JQuaternion3D(0.0, 0.0, 0.0, 0.0) &&
        module_b->getQuaternion() != JQuaternion3D(0.0, 0.0, 0.0, 0.0) &&
        getAngle(module_a->getQuaternion(), module_b->getQuaternion()) > precision.ccal) {
 
      DEBUG("  * different quaternion calibration: " 
            << module_a->getQuaternion() << " (A), "
            << module_b->getQuaternion() << " (B) "
            << ", B - A " << getAngle(module_b->getQuaternion(), module_a->getQuaternion()) << endl);
        
      is_equal = false;
    }
 
    // average DOM positions
 
    if (getDistance(module_a->getPosition(), module_b->getPosition()) > precision.pcal) {
 
      DEBUG("  * different position: "
            << module_a->getPosition() << " (A), " 
            << module_b->getPosition() << " (B)"
            << ", B - A " << JPosition3D(module_b->getPosition() - module_a->getPosition()) << endl);
 
      is_equal = false;
    }
 
    // number of PMTs
 
    if (module_a->size() != module_b->size()) {
 
      DEBUG("  * different number of PMTs: "
            << module_a->size() << " (A), " 
            << module_b->size() << " (B)" << endl);
 
      is_equal = false;
    }
           
    // average t0
 
    if (!module_a->empty() &&
        !module_b->empty()) { 
 
      const JQuantile q = getQuantile(*module_a, *module_b);
 
      if (fabs(q.getMean()) > precision.tcal) {
 
        DEBUG("  * different average t0: "
              << ", B - A " << q.getMean()
              << endl);
      
        is_equal = false;
      }
    }
 
    // status
 
    if (module_a->getStatus(precision.scal) != module_b->getStatus(precision.scal)) {
 
      DEBUG("  * different status module " << module_a->getID() << ": "
            << module_a->getStatus() << " (A), " 
            << module_b->getStatus() << " (B)"   
            << endl);
        
      is_equal = false;
    }
 
    // comparing by PMT
 
    // identifier
 
    for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
 
      const JPMT& pmt_a = module_a->getPMT(pmt);
      const JPMT& pmt_b = module_b->getPMT(pmt);
 
      if (pmt_a.getID() != pmt_b.getID()) {
 
        DEBUG("  * different identifier PMT " << setw(2) << pmt << ": "
              << setw(8) << pmt_a.getID() << " (A" << FILL(2,'0') << pmt << "), " << FILL()
              << setw(8) << pmt_b.getID() << " (B" << FILL(2,'0') << pmt << ")"   << FILL()
              << ", B - A " << pmt_b.getID() - pmt_a.getID()  
              << endl);
        
        is_equal = false;
      }
    }
    
    // t0
 
    for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
 
      const JPMT& pmt_a = module_a->getPMT(pmt);
      const JPMT& pmt_b = module_b->getPMT(pmt);
 
      if (fabs(pmt_a.getT0() - pmt_b.getT0()) > precision.tcal) {
 
        DEBUG("  * different T0 PMT " << setw(2) << pmt << ": "
              << FIXED(12,3) << pmt_a.getT0() << " (A" << FILL(2,'0') << pmt << "), " << FILL()
              << FIXED(12,3) << pmt_b.getT0() << " (B" << FILL(2,'0') << pmt << ")"   << FILL()
              << ", B - A " << pmt_b.getT0() - pmt_a.getT0()  
              << endl);
        
        is_equal = false;
      }
    }
 
    // positions
 
    for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
 
      const JPMT& pmt_a = module_a->getPMT(pmt);
      const JPMT& pmt_b = module_b->getPMT(pmt);
 
      // PMT positions
 
      if (pmt_a.getPosition().getDistance(pmt_b.getPosition()) > precision.pcal) {
 
        DEBUG("  * different PMT position: "
              << pmt_a.getPosition() << " (A" << FILL(2,'0') << pmt << "), " << FILL()
              << pmt_b.getPosition() << " (B" << FILL(2,'0') << pmt << ")"   << FILL()
              << ", B - A " << JPosition3D(pmt_b.getPosition() - pmt_a.getPosition()) << endl);
        
        is_equal = false;
      }
    }
    
    // orientations
 
    for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
 
      const JPMT& pmt_a = module_a->getPMT(pmt);
      const JPMT& pmt_b = module_b->getPMT(pmt);
 
      // PMT positions
 
      const double dot = pmt_a.getDirection().getDot(pmt_b.getDirection());
      
      if ((1.0 - dot) > precision.rcal) {
 
        DEBUG("  * different PMT direction: "
              << pmt_a.getDirection() << " (A" << FILL(2,'0') << pmt << "), " << FILL()
              << pmt_b.getDirection() << " (B" << FILL(2,'0') << pmt << ")"   << FILL()
              << ", B - A " << JPosition3D(JPosition3D(pmt_b.getDirection()) - JPosition3D(pmt_a.getDirection())) << ' ' << dot << endl);
        
        is_equal = false;
      }
    }
 
    // status
 
    for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
 
      const JPMT& pmt_a = module_a->getPMT(pmt);
      const JPMT& pmt_b = module_b->getPMT(pmt);
 
      if (pmt_a.getStatus(precision.scal) != pmt_b.getStatus(precision.scal)) {
 
        DEBUG("  * different status PMT " << setw(2) << pmt << ": "
              << pmt_a.getStatus() << " (A" << FILL(2,'0') << pmt << "), " << FILL()
              << pmt_b.getStatus() << " (B" << FILL(2,'0') << pmt << ")"   << FILL()
              << endl);
        
        is_equal = false;
      }
    }
  }    
  DEBUG(endl);
 
 
  if (outputFile != "") {
 
    set<int> string;
    set<int> floor;
    
    for (JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
      string.insert(module->getString());
      floor .insert(module->getFloor ());
    }
    
    for (JDetector::iterator module = detector_b.begin(); module != detector_b.end(); ++module) {
      string.insert(module->getString());
      floor .insert(module->getFloor ());
    }
 
    JManager<int, TH1D> H1(new TH1D("M[%]", NULL, numberOfPMTs, -0.5, numberOfPMTs - 0.5));
 
    TH2D M2("M2", NULL,
            string.size(), -0.5, string.size() - 0.5,
            floor .size(), -0.5, floor .size() - 0.5);
 
    for (set<int>::const_iterator i = string.begin(); i != string.end(); ++i) {
      M2.GetXaxis()->SetBinLabel(distance(string.begin(), i) + 1, MAKE_CSTRING(*i));
    }
 
    for (set<int>::const_iterator i = floor.begin();  i != floor.end();  ++i) {
      M2.GetYaxis()->SetBinLabel(distance(floor.begin(), i) + 1,  MAKE_CSTRING(*i));
    }
 
    TH2D* X2  = (TH2D*) M2.Clone("X2" );
    TH2D* Y2  = (TH2D*) M2.Clone("Y2" );
    TH2D* Z2  = (TH2D*) M2.Clone("Z2" );
    TH2D* T2  = (TH2D*) M2.Clone("T2" );
    TH2D* RMS = (TH2D*) M2.Clone("RMS");
    TH2D* R2  = (TH2D*) M2.Clone("R2" );
    TH2D* QA  = (TH2D*) M2.Clone("QA" );
    TH2D* QB  = (TH2D*) M2.Clone("QB" );
    TH2D* QC  = (TH2D*) M2.Clone("QC" );
    TH2D* QD  = (TH2D*) M2.Clone("QD" );
    TH2D* Q2  = (TH2D*) M2.Clone("Q2" );
 
    for (JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
      if (!module_router_b.hasModule(module->getID()) ) {
        M2.Fill(module->getString(), module->getFloor(), -1.0);
      }
    }
    
    for (JDetector::iterator module = detector_b.begin(); module != detector_b.end(); ++module) {
      if (!module_router_a.hasModule(module->getID()) ) {
        M2.Fill(module->getString(), module->getFloor(), +1.0);
      }
    }
 
 
    for (JDetector::iterator module_a = detector_a.begin(); module_a != detector_a.end(); ++module_a) {
 
      if (!module_router_b.hasModule(module_a->getID())) {
        continue;
      }
 
      const JModule* module_b = &module_router_b.getModule(module_a->getID());
 
      for (size_t i = 0; i != numberOfPMTs; ++i) {
 
        if (module_a->getFloor() != 0) {
          H1[module_a->getID()]->SetBinContent(i + 1, module_a->getPMT(i).getT0() - module_b->getPMT(i).getT0());
        }
      }
 
      X2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), module_a->getX() - module_b->getX() + numeric_limits<double>::min());
      Y2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), module_a->getY() - module_b->getY() + numeric_limits<double>::min());
      Z2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), module_a->getZ() - module_b->getZ() + numeric_limits<double>::min());
 
      if (module_a->getFloor() != 0 &&
          module_b->getFloor() != 0) {
 
        const JQuaternion3D Q = getRotation(*module_a, *module_b);
        const JQuantile     q = getQuantile(*module_a, *module_b);
 
        const JQuaternion3D::decomposition q1(Q, JVector3Z_t);
 
        const double phi = (JVector3Z_t.getDot(q1.twist) >= 0.0 ? +1.0 : -1.0) * q1.twist.getAngle();
 
        R2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), phi);
        QA ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), Q.getA());
        QB ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), Q.getB());
        QC ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), Q.getC());
        QD ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), Q.getD());
        Q2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), Q.getAngle());
        T2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), q.getMean());
        RMS->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), q.getSTDev());
      }
    }
 
 
    TFile out(outputFile.c_str(), "recreate");    
 
    for (TH2D* h2 : { &M2, X2, Y2, Z2, T2, RMS, R2, QA, QB, QC, QD, Q2 }) {
      out << *h2;
    }
 
    out << H1;
 
    out.Write();
    out.Close();
  }
 
  ASSERT(is_equal);
 
  return 0;
}