GammaLib

#ifndef ROOT_Tools_Region

#define ROOT_Tools_Region

#include <string>

#include <vector>

#include <TPolyLine.h>

#include <TH1D.h>

#include <TH2D.h>

#include <stash/Coordinate.hh>

#include <stash/Lambda.hh>

namespace Sash {

  class Time;

}

namespace Tools

{

  class Circle;

  class Box;

  /**

   * \brief Base class for all regions.

   * 

   * This is a base class that defines a physical region on the sky,

   * with methods to test if a point is inside (Region::Contains(

   * point )) or if regions overlap (Region::Overlaps(region )). The

   * base-class is abstract, and defines a set of virtual functions

   * that may be implemented by sub-classes, such as Tools::Box, or Tools::Circle.

   *

   * Regions can be:

   *  - rotated (about an arbitrary point) via SetRotation()

   *  - inverted (rejecting points inside and accepting points outside)  via SetInverted()

   *  - grouped into Tools::Compound regions (using inverted regions, you may construct complex shapes)

   * 

   * \note Since Tools::Region defines an internal rotation transformation

   * (about an arbitrary coordinate), sub-classes do not need to deal

   * with rotations directly, and can assume that the region is fixed at the

   * center position and aligned with the X-Y plane. 

   *

   * \note that this means that internally, fCentre can be assumed to

   * be unrotated, and should be used in methods like IsInside(),

   * while GetCentre() returns the true region center on the sky

   * (after applying the rotation transformation).

   *

   * \ingroup utilities

   *

   * \author Heidelberg team + CEA team

   */

  class Region : public TNamed, public TAttLine, public TAttFill

  {

  public:

    typedef std::vector<const Tools::Region*> RegionList;

    typedef std::vector<const Tools::Region*>::const_iterator RegionListIterator;

    typedef std::vector< std::vector< Stash::Coordinate > > VertexList;

    Region(const Stash::Coordinate& cor =

           Stash::Coordinate(0,0,Stash::System::GetDummySystem()),

           std::string name = "Region", 

           std::string type = "",

           Stash::Lambda phi = Stash::Lambda(0,Stash::Angle::Degrees));

    virtual ~Region(){;}

    double GetDistanceToCentre(const Stash::Coordinate& cor) const;

    /**

     * \brief Get the true center (after rotation is applied) of the region.

     *

     * This is the inverse-transform of fCentre.

     */

    Stash::Coordinate GetTrueCentre() const {return InverseTransform(fCentre);}

    /**

     * \brief Get center of the region (before rotation is applied, as given during construction).

     */

    Stash::Coordinate GetCentre() const {return fCentre;}

    /**

     * \brief Get coordinate system of the region center.

     */

    const Stash::SystemRef& GetSystem() const { return fCentre.GetSystem(); }

    /**

     * \brief Move the region to another part of the sky.

     */

    virtual void MoveTo(Stash::Coordinate newcen) {

      // generally move the center (doesn't work the same for some

      // region subtypes though, e.g. polygons where the center is not

      // strictly relavant). For those, it must be re-implemented

      fCentre = newcen.GetCoordinate(*fCentre.GetSystem()); 

    }

    bool Overlaps(const Tools::Region& reg) const;

    /**

     * \brief String representation of region type.

     */

    std::string GetType() { return fType; }

    double GetEquivalentCircleRadius() const;

    /**

     * \brief accuracy is used in calculations that sample points over

     * the region (e.g. to calculate overlaps or zenith distributions,

     * etc). Increasing it speeds things up at the expense of some accuracy.

     */

    void SetAccuracy( double maxerror ) { fAccuracy= maxerror;}

    Stash::Lambda GetRotationAngle() const { return fRotationAngle; }

    Stash::Coordinate GetRotationPoint() const { return fRotationPoint; }

    void SetRotation( Stash::Lambda phi=Stash::Lambda(0,Stash::Angle::Degrees), 

                      Stash::Coordinate point=Stash::Coordinate(0,0,Stash::System::GetDummySystem()) );

    /**

     * Control whether region is inverted or not (inverted means

     * points outside will return true Contains())

     */

    void SetInverted(bool inv=true) { fIsInverted=inv; } 

    bool IsInverted() const { return fIsInverted;}

