GammaLib

#include <algorithm>

#include <iomanip>

#define INFO_OUT_TAG "ReflectedBgMaker> " 

#define DEBUG 0

#include <utilities/debugging.hh>

#include <utilities/ConfigHandler.hh>

#include <crash/RotatedSystem.hh>

#include <crash/EquatorSystem.hh>

#include <sash/RunHeader.hh>

#include "ReflectedBgMaker.hh"

using std::cout;

using std::endl;

/**

 * Constructor

 *

 */

Background::ReflectedBgMaker::

ReflectedBgMaker(std::string name, Utilities::ConfigHandler& conf,

                 std::string input, bool init_now) :

  BgMaker(name.c_str(), conf, input, init_now)

{

  ConstructMaker(conf);

}

Background::ReflectedBgMaker::

ReflectedBgMaker(Utilities::ConfigHandler& conf,

                 std::string input, bool init_now) :

  BgMaker("ReflectedBgMaker", conf, input, init_now)

{

  ConstructMaker(conf);

}

void Background::ReflectedBgMaker::ConstructMaker(Utilities::ConfigHandler& conf)

{

  DEBUG_OUT << "Hello" << endl;

  fMaxNumberOffRegions = 100;

  fMinNumberOffRegions = 1;

  fMinOnOffDist = 0.1;

  fNumberOffRegions = 0;

  fRotationAngleDeg = 0.;

  fCurrentOffset = 0;

  fMaxOffRegionDist = 10.;

  fUseCircularOff = false; 

  std::string Name = "ReflectedBgMaker";

  //Read config options from the confighandler

  conf.SetSectionHelpText(Name, "Estimates background from a "

                          "series of on-region-shaped regions "

                          "that have the same offset from the observation "

                          "position as the test position. This BgMaker "

                          "is useful for making spectra, but is not well "

                          "suited for 2-D images." );

  conf.GetValue(Name, "RotationAngle", 

                fRotationAngleDeg, 

                "Rotation angle of the initial OFF "

                "region from the test position (deg)" );

  conf.GetValue(Name, "MaxNumberOffRegions", 

                fMaxNumberOffRegions,

                "Maximum allowed number of OFF regions");

  conf.GetValue(Name, "MinNumberOffRegions",

                fMinNumberOffRegions,

                "Minimum required number of OFF regions");

  conf.GetValue(Name, "MaxDistanceOffRegions", 

                fMaxOffRegionDist,

                "Maximum allowed distance between OFF regions and target");

  conf.GetValue(Name, "MinOnOffDistance", 

                fMinOnOffDist,

                "Minimum angular distance between the on region and "

                "any off region");

  fUseCircularOff = conf.GetFlag(Name, "UseCircularOff",

                                 "Use circular off regions "

                                 "(for backwards compatibility)" );

  // Make sure the exposure is weighted by area of the on and off regions

  ScaleExposureByRegionArea( true );

}

Background::ReflectedBgMaker::~ReflectedBgMaker() {

  DEBUG_OUT << "Goodbye!" << endl;

}

/**

 * Constructs OFF regions at the beginning of each run.

 * Calls CreateOffRegions().

 */

void  

Background::ReflectedBgMaker::

BeginRun() {

  // Look up the observation and test positions:

  // ON region gets automatically constructed by BgMaker::InitRun in

  // the Current BgStats object

  fNumberOffRegions = 0;

  const Tools::Compound& onregion ( GetBgStats()->GetOnRegion() );

  const Sash::RunHeader *runheader = fHess->Get<Sash::RunHeader>();

  if (!runheader) {

    WARN_OUT << "No runheader found" << endl;

    throw std::runtime_error( "No runheader was found!" );

  }

  const Stash::Coordinate &obspos_radec = runheader->GetObservationPosition().

    GetCoordinate(*(fHess->GetRADecJ2000System()));

  const Stash::Coordinate &testpos_radec = onregion.GetCentre().

    GetCoordinate(*(fHess->GetRADecJ2000System()));

  DEBUG_OUT << "Obspos =" << obspos_radec <<endl;

  DEBUG_OUT << "Testpos=" << testpos_radec <<endl;

  fCurrentOffset = testpos_radec.GetAngularDistance( obspos_radec ).

    GetDegrees();

  // Construct OFF regions

  Tools::Circle* circle = IsCircle(&onregion);

  if (circle)

    CreateOffRegionsForOnCircle(circle, obspos_radec);

  else

    CreateOffRegions(&onregion, obspos_radec);

  // Check that enough were found

  if ( GetBgStats()->GetOffRegion().GetNumRegions() < fMinNumberOffRegions) {

    WARN_OUT << "Run will be skipped: Too few off regions! (found "

             << GetBgStats()->GetOffRegion().GetNumRegions()

             << " / required " << fMinNumberOffRegions << " )" << std::endl;

    GetBgStats()->GetSkippedFlag() = true;   

  }

}

