GammaLib

#include <iostream>

#include <string>

#include <TPad.h>

#include <TBox.h>

#include <TLine.h>

#include <TRandom.h>

#define INFO_OUT_TAG "Region> "

#define DEBUG 0

#include <utilities/debugging.hh>

#include <sash/HESSArray.hh>

#include <sash/TimeDiff.hh>

#include <sash/Time.hh>

#include <stash/constants.h>

#include "Region.hh"

using namespace std;

/** 

 */

Tools::Region::Region(const Stash::Coordinate& center,string name,

                      string type, Stash::Lambda phi)

  : TNamed(name.c_str(),""),

    TAttLine(), TAttFill(),

    fCentre(center), 

    fAccuracy(1.0),

    fRotationPoint(center),

    fRotationAngle(phi),

    fType(type),

    fIsInverted(false)

{

}

/**

 * Currently Applies a rotation to the given point (for passively

 * rotating the region - this should be called in Contains() after

 * converting to the RegionSystem to get the correct point)

 */

Stash::Coordinate 

Tools::Region::

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

  if (fRotationPoint.GetSystem() == Stash::System::GetDummySystem() )

    return src;

  Stash::Coordinate tsrc = 

    src.GetCoordinate( *fRotationPoint.GetSystem() ); // eventually should be in RegionSystem

  Stash::DirectionVector dv = tsrc.GetDirectionVector();

  Stash::DirectionVector cendv =  fRotationPoint.GetDirectionVector();

  Stash::DirectionVector newdv = dv.RotateAroundPoint( cendv, fRotationAngle);

  return Stash::Coordinate(newdv, fRotationPoint.GetSystem());

}

/**

 * Undo the internal rotation transform of the region (useful for

 * getting back to proper "Sky" coordinates, rather than the internal

 * de-rotated coordinates of the region.  It should be used in

 * AsVertexList() to undo the internal transform.

 */

Stash::Coordinate 

Tools::Region::

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

  Stash::Coordinate tsrc = 

    src.GetCoordinate( *fRotationPoint.GetSystem() ); // eventually should be in RegionSystem

  Stash::DirectionVector dv = tsrc.GetDirectionVector();

  Stash::DirectionVector cendv =  fRotationPoint.GetDirectionVector();

  Stash::DirectionVector newdv = dv.RotateAroundPoint( cendv, fRotationAngle.GetNegative() );

  return Stash::Coordinate(newdv, GetSystem());

}

/**

 * Set the rotation of this region. Each region can be rotated about

 * an arbitrary point (if no point is specified, the region center is

 * used).  The rotation angle is relative to North in the Region's

 * cordinate System, and increases counter-clockwise (the is the same

 * definition of the "position angle" in astronomy).

 *

 * \code

 *     // Example: rotate a region at the test position about the observation position

 *     Tools::Box *box = new Tools::Box( testpos, wid, len );

 *     box->SetRotation( Stash::Beta(45,Stash::Angle::Degrees), obspos );

 * \endcode

 *

 * If you want to apply successive rotations (e.g. rotate first about

 * one point, and then again about another), you should first rotate

 * the region and then add it to a Tools::Compound (which contains its

 * own rotation transform), and apply the second rotation to that.

 */

void 

Tools::Region::

SetRotation( Stash::Lambda phi, Stash::Coordinate point ) {

  if (point == Stash::Coordinate(0,0,Stash::System::GetDummySystem())) {

    DEBUG_OUT << GetName() << ": USING CENTER POINT FOR ROTATION POINT" << std::endl;

    point = GetCentre();

  }

  // don't rotate if region centre is uninitialised

  // exception: zero angle (needed by Compound::Reset())

  if (point == Stash::Coordinate(0,0,Stash::System::GetDummySystem())) {

    if( !(phi == Stash::Lambda(0)) ) {

      WARN_OUT << GetName() << ": Region Center is in dummy system. Can't rotate."

               << " If this is a Compound Region, you must add a region first to define"

               << " the proper coordinate system. "

               << std::endl;

      return;

    }

  }

  else

    fRotationPoint = point.GetCoordinate( *GetSystem() );

  fRotationAngle = phi;

  DEBUG_OUT << "Rotation set to: "<< fRotationAngle << " about "<<fRotationPoint

            << std::endl;

}

/**

 * \returns distance of the given coordinate to the true center-point of the region.

 */

double Tools::Region::GetDistanceToCentre(const Stash::Coordinate& cor) const

{

  Stash::Coordinate dir = cor.GetCoordinate(*GetSystem());

  return GetTrueCentre().GetAngularDistance(dir).GetDegrees();

}

/**

 * Loops over the region and calculates an approximate center of

 * gravity.  This can be re-implemented in sub-classes where the

 * center of gravity is known mathematically, to increase speed.

 */

Stash::Coordinate Tools::Region::GetCentreOfGravity() const