    /**

     * \brief Grow the boundary of the region by a given amount of space angle.

     *

     * \note Has to be implemented by subclasses.

     */

    virtual Tools::Region* GrowBoundary(const double degrees,

                                        const std::string name="_grown") const =0;

    bool Contains(const Stash::Coordinate& dir) const;

    virtual bool IsCompletelyContained(const Tools::Region& reg) const;

    virtual Stash::Coordinate GetCentreOfGravity() const;

    // These Fill functions are used by the BgStats class to fill the

    // expected distributions of events in zenith/offset/azimuth.

    virtual void FillOffsetDistribution(const Stash::Coordinate& pos,

                                        TH1 &, double psicut=-1 ) const;

    virtual void FillZenithOrAzimuthDistribution( TH1 &, const Sash::Time, 

                                                  const Sash::Time, bool fillZenith=true) const;

    virtual void FillZenithOffsetDistribution(const Stash::Coordinate& pos,

                                              TH2 &, const Sash::Time, const Sash::Time, 

                                              double psicut=-1 ) const;

    virtual void Fill2AxesDistanceDistribution(const Stash::Coordinate& pos,

                                               TH2 &, double psicut=-1 ) const;

    virtual double GetArea() const =0;  ///< Returns area of region in square degrees

    virtual double GetBound2() const=0; ///< Return width of bounding box

    virtual double GetBound1() const=0; ///< Return length of bounding box

    virtual double GetRadialBound()const=0; ///< Return radius of circle containing reg

    Tools::Box* GetBoundingBox() const;

    Tools::Circle* GetBoundingCircle() const;

    double GetMonteCarloArea(int samples_per_square_deg=1e6) const;

    /**

     * \brief Return FITS region string representation of this region.

     *

     * Subclasses should implement this such that it returns a valid

     * FITS region string. Used by the output streamer and Print()

     */

    virtual std::string AsFITSRegionString() const {

      return "# "+ std::string(GetName())+ ": AsFITSRegionString() is not "

        "implemented for this region type ";

    }

    /**

     * \brief Return VertexList representation of this region.

     *

     * To be implemented by subclasses: shold return a

     * representation/approximatin of the region as a continuous list

     * of coordinates that follows the boundaries of the region (this

     * is used for drawing the regions in the correct transformation)

     *

     * The vertex list returned is in the native system of the region. 

     */

    virtual VertexList AsVertexList() const{

      VertexList blanklist;

      return blanklist;

    };

  protected:

    Stash::Coordinate fCentre; ///< Centre of Region

    std::string GetFITSSystemName() const;

    double fAccuracy; ///< accuracy of grid or monte-carlo operations 1.0=normal, larger is more accurate (more steps used in loops)

    Stash::Coordinate Transform( const Stash::Coordinate &src ) const;

    Stash::Coordinate InverseTransform( const Stash::Coordinate &src ) const;

    Stash::Coordinate GetPointAtOffsetFromCenter( double offsX, double offsY ) const ;

    void GetOffsetOfPointFromCenter( Stash::Coordinate point, double &offsX, double &offsY, bool untransformed=false ) const ;

    Stash::Coordinate GetPointAtRadius( Stash::Beta rad, Stash::Lambda phi ) const;

    /**

     * Returns the optimal number of spatial steps for a region of a

     * given size, taking into accunt the accuracy leve; \param

     * radius: approximate radius of region in deg

     */

    int GetNumSpaceSteps(double radius) const { 

      return 100*2.0*radius*fAccuracy;

    }

    /**

     * Returns the optimal number of spatial steps for looping over a

     * time range

     */

    int GetNumTimeSteps() const {

      return 20*fAccuracy;

    }

    /**

     * \brief Check if dir is inside this region.

     *

     * Must be implemented by the sub-classes. This function should

     * return true if the point is inside the *un-rotated* region (the

     * rotation is taken into account automatically by Contains()

     */

    virtual bool IsInside( const Stash::Coordinate& dir ) const = 0;    

  private:

    Stash::Coordinate fRotationPoint; ///< Center of rotation for Transform()

    Stash::Lambda fRotationAngle; ///< Rotation Angle for Transform()

    std::string fType; ///< type name

    bool fIsInverted; ///< is this a negative region? (Contains() gives points outside)

