#ifndef TP0DGeom_hxx_seen
#define TP0DGeom_hxx_seen

#include <vector>
#include <map>

#include <HEPUnits.hxx>
#include <IDelaunay2D.hxx>

#include "EoaGeomInfo.hxx"
#include "IScintBarGeom.hxx"

#include "IGeomBase.hxx"

namespace COMET {
    OA_EXCEPTION(ENoSuchBar,EoaGeomInfo);
    class IP0DGeom;
    class IP0DGeomVisitor;
    class IGeomInfo;
};


/// Summary of the P0D (Pizero Detector) Geometry.  This contains useful
/// information about the P0D geometry which is derived from the ROOT
/// TGeoManager geometry description of the off-axis detector.  This must
/// accessed through IGeomInfo::P0D().  IP0DGeom is a "singleton" object
/// that is managed by the IGeomInfo class (which really is a singleton).
///
/// Most of the information is provided through the GetBar(IGeometryId) method
/// which returns a summary of the bar based on its geometry id.  In
/// particular, this is the best way to get the super-P0Dule, P0Dule, Layer,
/// and bar index of a bar.  For example, to translate a geometry ID into a
/// bar P0Dule you would do:
///
/// \code
/// int p0dule = IGeomInfo::P0D().GetBar(geomid).GetP0Dule();
/// \endcode
///
class COMET::IP0DGeom: public COMET::IGeomBase {
    friend class IP0DGeomVisitor;
    friend class IGeomInfo;

public:
    /// A class with information about each bar.  
    class Bar : public IScintBarGeom, public COMET::IMeshPoint {
    public:
        typedef std::vector<COMET::IGeometryId> Neighbors;

        Bar(COMET::IGeometryId geomId, int orientation, double x, double y);
        virtual ~Bar();
        
        /// Get the Super P0Dule containing this bar.  The
        /// numbering starts at the upstream end and starts from
        /// zero.  0) USECal, 1) USWT, 2) CWT, 3) CECal.
        int GetSuperP0Dule() const;

        /// Get the P0Dule number counting from the upstream end of the P0D [0
        /// to 39].
        int GetP0Dule() const;

        /// Return the Layer for this bar [0 is X, 1 is Y].
        int GetLayer() const;

        /// Return the bar number [0 to 125 for X, 0 to 133 for Y].
        int GetNumber() const;

        /// Internal methods required by the COMET::IMeshPoint.
        virtual double X(void) const {return fX;}
        virtual double Y(void) const {return fY;}

    private:
        /// The X position of this bar that is used to find the nearest
        /// neighbors.
        double fX;

        /// The Y position of this bar that is used to find the nearest
        /// neighbors.
        double fY;
    };

    typedef std::map<COMET::IGeometryId,COMET::IP0DGeom::Bar*> BarMap;

public:
    IP0DGeom();
    ~IP0DGeom();

    /// Find the index of the active plane in the detector.  This method
    /// returns the layer index, but should only be used when you don't have a
    /// bar geometry id (for instance, you don't have a hit).  If you have the
    /// hit, then you should use the GetBar(IGeometryId) method using
    /// IHit::GetGeomId().  Return a negative number if the Z is not on a
    /// plane.  The active planes include planes that are part of the upstream
    /// ecal, the water target, the carbon target, and the central ecal
    /// regions.
    int ActivePlane(double Z) const;

    /// Return the Z position of the indexed active plane.  The indices are
    /// defined by the return result of ActivePlane.
    double ActivePlaneZ(int index) const;

    /// Return the number of active planes in the detector.
    int ActivePlaneCount(void) const;

    /// Find the index of an active X plane in the detector.  This method
    /// returns the layer index, but should only be used when you don't have a
    /// bar geometry id (for instance, you don't have a hit).  If you have the
    /// hit, then you should use the GetBar(IGeometryId) method using
    /// IHit::GetGeomId().  Return a negative number if the Z is not on an X
    /// plane.  The active planes include planes that are part of the upstream
    /// ecal, the water target, the carbon target, and the central ecal
    /// regions.
    int ActiveXPlane(double Z) const;

