JNuisance.hh

Go to the documentation of this file.

#ifndef __JASTRONOMY__JNUISANCE__
#define __JASTRONOMY__JNUISANCE__
 
#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 "TString.h"
 
#include "JLang/JException.hh"
 
 
/**
 * \file
 *
 * Nuisance parameter.
 * \author mdejong
 */
namespace JASTRONOMY {}
namespace JPP { using namespace JASTRONOMY; }
 
namespace JASTRONOMY {
 
  using JLANG::JParseError;
  using JLANG::JFileOpenException;
  using JLANG::JValueOutOfRange;
 
 
  /**
   * Interface for nuisance parameter.
   */
  struct JNuisance
  {
    /**
     * Central value.
     */
    static constexpr double  FACTOR  =  1.0;
 
    /**
     * Virtual destructor.
     */
    virtual ~JNuisance()
    {}
 
 
    /**
     * Read nuisance from input stream.
     *
     * \param  in                 input stream
     * \return                    input stream
     */
    virtual std::istream& read(std::istream& in) = 0;
 
 
    /**
     * Write nuisance to output stream.
     *
     * \param  out                output stream
     * \return                    output stream
     */
    virtual std::ostream& write(std::ostream& out) const = 0;
 
 
    /**
     * Get value.
     *
     * \return                    value
     */
    virtual double get() const = 0;
 
 
    /**
     * Read nuisance from input stream.
     *
     * \param  in                 input stream
     * \param  nuisance           nuisance
     * \return                    input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JNuisance& nuisance)
    {
      return nuisance.read(in);
    }
 
 
    /**
     * Write nuisance to output stream.
     *
     * \param  out                output stream
     * \param  nuisance           nuisance
     * \return                    output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JNuisance& nuisance)
    {
      return nuisance.write(out);
    }
  };
  
 
  /**
   * Implementation of fixed nuisance.
   */
  struct JNuisanceFixed :
    public JNuisance
  {
    /**
     * Read nuisance from input stream.
     *
     * \param  in                 input stream
     * \return                    input stream
     */
    virtual std::istream& read(std::istream& in) override
    {
      return in;
    }
 
 
    /**
     * Write nuisance to output stream.
     *
     * \param  out                output stream
     * \return                    output stream
     */
    virtual std::ostream& write(std::ostream& out) const override
    {
      return out;
    }
 
 
    /**
     * Get value.
     *
     * \return                    value
     */
    double get() const override
    {
      return FACTOR;
    }
  };
 
 
  /**
   * Implementation of Gaussian nuisance.
   */
  struct JNuisanceGauss :
    public JNuisance
  {
    /**
     * Read nuisance from input stream.
     *
     * \param  in                 input stream
     * \return                    input stream
     */
    virtual std::istream& read(std::istream& in) override
    {
      return in >> this->sigma;
    }
 
 
    /**
     * Write nuisance to output stream.
     *
     * \param  out                output stream
     * \return                    output stream
     */
    virtual std::ostream& write(std::ostream& out) const override
    {
      return out << this->sigma;
    }
 
 
    /**
     * Get value.
     *
     * \return                    value
     */
    double get() const override
    {
      return gRandom->Gaus(FACTOR, sigma);
    }
 
 
    double sigma;
  };
  
 
  /**
   * Implementation of nuisance based on a formula.
   */
  struct JNuisanceTF1 :
    public JNuisance
  {
    /**
     * Read nuisance 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 nuisance 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 value.
     *
     * \return                    value
     */
    double get() const override
    {
      return f1->GetRandom(gRandom);
    }
 
 
    std::shared_ptr<TF1> f1;
  };
 
 
  /**
   * Implementation of nuisance based on histogram.
   */
  struct JNuisanceTH1 :
    public JNuisance
  {
    /**
     * Virtual destructor.
     */
    virtual ~JNuisanceTH1()
    {
      if (in != NULL) {
        in->Close();
      }
    }
 
 
    /**
     * Read nuisance 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 nuisance 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 value.
     *
     * \return                    value
     */
    double get() const override
    {
      return h1->GetRandom(gRandom);
    }
 
 
    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;
  };
 
 
  /**
   * Type definition of generic nuisance.
   */
  typedef  std::shared_ptr<JNuisance>   nuisance_type;
 
 
  /**
   * Helper for nuisance I/O.
   */
  struct JNuisance_t :
    public nuisance_type
  {
    typedef  std::map<std::string, nuisance_type>  dictionary_type;
 
    /**
     * Default constructor.
     */
    JNuisance_t() :
      nuisance_type(std::make_shared<JNuisanceFixed>()),
      key("fixed")
    {
      dictionary[key] = static_cast<nuisance_type&>(*this);
 
      dictionary["Gauss"]    .reset(new JNuisanceGauss());
      dictionary["Formula"]  .reset(new JNuisanceTF1());
      dictionary["Histogram"].reset(new JNuisanceTH1());
    }
 
 
    /**
     * Get dictionary.
     *
     * \return                    dictionary
     */
    dictionary_type& getDictionary()
    {
      return dictionary;
    }
 
 
    /**
     * Read nuisance from input stream.
     *
     * \param  in                 input stream
     * \param  object             nuisance
     * \return                    input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JNuisance_t& object)
    {
      using namespace std;
 
      string key;
 
      if (in >> key) {
 
        if (object.dictionary.count(key)) {
 
          static_cast<nuisance_type&>(object) = object.dictionary[key];
 
          object.key = key;
 
          object->read(in);
 
        } else {
 
          THROW(JParseError, "Invalid key " << key);
        }
      }
 
      return in;
    }
 
 
    /**
     * Write nuisance to output stream.
     *
     * \param  out                output stream
     * \param  object             nuisance
     * \return                    output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JNuisance_t& object)
    {
      out << object.key << ' ';
 
      if (static_cast<const nuisance_type&>(object)) {
        object->write(out);
      }
 
      return out;
    }
 
  protected:
    std::string     key;
    dictionary_type dictionary;
  };
}
 
#endif