{

  double n=0;

  double suml=0.0, sumb=0.0;

  double radius = GetRadialBound();

  int steps = GetNumSpaceSteps(radius);

  DEBUG_OUT << "Steps=" << steps<< endl;

  double rl = radius/GetTrueCentre().GetBeta().Cos();

  for (int i=0; i<steps; ++i) {

    Stash::Lambda l = GetTrueCentre().GetLambda() 

      + Stash::Lambda(2.*rl*(i-steps*0.5)/(steps*0.9),Stash::Angle::Degrees);

    for (int j=0; j<steps; ++j) {

      Stash::Beta b = GetTrueCentre().GetBeta() 

        + Stash::Beta(2.*radius*(j-steps*0.5)/(steps*0.9),

                      Stash::Angle::Degrees);

      Stash::Coordinate test(l,b,GetSystem());

      if (Contains(test)) {

        sumb += b.GetDegrees();

        suml += l.GetDegrees();

        n += 1.0;

      }

    }

  }

  Stash::Coordinate cog(Stash::Lambda(suml/n,Stash::Angle::Degrees),

                        Stash::Beta(sumb/n,Stash::Angle::Degrees),

                        GetSystem());

  return cog;

}

/**

 * Returns the radius of a circle with the same area as the region

 */

double

Tools::Region::

GetEquivalentCircleRadius() const {

  return sqrt( GetArea()/M_PI );

}

/**

 * Fills a histogram of the Offset distribution from the given

 * reference position

 *

 * \param psicut - will be ignored, if set to -1.

 */

void

Tools::Region::FillOffsetDistribution(const Stash::Coordinate& refpos,

                                        TH1 &hist, double psicut ) const

{

  Stash::Coordinate pos = refpos.GetCoordinate(*GetSystem());

  double radius = GetRadialBound();

  double rl = radius/GetTrueCentre().GetBeta().Cos();

  const int steps = GetNumSpaceSteps(radius);

  DEBUG_OUT << "pos: " << pos

            << " refpos: " << refpos

            << " steps: "<< steps

            << endl;

  for (int i=0; i<steps; ++i) {

    Stash::Lambda l = GetTrueCentre().GetLambda() 

      + Stash::Lambda(2.*rl*(i-steps*0.5)/(steps*0.9),Stash::Angle::Degrees);

    for (int j=0; j<steps; ++j) {

      Stash::Beta b = GetTrueCentre().GetBeta() 

        + Stash::Beta(2.*radius*(j-steps*0.5)/(steps*0.9),

                      Stash::Angle::Degrees);

      Stash::Coordinate test(l,b,GetSystem());

      if (Contains(test)) {

        double d = pos.GetAngularDistance(test).GetDegrees();

        // ignore psi, if it is negativ

        if (psicut < 0. || d <= psicut) hist.Fill(d);

      }

    }

  }

  if( hist.Integral() == 0 ){

    WARN_OUT << "Offset histogram empty!" << std::endl;

  }

}

/**

 * \returns a string representing the coordinate system for use with AsFITSRegionString()

 */

string

Tools::Region::

GetFITSSystemName() const {

  string name = GetSystem().GetSystem()->GetName();

  if (name == "GalacticSystem")

    return "galactic";

  else if (name == "RADecJ2000System")

    return "fk5";

  return "unknown";

}

/**

 * \brief Check if dir is contained in this region.

 */

bool Tools::Region::Contains(const Stash::Coordinate& dir) const {

  Stash::Coordinate dir_trans = Transform(dir);

  if (fIsInverted)

    return !IsInside(dir_trans); // call the sub-class's virtual method

  else 

    return IsInside(dir_trans);

} 

/**

 * \brief Check if another region is completely contained in this region.

 */

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

{

  double radius = reg.GetRadialBound();

  double rl = radius/reg.GetTrueCentre().GetBeta().Cos();

  int steps = GetNumSpaceSteps(radius);

  for (int i=0; i<steps; ++i) {

    Stash::Lambda l = reg.GetTrueCentre().GetLambda() 

      + Stash::Lambda(2.*rl*(i-steps*0.5)/(steps*0.9),Stash::Angle::Degrees);

    for (int j=0; j<steps; ++j) {

      Stash::Beta b = reg.GetTrueCentre().GetBeta() 

        + Stash::Beta(2.*radius*(j-steps*0.5)/(steps*0.9),

                      Stash::Angle::Degrees);

      Stash::Coordinate test(l,b,GetSystem());

      if (reg.Contains(test) && (!Contains(test))) return false;

    }

  }

  return true;

}

/**

 * Check whether the given region overlaps this one

 */

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

{

  if( reg.GetSystem() == Stash::System::GetDummySystem() )

    return false;

  double radius = GetRadialBound();

  int steps = GetNumSpaceSteps(radius);

  double rl = radius/GetTrueCentre().GetBeta().Cos();

  // make sure to use the transformed region centre here

  // otherwise Overlaps() returns wrong results for rotated regions

  for (int i=0; i<steps; ++i) {

    Stash::Lambda l = GetTrueCentre().GetLambda() + 

      Stash::Lambda(2.*rl*(i-steps*0.5)/(steps*0.9),Stash::Angle::Degrees);

    for (int j=0; j<steps; ++j) {

      Stash::Beta b = GetTrueCentre().GetBeta() + 

        Stash::Beta(2.*radius*(j-steps*0.5)/(steps*0.9),Stash::Angle::Degrees);

      Stash::Coordinate test(l,b,GetSystem());

      if (reg.Contains(test) && (Contains(test))) return true;

    }

  }

  return false;

}

