SensorInfoBase.h

/**************************************************************************
 * BASF2 (Belle Analysis Framework 2)                                     *
 * Copyright(C) 2010 - Belle II Collaboration                             *
 *                                                                        *
 * Author: The Belle II Collaboration                                     *
 * Contributors: Martin Ritter                                            *
 *                                                                        *
 * This software is provided "as is" without any warranty.                *
 **************************************************************************/
 
#ifndef VXD_SENSORINFO_H
#define VXD_SENSORINFO_H
 
#include <vxd/dataobjects/VxdID.h>
#include <float.h>
 
#include <TGeoMatrix.h>
#include <TVector3.h>
 
namespace Belle2 {
  namespace VXD {
 
    class SensorInfoBase {
    public:
      enum SensorType {
        PXD = 0,  
        SVD = 1,  
        TEL = 2,  
        VXD = -1, 
      };
 
      SensorInfoBase(SensorType type, VxdID id, double width, double length, double thickness,
                     int uCells, int vCells, double width2 = -1, double splitLength = -1, int vCells2 = 0):
        m_type(type), m_id(id), m_width(width), m_length(length), m_thickness(thickness),
        m_deltaWidth(0), m_splitLength(0), m_uCells(uCells), m_vCells(vCells), m_vCells2(vCells2)
      {
        if (width2 > 0) m_deltaWidth = width2 - width;
        if (splitLength > 0) m_splitLength = splitLength / length;
      }
      virtual ~SensorInfoBase() {}
 
      SensorType getType() const { return m_type; }
      VxdID getID() const { return m_id; }
 
      double getWidth(double v = 0) const
      {
        if (m_deltaWidth == 0) return m_width;
        return m_width + (v / m_length + 0.5) * m_deltaWidth;
      }
 
      double getBackwardWidth() const
      {
        return getWidth(-0.5 * m_length);
      }
 
      double getForwardWidth() const
      {
        return getWidth(0.5 * m_length);
      }
 
      double getLength() const { return m_length; }
 
      double getThickness() const { return m_thickness; }
 
      double getUSize(double v = 0) const { return getWidth(v); }
 
      double getVSize() const { return getLength(); }
 
      double getWSize() const { return getThickness(); }
 
      double getUPitch(double v = 0) const { return getWidth(v) / m_uCells; }
 
      double getVPitch(double v = 0) const
      {
        if (m_splitLength <= 0) return m_length / m_vCells;
        if (v / m_length + 0.5 >= m_splitLength) return m_length * (1 - m_splitLength) / m_vCells2;
        return m_length * m_splitLength / m_vCells;
      }
 
      int getVPitchID(double v = 0) const
      {
        if (m_splitLength <= 0) return 0;
        if (v / m_length + 0.5 >= m_splitLength) return 0;
        return 1;
      }
 
      double getUCellPosition(int uID, int vID = -1) const
      {
        if (m_deltaWidth == 0) return ((uID + 0.5) / m_uCells - 0.5) * m_width;
        double v = 0;
        if (vID >= 0) v = getVCellPosition(vID);
        return ((uID + 0.5) / m_uCells - 0.5) * getWidth(v);
      }
 
      double getVCellPosition(int vID) const
      {
        if (m_splitLength <= 0) return ((vID + 0.5) / m_vCells - 0.5) * m_length;
        if (vID >= m_vCells) return ((vID - m_vCells + 0.5) / m_vCells2 * (1 - m_splitLength) - 0.5 + m_splitLength) * m_length;
        return ((vID + 0.5) / m_vCells * m_splitLength - 0.5) * m_length;
      }
 
      int getUCellID(double u, double v = 0, bool clamp = false) const
      {
        if (clamp) return std::min(getUCells() - 1, std::max(0, getUCellID(u, v, false)));
        return static_cast<int>((u / getWidth(v) + 0.5) * m_uCells);
      }
 
      int getVCellID(double v, bool clamp = false) const
      {
        if (clamp) return std::min(getVCells() - 1, std::max(0, getVCellID(v, false)));
        double nv = v / m_length + 0.5;
        if (m_splitLength <= 0) return static_cast<int>(nv * m_vCells);
        if (nv >= m_splitLength) return static_cast<int>((nv - m_splitLength) / (1 - m_splitLength) * m_vCells2) + m_vCells;
        return static_cast<int>(nv / m_splitLength * m_vCells);
      }
 
