JMorphology.hh

Go to the documentation of this file.

#ifndef __JASTRONOMY__JMORPHOLOGY__
#define __JASTRONOMY__JMORPHOLOGY__
 
#include <string>
#include <istream>
#include <ostream>
#include <memory>
#include <map>
 
#include "TRandom3.h"
#include "TFile.h"
#include "TH2D.h"
 
#include "JGeometry3D/JAngle3D.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JGeometry3D/JRotation3D.hh"
 
#include "JLang/JTitle.hh"
#include "JLang/JException.hh"
 
#include "JAstronomy/JAstronomy.hh"
 
 
/**
 * \file
 * 
 * Tools for simulation of source morphology.
 * \author mdejong
 */
namespace JASTRONOMY {}
namespace JPP { using namespace JASTRONOMY; }
 
namespace JASTRONOMY {
 
  using JLANG::JTitle;
  using JLANG::JParseError;
  using JLANG::JFileOpenException;
  using JLANG::JValueOutOfRange;
 
 
  /**
   * Interface for source morphology simulation.
   * Note that all input angles are in degrees.
   */
  struct JMorphology :
    public JSourceLocation
  {
    /**
     * Virtual destructor.
     */
    virtual ~JMorphology()
    {}
 
 
    /**
     * Read source location in degrees from input stream.
     *
     * \param  in                 input stream
     * \param  source             source
     * \return                    input stream
     */
    static std::istream& read(std::istream& in, JSourceLocation& source)
    {
      angle_type_deg angle;
 
      in >> angle;
 
      source.set(angle);
 
      return in;
    }
 
 
    /**
     * Write source location in degrees to output stream.
     *
     * \param  out                output stream
     * \param  source             source
     * \return                    output stream
     */
    static std::ostream& write(std::ostream& out, const JSourceLocation& source)
    {
      return out << angle_type_deg(source);
    }
 
 
    /**
     * Read source morphology from input stream.
     *
     * \param  in                 input stream
     * \return                    input stream
     */
    virtual std::istream& read(std::istream& in) = 0;
    
 
    /**
     * Write source morphology to output stream.
     *
     * \param  out                output stream
     * \return                    output stream
     */
    virtual std::ostream& write(std::ostream& out) const = 0;
 
 
    /**
     * Get location on sky.
     *
     * \return                    location on sky
     */
    virtual JSourceLocation get() const = 0;
  };
 
 
  /**
   * Implementation of point source morphology.
   * Note that all input angles are in degrees.
   */
  struct JMorphologyPoint :
    public JMorphology
  {
    /**
     * Read source morphology from input stream.
     *
     * \param  in                 input stream
     * \return                    input stream
     */
    virtual std::istream& read(std::istream& in) override
    {
      return JMorphology::read(in, *this);
    }
    
 
    /**
     * Write source morphology to output stream.
     *
     * \param  out                output stream
     * \return                    output stream
     */
    virtual std::ostream& write(std::ostream& out) const override
    {
      return JMorphology::write(out, *this);
    }
    
 
    /**
     * Get location on sky.
     *
     * \return                    location on sky
     */
    virtual JSourceLocation get() const override
    {
      using namespace JPP;
 
      JDirection3D u(0.0, 0.0, 1.0);
   
      const JRotation3D Rs(this->getSourceLocation());
 
      u.rotate_back(Rs);
 
      return JSourceLocation(u);
    }
  };
 
 
  /**
   * Implementation of Gaussian source morphology.
   * Note that all input angles are in degrees.
   */
  struct JMorphologyGauss :
    public JMorphology
  {
    /**
     * Read source morphology from input stream.
     *
     * \param  in                 input stream
     * \return                    input stream
     */
    virtual std::istream& read(std::istream& in) override
    {
      return JMorphology::read(in, *this)
        >> this->sigma_deg;
    }
    
 
    /**
     * Write source morphology to output stream.
     *
     * \param  out                output stream
     * \return                    output stream
     */
    virtual std::ostream& write(std::ostream& out) const override
    {
      return JMorphology::write(out, *this)
        << ' ' << this->sigma_deg;
    }
    
 
    /**
     * Get location on sky.
     *
     * \return                    location on sky
     */
    virtual JSourceLocation get() const override
    {
      using namespace JPP;
 
      JDirection3D u(0.0, 0.0, 1.0);
   
      if (sigma_deg > 0.0) {
 
        const double theta = gRandom->Gaus(0.0, getRadians(sigma_deg));
        const double phi   = gRandom->Uniform(-PI, +PI);
 
        const JRotation3D R(JAngle3D(theta, phi));
 
        u.rotate_back(R);
      }
 
      const JRotation3D Rs(this->getSourceLocation());
 
      u.rotate_back(Rs);
 
      return JSourceLocation(u);
    }
 
 
    double sigma_deg;
  };
 
 
  /**
   * Implementation of 2D-Gaussian source morphology.
   * Note that all input angles are in degrees.
   */
  struct JMorphologyGauss2D :
    public JMorphology
  {
    /**
     * Read source morphology from input stream.
     *
     * \param  in                 input stream
     * \return                    input stream
     */
    virtual std::istream& read(std::istream& in) override
    {
      return JMorphology::read(in, *this)
        >> this->sigmaX_deg >> this->sigmaY_deg;
    }
    
 
    /**
     * Write source morphology to output stream.
     *
     * \param  out                output stream
     * \return                    output stream
     */
    virtual std::ostream& write(std::ostream& out) const override
    {
      return JMorphology::write(out, *this)
        << ' ' << this->sigmaX_deg << ' ' << this->sigmaY_deg;
    }
    
 
    /**
     * Get location on sky.
     *
     * \return                    location on sky
     */
    virtual JSourceLocation get() const override
    {
      using namespace JPP;
 
      JDirection3D u(0.0, 0.0, 1.0);
   
      if (sigmaX_deg > 0.0 ||
          sigmaY_deg > 0.0) {
 
        const double x = sin(gRandom->Gaus(0.0, getRadians(sigmaX_deg)));
        const double y = sin(gRandom->Gaus(0.0, getRadians(sigmaY_deg)));
 
        u = JDirection3D(x, y, sqrt(1.0 - x*x - y*y));
      }
 
      const JRotation3D Rs(this->getSourceLocation());
 
      u.rotate_back(Rs);
 
      return JSourceLocation(u);
    }
 
 
    double sigmaX_deg;
    double sigmaY_deg;
  };
 
 
  /**
   * Implementation of binary source morphology.
   * Note that all input angles are in degrees.
   */
  struct JMorphologyBinary :
    public JMorphology
  {
    /**
     * Read source morphology from input stream.
     *
     * \param  in                 input stream
     * \return                    input stream
     */
    virtual std::istream& read(std::istream& in) override
    {
      return JMorphology::read(in, *this)
        >> this->x1_deg >> this->y1_deg >> this->w1 >> this->s1_deg
        >> this->x2_deg >> this->y2_deg >> this->w2 >> this->s2_deg;
    }
    
 
    /**
     * Write source morphology to output stream.
     *
     * \param  out                output stream
     * \return                    output stream
     */
    virtual std::ostream& write(std::ostream& out) const override
    {
      return JMorphology::write(out, *this)
        << ' ' << this->x1_deg << ' ' << this->y1_deg << ' ' << this->w1 << ' ' << this->s1_deg
        << ' ' << this->x2_deg << ' ' << this->y2_deg << ' ' << this->w2 << ' ' << this->s2_deg;
    }
    
 
    /**
     * Get location on sky.
     *
     * \return                    location on sky
     */
    virtual JSourceLocation get() const override
    {
      using namespace JPP;
 
      const double w = gRandom->Uniform(0.0, w1 + w2);
 
      const double x = sin(getRadians(w < w1 ? x1_deg : x2_deg));
      const double y = sin(getRadians(w < w1 ? y1_deg : y2_deg));
      const double s =     getRadians(w < w1 ? s1_deg : s2_deg);
 
      JDirection3D v(x, y, sqrt(1.0 - x*x - y*y));
 
      JDirection3D u(0.0, 0.0, 1.0);
 
      if (s > 0.0) {
 
        const double theta = gRandom->Gaus(0.0, s);
        const double phi   = gRandom->Uniform(-PI, +PI);
 
        const JRotation3D R(JAngle3D(theta, phi));
 
        u.rotate_back(R);
      }
      {
        const JRotation3D R(v);
 
        u.rotate_back(R);
      }
 
      const JRotation3D Rs(this->getSourceLocation());
 
      u.rotate_back(Rs);
 
      return JSourceLocation(u);
    }
 
 
    double x1_deg, y1_deg, w1, s1_deg;
    double x2_deg, y2_deg, w2, s2_deg;
  };
 
 
  /**
   * Implementation of histogram source morphology.
   * Note that all input angles are in degrees.
   */
  struct JMorphologyHistogram :
    public JMorphology
  {
    /**
     * Virtual destructor.
     */
    virtual ~JMorphologyHistogram()
    {
      if (in != NULL) {
        in->Close();
      }
    }
 
 
    /**
     * Read source morphology from input stream.
     *
     * \param  in                 input stream
     * \return                    input stream
     */
    virtual std::istream& read(std::istream& in) override
    {
      JMorphology::read(in, *this)
        >> this->filename
        >> this->histname;
 
      load();
 
      return in;
    }
    
 
    /**
     * Write source morphology to output stream.
     *
     * \param  out                output stream
     * \return                    output stream
     */
    virtual std::ostream& write(std::ostream& out) const override
    {
      return JMorphology::write(out, *this)
        << ' ' << this->filename
        << ' ' << this->histname;
    }
    
 
    /**
     * Get location on sky.
     *
     * \return                    location on sky
     */
    virtual JSourceLocation get() const override
    {
      using namespace JPP;
 
      Double_t x;
      Double_t y;
 
      h2->GetRandom2(x, y, gRandom);
 
      x = sin(getRadians(x));
      y = sin(getRadians(y));
        
      JDirection3D u(x, y, sqrt(1.0 - x*x - y*y));
 
      const JRotation3D Rs(this->getSourceLocation());
 
      u.rotate_back(Rs);
 
      return JSourceLocation(u);
    }
 
    std::string filename;
    std::string histname;
 
  protected:
    /**
     * Load histogram from file.
     */
    void load()
    {
      in = TFile::Open(filename.c_str(), "exist");
 
      if (in == NULL || !in->IsOpen()) {
        THROW(JFileOpenException, "File: " << filename << " not opened.");
      }
 
      h2 = dynamic_cast<TH2D*>(in->Get(histname.c_str()));
 
      if (h2 == NULL) {
        THROW(JValueOutOfRange, "Histogram: " << histname << " not found.");
      }
    }
 
    mutable TFile* in = NULL;
    mutable TH2D*  h2 = NULL;
  };
 
 
  /**
   * Type definition of generic morphology.
   */
  typedef  std::shared_ptr<JMorphology>    morphology_type;
 
 
  /**
   * Helper for source morphology I/O.
   */
  struct JMorphology_t :
    public JTitle,
    public morphology_type
  {
    typedef  std::map<std::string, morphology_type> dictionary_type;
 
    /**
     * Default constructor.
     */
    JMorphology_t()
    {
      dictionary["Point"]    .reset(new JMorphologyPoint());
      dictionary["Gauss"]    .reset(new JMorphologyGauss());
      dictionary["Gauss2D"]  .reset(new JMorphologyGauss2D());
      dictionary["Binary"]   .reset(new JMorphologyBinary());
      dictionary["Histogram"].reset(new JMorphologyHistogram());
    }
 
 
    /**
     * Get dictionary.
     *
     * \return                    dictionary
     */
    const dictionary_type& getDictionary() const
    {
      return dictionary;
    }
 
 
    /**
     * Put source morphology at given key.
     *
     * \param  key                key
     * \param  morphology         morphology
     */
    void put(const std::string& key, JMorphology* morphology)
    {
      return dictionary[key].reset(morphology);
    }
 
 
    /**
     * Read source morphology from input stream.
     *
     * \param  in                 input stream
     * \param  object             source morphology
     * \return                    input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JMorphology_t& object)
    {
      using namespace std;
 
      string key;
 
      if (in >> key) {
 
        if (object.dictionary.count(key)) {
 
          static_cast<shared_ptr<JMorphology>&>(object) = object.dictionary[key];
 
          object.setTitle(key);
 
          object->read(in);
          
        } else {
 
          THROW(JParseError, "Invalid key " << key);
        }
      }
      
      return in;
    }
 
 
    /**
     * Write source morphology to output stream.
     *
     * \param  out                output stream
     * \param  object             source morphology
     * \return                    output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JMorphology_t& object)
    {
      out << object.getTitle() << ' ';
 
      if (static_cast<const std::shared_ptr<JMorphology>&>(object)) {   
        object->write(out);
      }
 
      return out;
    }
 
  protected:
    dictionary_type  dictionary;
  };
}
 
#endif