/**

 * Gives an approximation of the area using Monte-Carlo integration.

 * This is useful when you have a Compound region with overlapping

 * sub-regions, which would not be properly taken into account in

 * GetArea(), which simply sums the area of all sub-regions.  This

 * method is much slower than GetArea, however.

 *

 * \param samples_per_square_deg: approximate sample density to use

 *

 * \todo: this can be optimized to use the region's bounding box, but

 * bounding-boxes needs to be well-tested first.

 *

 * \note this function is used for testing purposes when debugging

 * regions and the Contains() function. It's not intended for

 * production use. Generally one should use the GetArea() function,

 * which does it analytically.

 */

double 

Tools::Region::

GetMonteCarloArea(int samples_per_square_deg) const {

  // here we use a bounding box that is the size of the containment

  // radius (just to be safe). Really the Bounding box of the region

  // could be used, but in some cases that hasn't been tested.

  double radius = GetRadialBound();

  Tools::Box box( GetTrueCentre(), radius*2, radius*2 );

  double thrownarea = box.GetArea();

  int nsamples = samples_per_square_deg * thrownarea;

  DEBUG_OUT << "Thrown Area: " << thrownarea << endl;

  DEBUG_OUT << "    Samples: " << nsamples << endl;

  // generate random samples over the box (defined as [-R,R] from the region center)

  int count = 0;

  for (int ii=0; ii < nsamples; ii++) {

    if (this->Contains( box.GetRandomSample() ))

        count++;

  }

  DEBUG_OUT << " Contained:" << count << endl;

  return thrownarea*count/double(nsamples);

}

/**

 * Generates the zenith or azimuth angle distribution of a given region over the

 * specified time range.

 *

 * \param hist - histogram to fill 

 * \param starttime - starting time of the run

 * \param endtime - ending time of the run

 * \param fillZenith - decide whether the zenith (true, default) or azimuth (false) angle will be filled

 *

 * \todo Document the sampling method implemented here!

 *

 */

void Tools::Region::FillZenithOrAzimuthDistribution( TH1 &hist, 

                                                     const Sash::Time starttime,

                                                     const Sash::Time endtime,

                                                     bool fillZenith ) const

{

  if( fillZenith ){

    DEBUG_OUT << "Filling zenith distribution..." << std::endl;

  }

  else{

    DEBUG_OUT << "Filling azimuth distribution..." << std::endl;

  }

  Sash::HESSArray *hess = &Sash::HESSArray::GetHESSArray();

  const Stash::System* altAzSys = hess->GetAltAzSystem().GetSystem();

  double radius = GetRadialBound();

  double rl = radius/GetTrueCentre().GetBeta().Cos();

  int tsteps = GetNumTimeSteps();

  int steps = GetNumSpaceSteps(radius);

  Sash::TimeDiff ttotal(endtime-starttime);

  double totsec = (double) ttotal;

  Sash::TimeDiff tdelta( totsec/((double) tsteps) );

  // store the HESS time for later resetting

  Sash::Time oldtime = hess->GetCurrentTime();

  for ( int t = 0; t<tsteps; t++) {

    Sash::Time tnow = starttime + Sash::TimeDiff(t*(double)tdelta);

    hess->SetCurrentTime( tnow );

    for (int i=0; i<steps; ++i) {

      Stash::Lambda l = GetTrueCentre().GetLambda() 

        + Stash::Lambda(2.*rl*(i-steps*0.5)/(steps*0.9),Stash::Angle::Degrees);

      for (int j=0; j<steps; ++j) {

        Stash::Beta b = GetTrueCentre().GetBeta() 

          + Stash::Beta(2.*radius*(j-steps*0.5)/(steps*0.9),

                        Stash::Angle::Degrees);

        Stash::Coordinate test(l,b,GetSystem());

        if (Contains(test)) {

          if( fillZenith ){

            // zenith

            double z = 90. - test.GetCoordinate(*altAzSys).GetBeta().GetDegrees();

            hist.Fill(z);}

          else{

            // azimuth

            double a = test.GetCoordinate(*altAzSys).GetLambda().GetDegrees();

            hist.Fill(a);

          }

        }

      }

    }

  }

  if (hist.GetEntries()<5) {

    WARN_OUT << "less than 5 entries used to evaluate distribution "

             << "for region " << GetName() << endl;

  }

  hess->SetCurrentTime( oldtime );

}

/**

 * Generates the zenith angle / offset to obs. pos 2D - distribution

 * of a given region over the specified time range.

 *

 * \param hist - histogram to fill 

 * \param starttime - starting time of the run

 * \param endtime - ending time of the run

 * \param psicut - will be ignored, if set to -1.

 *

 */

void 

Tools::Region::

FillZenithOffsetDistribution(const Stash::Coordinate& refpos,

                             TH2 &hist,

                             const Sash::Time starttime,

                             const Sash::Time endtime,

                             double psicut ) const