/**

 * Calculates exposure(On,Off) at the end of the run.

 * Sets exposureOn to alpha*NOff (with alpha=areaOn/areaOff)

 * and exposureOff to NOff

 */

void Background::ReflectedBgMaker::EndRun()

{

  DEBUG_OUT << std::endl;

  double alpha = GetBgStats()->GetOnRegion().GetArea()

    / GetBgStats()->GetOffRegion().GetArea();

  double nOff = GetBgStats()->GetRateStats().nOff;

  GetBgStats()->GetRateStats().exposureOff = nOff;

  GetBgStats()->GetRateStats().exposureOn  = alpha * nOff;

  DEBUG_OUT << "on region area:  " << GetBgStats()->GetOnRegion().GetArea() << endl;

  DEBUG_OUT << "off region area: " << GetBgStats()->GetOffRegion().GetArea() << endl;

  DEBUG_OUT << "-> alpha:        " << alpha << endl;

  DEBUG_OUT << "exposure off:    " << GetBgStats()->GetRateStats().exposureOff << endl;

  DEBUG_OUT << "-> exposure on:  " << GetBgStats()->GetRateStats().exposureOn << endl;

}

/*

 * Fill On events into OnMap for debug purposes

 */

bool Background::ReflectedBgMaker::IsOnMapEvent()

{

  DEBUG_OUT << std::endl;

  if(Background::BgMaker::IsOnEvent()) return true;

  else return false;

}

/*

 * Fill Off events into OffMap for debug purposes

 */

bool Background::ReflectedBgMaker::IsOffMapEvent()

{

  DEBUG_OUT << std::endl;

  if(Background::BgMaker::IsOffEvent()) return true;

  else return false;

}

/*

 * Leave OnExposureMap empty

 */

void Background::ReflectedBgMaker::FillOnExposure()

{

  DEBUG_OUT << std::endl;

}

/*

 * Leave OffExposureMap empty

 */

void Background::ReflectedBgMaker::FillOffExposure()

{

  DEBUG_OUT << std::endl;

}

/**

 * \brief Check if a region is a circle and returns a pointer to it.

 *

 * This is needed because circular on regions are created as

 * Tools::Compound objects containing one Tools::Circle,

 * and Tools::Compound is not a superclass of Tools::Circle.

 *

 * In that case a pointer to that Tools::Circle in the compound is returned.

 */

Tools::Circle* Background::ReflectedBgMaker::IsCircle(const Tools::Region* region) {

  // Is the type Tools::Circle?

  DEBUG_OUT << "region type: " << typeid(*region).name() << endl;

  if (typeid(*region) == typeid(Tools::Circle))

    return (Tools::Circle*)region;

  // Is the type Tools::Compound, but it contains exactly one circle?

  if (typeid(*region) == typeid(Tools::Compound)) {

    Tools::Compound* compound = (Tools::Compound*)region;

    DEBUG_OUT << "compound type: " << typeid(*compound).name() << endl;

    DEBUG_OUT << "NumRegions: " << compound->GetNumRegions() << endl;

    if (compound->GetNumRegions() == 1) {

      const Tools::Region* first_region = compound->GetRegions()[0];

      DEBUG_OUT << "First region type: " << typeid(*first_region).name() << endl;

      if (typeid(*first_region) == typeid(Tools::Circle))

        return (Tools::Circle*)first_region;

    }

  }

  return 0;

}