    /// Return the Z position of the indexed active X plane.  The indices are
    /// defined by the return result of ActiveXPlane.
    double ActiveXPlaneZ(int index) const;

    /// Return the number of active X planes in the detector.
    int ActiveXPlaneCount(void) const;

    /// Find the index of an active Y plane in the detector.  This method
    /// returns the layer index, but should only be used when you don't have a
    /// bar geometry id (for instance, you don't have a hit).  If you have the
    /// hit, then you should use the GetBar(IGeometryId) method using
    /// IHit::GetGeomId().  Return a negative number if the Z is not on a Y
    /// plane. The active planes include planes that are part of the upstream
    /// ecal, the water target, the carbon target, and the central ecal
    /// regions.
    int ActiveYPlane(double Z) const;

    /// Return the Z position of the indexed active Y plane.  The indices are
    /// defined by the return result of ActiveXPlane.
    double ActiveYPlaneZ(int index) const;

    /// Return the number of active Y planes in the detector.
    int ActiveYPlaneCount(void) const;

    /// The minimum corner for the scintillator in the P0D.
    const TVector3& ActiveMin(void) const {return fActiveMin;}
    
    /// The minimum corner for the scintillator in the P0D.
    const TVector3& ActiveMax(void) const {return fActiveMax;}

    /// The minimum corner for the upstream ECal and Preshower in the P0D.
    const TVector3& USECalMin(void) const {return fUSECalMin;}
    
    /// The maximum corner for the upstream ECal and Preshower in the P0D.
    const TVector3& USECalMax(void) const {return fUSECalMax;}

    /// The minimum corner for the central ECal and Preshower in the P0D.
    const TVector3& CECalMin(void) const {return fCECalMin;}
    
    /// The maximum corner for the central ECal and Preshower in the P0D.
    const TVector3& CECalMax(void) const {return fCECalMax;}

    /// The minimum corner for the water target region in the P0D.
    const TVector3& TargetMin(void) const {return fWaterMin;}
    
    /// The minimum corner for the water target region in the P0D.
    const TVector3& TargetMax(void) const {return fWaterMax;}

    /// \deprecated The minimum corner for the water target region in the P0D.
    const TVector3& WaterMin(void) const {return fWaterMin;}
    
    /// \deprecated The minimum corner for the water target region in the P0D.
    const TVector3& WaterMax(void) const {return fWaterMax;}

    /// \deprecated The minimum corner for the carbon target region in the
    /// P0D.
    const TVector3& CarbonMin(void) const {return fCarbonMin;}
    
    /// \deprecated The minimum corner for the carbon target region in the
    /// P0D.
    const TVector3& CarbonMax(void) const {return fCarbonMax;}

    /// Find the distance to the nearest water target.  This is the signed
    /// distance along the Z axis.
    double WaterDistance(double Z) const;

    /// @{ Get the neighbors and other information for a particular bar.  This
    /// will throw a ENoSuchBar exception if the geometry id doesn't
    /// correspond to a P0D bar.
    const Bar& GetBar(COMET::IGeometryId geomId) const;
    const Bar& GetBar(const COMET::IHit& hit) const;
    /// @}

    /// Get all the P0D Bars.
    const BarMap& GetBars() const;

    /// Get a vector of the node ID's for the X bars
    const std::vector<COMET::IGeometryId>& GetXBars(void) const {return fXBars;}

    /// Get a vector of the node ID's for the Y bars
    const std::vector<COMET::IGeometryId>& GetYBars(void) const {return fYBars;}

    /// Get the distance to the edge of the fiducial volume along the Z axis.
    double WaterFiducialDistance(double Z) const;
    