{

  DEBUG_OUT << "refpos: " << refpos << std::endl;

  DEBUG_OUT << "starttime: " << starttime << std::endl;

  DEBUG_OUT << "endtime: " << endtime << std::endl;

  DEBUG_OUT << "psicut: " << psicut << std::endl;

  Sash::HESSArray *hess = &Sash::HESSArray::GetHESSArray();

  const Stash::System* altAzSys = hess->GetAltAzSystem().GetSystem();

  Stash::Coordinate pos = refpos.GetCoordinate(*(GetSystem()));

  double radius = GetRadialBound();

  double rl = radius/GetTrueCentre().GetBeta().Cos();

  int steps = GetNumSpaceSteps(radius);

  int tsteps = GetNumTimeSteps();

  DEBUG_OUT << "SPACE STEPS:" << steps << endl;

  DEBUG_OUT << " TIME STEPS:" << tsteps << endl;

  Sash::TimeDiff ttotal(endtime-starttime);

  double totsec = (double) ttotal;

  Sash::TimeDiff tdelta( totsec/((double) tsteps) );

  // store the HESS time for later resetting

  Sash::Time oldtime = hess->GetCurrentTime();

  DEBUG_OUT << "pos: " << pos << std::endl;

  DEBUG_OUT << "radius: " << radius << std::endl;

  DEBUG_OUT << "rl: " << rl << std::endl;

  DEBUG_OUT << "tdelta: " << tdelta << std::endl;

  for ( int t = 0; t<tsteps; t++) {

    Sash::Time tnow = starttime + Sash::TimeDiff(t*(double)tdelta);

    hess->SetCurrentTime( tnow );

    DEBUG_OUT_L(10) << "tnow: " << tnow << std::endl;

    for (int i=0; i<steps; ++i) {

      Stash::Lambda l = GetTrueCentre().GetLambda() 

        + Stash::Lambda(2.*rl*(i-steps*0.5)/(steps*0.9),Stash::Angle::Degrees);

      for (int j=0; j<steps; ++j) {

        Stash::Beta b = GetTrueCentre().GetBeta() 

          + Stash::Beta(2.*radius*(j-steps*0.5)/(steps*0.9),

                        Stash::Angle::Degrees);

        Stash::Coordinate test(l,b,GetSystem());

        DEBUG_OUT_L(10) << "test: " << test << std::endl;

        if (Contains(test)) {

          double z = 90. - test.GetCoordinate(*altAzSys).GetBeta().GetDegrees();

          double d = pos.GetAngularDistance(test).GetDegrees();

          DEBUG_OUT_L(10) << "z = " << z << ", d = " << d << std::endl;

          // ignore psi, if it is negativ

          if (psicut < 0. || d <= psicut) {

            DEBUG_OUT_L(10) << "Filling hist" << std::endl;

            hist.Fill(z,d);

          }

        }

      }

    }

  }

  hess->SetCurrentTime( oldtime );

}

/**

 * Fills histogram off 2d offsets in the sky from 

 *

 * \param psicut - will be ignored, if set to -1.

 */

void Tools::Region::Fill2AxesDistanceDistribution(const Stash::Coordinate& 

                                                  refpos,

                                                  TH2 &hist2, 

                                                  double psicut ) const

{

  Stash::Coordinate pos = refpos.GetCoordinate(*(GetCentre().GetSystem() ) );

  DEBUG_OUT_L(2) << "pos " << pos

                 << "refpos " << refpos

                 << std::endl;

  double radius = GetRadialBound();

  int steps = GetNumSpaceSteps(radius);

  double rl = radius/GetTrueCentre().GetBeta().Cos();

  double refL = pos.GetLambda().GetDegrees();

  double refB = pos.GetBeta().GetDegrees();

  double cenL = GetTrueCentre().GetLambda().GetDegrees();

  double cenB = GetTrueCentre().GetBeta().GetDegrees();

  if (refL < -180.) refL += 360.;

  if (refL > 180.) refL -= 360.;

  if (refB < -180.) refB += 360.;

  if (refB > 180.) refB -= 360.;

  if (cenL < -180.) cenL += 360.;

  if (cenL > 180.) cenL -= 360.;

  if (cenB < -180.) cenB += 360.;

  if (cenB > 180.) cenB -= 360.;

  double DL = cenL - refL;

  double DB = cenB - refB;

  for (int i=0; i<steps; ++i) {

    double deltL = 2.*rl*(i-steps*0.5)/(steps*0.9);

    Stash::Lambda l = GetTrueCentre().GetLambda() 

      + Stash::Lambda(deltL,Stash::Angle::Degrees);

    for (int j=0; j<steps; ++j) {

      double deltB = 2.*radius*(j-steps*0.5)/(steps*0.9);

      Stash::Beta b = GetTrueCentre().GetBeta() 

        + Stash::Beta(deltB, Stash::Angle::Degrees);

      Stash::Coordinate test(l,b,GetCentre().GetSystem());

      if (Contains(test)) {

        double d = pos.GetAngularDistance(test).GetDegrees();

        // ignore psi, if it is negative

        if (psicut < 0. || d <= psicut) {

          hist2.Fill(deltL+DL, deltB+DB);

        }

      }

    }

  }

}

/**

 * Returns a Tools::Box that bounds the region

 */

Tools::Box*

Tools::Region::

GetBoundingBox() const {

  return new Tools::Box(GetTrueCentre(), GetBound1(), GetBound2(), 

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

                        string(GetName())+"_bbox");

}

Tools::Circle*

Tools::Region::

GetBoundingCircle() const {

  return new Tools::Circle(GetTrueCentre(), GetRadialBound(), string(GetName())+"_bcirc");

}

/////////////////////////////////////////////////////////////////////////////

// CircleSegment ////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////////////////

Tools::CircleSegment::CircleSegment(const Stash::Coordinate& cor,

                                    double r1,double r2,

                                    double phi1, double phi2,

                                    string name):

  Region(cor,name,"SEGMENT", Stash::Lambda(0)),

  fR1(r1),fR2(r2),

  fPhi1(phi1), fPhi2(phi2), fLine(0)