/**

 * \brief Create off regions on a circle around the observation position.

 *

 * OFF regions are found by incrementally rotating the ON region.

 * For each OFF region candidate, it is checked whether it overlaps with

 * an already existing OFF region or an exclusion region. To exclude any

 * potential contamination of OFF regions by signal, an additional safety

 * margin of fMinOnOffDist is required between the ON region and

 * any OFF region.

 *

 * OFF regions are placed in the so-called TangentialRADecSystem, which

 * is a Cartesian projection of the RADec system at the observation position.

 * This approach avoids a potential bias of the OFF area due to projection effects.

 *

 * The OFF regions are added to Offregions compound in CurrentBgStats

 *

 * \param onregion - Compound ON region

 * \param obspos   - Observation position in Tangential system

 *

 * \todo Make sure this also works for Tools::Compund ON regions

 */

bool

Background::ReflectedBgMaker::

CreateOffRegions(const Tools::Region* onregion,

                 const Stash::Coordinate &obspos_radec)

{

  INFO_OUT << "On region is not a circle. "

    "Creating off regions using the general algorithm." << endl;

  fNumberOffRegions = 0;

  // Skip run if the ON region contains the observation position

  if( onregion->Contains(obspos_radec) ) {

    WARN_OUT << "ON region contains the observation position." << std::endl;

    return false;

  }

  // test position

  const Stash::Coordinate testpos_radec = onregion->GetCentre().

    GetCoordinate(*fHess->GetRADecJ2000System());

  // Enlarge a copy of the on region by fMinOnOffDist. This copy is then

  // used to check for overlap with potential off regions, and thus

  // guarantees a proper spacing between the on region and any off

  Tools::Region *grownOn

    = onregion->GrowBoundary(fMinOnOffDist);

  DEBUG_OUT << "on region: " << *onregion << std::endl;

  DEBUG_OUT << "grown on region: " << *grownOn << std::endl;

  // Create off region only once, for reasons of speed,

  // and later only modify rotation angle and position

  Tools::Region* tmpOff = 0;

  if( fUseCircularOff ) {

    double off_radius = onregion->GetEquivalentCircleRadius();

    DEBUG_OUT << "using circular off regions with radius radius " 

              << off_radius << " deg" << std::endl;

    tmpOff = new Tools::Circle( testpos_radec, off_radius, "Off");

  }

  else {

    DEBUG_OUT << "using on-shaped off regions" << std::endl;

    tmpOff = (Tools::Region*)onregion->Clone("Off");

  }

  DEBUG_OUT << *tmpOff << std::endl;

  // Brute-force off region finding.

  // Rotate on region around observation position.

  // step size is adapted according to the ratio of equivalent circle radius

  // and offset.

  //

  // Slow, but reliable for even strangly shaped regions.

  double radius = onregion->GetEquivalentCircleRadius();

  double offset = testpos_radec.GetAngularDistance(obspos_radec).GetDegrees();

  double stepsize = asin(radius/offset) * 180.0 / M_PI / 5.;

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

  DEBUG_OUT << "offset:                   " << offset << std::endl;

  DEBUG_OUT << "stepsize:                 " << stepsize << std::endl;

  if( stepsize != stepsize ) {

    WARN_OUT << "Offset too small (" << offset 

             << "deg), no Off regions could be found." << std::endl;

    return false;

  }

  Stash::Lambda rotAngle(fRotationAngleDeg, Stash::Angle::Degrees);

  Stash::SmallAngle dRot(stepsize, Stash::Angle::Degrees);

  int numSteps = 360. / stepsize;

  do {

    tmpOff->SetRotation( rotAngle, obspos_radec );

    DEBUG_OUT << "Trying to find off region " << *tmpOff

              << " at rotation angle " 

              << rotAngle << std::endl;

    // add the new region only if it does not overlap with any other

    // exclusion, on, or off region. Leave fMinOnOffDist space

    // between the on and the off regions (use the grown on region)

    if( !(tmpOff->Overlaps(GetBgStats()->GetOffRegion()) )) {

      if( !(tmpOff->Overlaps(GetBgStats()->GetOnExclusionRegion()) )) {

        if( !tmpOff->Overlaps( *grownOn ) ) {

          // add off region to bgstats

          std::ostringstream offName;

          offName << std::fixed << std::setprecision(1)

                  << rotAngle.GetDegrees() << "deg";

          Tools::Region* newOff 

            = (Tools::Region*)tmpOff->Clone( offName.str().c_str() );

          GetBgStats()->GetOffRegion().AddRegion( newOff );

          fNumberOffRegions++;

          INFO_OUT << "Added new off region at "<< rotAngle.GetDegrees() << " deg from North at "

                   << newOff->GetTrueCentre().GetCoordinate(*fHess->GetRADecJ2000System()) << std::endl;

          DEBUG_OUT << "offset from observation position: " 

                    << newOff->GetTrueCentre().GetCoordinate(*fHess->GetRADecJ2000System())

            .GetAngularDistance( obspos_radec ).GetDegrees() << " deg" << std::endl;

          if( fNumberOffRegions == fMaxNumberOffRegions ) {

            INFO_OUT << "Reached maximum allowed number of OFF regions ("

                     << fMaxNumberOffRegions << "). Stopping creating new "

                     << "ones at your request." << std::endl;

            break;

          }

        }

      }

    }

    rotAngle += dRot;

  } while( --numSteps != 0 );

  delete grownOn;

  delete tmpOff;

  if ( fNumberOffRegions == 0) {

    WARN_OUT << "No Off regions found." << std::endl;

    return false;

  } 

  INFO_OUT << "Added "<< fNumberOffRegions

           << " OFF regions at offset " << fCurrentOffset << endl;

  return true;

}

