JBallarat.hh

Go to the documentation of this file.

#ifndef __JCOMPASS__JBALLARAT__
#define __JCOMPASS__JBALLARAT__

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JDetectorAddressMapToolkit.hh"

#include "JCompass/JEvt.hh"
#include "JCompass/JEvtToolkit.hh"

#include "JLang/JObjectIterator.hh"
#include "JLang/JException.hh"
#include "JLang/JComparator.hh"

#include "JGeometry3D/JQuaternion3D.hh"
#include "JGeometry3D/JEigen3D.hh"

#include "JTools/JElement.hh"
#include "JTools/JCollection.hh"
#include "JTools/JPolfit.hh"
#include "JTools/JHashMap.hh"


namespace JCOMPASS {}
namespace JPP { using namespace JCOMPASS; }

namespace JCOMPASS {

  using JLANG::JObjectIterator;
  using JLANG::JValueOutOfRange;
  using JDETECTOR::JDetector;


  /**
   * Dynamic orientation calibration.
   *
   * Note that the quaternion data of the internal detector are used to track the consecutive rotations.
   */
  struct JBallarat {

    enum {
      NUMBER_OF_POINTS  = 20,    //!< number of points  for interpolation
      NUMBER_OF_DEGREES =  1     //!< number of degrees for interpolation
    };

    typedef JTOOLS::JElement2D<double, JGEOMETRY3D::JQuaternion3D>          element_type;
    typedef JTOOLS::JPolfitFunction1D<NUMBER_OF_POINTS,
                                      NUMBER_OF_DEGREES,
                                      element_type, JTOOLS::JCollection>    function_type;
    typedef typename function_type::collection_type::container_type         container_type;
  
    typedef JTOOLS::JHashMap<int, function_type>                            data_type;

    typedef data_type::const_iterator                                       const_iterator;
    typedef data_type::const_reverse_iterator                               const_reverse_iterator;


    /**
     * Constructor.
     *
     * \param  detector                 detector
     * \param  Tmax_s                   applicability period [s]
     */
    JBallarat(const JDetector& detector,
              const double     Tmax_s) :
      detector(detector),
      Tmax_s(Tmax_s),
      t0_s(std::numeric_limits<double>::lowest())
    {
      using namespace JPP;

      if (hasDetectorAddressMap(this->detector.getID())) {

        const JDetectorAddressMap& demo = getDetectorAddressMap(detector.getID());

        for (JDetector::iterator module = this->detector.begin(); module != this->detector.end(); ++module) {
          module->setQuaternion(getRotation(getModule(demo.get(module->getID())), *module));
        }

      } else {

        THROW(JValueOutOfRange, "No detector address map for detector identier " << detector.getID());
      }
    }


    /**
     * Load calibration data.
     *
     * \param  input                    detector calibration data
     */
    void load(JObjectIterator<JCOMPASS::JOrientation>& input)
    {
      using namespace JPP;

      t0_s = std::numeric_limits<double>::lowest();

      while (input.hasNext()) {

        const JOrientation* orientation = input.next();

        static_cast<container_type&>(calibration[orientation->id]).push_back(element_type(orientation->t, getQuaternion(*orientation)));
      }

      for (data_type::iterator i = calibration.begin(); i != calibration.end(); ++i) {
        i->second.sort();
        i->second.compile();
      }
    }


    /**
     * Check validity of calibration.
     *
     * \return                          true if valid; else false
     */
    bool is_valid() const
    {
      for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
        if (!calibration.has(module->getID())) {
          return false;
        }
      }

      return true;
    }


    /**
     * Get minimal abscissa.
     *
     * \return                           minimal abscissa
     */
    double getXmin() const
    {
      double xmin = std::numeric_limits<double>::max();

      for (const_iterator i = this->begin(); i != this->end(); ++i) {
        if (!i->second.empty() && i->second.getXmin() < xmin) {
          xmin = i->second.getXmin();
        }
      }

      return xmin;
    }


    /**
     * Get maximal abscissa.
     *
     * \return                           maximal abscissa
     */
    double getXmax() const
    {
      double xmax = std::numeric_limits<double>::lowest();

      for (const_iterator i = this->begin(); i != this->end(); ++i) {
        if (!i->second.empty() && i->second.getXmax() > xmax) {
          xmax = i->second.getXmax();
        }
      }

      return xmax;
    }


    const_iterator         begin()  const { return calibration.begin(); }    //!< begin of calibration data
    const_iterator         end()    const { return calibration.end(); }      //!< end   of calibration data
    const_reverse_iterator rbegin() const { return calibration.rbegin(); }   //!< begin of reverse of calibration data
    const_reverse_iterator rend()   const { return calibration.rend(); }     //!< begin of reverse of calibration data


    /**
     * Get detector calibrated at given time.
     *
     * \param  t1_s                     time [s]
     * \return                          detector
     */
    const JDetector& operator()(const double t1_s)
    {
      using namespace std;
      using namespace JPP;

      if (!calibration.empty()) {

        if (fabs(t1_s - t0_s) > Tmax_s) {

          for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
          
            if (calibration.has(module->getID())) {

              const function_type& f1 = calibration.get(module->getID());

              if (!f1.empty()) {

                if (t1_s >= f1.getXmin() && t1_s <= f1.getXmax()) {

                  JQuaternion3D Q0 = module->getQuaternion();
                  JQuaternion3D Q1 = f1(t1_s);

                  module->rotate(Q1 * Q0.conjugate());
                  module->setQuaternion(Q1);
                }
              }
            }
          }

          t0_s = t1_s;
        }
      }

      return detector;
    }


  protected:
    data_type  calibration;
    JDetector  detector;
    double     Tmax_s;

  private:
    double     t0_s;
  };
}

#endif