{

  if( fPhi1 < 0 ) {

    WARN_OUT << "please define phi1 within the boundaries [0;360]!" << std::endl;

    fPhi1 += 360.;

  }

  else if( fPhi1 > 360. ) {

    WARN_OUT << "please define phi1 within the boundaries [0;360]!" << std::endl;

    fPhi1 -= 360.;

  }

  if( fPhi2 < 0 ) {

    WARN_OUT << "please define phi1 within the boundaries [0;360]!" << std::endl;

    fPhi2 += 360.;

  }

  else if( fPhi2 > 360. ) {

    WARN_OUT << "please define phi1 within the boundaries [0;360]!" << std::endl;

    fPhi2 -= 360.;

  }

}

Tools::CircleSegment::~CircleSegment()

{

  if(fLine) delete fLine;

}

double Tools::CircleSegment::GetArea() const {

  double area = 0;

  if( fabs(fPhi2-fPhi1) > 1e-10 ) {

    bool overzero = (fPhi2 < fPhi1) ? (true) : (false);

    if( overzero )

      area = (360. - fPhi1 + fPhi2) / 360.;

    else

      area = (fPhi2 - fPhi1) / 360.;

    area *= M_PI * (pow(fR2,2.) - pow(fR1,2.)); 

  }

  return area;

}

bool Tools::CircleSegment::IsInside(const Stash::Coordinate& dir) const

{

  double r = GetCentre().GetAngularDistance(dir).GetDegrees();

  DEBUG_OUT_L(2) << "Pos: (" << dir.GetLambda().GetStringInHours() << "," 

                 << dir.GetBeta() << "), Centre: (" 

                 << fCentre.GetLambda().GetStringInHours() 

                 << "," << GetCentre().GetBeta() << "), r " << r << " [" << fR1 

                 << "," << fR2 << "]"  << std::endl;

  if ((r>fR2) || (r<fR1)) {

    DEBUG_OUT_L(2) << "--> Not Contained" << endl;

    return false;

  }

  double dl = dir.GetLambda().GetDegrees() - fCentre.GetLambda().GetDegrees();

  double sign = dl/fabs(dl+1e-10); // need epsilon to avoid NaN, which

                                   // get marked as Contained!

  Stash::Coordinate ldl = dir + Stash::Lambda(dl,Stash::Angle::Degrees);

  double dx = sign * dir.GetAngularDistance(ldl).GetDegrees();

  double dy = (dir.GetBeta() - fCentre.GetBeta()).GetDegrees();

  Stash::Lambda ang(atan2(dx,dy), Stash::Angle::Radians);

  double phi = ang.GetDegrees();

  DEBUG_OUT_L(2) << " phi " << phi << " [" << fPhi1 << "," << fPhi2 << "]" 

                 << std::endl;

  bool overzero = (fPhi2 < fPhi1) ? (true) : (false);

  if( overzero ) {

    if( phi > fPhi2 && phi < fPhi1 )

      return false;

  }

  else {

    if( phi < fPhi1 || phi > fPhi2 )

      return false;

  }

  return true;

}

/**

 * Grow the boundary of the CircleSegment by a given amount of space angle.

 *

 * The CircleSegment is grown both towards its inner and outer boundary, and 

 * the azmuthal coverage is enlarged as well.

 *

 */

Tools::CircleSegment* 

Tools::CircleSegment::

GrowBoundary(const double degrees,

             const string name) const

{

  string grownName(GetName());

  grownName += name;

  Tools::CircleSegment* grownReg 

    = (Tools::CircleSegment*)this->Clone(grownName.c_str());

  if( degrees > 0 ) {

    // grow radially

    grownReg->SetR1( ((fR1-degrees) < 0) ? (0) : (fR1-degrees) );

    grownReg->SetR2( fR2 + degrees );

    // grow azimuthally (into full ring if necessary)

    double newPhi1, newPhi2;

    newPhi1 = fPhi1 - degrees 

      / ((fR1 + fR2) / 2 ) * 180./TMath::Pi();

    newPhi2 = fPhi2 + degrees 

      / ((fR1 + fR2) / 2 ) * 180./TMath::Pi();

    // bring angles back into boundaries

    // both angles crossed zero, so full circle

    if( newPhi1 < 0 && newPhi2 > 360. ) {

      newPhi1 = 0;

      newPhi2 = 360.;

    }

    // ordinary back-to-boundary

    else { 

      if( newPhi1 < 0 ) {

        newPhi1 += 360.;

      }

      if( newPhi2 > 360.) {

        newPhi2 -= 360.;

      }

    }

    // check if we have grown into a full ring

    bool wasOverZero = (fPhi2 < fPhi1) ? (true) : (false);

    bool isOverZero  = (newPhi2 < newPhi1) ? (true) : (false);

    if( wasOverZero && newPhi2 > newPhi1 ) {

      newPhi1 = 0.;

      newPhi2 = 360.;

    }

    else if( isOverZero && newPhi2 > newPhi1 ) {

      newPhi1 = 0.;

      newPhi2 = 360.;

    }

    grownReg->SetPhi1( newPhi1 );

    grownReg->SetPhi2( newPhi2 );

  }

  return grownReg;

}

Tools::Region::VertexList 