/**

 * \brief Create off regions for a circular on region.

 *

 * The off regions are stored in CurrentBgStats.

 *

 * \param onregion - Compound ON region

 * \param obspos   - Observation position in Tangential system

 *

 * \note This is the common case that warrants a fast implementation,

 * the general algorithm is in CreateOffRegions().

 */

bool

Background::ReflectedBgMaker::

CreateOffRegionsForOnCircle(const Tools::Circle* onregion,

                            const Stash::Coordinate &obspos_radec)

{

  INFO_OUT << "On region is a circle. "

    "Creating off regions using the fast algorithm." << endl;

  DEBUG_OUT << "on region: " << *onregion << std::endl;

  fNumberOffRegions = 0;

  // Skip run if the ON region contains the observation position

  if( onregion->Contains(obspos_radec) ) {

    WARN_OUT << "ON region contains the observation position." << std::endl;

    return false;

  }

  // test position

  const Stash::Coordinate onpos_radec = onregion->GetCentre().

    GetCoordinate(*fHess->GetRADecJ2000System());

  DEBUG_OUT << "onpos_radec: " << onpos_radec << std::endl;

  // Create off region only once, for reasons of speed,

  // and later only modify rotation angle and position

  Tools::Region* tmpOff = (Tools::Circle*)onregion->Clone("Off");

  DEBUG_OUT << "off region: " << *tmpOff << std::endl;

  DEBUG_OUT << "--- OnExclusionRegion START --- " << endl

            << GetBgStats()->GetOnExclusionRegion() << endl

            << "--- OnExclusionRegion END --- " << endl;

  // Compute the relevant numbers needed for computing positions

  // segment = angle a test circle spans on the circle of radius

  // offset around the observation position

  double radius = onregion->GetR2();

  double offset = onpos_radec.GetAngularDistance(obspos_radec).GetDegrees();

  double segment = 2 * asin(radius/offset) * 180.0 / M_PI;

  // \note This scheme of deciding on the number of test positions

  // was chosen to match the general algorithm in CreateOffRegions()

  int steps_per_segment = 10;

  double stepsize = segment / steps_per_segment;

  int n_test_positions = 360 / segment * steps_per_segment;

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

  DEBUG_OUT << "offset:   " << offset << std::endl;

  DEBUG_OUT << "segment: " << segment << std::endl;

  DEBUG_OUT << "steps_per_segment: " << steps_per_segment << std::endl;

  DEBUG_OUT << "stepsize: " << stepsize << std::endl;

  DEBUG_OUT << "n_test_positions: " << n_test_positions << std::endl;

  // Loop test positions from 0 to 360 deg

  double rot_angle = 0;

  do {

    // Rotate off region to current rot_angle

    Stash::Lambda rot_angle_lambda(rot_angle, Stash::Angle::Degrees);

    tmpOff->SetRotation( rot_angle_lambda, obspos_radec );

    // Compute the new position and offset to on center position

    const Stash::Coordinate tmpOffpos_radec = tmpOff->GetTrueCentre().

      GetCoordinate(*fHess->GetRADecJ2000System());

    double on_off_dist = tmpOffpos_radec.GetAngularDistance( onpos_radec ).GetDegrees();

    double off_radius = tmpOff->GetEquivalentCircleRadius();

    // The criterion is that the on and off circle must have a separation

    // of at least fMinOnOffDist, i.e. that the centers of the on and off

    // circles have a separation < allowed_dist defined like this:

    double allowed_dist = fMinOnOffDist + radius + off_radius;

    DEBUG_OUT << "rot_angle " << rot_angle << std::endl;

    DEBUG_OUT << "rot_angle_lambda " << rot_angle_lambda << std::endl;

    DEBUG_OUT << "tmpOff " << *tmpOff << std::endl;

    DEBUG_OUT << "tmpOff rotation " << tmpOff->GetRotationAngle() << std::endl;

    DEBUG_OUT << "tmpOffpos_radec " << tmpOffpos_radec << std::endl;

    DEBUG_OUT << "onpos_radec " << onpos_radec << std::endl;

    DEBUG_OUT << "on_off_dist " << on_off_dist << std::endl;

    DEBUG_OUT << "off_radius " << off_radius << std::endl;

    DEBUG_OUT << "allowed_dist " << allowed_dist << std::endl;

    // Skip to next test position if off region too close to on region

    if (on_off_dist < allowed_dist) {

      DEBUG_OUT << "Skipping rot_angle " << rot_angle << " because on_off_dist " << on_off_dist

                << " is less than " << fMinOnOffDist << std::endl;

      rot_angle += stepsize;

      continue;

    }

    // Skip to next test position if off region overlaps with exclusion region

    if (tmpOff->Overlaps(GetBgStats()->GetOnExclusionRegion())) {

      DEBUG_OUT << "Skipping rot_angle " << rot_angle

                << " because it overlaps with the exclusion region." << std::endl;

      rot_angle += stepsize;

      continue;

    }

    // Now we know this is a valid off region position.

    // We create it and then move on by one segment length

    // to avoid placing new off regions over existing ones.

    DEBUG_OUT << "Accepting rot_angle " << rot_angle << std::endl;

    // add off region to bgstats

    std::ostringstream offName;

    offName << std::fixed << std::setprecision(1) << rot_angle << "deg";

    Tools::Region* newOff = (Tools::Region*)tmpOff->Clone( offName.str().c_str() );

    GetBgStats()->GetOffRegion().AddRegion( newOff );

    INFO_OUT << "Added new off region at "<< rot_angle << " deg from North at "

             << newOff->GetTrueCentre().GetCoordinate(*fHess->GetRADecJ2000System()) << std::endl;

    DEBUG_OUT << "offset from observation position: "

              << newOff->GetTrueCentre().GetCoordinate(*fHess->GetRADecJ2000System())

      .GetAngularDistance( obspos_radec ).GetDegrees() << " deg" << std::endl;

    if( fNumberOffRegions == fMaxNumberOffRegions ) {

      INFO_OUT << "Reached maximum allowed number of OFF regions ("

               << fMaxNumberOffRegions << "). Stopping creating new "

               << "ones at your request." << std::endl;

      break;

    }

    fNumberOffRegions++;

    rot_angle += segment;

  } while( rot_angle < 360 );

  delete tmpOff;

  if (fNumberOffRegions == 0) {

    WARN_OUT << "No Off regions found." << std::endl;

    return false;

  }

  INFO_OUT << "Added "<< fNumberOffRegions << " OFF regions at offset " << fCurrentOffset << endl;

  return true;

}