JResolution.hh

Go to the documentation of this file.

#ifndef __JASTRONOMY__JRESOLUTION__
#define __JASTRONOMY__JRESOLUTION__
 
#include <string>
#include <istream>
#include <ostream>
#include <iomanip>
#include <memory>
#include <map>
 
#include "TRandom3.h"
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TF1.h"
 
#include "JGeometry3D/JAngle3D.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JGeometry3D/JRotation3D.hh"
 
#include "JLang/JClonable.hh"
#include "JLang/JTitle.hh"
#include "JLang/JException.hh"
 
#include "JAstronomy/JAstronomy.hh"
 
 
/**
 * \file
 *
 * Tools for simulation of detector resolution.
 * \author mdejong
 */
namespace JASTRONOMY {}
namespace JPP { using namespace JASTRONOMY; }
 
namespace JASTRONOMY {
 
  using JLANG::JClonable;
  using JLANG::JTitle;
  using JLANG::JParseError;
  using JLANG::JFileOpenException;
  using JLANG::JValueOutOfRange;
  using JGEOMETRY3D::JDirection3D;
 
 
  /**
   * Interface for detector resolution simulation.
   * Note that all input angles are in degrees.
   */
  struct JResolution :
    public JClonable<JResolution>,
    public JTitle
  {
    /**
     * Virtual destructor.
     */
    virtual ~JResolution()
    {}
 
 
    /**
     * Read detector resolution from input stream.
     *
     * \param  in                 input stream
     * \return                    input stream
     */
    virtual std::istream& read(std::istream& in) = 0;
 
 
    /**
     * Write detector resolution to output stream.
     *
     * \param  out                output stream
     * \return                    output stream
     */
    virtual std::ostream& write(std::ostream& out) const = 0;
 
 
    /**
     * Get direction.
     *
     * \return                    direction
     */
    virtual JDirection3D get() const = 0;
  };
  
 
  /**
   * Implementation of Gaussian detector resolution.
   * Note that all input angles are in degrees.
   */
  struct JResolutionGauss :
    public JClonable<JResolution, JResolutionGauss>
  {
    /**
     * Read detector resolution from input stream.
     *
     * \param  in                 input stream
     * \return                    input stream
     */
    virtual std::istream& read(std::istream& in) override
    {
      return in >> this->sigma_deg;
    }
 
 
    /**
     * Write detector resolution to output stream.
     *
     * \param  out                output stream
     * \return                    output stream
     */
    virtual std::ostream& write(std::ostream& out) const override
    {
      return out << this->sigma_deg;
    }
 
 
    /**
     * Get direction.
     *
     * \return                    direction
     */
    JDirection3D 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);
      }
 
      return u;
    }
 
 
    double sigma_deg;    
  };
  
 
  /**
   * Implementation of detector resolution based on a formula.
   * Note that all input angles are in degrees.
   */
  struct JResolutionTF1 :
    public JClonable<JResolution, JResolutionTF1>
  {
    /**
     * Read detector resolution from input stream.
     *
     * \param  in                 input stream
     * \return                    input stream
     */
    virtual std::istream& read(std::istream& in) override
    {
      TString buffer;
 
      if (buffer.ReadLine(in)) {
      
        f1.reset(new TF1("f1", buffer));
 
        if (!f1->IsValid()) {
          THROW(JParseError, "Invalid function " << buffer);
        }
      }
      
      return in;
    }
 
 
    /**
     * Write detector resolution to output stream.
     *
     * \param  out                output stream
     * \return                    output stream
     */
    virtual std::ostream& write(std::ostream& out) const override
    {
      if (f1)
        return out << f1->GetExpFormula();
      else
        return out;
    }
 
 
    /**
     * Get direction.
     *
     * \return                    direction
     */
    JDirection3D get() const override
    {
      using namespace JPP;
 
      JDirection3D u(0.0, 0.0, 1.0);
 
      if (f1) {
 
        const double theta = getRadians(f1->GetRandom(gRandom));
        const double phi   = gRandom->Uniform(-PI, +PI);
 
        const JRotation3D R(JAngle3D(theta, phi));
 
        u.rotate_back(R);
      }
 
      return u;
    }
 
 
    std::shared_ptr<TF1> f1;
  };
 
 
  /**
   * Implementation of detector resolution based on histogram.
   * Note that all input angles are in degrees.
   */
  struct JResolutionTH1 :
    public JClonable<JResolution, JResolutionTH1>
  {
    /**
     * Virtual destructor.
     */
    virtual ~JResolutionTH1()
    {
      if (in != NULL) {
        in->Close();
      }
    }
 
 
    /**
     * Read detector resolution from input stream.
     *
     * \param  in                 input stream
     * \return                    input stream
     */
    virtual std::istream& read(std::istream& in) override
    {
      in >> this->filename
         >> this->histname;
 
      load();
 
      return in;
    }
 
 
    /**
     * Write detector resolution to output stream.
     *
     * \param  out                output stream
     * \return                    output stream
     */
    virtual std::ostream& write(std::ostream& out) const override
    {
      return out << this->filename << ' '
                 << this->histname;
    }
 
 
    /**
     * Get direction.
     *
     * \return                    direction
     */
    JDirection3D get() const override
    {
      using namespace JPP;
 
      JDirection3D u(0.0, 0.0, 1.0);
 
      if (h1 != NULL) {
 
        const double theta = getRadians(h1->GetRandom(gRandom));
        const double phi   = gRandom->Uniform(-PI, +PI);
 
        const JRotation3D R(JAngle3D(theta, phi));
 
        u.rotate_back(R);
      }
 
      return 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.");
      }
 
      h1 = dynamic_cast<TH1D*>(in->Get(histname.c_str()));
 
      if (h1 == NULL) {
        THROW(JValueOutOfRange, "Histogram: " << histname << " not found.");
      }
    }
 
    mutable TFile* in = NULL;
    mutable TH1D*  h1 = NULL;
  };
 
 
  /**
   * Implementation of detector resolution based on histogram filled according <tt>log10(alpha)</tt>.
   * Note that all input angles are in degrees.
   */
  struct JResolutionLog :
    public JResolutionTH1
  {
    /**
     * Get direction.
     *
     * \return                    direction
     */
    JDirection3D get() const override
    {
      using namespace JPP;
 
      JDirection3D u(0.0, 0.0, 1.0);
 
      if (h1 != NULL) {
 
        const double theta = getRadians(pow(10.0,h1->GetRandom(gRandom)));
        const double phi   = gRandom->Uniform(-PI, +PI);
 
        const JRotation3D R(JAngle3D(theta, phi));
 
        u.rotate_back(R);
      }
 
      return u;
    }
  };
 
 
  /**
   * Type definition of generic resolution.
   */
  typedef  std::shared_ptr<JResolution>    resolution_type;
 
 
  /**
   * Helper data structure for detector resolution I/O.
   */
  struct JResolutionHelper :
    public std::map<std::string, resolution_type>
  {
    /**
     * Default constructor.
     */
    JResolutionHelper()
    {
      (*this)["Gauss"]      .reset(new JResolutionGauss());
      (*this)["Formula"]    .reset(new JResolutionTF1());
      (*this)["Histogram"]  .reset(new JResolutionTH1());
      (*this)["Logarithmic"].reset(new JResolutionLog());
    }
 
 
    /**
     * Parse detector resolution from input stream.
     *
     * \param  in                 input stream
     * \param  object             detector resolution
     * \return                    input stream
     */
    inline std::istream& operator()(std::istream& in, resolution_type& object) const
    {
      std::string key;
 
      if (in >> key) {
 
        const_iterator p = this->find(key);
 
        if (p != this->end()) {
 
          object.reset(p->second->clone());
 
          object->setTitle(key);
          object->read(in);
 
        } else {
 
          THROW(JParseError, "Invalid key " << key);
        }
      }
 
      return in;
    }
  };
 
 
  /**
   * Helper object for resolution I/O.
   */
  static const JResolutionHelper resolution_helper;
 
 
  /**
   * Read resolution from input stream.
   *
   * \param  in                 input stream
   * \param  object             resolution
   * \return                    input stream
   */
  inline std::istream& operator>>(std::istream& in, resolution_type& object)
  {
    return resolution_helper(in, object);
  }
 
 
  /**
   * Write resolution to output stream.
   *
   * \param  out                output stream
   * \param  object             resolution
   * \return                    output stream
   */
  inline std::ostream& operator<<(std::ostream& out, const resolution_type& object)
  {
    if (object) {
      out << object->getTitle() << ' ';
      object->write(out);
    }
 
    return out;
  }
}
 
#endif