Tools::CircleSegment::

AsVertexList() const {

  Tools::Region::VertexList vlist;

  vlist.resize(1);

  const int NPOINTS = 50;

  Stash::Lambda step(0);

  bool overzero = (fPhi2 < fPhi1) ? (true) : (false);

  if( overzero )

    step = Stash::Lambda((360.-fPhi1+fPhi2)/NPOINTS, Stash::Angle::Degrees);

  else

    step = Stash::Lambda((fPhi2-fPhi1)/NPOINTS, Stash::Angle::Degrees);

  DEBUG_OUT << "STEP: " << step.GetDegrees() << " for "<< AsFITSRegionString() << endl;

  Stash::Beta r1(fR1, Stash::Angle::Degrees);

  Stash::Beta r2(fR2, Stash::Angle::Degrees);

  Stash::Lambda phi1(fPhi1, Stash::Angle::Degrees);

  Stash::Lambda phi2(fPhi2, Stash::Angle::Degrees);

  // Inner arc:

  if (fR1 > 1e-10) {

    Stash::Lambda phi = phi1;

    for (int ii=0; ii<NPOINTS; ii++) {

      Stash::Coordinate point = InverseTransform( GetPointAtRadius(r1, phi) );

      vlist[0].push_back( point );

      DEBUG_OUT_L(2) << "INNER Added point " << point << endl;

      phi+=step;

    }

  }

  // first connector:

  Stash::Coordinate p1 = InverseTransform( GetPointAtRadius( r1,phi2 ) );

  Stash::Coordinate p2 = InverseTransform( GetPointAtRadius( r2,phi2 ) );

  vlist[0].push_back( p1 );

  vlist[0].push_back( p2 );

  // Outer arc:

  if (fR2 > 1e-10) {

    Stash::Lambda phi = phi2;

    for (int ii=0; ii<NPOINTS; ii++) {

      Stash::Coordinate point = InverseTransform( GetPointAtRadius(r2, phi) );

      vlist[0].push_back( point );

      DEBUG_OUT_L(2) << "OUTTER Added point " << point << endl;

      phi-=step;

    }

  }

  // second connector:

  Stash::Coordinate p3 = InverseTransform( GetPointAtRadius( r2,phi1 ) );

  Stash::Coordinate p4 = InverseTransform( GetPointAtRadius( r1,phi1 ) );

  vlist[0].push_back( p3 );

  vlist[0].push_back( p4 );

  return vlist;

}

/**

 * \returns a point on a circle with given radius from the Centre

 * position. Used by CircleSegment::AsVertexList.

 *

 * \note rad does need to be a Stash::Beta (not Lambda) in order for

 * the "RA factor" to be taken into account properly)

 */

Stash::Coordinate

Tools::Region::

GetPointAtRadius( Stash::Beta rad, Stash::Lambda phi ) const {

  Stash::Coordinate point = fCentre + rad;

  Stash::DirectionVector dv = point.GetDirectionVector();

  Stash::DirectionVector cendv =  fCentre.GetDirectionVector();

  Stash::DirectionVector newdv = dv.RotateAroundPoint( cendv, phi.GetNegative() );

  return Stash::Coordinate(newdv, GetSystem());

}

/**

 * Returns a point on the line between A and B, where t is the parameter in range [0,1]

 */

// Stash::Coordinate 

// Tools::Region::

// GetPointOnLine( Stash::Coordinate A, Stash::Coordinate B, double t) {

//   //  return (1.0-t)*A + t*B;

// }

string

Tools::CircleSegment::

AsFITSRegionString() const {

  char rep[512];

  Stash::Coordinate realcen = GetTrueCentre();

  string sys = GetFITSSystemName();

  sprintf( rep, "%s; panda(%.4fd,%.4fd,%.4fd,%.4fd,1,%.4fd,%.4fd,1) # text={%s} ",

           sys.c_str(), 

           realcen.GetLambda().GetDegrees(), 

           realcen.GetBeta().GetDegrees(),

           GetPhi1() + GetRotationAngle().GetDegrees(),  

           GetPhi2() + GetRotationAngle().GetDegrees(),

           GetR1() / realcen.GetBeta().Cos(), GetR2() / realcen.GetBeta().Cos(),

           GetName()

           );

  return string(rep);

}

bool 

Tools::Ring::

IsInside(const Stash::Coordinate& dir) const {

  double rad=0;

  rad = fCentre.GetAngularDistance( dir ).GetDegrees();

  if (rad < GetR1() || rad > GetR2()) {

    return false;

  }

  else return true;

}

string

Tools::Ring::

AsFITSRegionString() const {

  char rep[256];

  Stash::Coordinate realcen = GetTrueCentre();

  string sys = GetFITSSystemName();

  sprintf( rep, "%s; annulus(%.4f,%.4f,%.4f,%.4f) # text={%s}", sys.c_str(), 

           realcen.GetLambda().GetDegrees(), realcen.GetBeta().GetDegrees(),

           GetR1(), GetR2(), GetName());

  return string(rep);

}

bool 

Tools::Circle::

IsInside(const Stash::Coordinate& dir) const {

  double rad=0;

  rad = fCentre.GetAngularDistance( dir ).GetDegrees();

  if (rad > GetR2()) {

    return false;

  }

  else return true;

}

string

Tools::Circle::

