JMonteCarloDetector.hh

Go to the documentation of this file.

#ifndef __JDETECTOR__JMONTECARLODETECTOR__
#define __JDETECTOR__JMONTECARLODETECTOR__

#include <string>
#include <istream>
#include <vector>
#include <cmath>
#include <algorithm>

#include "JMath/JConstants.hh"
#include "JLang/JObjectID.hh"
#include "JDetector/JDetector.hh"
#include "JGeometry3D/JEulerAngle3D.hh"
#include "JGeometry3D/JEulerMatrix3D.hh"
#include "JGeometry3D/JAngle3D.hh"


/**
 * \author mdejong
 */

namespace JDETECTOR {}
namespace JPP { using namespace JDETECTOR; }

namespace JDETECTOR {

  using JMATH::PI;
  using JLANG::JObjectID;
  using JLANG::getUndefinedObjectID;
  using JGEOMETRY3D::JEulerAngle3D;
  using JGEOMETRY3D::JEulerMatrix3D;
  using JGEOMETRY3D::JAngle3D;
  using JGEOMETRY3D::JDirection3D;
  using JGEOMETRY3D::JAxis3D;


  static const int LED_BEACON_PMT_TYPE = 2;          //!< PMT type of LED beacon


  /**
   * Monte Carlo detector (i.e.\ \".det\" file).
   */
  class JMonteCarloDetector : 
    public JDetector
  {
  private:
    /**
     * Auxiliary class for OM 
     */
    class OM : 
      public JObjectID 
    {
    public:
      /**
       * Constructor.
       *
       * \param  in        input stream
       */
      OM(std::istream& in)
      {
        in >> static_cast<JObjectID&>(*this) >> pmtType >> serialNumber >> address;
      }

      int pmtType;        //!< PMT type
      int serialNumber;   //!< PMT serial number
      int address;        //!< address
    };


    /**
     * Auxiliary class for OM cluster parameters
     */
    class OM_cluster : 
      public JObjectID, 
      public std::vector<int> 
    {
    public:
      /**
       * Constructor.
       *
       * \param  in        input stream
       */
      OM_cluster(std::istream& in)
      {
        int n;
          
        in >> static_cast<JObjectID&>(*this) >> string_id >> z >> n;
          
        this->resize(n);

        iterator p = this->begin();

        for ( ; n != 0; --n, ++p)
          in >> *p;
      }


      int    string_id;   //!< string identifier
      double z;           //!< z position
    };


    /**
     * Auxiliary class for OM cluster data
     */
    class OM_cluster_data : 
      public JObjectID
    {
    public:
      /**
       * Constructor.
       *
       * \param  in        input stream
       */
      OM_cluster_data(std::istream& in)
      {
        in >> static_cast<JObjectID&>(*this) >> x >> y >> z >> theta >> phi >> psi;
      }


      double x;           //!< x position
      double y;           //!< y position
      double z;           //!< z position
      double theta;       //!< Euler angle
      double phi;         //!< Euler angle
      double psi;         //!< Euler angle
    };

    /**
     * Auxiliary class for OM position
     */
    class OM_position : 
      public JObjectID 
    {
    public:
      /**
       * Constructor.
       *
       * \param  in        input stream
       */
      OM_position(std::istream& in)
      {
        in >> static_cast<JObjectID&>(*this) >> x >> y >> z >> theta >> phi;
      }


      double x;           //!< x position
      double y;           //!< y position
      double z;           //!< z position
      double theta;       //!< zenit angle of orientation
      double phi;         //!< azimuth angle of orientation
    };


    /**
     * Auxiliary class for string parameters
     */
    class String : 
      public JObjectID 
    {
    public:
      /**
       * Constructor.
       *
       * \param  in        input stream
       */
      String(std::istream& in)
      {
        in >> static_cast<JObjectID&>(*this) >> x >> y >> z >> tOffset >> twist >> twistRate;
      }


      double x;           //!< x position
      double y;           //!< y position
      double z;           //!< z position
      double tOffset;     //!< time offset
      double twist;       //!< twist offset
      double twistRate;   //!< twist rate
    };


    /**
     * Auxiliary class for LCM logic parameters
     */
    class LCM_logic : 
      public JObjectID, 
      public std::vector<int> 
    {
    public:
      /**
       * Constructor.
       *
       * \param  in        input stream
       */
      LCM_logic(std::istream& in)
      {
        int n;
          
        in >> static_cast<JObjectID&>(*this) >> string_id >> lcm_id >> n;

        this->resize(n);

        iterator p = this->begin();
          
        for ( ; n != 0; --n, ++p)
          in >> *p;
      }

      int string_id;      //!< string identifier
      int lcm_id;         //!< LCM identifier
    };


    /**
     * Auxiliary class for LCM logic parameters
     */
    class LCM_reverse_logic : 
      public JObjectID 
    {
    public:
      /**
       * Default constructor.
       */
      LCM_reverse_logic() :
        JObjectID(),
        lcm_id(-1),
        pmt_id(-1)
      {}


      /**
       * Constructor.
       *
       * \param  __id      identifier
       * \param  __lcm_id  LCM identifier
       * \param  __pmt_id  PMT identifier
       */
      LCM_reverse_logic(const int __id, 
                        const int __lcm_id, 
                        const int __pmt_id) :
        JObjectID(__id),
        lcm_id   (__lcm_id),
        pmt_id   (__pmt_id)
      {}


      int lcm_id;         //!< LCM identifier
      int pmt_id;         //!< PMT identifier
    };


  public:
    /**
     * Default constructor.
     */
    JMonteCarloDetector() :
      JDetector(),
      use_logic(false)
    {}


    /**
     * Constructor.
     *
     * \param  useLogic    apply logic data in file
     */
    JMonteCarloDetector(const bool useLogic) :
      JDetector(),
      use_logic(useLogic)
    {}



    /**
     * Read detector from input.
     *
     * \param  in         input stream
     * \param  detector   detector
     * \return            input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JMonteCarloDetector& detector)
    {
      using namespace std;

      detector.clear();


      vector<OM>              OMaddress;
      vector<OM_cluster>      OMcluster;
      vector<OM_position>     OMposition;
      vector<String>          DetectorString;
      vector<OM_cluster_data> OMclusterdata;
      vector<LCM_logic>       LCMlogic;


      for (string key, buffer; in >> key; ) {

        if      (key == "OM:")

          OMaddress.push_back(OM(in));

        else if (key == "OM_cluster:")

          OMcluster.push_back(OM_cluster(in));

        else if (key == "OM_position:")

          OMposition.push_back(OM_position(in));

        else if (key == "string:")

          DetectorString.push_back(String(in));

        else if (key == "OM_cluster_data:")

          OMclusterdata.push_back(OM_cluster_data(in));

        else if (key == "LCM_logic:")

          LCMlogic.push_back(LCM_logic(in));

        else

          getline(in, buffer);
      }


      // invert LCM_Logic data for fast lookup

      vector<LCM_reverse_logic> LCMreverselogic;

      for (vector<LCM_logic>::const_iterator i  = LCMlogic.begin(); i != LCMlogic.end(); ++i) {
        
        const int lcm_id = i->getID();

        for (LCM_logic::const_iterator j = i->begin(); j != i->end(); ++j) {

          const int id     = *j;
          const int pmt_id = distance(i->begin(),j) + 1;

          LCMreverselogic.push_back(LCM_reverse_logic(id, lcm_id, pmt_id));
        }
      }


      // sort data according to identifier

      sort(OMaddress      .begin(), OMaddress      .end());
      sort(OMcluster      .begin(), OMcluster      .end());
      sort(OMposition     .begin(), OMposition     .end());
      sort(DetectorString .begin(), DetectorString .end());
      sort(LCMreverselogic.begin(), LCMreverselogic.end());


      vector<OM_cluster>     ::const_iterator   OMclusterIterator;
      vector<String>         ::const_iterator   StringIterator;
      vector<int>            ::const_iterator   addressIterator;
      vector<OM>             ::const_iterator   OMIterator;
      vector<OM_position>    ::const_iterator   OMpositionIterator;
      vector<OM_cluster_data>::const_iterator   OMclusterdataIterator;


      for (OMclusterIterator  = OMcluster.begin(); 
           OMclusterIterator != OMcluster.end(); 
           OMclusterIterator++) {
        
        StringIterator = find(DetectorString.begin(),
                              DetectorString.end(),
                              OMclusterIterator->string_id);
        
        if (StringIterator != DetectorString.end()) {
          
          const JLocation location(StringIterator->getID(), -1);

          for (addressIterator  = OMclusterIterator->begin(); 
               addressIterator != OMclusterIterator->end(); 
               addressIterator++) {
            
            OMIterator = find(OMaddress.begin(),
                              OMaddress.end(),
                              *addressIterator);
            
            if (OMIterator != OMaddress.end()) {
              
              OMpositionIterator = find(OMposition.begin(),
                                        OMposition.end(),
                                        OMIterator->address);
              
              if (OMpositionIterator != OMposition.end()) {
                
                // rotate string
                
                double ct = cos(StringIterator->twist * PI / 180.0);
                double st = sin(StringIterator->twist * PI / 180.0);
                
                double x  =  ct * OMpositionIterator->x - st * OMpositionIterator->y;
                double y  =  st * OMpositionIterator->x + ct * OMpositionIterator->y;

                JPosition3D  pos(StringIterator->x + x,
                                 StringIterator->y + y,
                                 StringIterator->z + OMpositionIterator->z + OMclusterIterator->z);
                
                JDirection3D dir(JAngle3D(OMpositionIterator->theta, 
                                          StringIterator->twist * PI / 180.0 + OMpositionIterator->phi));
                
                int lcm_id  = OMclusterIterator->getID();
                int pmt_id  = OMIterator->address;
                
                if (detector.use_logic) {
                  
                  vector<LCM_reverse_logic>::const_iterator p = find(LCMreverselogic.begin(), 
                                                                     LCMreverselogic.end(), 
                                                                     OMIterator->getID());
                  
                  if (p != LCMreverselogic.end()) {
                    lcm_id = p->lcm_id;
                    pmt_id = p->pmt_id;
                  }
                }
                
                detector.put(lcm_id, pmt_id - 1, location, JPMT(OMIterator->getID(), JAxis3D(pos, dir)));
              }
            }
          }
        }
      }


      // 'event' data

      for (OMclusterdataIterator  = OMclusterdata.begin(); 
           OMclusterdataIterator != OMclusterdata.end(); 
           OMclusterdataIterator++) {
          
        const JPosition3D    P(OMclusterdataIterator->x,
                               OMclusterdataIterator->y,
                               OMclusterdataIterator->z);
        
        const JEulerMatrix3D R(JEulerAngle3D(OMclusterdataIterator->theta,
                                             OMclusterdataIterator->phi,
                                             OMclusterdataIterator->psi));
        
        OMclusterIterator = find(OMcluster.begin(),
                                 OMcluster.end(),
                                 OMclusterdataIterator->getID());

        if (OMclusterIterator != OMcluster.end()) {

          const JLocation location(OMclusterIterator->string_id, -1);

          for (addressIterator  = OMclusterIterator->begin(); 
               addressIterator != OMclusterIterator->end(); 
               addressIterator++) {

            OMIterator = find(OMaddress.begin(),
                              OMaddress.end(),
                              *addressIterator);
              
            if (OMIterator != OMaddress.end()) {
                
              OMpositionIterator = find(OMposition.begin(),
                                        OMposition.end(),
                                        OMIterator->address);
                
              if (OMpositionIterator != OMposition.end()) {
                  
                JPosition3D  pos(OMpositionIterator->x,
                                 OMpositionIterator->y,
                                 OMpositionIterator->z);
                  
                JDirection3D dir(JAngle3D(OMpositionIterator->theta, 
                                          OMpositionIterator->phi));
                                  
                pos.transform(R);
                dir.transform(R);
                  
                pos += P;
                  
                int lcm_id  = OMclusterIterator->getID();
                int pmt_id  = OMIterator->address;

                if (detector.use_logic) {

                  vector<LCM_reverse_logic>::const_iterator p = find(LCMreverselogic.begin(), 
                                                                     LCMreverselogic.end(), 
                                                                     OMIterator->getID());
                  
                  if (p != LCMreverselogic.end()) {
                    lcm_id = p->lcm_id;
                    pmt_id = p->pmt_id;
                  }
                }

                detector.put(lcm_id, pmt_id - 1, location, JPMT(OMIterator->getID(), JAxis3D(pos, dir)));
              }
            }
          }
        }
      }

      // set floor identifier

      sort(detector.begin(), detector.end(), JMonteCarloDetector::compare);

      int string = -1;
      int floor  = -1;
      
      for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
        
        if (module->getString() != string) {
          string = module->getString();
          floor  = 0;
        }

        module->setLocation(JLocation(string,++floor));
      }
      
      return in;
    }


    /**
     * Set usage of logic.
     *
     * \param  useLogic    apply logic data in file
     */
    void setUseLogic(const bool useLogic) 
    {
      use_logic = useLogic;
    }


  protected:
    /**
     * Set PMT.
     * 
     * \param  moduleID    module identifier
     * \param  pmtAddress  PMT address
     * \param  location    module location
     * \param  pmt         JPMT
     * \return             this JDetector
     */
    JDetector& put(const int        moduleID,
                   const int        pmtAddress,
                   const JLocation& location,
                   const JPMT&      pmt)
    {
      JDetector::iterator p = std::lower_bound(this->begin(), this->end(), moduleID);

      if (p == this->end() || p->getID() != moduleID) {
        
        JModule module(moduleID, location);
        
        module.resize(pmtAddress + 1);

        module[pmtAddress] = pmt;

        module.setPosition(pmt.getPosition());
        
        this->insert(p, module);
        
      } else {
        
        // average position
        
        int n = 0;

        for (JModule::const_iterator i = p->begin(); i != p->end(); ++i) {
          if (i->getID() != getUndefinedObjectID().getID()) {
            ++n;
          }
        }

        if (pmtAddress + 1 > (int) p->size()) {
          p->resize(pmtAddress + 1);
        }

        JPosition3D& P = static_cast<JPosition3D&>(*p);

        P.mul(n);

        if (p->at(pmtAddress).getID() != getUndefinedObjectID().getID())
          P.sub(p->at(pmtAddress).getPosition());
        else
          ++n;

        P.add(pmt.getPosition());
        P.div(n);

        p->at(pmtAddress) = pmt;
      }

      return *this;
    }


    /**
     * Binary search method.
     *
     * \param  __begin     begin of data
     * \param  __end       end   of data
     * \param  id          target value
     * \return             pointer to corresponding element 
     */
    template<class T>
    static T find(T __begin, T __end, const int id) 
    {
      T i = std::lower_bound(__begin, __end, id);

      if (i != __end && *i != id)
        return __end;
      else
        return i;
    }


    /**
     * Module comparator.\n
     * The comparison is based on:
     *    -# string identifier;
     *    -# z-position.
     *
     * \param  first       first  module
     * \param  second      second module
     * \return             true if first module before second; else false
     */
    static bool compare(const JModule& first, const JModule& second)
    {
      if (first.getString() == second.getString())
        return first.getZ()      < second.getZ();
      else
        return first.getString() < second.getString();
    }

  private:
    bool use_logic;
  };
}

#endif