    /// Get the distance to the edge of the fiducial volume.  This water
    /// fiducial distance takes 3 optional parameters.  The first is the Z
    /// position, if this is excluded then the calculated value is how far the
    /// point is from the horizontal and vertical boundary of the fiducial
    /// volume.  The second (xyFiducial) is the size of the horizontal and
    /// vertical fiducial cut.  The third optional parameter is the size of
    /// the fiducial cut along the Z axis.  The fiducial volume is found
    /// relative to the maximum X, Y, or Z of a scintillator bar in the water
    /// target.  For example, the upstream end of the fiducial volume is
    /// relative to the center of the most upstream layer (the center of the X
    /// layer in the first P0Dule in the upstream water target super-P0Dule).
    /// The downstream end of the fiducial volum is relative to the most
    /// downstream layer (the center of the Y layer in the last P0Dule of the
    /// central water target super-P0Dule).
    double WaterFiducialDistance(double X, double Y, 
                                 double Z=-100*unit::m,
                                 double xyFiducial=25*unit::cm,
                                 double zFiducial=0*unit::cm) const;

    /// Get the distance to the edge of the P0D along the Z axis.
    double EdgeDistance(double Z) const;
    
    /// Get the distance to the edge of the P0D.
    double EdgeDistance(double X, double Y, double Z=-100*unit::m) const;

    /// Get the average water depth for the P0D.
    double TargetDepth() const;

    /// Get the depth of water for a particular target module.  If side is
    /// zero then the average depth for the layer is returned.  If side is
    /// positive (negative), then the depth for the positive (negative) X side
    /// is returned.  If there is an error, a negative depth will be returned.
    double TargetDepth(COMET::IGeometryId target, int side=0) const;

    void ls(const char* opt);
    
private:
    /// Find the index of a Z value in a vector.
    int FindZPlane(std::vector<double> planes, double z, double toler) const;

    void CompressPlanes(std::vector<double>& planes);

    /// The list of Z coordinates for the scintillator planes in the P0D.
    std::vector<double> fActiveZ;

    /// The list of Z coordinates for the XZ scintillator planes in the P0D.
    std::vector<double> fActiveXZ;

    /// The list of Z coordinates for the YZ scintillator planes in the P0D.
    std::vector<double> fActiveYZ;

    /// The list of Z coordinates for the water target planes in the P0D.
    std::vector<double> fWaterZ;

    /// The minimum corner of the scintillator in the P0D.
    TVector3 fActiveMin;

    /// The maximum corner of the scintillator in the P0D.
    TVector3 fActiveMax;

    /// The minimum corner of the scintillator in the Upstream ecal and
    /// preshower.
    TVector3 fUSECalMin;

    /// The maximum corner of the scintillator in the Upstream ecal and
    /// preshower.
    TVector3 fUSECalMax;

    /// The minimum corner of the scintillator in the Central ecal and
    /// preshower.
    TVector3 fCECalMin;

    /// The maximum corner of the scintillator in the Central ecal and
    /// preshower.
    TVector3 fCECalMax;

    /// The minimum corner of the scintillator in the water target region.
    TVector3 fWaterMin;

    /// The minimum corner of the scintillator in the water target region.
    TVector3 fWaterMax;

    /// The minimum corner of the scintillator in the carbon target region.
    TVector3 fCarbonMin;

    /// The minimum corner of the scintillator in the carbon target region.
    TVector3 fCarbonMax;

    // A method called by IGeomInfo when the geometry information needs to be
    // filled.
    void Fill();

    // Clear the current geometry information
    void Clear();

    // The triangulation to determine the neighbors.
    COMET::IDelaunay2D *fNeighborsXZ;
    COMET::IDelaunay2D *fNeighborsYZ;

    // The map of geometry node id to the bar neighbor information.
    BarMap fBarMap;

    // A vector of all the X bars.
    std::vector<COMET::IGeometryId> fXBars;

    // A vector of all the Y bars.
    std::vector<COMET::IGeometryId> fYBars;
};
#endif