JPMTParametersMap.hh

Go to the documentation of this file.

#ifndef __JDETECTOR__JPMTPARAMETERSMAP__
#define __JDETECTOR__JPMTPARAMETERSMAP__
 
#include <istream>
#include <ostream>
#include <sstream>
#include <fstream>
#include <map>
 
#include "JDetector/JPMTIdentifier.hh"
#include "JDetector/JPMTParameters.hh"
#include "JDetector/JPMTParametersToolkit.hh"
#include "JLang/JEquationParameters.hh"
#include "JLang/JException.hh"
#include "JLang/JThrow.hh"
#include "JLang/JStringStream.hh"
#include "JLang/JObjectStreamIO.hh"
#include "Jeep/JComment.hh"
#include "Jeep/JProperties.hh"
#include "JSystem/JStat.hh"
 
 
/**
 * \author mdejong
 * \file
 * Map of PMT parameters.
 */
 
namespace JDETECTOR {}
namespace JPP { using namespace JDETECTOR; }
 
namespace JDETECTOR {
 
  using JLANG::JEquationParameters;
  using JLANG::JObjectStreamIO;
  using JLANG::JFileOpenException;
  using JLANG::JFileReadException;
  using JLANG::JValueOutOfRange;
  using JLANG::JException;
  using JLANG::JThrow;
  using JEEP::JComment;
 
 
  /**
   * Type definition of map PMT identifier to PMT parameters.
   */
  typedef std::map<JPMTIdentifier, JPMTParameters>  JPMTParametersMap_t;
 
 
  /**
   * Auxiliary class for map of PMT parameters.
   *
   * Input syntax as follows:
   * <pre>
   *            QE=<value>
   *            mu=<value>
   *
   *            %.<key>=<value>
   *
   *            PMT=<module> <pmt> <PMT parameters>
   *
   *            pmt=<module> <pmt> <key>=<value>[, <key>=<value>]
   *
   *            file=<file name>
   *
   * </pre>
   * where 
   *  - QE   refers to the global quantum efficiency;
   *  - mu   refers to the expectation value for number of photo-electrons given two-fold (or higher) coincidence rate;
   *  - '\%'        to the default values;
   *  - <module>  to the module identifier;
   *  - <pmt>     to the readout channel;
   *  - <key>     to one of the data members of the JDETECTOR::JPMTParameters data structure; and 
   *  - <value>   to the corresponding value.
   *
   * The global QE is applied as a multiplicative factor to the default and PMT specific QEs.
   *
   * The default values are provided by the method getDefaultPMTParameters and
   * can also be modified using method setDefaultPMTParameters.
   * If no optional values are given, the default values are maintained.
   *
   * Finally, <PMT parameters> corresponds to the list of PMT parameter values that should
   * be compatible with the I/O methods of the JDETECTOR::JPMTParameters data structure.
   *
   * Note that multiple input sequences should be seperated by a semicolumn ';' or be terminated by a newline.
   * When multiple input sequences appear for the same target, the last sequence prevails.
   */
  class JPMTParametersMap :
    public JPMTParametersMap_t,
    public JObjectStreamIO<JPMTParametersMap>,
    public JThrow<JPMTParametersMap>
  {
  public:
    /**
     * Constructor.
     *
     * \param  parameters              PMT parameters
     */
    JPMTParametersMap(const JPMTParameters& parameters = JPMTParameters()) :
      JPMTParametersMap_t(),
      QE(1.0),
      mu(0.0)
    {
      setDefaultPMTParameters(parameters);
    }
 
 
    /**
     * Constructor.
     *
     * \param  file_name               file name
     */
    JPMTParametersMap(const char* const file_name) : 
      JPMTParametersMap_t(),
      QE(1.0),
      mu(0.0)
    {
      using namespace std;
 
      ifstream in(file_name);
 
      in >> *this;
 
      in.close();
    }
 
 
    /**
     * Get PMT parameters.\n
     * This method returns the default PMT parameters if the parameters 
     * corresponding to the given PMT identifier have not been defined. 
     *
     * Note that the value of QE as part of the return value has been scaled
     * with the global QE of the PMT parameters map.
     *
     * \param  id                      PMT identifier
     * \return                         PMT parameters
     */
    const JPMTParameters& getPMTParameters(const JPMTIdentifier& id) const
    {
      static JPMTParameters parameters;
      
      JPMTParametersMap_t::const_iterator i = find(id);
 
      if (i != end())
        parameters = i->second;
      else
        parameters = getDefaultPMTParameters();
 
      parameters.QE *= this->QE;
      
      return parameters;
    }
 
 
    /**
     * Get QE of given PMT.
     *
     * \param  id                      PMT identifier
     * \return                         QE
     */
    double getQE(const JPMTIdentifier& id) const
    {
      return getPMTParameters(id).QE;
    }
 
 
    /**
     * Get ratio of hit probabilities of given PMT.
     *
     * \param  id                      PMT identifier
     * \return                         ratio
     */
    double getHitProbability(const JPMTIdentifier& id) const
    {
      using JDETECTOR::getHitProbability;
 
      return getHitProbability(getQE(id), this->mu);
    }
 
 
    /**
     * Get default PMT parameters.
     *
     * \return                         PMT parameters
     */
    const JPMTParameters& getDefaultPMTParameters() const
    {
      return defaultPMTParameters;
    }
 
 
    /**
     * Set default PMT parameters.
     *
     * \param  parameters              PMT parameters
     */
    void setDefaultPMTParameters(const JPMTParameters& parameters)
    {
      defaultPMTParameters = parameters;
    }
 
 
    /**
     * Get global QE.
     *
     * \return                         QE
     */
    double getQE() const
    {
      return this->QE;
    }
 
 
    /**
     * Set global QE.
     *
     * \param  QE                      QE
     */
    void setQE(const double QE)
    {
      this->QE = QE;
    }
 
 
    /**
     * Get expectation value for number of photo-electrons given two-fold (or higher) coincidence rate.
     *
     * \return                         expectation value
     */
    double getMu() const
    {
      return this->mu;
    }
 
 
    /**
     * Check validity of PMT parameters.
     *
     * \return                         true if valid; else false
     */
    bool is_valid() const
    {
      try {
 
        if (this->QE < 0) {
          THROW(JException, "Invalid global QE " << QE);
        }
 
        if (this->mu < 0) {
          THROW(JException, "Invalid expectation value number of photo-electrons " << mu);
        }
        
        if (!defaultPMTParameters.is_valid()) {
          THROW(JException, "Invalid default PMT parameters");
        }
        
        for (const_iterator i = this->begin(); i != this->end(); ++i) {
          
          if (!i->second.is_valid()) {
            THROW(JException, "Invalid PMT parameters at " << i->first);
          }
 
          if (i->second.QE * this->QE > 1.0 ||
              i->second.QE * this->QE < 0.0) {
            THROW(JException, "Invalid PMT QE at " << i->first << ' ' << i->second.QE * this->QE);
          }
        }
 
        return true;
      }
      catch(const JException& error) {
 
        Throw(error);
 
        return false;
      }
    }
 
 
    /**
     * Convert the hit probabilities to QEs for given expectation value.
     *
     * The expectation value corresponds to the number of photo-electrons given 
     * the multiplicity range of coincidences as specified at JCalibrateK40.cc.\n
     * The survival probability is taken into account in the conversion 
     * (see method JDETECTOR::getSurvivalProbability).
     *
     * \param  mu                      expectation value
     */
    void convertHitProbabilityToQE(const double mu)
    {
      using namespace std;
      using JDETECTOR::getHitProbability;
      using JDETECTOR::getQE;
 
      this->mu = mu;
 
      const double Pmax = getMaximalHitProbability(this->mu);
 
      for (JPMTParametersMap::const_iterator i = this->begin(); i != this->end(); ++i) {
        if (i->second.QE > Pmax) {
          THROW(JValueOutOfRange, "Hit probability PMT " << i->first << ' ' << i->second.QE << " > maximum probability given expectation value " << this->mu << endl);
        }
      }
      
      for (JPMTParametersMap::iterator i = this->begin(); i != this->end(); ++i) {
        
        const double P = getSurvivalProbability(i->second);
        
        if (i->second.QE > 0.0 && P > 0.0)
          i->second.QE = getQE(i->second.QE, this->mu) / P; 
        else
          i->second.QE = 0.0;
      }
    }
 
 
    /**
     * Stream input
     *
     * \param  in                      input stream
     * \param  object                  PMT parameters map
     * \return                         input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JPMTParametersMap& object)
    {
      using namespace std;
      using namespace JPP;
      
      JStringStream is(in);
 
      if (getFileStatus(is.str().c_str())) {
        is.load();
      }
 
      JProperties properties(getEquationParameters(), 1);
        
      JPMTParametersMapHelper  helper(object, getEquationParameters());
      JProperties              demo = object.defaultPMTParameters.getProperties();
      
      properties["PMT"] = static_cast<JPMTParametersMap_t&>(object);
      properties["pmt"] = helper;
      properties["%"]   = demo;
      properties["QE"]  = object.QE;
      properties["mu"]  = object.mu;
 
      is >> object.comment;
      is >> properties;
      
      return in;
    }
 
 
    /**
     * Stream output
     *
     * \param  out                     output stream
     * \param  object                  PMT parameters map
     * \return                         output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JPMTParametersMap& object)
    {
      using namespace JLANG;
      
      JProperties properties(getEquationParameters());
 
      JProperties demo = object.defaultPMTParameters.getProperties();
 
      properties["PMT"] = static_cast<const JPMTParametersMap_t&>(object);
      properties["%"]   = demo;
 
      out << object.comment;
      out << setequation(getEquationParameters()) << JEquation::make_equation("QE", object.QE);
      out << setequation(getEquationParameters()) << JEquation::make_equation("mu", object.mu);
      out << properties;
 
      return out;
    }
 
 
    /**
     * Get equation parameters.
     *
     * \return                         equation parameters
     */
    static inline JEquationParameters& getEquationParameters()
    {
      static JEquationParameters parameters;
 
      return parameters;
    }
 
 
    /**
     * Set equation parameters.
     *
     * \param  equation                equation parameters
     */
    static inline void setEquationParameters(const JEquationParameters& equation)
    {
      getEquationParameters() = equation;
    }
 
 
    JComment comment;
 
  protected:
    /**
     * Default PMT parameters.
     */
    JPMTParameters defaultPMTParameters;
 
    /**
     * Global QE.
     */
    double QE;
 
 
    /**
     * Expectation value for number of photo-electrons given two-fold (or higher) coincidence rate.
     */
    double mu;
 
  private:
    /**
     * Auxiliary class for I/O of PMT parameters map.
     */
    class JPMTParametersMapHelper {
    public:
      /**
       * Constructor.
       *
       * \param  __object                PMT parameters map
       * \param  __equation              equation parameters     
       */
      JPMTParametersMapHelper(JPMTParametersMap&         __object,
                              const JEquationParameters& __equation) :
        object  (__object),
        equation(__equation)
      {}
 
 
      /**
       * Stream input PMT parameters map helper.
       *
       * \param  in                      input stream
       * \param  helper                  PMT parameters map helper
       * \return                         input stream
       */
      friend inline std::istream& operator>>(std::istream& in, JPMTParametersMapHelper& helper)
      {
        JPMTIdentifier id;
 
        if (in >> id) {
          helper.object[id].getProperties().read(in);
        }
 
        return in;
      }
 
 
      /**
       * Stream output PMT parameters map helper.
       *
       * \param  out                     output stream
       * \param  helper                  PMT parameters map helper
       * \return                         output stream
       */
      friend inline std::ostream& operator<<(std::ostream& out, const JPMTParametersMapHelper& helper)
      {
        for (JPMTParametersMap::const_iterator i = helper.object.begin(); i != helper.object.end(); ++i) {
 
          out << i->first;
 
          i->second.getProperties(helper.equation).write(out);
 
          out << std::endl;
        }
 
        return out;
      }
      
    private:
      JPMTParametersMap&         object;
      const JEquationParameters& equation;
    };
  };
}
 
#endif