JModule.hh

Go to the documentation of this file.

#ifndef __JDETECTOR__JMODULE__
#define __JDETECTOR__JMODULE__

#include <istream>
#include <ostream>
#include <vector>

#include "JDetector/JModuleIdentifier.hh"
#include "JDetector/JLocation.hh"
#include "JDetector/JPMT.hh"
#include "JDetector/JDetectorVersion.hh"
#include "Jeep/JStatus.hh"

#include "JGeometry3D/JPosition3D.hh"
#include "JGeometry3D/JRotation3D.hh"
#include "JGeometry3D/JQuaternion3D.hh"
#include "JGeometry3D/JTransformation3D.hh"

#include "JLang/JException.hh"
#include "JMath/JMatrix3S.hh"
#include "JMath/JMath.hh"

#include "JIO/JSerialisable.hh"


/**
 * \file
 *
 * Data structure for optical module.
 * \author mdejong
 */
namespace JDETECTOR {}
namespace JPP { using namespace JDETECTOR; }

namespace JDETECTOR {

  using JGEOMETRY3D::JPosition3D;
  using JGEOMETRY3D::JRotation3D;
  using JGEOMETRY3D::JQuaternion3D;
  using JGEOMETRY3D::JTransformation3D;
  using JGEOMETRY3D::JVector3D;
  using JGEOMETRY3D::JVersor3D;
  using JGEOMETRY3D::JAxis3D;
  using JEEP::JStatus;
  using JLANG::JValueOutOfRange;
  using JIO::JReader;
  using JIO::JWriter;