AsFITSRegionString() const {

  char rep[256];

  Stash::Coordinate realcen = GetTrueCentre();

  string sys = GetFITSSystemName();

  sprintf( rep, "%s; circle(%.4fd,%.4fd,%.4fd) # text={%s}", sys.c_str(), 

           realcen.GetLambda().GetDegrees(), realcen.GetBeta().GetDegrees(),

           GetR2()/realcen.GetBeta().Cos(), GetName());

  return string(rep);

}

/////////////////////////////////////////////////////////////////////////////

// BOX //////////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////////////////

Tools::Box::~Box()

{ }

Tools::Box::Box(const Stash::Coordinate& cor,

                double l1,double l2, Stash::Lambda phi,

                string name) :

  Region(cor,name,"BOX",phi),fL1(l1),fL2(l2)

{}

bool Tools::Box::IsInside(const Stash::Coordinate& dir) const

{

  double dist = fCentre.GetAngularDistance(dir).GetDegrees();

  if(dist > fL1 && dist > fL2) {

    DEBUG_OUT_L(2) << "--> Not Contained" << endl;

    return false;

  }

  double dl=0,db=0;

  GetOffsetOfPointFromCenter( dir, dl, db );

  if (dl < fL1/2.0 && db <fL2/2.0)

    return true;

  return false;

}

double 

Tools::Box::

GetBound1() const

{ 

  double dx,dy;

  Tools::Region::VertexList vlist = AsVertexList();

  std::vector< Stash::Coordinate >::iterator ivert;

  double max = -1000;

  for (ivert=vlist[0].begin(); ivert!=vlist[0].end(); ivert++) {

    GetOffsetOfPointFromCenter( *ivert, dx, dy, false );

    if (fabs(dx) >max) max=fabs(dx);

  }

  return 2.0*max;

}

double 

Tools::Box::

GetBound2() const

{ 

  double dx,dy;

  Tools::Region::VertexList vlist = AsVertexList();

  std::vector< Stash::Coordinate >::iterator ivert;

  double max = -1000;

  for (ivert=vlist[0].begin(); ivert!=vlist[0].end(); ivert++) {

    GetOffsetOfPointFromCenter( *ivert, dx, dy, false );

    if (fabs(dy) >max) max=fabs(dy);

  }

  return 2.0*max;

}

/**

 * Grow the boundary of the box by a given amount of space angle.

 *

 * Grows both the width and length of the box.

 */

Tools::Box*

Tools::Box::

GrowBoundary(const double degrees, const string name) const

{

  string grownName(GetName());

  grownName += name;

  Tools::Box* grownReg = (Tools::Box*)this->Clone(grownName.c_str());

  if( degrees > 0 ) {

    grownReg->SetL1( fL1 + 2*degrees );

    grownReg->SetL2( fL2 + 2*degrees );

  }

  return grownReg;

}

string

Tools::Box::

AsFITSRegionString() const {

  char rep[256];

  Stash::Coordinate realcen = GetTrueCentre();

  string sys = GetFITSSystemName();

  sprintf( rep, "%s; box(%.4fd,%.4fd,%.4fd,%4fd,%4fd) # text={%s}",  sys.c_str(), 

           realcen.GetLambda().GetDegrees(), realcen.GetBeta().GetDegrees(),

           fL1, fL2, GetRotationAngle().GetDegrees(), GetName()

           );

  return string(rep);

}

/**

 * \todo: also step along each axis and insert multiple points (for

 * better plotting on distorted projections)

 */

Tools::Region::VertexList 

Tools::Box::

AsVertexList() const {

  Tools::Region::VertexList vlist;

  vlist.resize(1);

  Stash::Coordinate p1 = GetPointAtOffsetFromCenter( fL1/2.0, fL2/2.0 );

  Stash::Coordinate p2 = GetPointAtOffsetFromCenter( fL1/2.0, -fL2/2.0 );

  Stash::Coordinate p3 = GetPointAtOffsetFromCenter( -fL1/2.0, -fL2/2.0 );

  Stash::Coordinate p4 = GetPointAtOffsetFromCenter( -fL1/2.0, fL2/2.0 );

  vlist[0].push_back( InverseTransform(p1) );

  vlist[0].push_back( InverseTransform(p2) );

  vlist[0].push_back( InverseTransform(p3) );

  vlist[0].push_back( InverseTransform(p4) );

  vlist[0].push_back( vlist[0][0] ); // close the curve

  return vlist;

}

/**

 * \returns a point at a given offset in "sky degrees" (e.g. not RA/Dec degrees).

 */

Stash::Coordinate 

Tools::Region::

GetPointAtOffsetFromCenter( double offsX, double offsY ) const {

  // could do this simpler with 3D direction vectors instead of a cos(dec)

  // factor, but this works fine (not easy to do it with how Stash

  // is set up though)

  Stash::Coordinate cen = fCentre;

  double dist = sqrt(pow(offsX,2) + pow(offsY,2));

  Stash::Beta db( offsY, Stash::Angle::Degrees );

  Stash::Lambda dl( offsX/cen.GetBeta().Cos(), Stash::Angle::Degrees );

  Stash::Coordinate newcoord = cen + dl + db;

  DEBUG_OUT << "Distance Error: " 

            <<newcoord.GetAngularDistance(fCentre).GetDegrees() - dist << endl;

  double dx,dy;

  GetOffsetOfPointFromCenter( newcoord, dx, dy, false );

  DEBUG_OUT << "error="<< fabs(offsX)-fabs(dx) << " offsX=" << offsX << " dx="<<dx << endl;

  DEBUG_OUT << "error="<< fabs(offsY)-fabs(dy) << " offsY=" << offsY << " dy="<<dx << endl;

  return newcoord;

}