      int getUCells() const { return m_uCells; }
      int getVCells() const { return m_vCells + m_vCells2; }
      int getVCells2() const { return m_vCells2; }
 
      bool inside(double u, double v, double uTolerance = DBL_EPSILON, double vTolerance = DBL_EPSILON) const
      {
        double nu = u / (getWidth(v) + 2 * uTolerance) + 0.5;
        double nv = v / (getLength() + 2 * vTolerance) + 0.5;
        return 0 <= nu && nu <= 1 && 0 <= nv && nv <= 1;
      }
 
      bool inside(const TVector3& local) const
      {
        double nw = local.z() / getThickness() + 0.5;
        return inside(local.x(), local.y()) && 0 <= nw && nw <= 1;
      }
 
      void forceInside(double& u, double& v) const
      {
        double length = getLength() / 2.0;
        v = std::min(length, std::max(-length, v));
        double width = getWidth(v) / 2.0;
        u = std::min(width, std::max(-width, u));
      }
 
      void forceInside(TVector3& local) const;
 
      TVector3 pointToGlobal(const TVector3& local, bool reco = false) const;
 
      TVector3 vectorToGlobal(const TVector3& local, bool reco = false) const;
 
      TVector3 pointToLocal(const TVector3& global, bool reco = false) const;
 
      TVector3 vectorToLocal(const TVector3& global, bool reco = false) const;
 
      void setTransformation(const TGeoHMatrix& transform, bool reco = false)
      {
        if (reco) m_recoTransform = transform;
        else m_transform = transform;
      }
 
      const TGeoHMatrix& getTransformation(bool reco = false) const
      {
        if (reco) return m_recoTransform;
        else return m_transform;
      }
 
      void setSurfaceParameters(const std::vector<double>& planarParameters)
      {
        m_surfaceDeformationParameters = planarParameters;
      }
 
      const std::vector<double>& getSurfaceParameters() const
      {
        return m_surfaceDeformationParameters;
      }
 
    protected:
      SensorType m_type;
      unsigned short m_id;
      double m_width;
      double m_length;
      double m_thickness;
      double m_deltaWidth;
      double m_splitLength;
      int m_uCells;
      int m_vCells;
      int m_vCells2;
      TGeoHMatrix m_transform;
      TGeoHMatrix m_recoTransform;
      std::vector<double> m_surfaceDeformationParameters = std::vector<double>(12, 0.0);
    };
 
    inline void SensorInfoBase::forceInside(TVector3& local) const
    {
      double u = local.x();
      double v = local.y();
      double thickness = getThickness() / 2.0;
      forceInside(u, v);
      local.SetX(u);
      local.SetY(v);
      local.SetZ(std::min(thickness, std::max(-thickness, local.z())));
    }
 
    inline TVector3 SensorInfoBase::pointToGlobal(const TVector3& local, bool reco) const
    {
      double clocal[3];
      double cmaster[3];
      local.GetXYZ(clocal);
      if (reco) m_recoTransform.LocalToMaster(clocal, cmaster);
      else m_transform.LocalToMaster(clocal, cmaster);
      return TVector3(cmaster);
    }
 
    inline TVector3 SensorInfoBase::vectorToGlobal(const TVector3& local, bool reco) const
    {
      double clocal[3];
      double cmaster[3];
      local.GetXYZ(clocal);
      if (reco) m_recoTransform.LocalToMasterVect(clocal, cmaster);
      else m_transform.LocalToMasterVect(clocal, cmaster);
      return TVector3(cmaster);
    }
 
    inline TVector3 SensorInfoBase::pointToLocal(const TVector3& global, bool reco) const
    {
      double clocal[3];
      double cmaster[3];
      global.GetXYZ(cmaster);
      if (reco) m_recoTransform.MasterToLocal(cmaster, clocal);
      else m_transform.MasterToLocal(cmaster, clocal);
      return TVector3(clocal);
    }
 
    inline TVector3 SensorInfoBase::vectorToLocal(const TVector3& global, bool reco) const
    {
      double clocal[3];
      double cmaster[3];
      global.GetXYZ(cmaster);
      if (reco) m_recoTransform.MasterToLocalVect(cmaster, clocal);
      else m_transform.MasterToLocalVect(cmaster, clocal);
      return TVector3(clocal);
    }
  }
} //Belle2 namespace
#endif