  /**
   * Data structure for a composite optical module.
   *
   * A module consists of a set of PMTs. A JPMT object is used for each PMT.\n
   * The index of the PMT in the module corresponds to the readout channel (TDC).\n
   * The quaternion data and time offset correspond to the calibration of the compass and piezo sensor inside the module, respectively.\n
   * There are no PMTs and piezo sensor in the base module. The time offset then corresponds to the hydrophone.
   * Note that the positions of the PMTs are absolute in space (i.e.\ not relative to the position of the module).
   *
   * The I/O of the position, quaternion data and time offset of the module depends on the detector version.\n
   * The member method JModule::compile is used to set the position, quaternion data and time offset
   * for detector versions for which these are not defined.\n
   * In this, the position of the module is determined from the intersection point of the PMT axes.
   *
   * Note finally that the positions of the reference modules may not exactly be at the origin 
   * due to the finite accuracy of the PMT axes.\n
   */
  class JModule :
    public JModuleIdentifier,
    public JLocation,
    public JPosition3D,
    public JQuaternion3D,
    public JCalibration,
    public JStatus,
    public std::vector<JPMT>
  {
  public:

    using JStatus::has;
    using JStatus::set;
    using JStatus::reset;

    /**
     * Default constructor.
     */
    JModule() :
      JModuleIdentifier(),
      JLocation(),
      JPosition3D(),
      JQuaternion3D(),
      JCalibration(),
      JStatus(),
      std::vector<JPMT>()
    {}


    /**
     * Constructor.
     *
     * \param  id         identifier
     * \param  location   location
     */
    JModule(const int        id,
            const JLocation& location) :
      JModuleIdentifier(id),
      JLocation(location),
      JPosition3D(),
      JQuaternion3D(),
      JCalibration(),
      JStatus(),
      std::vector<JPMT>()
    {}


    /**
     * Get detector version.
     */
    static JDetectorVersion& getVersion()
    {
      static JDetectorVersion version;

      return version;
    }


    /**
     * Set detector version.
     *
     * \param  version    version
     */
    static void setVersion(const JVersion& version)
    {
      getVersion() = JDetectorVersion(version);
    }


    /**
     * Compare modules.
     *
     * The comparison only covers the orientations of the modules.
     *
     * \param  first          first  module
     * \param  second         second module
     * \param  precision      precision
     * \return                true if two modules are equal; else false
     */
    static inline bool compare(const JModule& first,
                               const JModule& second,
                               const double   precision = 1.0e-3)
    {
      if (first.size() == second.size()) {

        for (size_t i = 0; i != first.size(); ++i) {
          if (first[i].getDirection().getDot(second[i].getDirection()) < 1.0 - precision) {
            return false;
          }
        }

        return true;
      }

      return false;
    }



    /**
     * Get PMT.
     *
     * \param  index      readout channel (TDC)
     * \return            PMT at given index
     */
    const JPMT& getPMT(const int index) const
    {
      return at(index);
    }


    /**
     * Get PMT.
     *
     * \param  index      readout channel (TDC)
     * \return            PMT at given index
     */
    JPMT& getPMT(const int index)
    {
      return at(index);
    }


    /**
     * Set PMT.
     *
     * \param  index      readout channel (TDC)
     * \param  pmt        PMT
     */
    void setPMT(const int index, const JPMT& pmt)
    {
      if (index >= (int) size()) {
        resize(index + 1);
      }

      (*this)[index] = pmt;
    }


    /**
     * Get center of module based on crossing point of PMT axes.
     *
     * This method perform a fit of the crossing point of the PMT axes.\n
     * A general purpose implementation of such a fit is available in JFIT::JEstimator<JPoint3D>.
     *
     * \return            center
     */
    JVector3D getCenter() const
    {
      using namespace JPP;

      if (this->size() > 1u) {

        double x = 0;
        double y = 0;
        double z = 0;

        JMatrix3S V;

        for (const_iterator pmt = this->begin(); pmt != this->end(); ++pmt) {

          const double xx = 1.0  -  pmt->getDX() * pmt->getDX();
          const double yy = 1.0  -  pmt->getDY() * pmt->getDY();
          const double zz = 1.0  -  pmt->getDZ() * pmt->getDZ();

          const double xy = -pmt->getDX() * pmt->getDY();
          const double xz = -pmt->getDX() * pmt->getDZ();
          const double yz = -pmt->getDY() * pmt->getDZ();

          V.a00 += xx;
          V.a01 += xy;
          V.a02 += xz;

          V.a11 += yy;
          V.a12 += yz;

          V.a22 += zz;

          x  +=  xx * pmt->getX()  + xy * pmt->getY()  +  xz * pmt->getZ();
          y  +=  xy * pmt->getX()  + yy * pmt->getY()  +  yz * pmt->getZ();
          z  +=  xz * pmt->getX()  + yz * pmt->getY()  +  zz * pmt->getZ();
        }

        V.a10 = V.a01;
        V.a20 = V.a02;
        V.a21 = V.a12;

        V.invert();

        return JVector3D(V.a00 * x  +  V.a01 * y  +  V.a02 * z,
                         V.a10 * x  +  V.a11 * y  +  V.a12 * z,
                         V.a20 * x  +  V.a21 * y  +  V.a22 * z);

      } else {
        throw JValueOutOfRange("JModule::getCenter(): Not enough PMTs.");
      }
    }


    /**
     * Compile module data.
     *
     * For detector versions before JDetectorVersion::V4, 
     * the position,
     * quaternion data and
     * time offset of the module should be set.\n
     * The position is set to the intersection point of the PMT axes (or their average position if this is not possible).\n
     * The quaternion data are maintained.\n
     * For an optical module (i.e.\ floor > 0),
     * the time offset is set to the average time offset of the PMTs and
     * for a base module (i.e.\ floor = 0),
     * the time offset is set to zero.\n
     * These time offsets should be corrected for delay in the piezo sensor and hydrophone, respectively.
     */
    void compile()
    {
      using namespace std;
      using namespace JPP;

      if (!this->empty()) {

        JPosition3D& pos = this->getPosition();

        try {
          pos = this->getCenter();
        }
        catch(const exception&) {

          pos = JPosition3D(0.0, 0.0, 0.0);

          for (iterator i = this->begin(); i != this->end(); ++i) {
            pos.add(i->getPosition());
          }

          pos.div(size());
        }
      }

      this->setCalibration(getAverage(make_array(this->begin(), this->end(), &JModule::getT0), 0.0));
    }


    /**
     * Rotate module.
     *
     * \param  R          rotation matrix
     */
    void rotate(const JRotation3D& R)
    {
      JPosition3D::rotate(R);

      for (iterator i = this->begin(); i != this->end(); ++i) {
        i->rotate(R);
      }
    }


    /**
     * Rotate back module.
     *
     * \param  R          rotation matrix
     */
    void rotate_back(const JRotation3D& R)
    {
      JPosition3D::rotate_back(R);

      for (iterator i = this->begin(); i != this->end(); ++i) {
        i->rotate_back(R);
      }
    }


    /**
     * Transformation of geometry (see method JGEOMETRY3D::JPosition3D::transform(const JRotation3D&, const JVector3D&)).
     *
     * \param  R          rotation matrix
     * \param  pos        position of origin (after rotation) 
     */
    void transform(const JRotation3D& R,
                   const JVector3D&   pos)
    {
      JPosition3D::transform(R, pos);

      for (iterator i = this->begin(); i != this->end(); ++i) {
        i->transform(R, pos);
      }
    }


    /**
     * Transformation of geometry.
     *
     * \param  T          transformation
     */
    void transform(const JTransformation3D& T)
    {
      JPosition3D::transform(T.getRotation(), T.getPosition());

      for (iterator i = this->begin(); i != this->end(); ++i) {
        i->transform(T);
      }
    }


    /**
     * Rotate module.
     *
     * \param  Q          quaternion
     */
    void rotate(const JQuaternion3D& Q)
    {
      JPosition3D::rotate(Q);

      for (iterator i = this->begin(); i != this->end(); ++i) {
        i->rotate(Q);
      }
    }


    /**
     * Rotate back module.
     *
     * \param  Q          quaternion
     */
    void rotate_back(const JQuaternion3D& Q)
    {
      JPosition3D::rotate_back(Q);

      for (iterator i = this->begin(); i != this->end(); ++i) {
        i->rotate_back(Q);
      }
    }


    /**
     * Set position.
     *
     * \param  pos        position
     * \return            this module
     */
    JModule& set(const JVector3D& pos)
    {
      return add(pos - this->getPosition());
    }


    /**
     * Add position.
     *
     * \param  pos        position
     * \return            this module
     */
    JModule& add(const JVector3D& pos)
    {
      for (iterator i = begin(); i != end(); ++i) {
        i->add(pos);
      }

      JPosition3D::add(pos);

      return *this;
    }


    /**
     * Subtract position.
     *
     * \param  pos        position
     * \return            this module
     */
    JModule& sub(const JVector3D& pos)
    {
      for (iterator i = begin(); i != end(); ++i) {
        i->sub(pos);
      }

      JPosition3D::sub(pos);

      return *this;
    }


    /**
     * Set time offset.
     *
     * \param  t0         time offset [ns]
     * \return            this module
     */
    JModule& set(const double t0)
    {
      for (iterator i = begin(); i != end(); ++i) {
        i->setT0(t0);
      }

      return *this;
    }


    /**
     * Add time offset.
     *
     * \param  t0         time offset [ns]
     * \return            this module
     */
    JModule& add(const double t0)
    {
      for (iterator i = begin(); i != end(); ++i) {
        i->addT0(t0);
      }

      return *this;
    }


    /**
     * Subtract time offset.
     *
     * \param  t0         time offset [ns]
     * \return            this module
     */
    JModule& sub(const double t0)
    {
      for (iterator i = begin(); i != end(); ++i) {
        i->subT0(t0);
      }

      return *this;
    }


    /**
     * Add position.
     *
     * \param  pos        position
     * \return            this module
     */
    JModule& operator+=(const JVector3D& pos)
    {
      return this->add(pos);
    }


    /**
     * Subtract position.
     *
     * \param  pos        position
     * \return            this module
     */
    JModule& operator-=(const JVector3D& pos)
    {
      return this->sub(pos);
    }


    /**
     * Read module from input.
     *
     * \param  in         input stream
     * \param  module     module
     * \return            input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JModule& module)
    {
      module = JModule();

      in >> static_cast<JModuleIdentifier&>(module);
      in >> static_cast<JLocation&>        (module);

      if (getDetectorVersion(getVersion()) >= JDetectorVersion::V4) {
        in >> static_cast<JPosition3D&>  (module);
        in >> static_cast<JQuaternion3D&>(module);
        in >> static_cast<JCalibration&> (module);
      }

      if (getDetectorVersion(getVersion()) >= JDetectorVersion::V5) {
        in >> static_cast<JStatus&>      (module);
      }

      unsigned int n;

      in >> n;

      for (JPMT pmt; n != 0 && in >> pmt; --n) {
        module.push_back(pmt);
      }

      if (getDetectorVersion(getVersion()) <  JDetectorVersion::V4) {
        module.compile();
      }

      return in;
    }


    /**
     * Write module to output.
     *
     * \param  out        output stream
     * \param  module     module
     * \return            output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JModule& module)
    {
      using namespace std;

      out << setw(10);
      out << static_cast<const JModuleIdentifier&>(module);
      out << ' ';
      out << static_cast<const JLocation&>        (module);

      if (getDetectorVersion(getVersion()) >= JDetectorVersion::V4) {
        out << ' ';
        out << static_cast<const JPosition3D&>  (module);
        out << ' ';
        out << static_cast<const JQuaternion3D&>(module);
        out << ' ';
        out << static_cast<const JCalibration&> (module);
      }

      if (getDetectorVersion(getVersion()) >= JDetectorVersion::V5) {
        out << ' ';
        out << static_cast<const JStatus&>      (module);
      }

      out << ' ' << module.size() << endl;

      for (const_iterator i = module.begin(); i != module.end(); ++i) {
        out << ' ' << *i << endl;;
      }

      return out;
    }


    /**
     * Read module from input.
     *
     * \param  in       reader
     * \param  module   module
     * \return          rrreader
     */
    friend inline JReader& operator>>(JReader& in, JModule& module)
    {
      module = JModule();

      in >> static_cast<JModuleIdentifier&>(module);
      in >> static_cast<JLocation&>        (module);

      if (getDetectorVersion(getVersion()) >= JDetectorVersion::V4) {
        in >> static_cast<JPosition3D&>  (module);
        in >> static_cast<JQuaternion3D&>(module);
        in >> static_cast<JCalibration&> (module);
      }

      if (getDetectorVersion(getVersion()) >= JDetectorVersion::V5) {
        in >> static_cast<JStatus&>      (module);
      }

      int n;

      in >> n;

      module.resize(n);

      for (JModule::iterator out = module.begin(); n != 0; --n, ++out) {
        in >> *out;
      }

      if (getDetectorVersion(getVersion()) <  JDetectorVersion::V4) {
        module.compile();
      }

      return in;
    }


    /**
     * Write module to output.
     *
     * \param  out      writer
     * \param  module   module
     * \return          writer
     */
    friend inline JWriter& operator<<(JWriter& out, const JModule& module)
    {
      out << static_cast<const JModuleIdentifier&>(module);
      out << static_cast<const JLocation&>        (module);

      if (getDetectorVersion(getVersion()) >= JDetectorVersion::V4) {
        out << static_cast<const JPosition3D&>  (module);
        out << static_cast<const JQuaternion3D&>(module);
        out << static_cast<const JCalibration&> (module);
      }

      if (getDetectorVersion(getVersion()) >= JDetectorVersion::V5) {
        out << static_cast<const JStatus&>      (module);
      }

      int n = module.size();

      out << n;

      for (const_iterator i = module.begin(); i != module.end(); ++i) {
        out << *i;
      }

      return out;
    }
  };
}

#endif