/**

 * \returns offset in sky degrees of point from the region center (similar to GetAngularDistance(), but returning both directions)

 *

 * \param untransformed - return the offset from the untransformed center point

 */

void 

Tools::Region::

GetOffsetOfPointFromCenter( Stash::Coordinate point, double &offsX, double &offsY,

                            bool untransformed ) const  {

  Stash::Coordinate cen = fCentre;

  if (untransformed)

    cen = GetTrueCentre();

  Stash::Coordinate px( point.GetLambda(), cen.GetBeta(), cen.GetSystem() );

  Stash::Coordinate py( cen.GetLambda(), point.GetBeta(), cen.GetSystem());

  offsX = px.GetAngularDistance( cen ).GetDegrees();

  offsY = py.GetAngularDistance( cen ).GetDegrees();

}

Stash::Coordinate 

Tools::Box::

GetRandomSample() {

  double dx = gRandom->Uniform( -fL1/2.0, fL1/2.0 );

  double dy = gRandom->Uniform( -fL2/2.0, fL2/2.0 );

  return InverseTransform(GetPointAtOffsetFromCenter( dx, dy ));

}

/////////////////////////////////////////////////////////////////////////////

// Compound /////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////////////////

/**

 * Returns total area of compound region.

 *

 * \note this area is not correct if regions overlap!

 */

double 

Tools::Compound::

GetArea() const {

  RegionList::const_iterator region;

  double area=0.0;

  for (region = fRegionList.begin(); region != fRegionList.end(); region++) {

    area += (*region)->GetArea();

  }

  return area;

}

/**

 * Creates a clone of the compound, using ROOT's Clone() mechanism.

 * Additionally, clone all Regions inside the compound.

 */

Tools::Compound*

Tools::Compound::

Clone(const char* name) const

{

  Tools::Compound* clone = (Tools::Compound*)Tools::Region::Clone(name);

  clone->ClearRegionList();

  // clone each member of the compound

  RegionList::const_iterator it = fRegionList.begin();

  for( ; it != fRegionList.end(); ++it ) {

    std::string newName = (*it)->GetName();

    newName.append("_");

    newName += std::string(name);

    clone->AddRegion( (Tools::Region*)((*it)->Clone(newName.c_str())) );

  }

  return clone;

}

/**

 * \brief Move all subregions to the new position (mainting relative offsets).

 *

 * \note this function modifies the regions inside the compound! 

 */

void 

Tools::Compound::

MoveTo( Stash::Coordinate newpos ) {

  Stash::Lambda dl = newpos.GetLambda() - fCentre.GetLambda();

  Stash::Beta db = newpos.GetBeta() - fCentre.GetBeta();

  RegionList::const_iterator it = fRegionList.begin();

  for( ; it != fRegionList.end(); ++it ) {

    Stash::Coordinate newcen = (*it)->GetCentre();

    newcen += dl;

    newcen += db;

    const_cast<Tools::Region*>(*it)->MoveTo(newcen);

  }

  fCentre = newpos;

}

/**

 * Add a region to the list of regions which make up this compound region

 *

 * \note the center of the compound will be the center of the first

 * region added. For speed, it is not updated to a center of

 * gravity. In the case where you would like it to be so, simply call:

 * \code region->SetCentre( region->GetCentreOfGravity() ) \endcode.

 */

void Tools::Compound::AddRegion( const Tools::Region *reg ) {

  fRegionList.push_back( reg );

  // Ensure that the Center point of the region is not in the

  // DummySystem, if so, use the system of the first region added as

  // the Region system (rotations will be done in that system)

  if (fCentre == Stash::Coordinate(0,0,Stash::System::GetDummySystem())) {

    DEBUG_OUT << GetName() << ": Setting the System for this compound to: " 

              << reg->GetSystem()->GetName() << endl;

    fCentre = reg->GetCentre();

    SetRotation(GetRotationAngle()); // updates the rotation center 

  }

  // the following is nice, but slow for large compounds. Instead user

  // should explicitly call compound->SetCentre(

  // compound->GetCentreOfGravity() ); else { fCentre =

  // GetCentreOfGravity(); }

}

Tools::Region::VertexList Tools::Compound::AsVertexList() const

{

  VertexList list;

  RegionList::const_iterator region;

  for (region = fRegionList.begin(); region != fRegionList.end(); region++) {

    DEBUG_OUT_L(2) << "Getting vertex list for: " << (*region)->GetName() << endl;

    VertexList regVertexList = (*region)->AsVertexList();

    VertexList::iterator it = regVertexList.begin();

    for( ; it != regVertexList.end(); ++it ){

      vector<Stash::Coordinate> clist;

      vector<Stash::Coordinate>::iterator cc;

      for (cc=(*it).begin(); cc != (*it).end(); cc++)

        clist.push_back( InverseTransform( *cc ) ); // apply the compound's transform

      list.push_back( clist );

    }

  }

  return list;

}

double Tools::Compound::GetBound1() const {

  double maxl1=0;