    ClassDef(Tools::Region, 5);

  };

  /**

   * Defines a segment of an annulus.

   *

   * The annulus is defined by inner (r1) and outer (r2) radius. Radii are

   * in sky degrees.

   * The azimuthal coverage of the segment is given by two angles (restricted

   * in range as 0 <= phi <= 360 deg). 

   *

   * Note that a segment is always defined counter-clockwise, i.e. phi2 > phi1

   * in general. Zero crossing is possible, e.g. a half ring in the north is

   * defined by setting phi1=270, phi2=90 (phi2 < phi1 in this case).

   *

   * A full ring can be defined by setting phi1=0, phi2=360.

   */

  class CircleSegment: public Region

  {

  public:

    ~CircleSegment();

    CircleSegment(const Stash::Coordinate& cor =

                  Stash::Coordinate(0,0,Stash::System::GetDummySystem()),

                  double r1 = 0.,double r2 = 0., 

                  double phi1 = 0., double phi2 = 360.,

                  std::string name = "Segment");

    double GetR1() const {return fR1;};

    double GetR2() const {return fR2;};

    double GetPhi1() const {return fPhi1;};

    double GetPhi2() const {return fPhi2;};

    void SetR1(double val){fR1 = val;}

    void SetR2(double val){fR2 = val;}

    virtual void SetPhi1(double val){fPhi1=val;}

    virtual void SetPhi2(double val){fPhi2=val;}

    virtual Tools::CircleSegment* GrowBoundary(const double degrees,

                                               const std::string name="_grown") const;

    virtual double GetArea() const;

    virtual double GetBound2() const { return fR2*2.0; } 

    virtual double GetBound1() const { return fR2*2.0/Transform(fCentre).GetBeta().Cos(); }

    virtual double GetRadialBound() const { return fR2; }

    virtual std::string AsFITSRegionString() const;

    virtual VertexList AsVertexList() const;

  protected:

    virtual bool IsInside(const Stash::Coordinate& dir) const;

  private:

    double fR1; 

    double fR2; 

    double fPhi1;

    double fPhi2;

    mutable TPolyLine* fLine; //!

    ClassDef(Tools::CircleSegment, 3);

  };

  /**

   * A Ring with inner and outer radius. 

   * 

   * Same as a CircleSegment from phi=0-360, but with some speed optimizations.

   */

  class Ring : public CircleSegment {

  public:

    Ring(const Stash::Coordinate &pos 

         = Stash::Coordinate(0,0,Stash::System::GetDummySystem()),

         double r1=0, double r2=0, std::string name = "Ring") 

      : CircleSegment( pos, r1, r2, 0.0, 360.0, name ) {

    }

    void SetPhi1(double ) { 

      //      std::cout << "Ring: Can't set Phi1 for a Ring!" << std::endl;

    } 

    void SetPhi2(double ) { 

      //      std::cout << "Ring: Can't set Phi2 for a Ring!" << std::endl;

    } 

    virtual std::string AsFITSRegionString() const;

  protected:

    virtual bool IsInside( const Stash::Coordinate& dir ) const;

    ClassDef(Tools::Ring, 1);

  };

  /**

   * A circular region.

   *

   * Same as Ring, but with r1=0 for speed optimisations.

   */

  class Circle : public Ring {

  public:

    Circle(const Stash::Coordinate &pos 

           = Stash::Coordinate(0,0,Stash::System::GetDummySystem()),

         double radius=0, std::string name = "Ring") 

      : Ring( pos, 0.0, radius, name ) {

    }

    virtual std::string AsFITSRegionString() const;

  protected:

    virtual bool IsInside( const Stash::Coordinate& dir ) const;

    ClassDef(Tools::Circle, 1);

  };

  /**

   * Defines a rectangular region. L1 and L2 are the length and width

   * of the region.

   */

  class Box: public Region

  {

  public:

    ~Box();

     Box(const Stash::Coordinate& cor =

         Stash::Coordinate(0,0,Stash::System::GetDummySystem()),

         double l1 = 0.,double l2 = 0., 

         Stash::Lambda phi = Stash::Lambda(0,Stash::Angle::Degrees),

         std::string name = "Box");

     virtual VertexList AsVertexList() const;

