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 "JGeometry3D/JPosition3D.hh"
#include "JGeometry3D/JRotation3D.hh"
#include "JGeometry3D/JQuaternion3D.hh"
#include "JGeometry3D/JTransformation3D.hh"
#include "JFit/JPoint3DEstimator.hh"
#include "JIO/JSerialisable.hh"


/**
 * \file
 *
 * Data structure for a composite 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 JIO::JReader;
  using JIO::JWriter;


  /**
   * Data structure for a composite optical module.
   *
   * A module consists of a set of JPMT objects.
   * The index of the PMT object in the module container corresponds to the readout channel (TDC).
   * The positions of the PMTs are absolute in space (i.e.\ not relative to the position of the module).
   * The position of the module is not subject to I/O but derived from the positions of the PMTs
   * using method compile().
   */
  class JModule :
    public JModuleIdentifier,
    public JLocation,
    public JPosition3D,
    public std::vector<JPMT>
  {
  public:
    /**
     * Default constructor.
     */
    JModule() :
      JModuleIdentifier(),
      JLocation(),
      JPosition3D(),
      std::vector<JPMT>()
    {}


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


    /**
     * Get reference to unique instance of this class object.
     *
     * This method returns a module with the standard configuration of PMTs.
     *
     * \return                 reference to this class object
     */
    static const JModule& getInstance()
    {
      static JModule module;

      if (module.empty()) {

        module.push_back(JPMT( 1, JAxis3D(JVector3D(+0.000, +0.000, -0.200), JVersor3D(+0.000, +0.000, -1.000))));

        module.push_back(JPMT( 2, JAxis3D(JVector3D(+0.000, +0.105, -0.170), JVersor3D(+0.000, +0.527, -0.850))));
        module.push_back(JPMT( 3, JAxis3D(JVector3D(+0.091, +0.053, -0.170), JVersor3D(+0.456, +0.263, -0.850))));
        module.push_back(JPMT( 4, JAxis3D(JVector3D(+0.091, -0.053, -0.170), JVersor3D(+0.456, -0.263, -0.850))));
        module.push_back(JPMT( 5, JAxis3D(JVector3D(+0.000, -0.105, -0.170), JVersor3D(+0.000, -0.527, -0.850))));
        module.push_back(JPMT( 6, JAxis3D(JVector3D(-0.091, -0.053, -0.170), JVersor3D(-0.456, -0.263, -0.850))));
        module.push_back(JPMT( 7, JAxis3D(JVector3D(-0.091, +0.053, -0.170), JVersor3D(-0.456, +0.263, -0.850))));

        module.push_back(JPMT( 8, JAxis3D(JVector3D(+0.083, +0.144, -0.111), JVersor3D(+0.416, +0.720, -0.555))));
        module.push_back(JPMT( 9, JAxis3D(JVector3D(+0.166, +0.000, -0.111), JVersor3D(+0.832, +0.000, -0.555))));
        module.push_back(JPMT(10, JAxis3D(JVector3D(+0.083, -0.144, -0.111), JVersor3D(+0.416, -0.720, -0.555))));
        module.push_back(JPMT(11, JAxis3D(JVector3D(-0.083, -0.144, -0.111), JVersor3D(-0.416, -0.720, -0.555))));
        module.push_back(JPMT(12, JAxis3D(JVector3D(-0.166, +0.000, -0.111), JVersor3D(-0.832, +0.000, -0.555))));
        module.push_back(JPMT(13, JAxis3D(JVector3D(-0.083, +0.144, -0.111), JVersor3D(-0.416, +0.720, -0.555))));

        module.push_back(JPMT(14, JAxis3D(JVector3D(+0.000, +0.191, -0.059), JVersor3D(+0.000, +0.955, -0.295))));
        module.push_back(JPMT(15, JAxis3D(JVector3D(+0.165, +0.096, -0.059), JVersor3D(+0.827, +0.478, -0.295))));
        module.push_back(JPMT(16, JAxis3D(JVector3D(+0.165, -0.096, -0.059), JVersor3D(+0.827, -0.478, -0.295))));
        module.push_back(JPMT(17, JAxis3D(JVector3D(+0.000, -0.191, -0.059), JVersor3D(+0.000, -0.955, -0.295))));
        module.push_back(JPMT(18, JAxis3D(JVector3D(-0.165, -0.096, -0.059), JVersor3D(-0.827, -0.478, -0.295))));
        module.push_back(JPMT(19, JAxis3D(JVector3D(-0.165, +0.096, -0.059), JVersor3D(-0.827, +0.478, -0.295))));

        module.push_back(JPMT(20, JAxis3D(JVector3D(+0.096, +0.165, +0.059), JVersor3D(+0.478, +0.827, +0.295))));
        module.push_back(JPMT(21, JAxis3D(JVector3D(+0.191, +0.000, +0.059), JVersor3D(+0.955, +0.000, +0.295))));
        module.push_back(JPMT(22, JAxis3D(JVector3D(+0.096, -0.165, +0.059), JVersor3D(+0.478, -0.827, +0.295))));
        module.push_back(JPMT(23, JAxis3D(JVector3D(-0.096, -0.165, +0.059), JVersor3D(-0.478, -0.827, +0.295))));
        module.push_back(JPMT(24, JAxis3D(JVector3D(-0.191, +0.000, +0.059), JVersor3D(-0.955, +0.000, +0.295))));
        module.push_back(JPMT(25, JAxis3D(JVector3D(-0.096, +0.165, +0.059), JVersor3D(-0.478, +0.827, +0.295))));

        module.push_back(JPMT(26, JAxis3D(JVector3D(+0.000, +0.166, +0.111), JVersor3D(+0.000, +0.832, +0.555))));
        module.push_back(JPMT(27, JAxis3D(JVector3D(+0.144, +0.083, +0.111), JVersor3D(+0.720, +0.416, +0.555))));
        module.push_back(JPMT(28, JAxis3D(JVector3D(+0.144, -0.083, +0.111), JVersor3D(+0.720, -0.416, +0.555))));
        module.push_back(JPMT(29, JAxis3D(JVector3D(+0.000, -0.166, +0.111), JVersor3D(+0.000, -0.832, +0.555))));
        module.push_back(JPMT(30, JAxis3D(JVector3D(-0.144, -0.083, +0.111), JVersor3D(-0.720, -0.416, +0.555))));
        module.push_back(JPMT(31, JAxis3D(JVector3D(-0.144, +0.083, +0.111), JVersor3D(-0.720, +0.416, +0.555))));
      }

      return module;
    }


    /**
     * 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;
    }


    /**
     * Compile position of module from the positions and directions of the PMTs.
     */
    void compile()
    {
      using namespace std;
      using namespace JPP;

      if (!empty()) {

        JPosition3D& pos = getPosition();

        try {
          pos = JEstimator<JPoint3D>(this->begin(), this->end());
        }
        catch(const exception&) {

          pos = JPosition3D(0.0, 0.0, 0.0);

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

          pos.div(size());
        }
      }
    }


    /**
     * 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 - static_cast<JPosition3D&>(*this));
    }


    /**
     * 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.clear();

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

      unsigned int n;

      in >> n;

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

      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 << static_cast<const JModuleIdentifier&>(module);
      out << ' ';
      out << static_cast<const JLocation&>        (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.clear();

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

      int n;

      in >> n;

      for (JPMT pmt; n != 0; --n) {

        in >> pmt;

        module.push_back(pmt);
      }

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

      int n = module.size();

      out << n;

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

      return out;
    }
  };
}

#endif