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/JClonable.hh"
#include "JLang/JTitle.hh"
#include "JLang/JException.hh"
 
 
/**
 * \file
 *
 * Nuisance parameter.
 * \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;
 
 
  /**
   * Interface for nuisance parameter.
   */
  struct JNuisance :
    public JClonable<JNuisance>,
    public JTitle
  {
    /**
     * 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 JClonable<JNuisance, JNuisanceFixed>
  {
    /**
     * 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 JClonable<JNuisance, JNuisanceGauss>
  {
    /**
     * 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 JClonable<JNuisance, JNuisanceTF1>
  {
    /**
     * 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 JClonable<JNuisance, JNuisanceTH1>
  {
    /**
     * 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 data structure for nuisance I/O.
   */
  struct JNuisanceHelper :
    public std::map<std::string, nuisance_type>
  {
    /**
     * Default constructor.
     */
    JNuisanceHelper()
    {
      (*this)["fixed"]    .reset(new JNuisanceFixed());
      (*this)["Gauss"]    .reset(new JNuisanceGauss());
      (*this)["Formula"]  .reset(new JNuisanceTF1());
      (*this)["Histogram"].reset(new JNuisanceTH1());
    }
 
 
    /**
     * Parse nuisance from input stream.
     *
     * \param  in                 input stream
     * \param  object             nuisance
     * \return                    input stream
     */
    std::istream& operator()(std::istream& in, nuisance_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 nuisance I/O.
   */
  static const JNuisanceHelper nuisance_helper;
 
 
  /**
   * Read nuisance from input stream.
   *
   * \param  in                 input stream
   * \param  object             nuisance
   * \return                    input stream
   */
  inline std::istream& operator>>(std::istream& in, nuisance_type& object)
  {
    return nuisance_helper(in, object);
  }
 
 
  /**
   * Write nuisance to output stream.
   *
   * \param  out                output stream
   * \param  object             nuisance
   * \return                    output stream
   */
  inline std::ostream& operator<<(std::ostream& out, const nuisance_type& object)
  {
    if (object) {
      out << object->getTitle() << ' ' << *object.get();
    }
 
    return out;
  }
}
 
#endif