     double GetL1() const {return fL1;};

     double GetL2() const {return fL2;};

     void SetL1(double val){fL1 = val;};

     void SetL2(double val){fL2 = val;};

    virtual Tools::Box* GrowBoundary(const double degrees,

                                     const std::string name="_grown") const;

    virtual double GetArea() const { return fL1*fL2; }

    virtual double GetRadialBound() const { 

      return sqrt(fL1*fL1/4.0+fL2*fL2/4.0); 

    }    

    virtual double GetBound1() const;

    virtual double GetBound2() const;

    Stash::Coordinate GetRandomSample();

    virtual std::string AsFITSRegionString() const;

   protected:

    virtual bool IsInside(const Stash::Coordinate& dir) const;

  private:

    double fL1;        /// length in degrees

    double fL2;        /// width in degrees

    ClassDef(Tools::Box, 4);

  };

  /**

   * A region which is a union of multiple simple regions. Compounds

   * can be thought of as region "groups", and may contain secondary

   * Compounds as well. This allows a hierarchy of regions to be

   * defined. Like all Regions, each Compound can be rotated (which

   * will rotate all regions inside the compound as a single group).

   *

   * \todo implement a method to add a negative region AddNegativeRegion() to allow

   * constructive geometry.

   */

  class Compound : public Region 

  {

  public:

    ~Compound() {;}

    Compound( std::string name = "Compound",const Stash::Coordinate& cor=

              Stash::Coordinate(0,0,Stash::System::GetDummySystem()) )

      : Region(cor, name){;}

    Compound* Clone(const char* newname) const;

    void AddRegion( const Tools::Region *reg );

    void Reset();

    void ClearRegionList() { fRegionList.clear(); }

    virtual Tools::Compound* GrowBoundary(const double degrees,

                                          const std::string name="_grown") const;

    virtual VertexList AsVertexList() const;

    virtual bool IsCompletelyContained(const Tools::Region& reg) const;

    virtual void FillZenithOffsetDistribution(const Stash::Coordinate& pos,

                                              TH2 & , const Sash::Time, 

                                              const Sash::Time, 

                                              double psicut=-1.) const;

    virtual void FillZenithOrAzimuthDistribution( TH1 &, const Sash::Time, 

                                                  const Sash::Time,

                                                  bool fillZenith=true) const;

    virtual void FillOffsetDistribution(const Stash::Coordinate& pos, TH1&, double psicut=-1.) const;

    virtual double GetArea() const;

    virtual double GetBound2() const;

    virtual double GetBound1()const;

    virtual double GetRadialBound() const;

    virtual std::string AsFITSRegionString() const;

    RegionList GetRegions() const {return fRegionList;};

    int GetNumRegions() const { return fRegionList.size(); }

    bool IsEmptyRegion() const { return (GetNumRegions() == 0); }

    virtual void MoveTo( Stash::Coordinate newpos );

    /**

     * Move the reference position to the Center-of-gravity

     * (doesn't move the region itself, just the point from which

     * bounding boxes are calculated.) Use MoveTo() instead to move the region

     */

    void RecenterOnCOG() { fCentre = GetCentreOfGravity(); }

    void Recenter(Stash::Coordinate refpos) { fCentre = refpos; }

  protected:

    virtual bool IsInside( const Stash::Coordinate &dir ) const;

    RegionList fRegionList;

  private:

    ClassDef(Tools::Compound, 3);

  };

  /**

   * A Compound, where Contains() only returns points common to all

   * regions inside it.  

   *

   * \note that currently the AsVertexList() method

   * (and thus drawing on a SkyHist) draws all the sub-regions, not

   * just the intersection.

   */

  class CompoundIntersection : public Compound

  {

  public:

    CompoundIntersection(std::string name = "Compound",const Stash::Coordinate& cor=

                         Stash::Coordinate(0,0,Stash::System::GetDummySystem()))

      : Compound( name, cor ) { ; } 

    virtual ~CompoundIntersection(){;}

  protected:

    virtual bool IsInside( const Stash::Coordinate &dir ) const;

    ClassDef(Tools::CompoundIntersection, 0);

  };

}

std::ostream& operator<<( std::ostream &stream, 

                          const Tools::Region& );

#endif