JHead.hh

Go to the documentation of this file.

#ifndef __JAANET__JHEAD__
#define __JAANET__JHEAD__

#include <string>
#include <vector>
#include <cmath>
#include <ctype.h>

#include "TObject.h"

#include "km3net-dataformat/offline/Vec.hh"
#include "km3net-dataformat/offline/Head.hh"

#include "JLang/JException.hh"
#include "JLang/JEquationParameters.hh"
#include "JROOT/JRootClass.hh"


/**
 * \author mdejong
 */

namespace JAANET {}
namespace JPP { using namespace JAANET; }

namespace JAANET {

  using JLANG::JException;
  using JLANG::JEquationParameters;
  
  
  static const char AANET_TAG_SEPARATOR  =  '_';     //!< Separator for AAnet tag extension for multiple tags ("_<counter>").

  class JHead;               // forward declaration for copy methods


  /**
   * Copy header from <tt>from</tt> to <tt>to</tt>.
   *
   * \param  from            header
   * \param  to              header
   */
  extern void copy(const JHead& from, Head& to);


  /**
   * Copy header from <tt>from</tt> to <tt>to</tt>.
   *
   * \param  from            header
   * \param  to              header
   */
  extern void copy(const Head& from, JHead& to);


  /**
   * Get tag without aanet extension "_<counter>" for identical tags.
   *
   * \param  tag             tag
   * \return                 tag
   */
  inline std::string getTag(const std::string& tag)
  {
    using namespace std;
    
    const string::size_type pos = tag.find(AANET_TAG_SEPARATOR);
    
    if (pos != string::npos) {
      
      for (string::size_type i = pos + 1; i != tag.size(); ++i) {
        if (!isdigit(tag[i])) {
          return tag;
        }
      }
      
      return tag.substr(0, pos);
    }
    
    return tag;
  }


  /**
   * Start of run record.
   */
  struct start_run {
    /**
     * Default constructor.
     */
    start_run() : 
      run_id(0)
    {}
    
    int run_id;                ///< MC run number
    
    ClassDefNV(start_run,1);
  };
  

  /**
   * General purpose string class.
   */
  struct String {
    /**
     * Default constructor.
     */
    String() : 
      buffer()
    {}

    /**
     * Test match.
     *
     * \param  object               string
     * \return                      true if matches; else false
     */
    inline bool match(const String& object) const
    {
      return buffer == object.buffer;
    }

    std::string buffer;        ///< General purpose name

    ClassDefNV(String,1);
  };


  /**
   * Detector file.
   */
  struct detector :
    public String
  {
    ClassDefNV(detector,1);
  };


  /**
   * Muon descriptor file.
   */
  struct muon_desc_file :
    public String
  { 
    ClassDefNV(muon_desc_file,1);
  };


  /**
   * Target.
   */
  struct target :
    public String
  { 
    ClassDefNV(target,1);
  };


  /**
   * Neutrino cross section file.
   */
  struct XSecFile : 
    public String
  { 
    ClassDefNV(XSecFile,1);
  };


  /**
   * Drawing.
   */
  struct drawing :
    public String
  { 
    ClassDefNV(drawing,1);
  };


  /**
   * General purpose class of phase space generation.
   */
  struct cut {
    /**
     * Default constructor.
     */
    cut() :
      Emin(0),
      Emax(0),
      cosTmin(0),
      cosTmax(0)
    {}

    /**
     * Comparison.
     *
     * \param  object               cut
     * \return                      true if this cut less than given cut; else false
     */
    inline bool less(const cut& object) const
    {
      if (Emin    == object.Emin    &&
          Emax    == object.Emax)
        return cosTmax < object.cosTmin;
      else
        return Emax    < object.Emin;
    }

    /**
     * Test match.
     *
     * \param  object               cut
     * \return                      true if matches; else false
     */
    inline bool match(const cut& object) const
    {
      return (Emin    == object.Emin    &&
              Emax    == object.Emax    &&
              cosTmin == object.cosTmin &&
              cosTmax == object.cosTmax);
    }

    double Emin;               ///< Minimal energy [GeV]
    double Emax;               ///< Maximal energy [GeV]
    double cosTmin;            ///< Minimal cosine zenith angle
    double cosTmax;            ///< Maximal cosine zenith angle

    ClassDefNV(cut,1);
  };


  /**
   * Phase space of primary particle.
   */
  struct cut_primary : 
    public cut 
  { 
    ClassDefNV(cut_primary,1);
  };


  /**
   * Phase space of atmospheric muon generation.
   */
  struct cut_seamuon : 
    public cut 
  { 
    ClassDefNV(cut_seamuon,1);
  };


  /**
   * Phase space of incoming particle
   */
  struct cut_in :
    public cut
  { 
    ClassDefNV(cut_in,1);
  };


  /**
   * Phase space of incident neutrino.
   */
  struct cut_nu : 
    public cut 
  { 
    ClassDefNV(cut_nu,1);
  };


  /**
   * Description of Monte Carlo event generation applications.
   */
  struct generator {
    /**
     * Default constructor.
     */
    generator() :
      program(),
      version(),
      date(),
      time()
    {}

    /**
     * Test match.
     *
     * Note that only the name of program is matched.
     *
     * \param  object               generator
     * \return                      true if matches; else false
     */
    inline bool match(const generator& object) const
    {
      return program == object.program;
    }

    std::string program;       ///< program name
    std::string version;       ///< program version
    std::string date;          ///< processing date
    std::string time;          ///< processing time

    ClassDefNV(generator,1);
  };


  /**
   * Generator for neutrino interaction
   */
  struct physics : 
    public generator
  { 
    ClassDefNV(physics,1);
  };


  /**
   * Generator for simulation.
   */
  struct simul : 
    public generator
  { 
    ClassDefNV(simul,1);
  };

  
  /**
   * Neutrino energy spectrum.
   */
  struct spectrum {
    /**
     * Default constructor.
     */
    spectrum() : 
      alpha(0)
    {}

    /**
     * Test match.
     *
     * \param  object               spectrum
     * \return                      true if matches; else false
     */
    inline bool match(const spectrum& object) const
    {
      return (alpha == object.alpha);
    }

    double alpha;              ///< Energy spectrum: \f$ \Phi \propto E^{-\alpha} \f$

    ClassDefNV(spectrum,1);
  };


  /**
   * The cylinder used for photon tracking.
   */
  struct can {
    /**
     * Default constructor.
     */
    can() :
      zmin(0),
      zmax(0),
      r(0)
    {}

    /**
     * Test match.
     *
     * \param  object               can
     * \return                      true if matches; else false
     */
    inline bool match(const can& object) const
    {
      return (zmin == object.zmin &&
              zmax == object.zmax &&
              r    == object.r);
    }

    double zmin;               ///< Bottom [m]
    double zmax;               ///< Top    [m]
    double r;                  ///< Radius [m]

    ClassDefNV(can,1);
  };


  /**
   * The fixe cylinder used for photon tracking.
   */
  struct fixedcan {
    /**
     * Default constructor.
     */
    fixedcan() :
      xcenter(0),
      ycenter(0),
      zmin(0),
      zmax(0),
      radius(0)
    {}

    /**
     * Test match.
     *
     * \param  object               can
     * \return                      true if matches; else false
     */
    inline bool match(const fixedcan& object) const
    {
      return (xcenter == object.xcenter &&
              xcenter == object.ycenter &&
              zmin    == object.zmin    &&
              zmax    == object.zmax    &&
              radius  == object.radius);
    }

    double xcenter;            ///< x-center [m]
    double ycenter;            ///< y-center [m]
    double zmin;               ///< Bottom [m]
    double zmax;               ///< Top    [m]
    double radius;             ///< Radius [m]

    ClassDefNV(fixedcan,1);
  };


  /**
   * Neutrino vertex volume.
   */
  struct genvol {
    /**
     * Default constructor.
     */
    genvol() :
      zmin(0),
      zmax(0),
      r(0),
      volume(0),
      numberOfEvents(0)
    {}

    /**
     * Test match.
     *
     * \param  object               generation volume
     * \return                      true if matches; else false
     */
    inline bool match(const genvol& object) const
    {
      return (zmin == object.zmin &&
              zmax == object.zmax &&
              r    == object.r);
    }

    /**
     * Addition.
     *
     * \param  object               generation volume
     * \return                      this generation volume
     */
    inline genvol& add(const genvol& object)
    {
      numberOfEvents += object.numberOfEvents;

      return *this;
    }

    /**
     * Scale.
     *
     * \param  factor               1D scale factor
     * \param  z                    z position
     * \return                      this genvol
     */
    inline genvol& mul(const double factor,
                       const double z = 0.0)
    {
      zmin    = (zmin - z) * factor  +  z;
      zmax    = (zmax - z) * factor  +  z;
      r      *= factor;
      volume *= factor * factor * factor;
 
      return *this;
    }

    double zmin;               ///< Bottom [m]
    double zmax;               ///< Top    [m]
    double r;                  ///< Radius [m]
    double volume;             ///< Volume [m^3]
    double numberOfEvents;     ///< Number of events

    ClassDefNV(genvol,1);
  };


  /**
   * Coordinate origin.
   */
  struct coord_origin :
    public Vec
  {
    /**
     * Default constructor.
     */
    coord_origin() : 
      Vec()
    {}

    /**
     * Constructor.
     *
     * \param  x                    x position
     * \param  y                    y position
     * \param  z                    z position
     */
    coord_origin(const double x,
                 const double y,
                 const double z) :
      Vec(x,y,z)
    {}

    /**
     * Test match.
     *
     * \param  object               coordinate origin
     * \return                      true if matches; else false
     */
    inline bool match(const coord_origin& object) const
    {
      return (x == object.x &&
              y == object.y &&
              z == object.z);
    }

    ClassDefNV(coord_origin,1);
  };


  /**
   * Phase space for incident neutrino.
   */
  struct genhencut {
    /**
     * Default constructor.
     */
    genhencut() :
      gDir(0),
      Emin(0)
    {}

    double gDir;               ///< Angle
    double Emin;               ///< Minimal energy [GeV]

    ClassDefNV(genhencut,1);
  };


  /**
   * Normlisation of CORSIKA events.
   */
  struct norma {
    /**
     * Default constructor.
     */
    norma() :
      primaryFlux(0),
      numberOfPrimaries(0)
    {}

    /**
     * Test match.
     *
     * \param  object               normalisation
     * \return                      true if matches; else false
     */
    inline bool match(const norma& object) const
    {
      return (primaryFlux == object.primaryFlux);
    }

    /**
     * Addition.
     *
     * \param  object               normalisation
     * \return                      this normalisation
     */
    inline norma& add(const norma& object)
    {
      numberOfPrimaries += object.numberOfPrimaries;

      return *this;
    }

    double primaryFlux;        ///< Primary flux 
    double numberOfPrimaries;  ///< Number of primaries

    ClassDefNV(norma,1);
  };


  /**
   * Normalisation of MUPAGE events.
   */
  struct livetime {
    /**
     * Default constructor.
     */
    livetime() :
      numberOfSeconds(0),
      errorOfSeconds(0)
    {}

    /**
     * Test match.
     *
     * \param  object               live time
     * \return                      true if matches; else false
     */
    inline bool match(const livetime& object) const
    {
      return ((numberOfSeconds == 0.0 && object.numberOfSeconds == 0.0) ||
              (numberOfSeconds >  0.0 && object.numberOfSeconds >  0.0));
    }

    /**
     * Addition.
     *
     * \param  object               live time
     * \return                      this live time
     */
    inline livetime& add(const livetime& object)
    {
      numberOfSeconds += object.numberOfSeconds;
      errorOfSeconds   = sqrt(errorOfSeconds        * errorOfSeconds  +  
                              object.errorOfSeconds * object.errorOfSeconds);

      return *this;
    }

    /**
     * Scale.
     *
     * \param  factor               1D scale factor
     * \return                      this livetime
     */
    inline livetime& mul(const double factor)
    {
      numberOfSeconds /= (factor * factor);
      errorOfSeconds  /= (factor * factor);
      
      return *this;
    }

    double numberOfSeconds;    ///< Live time [s]
    double errorOfSeconds;     ///< Uncertainty on live time [s]

    ClassDefNV(livetime,1);
  };


  /**
   * Neutrino flux.
   */
  struct flux {
  public:
    /**
     * Default constructor.
     */
    flux() :
      type(0),
      key(),
      file_1(),
      file_2()
    {}

    int         type;          ///< Type
    std::string key;           ///< Key
    std::string file_1;        ///< File name
    std::string file_2;        ///< File name

    ClassDefNV(flux,1);
  };


  /**
   * The bottom of the sea.
   */
  struct seabottom {
    /**
     * Default constructor.
     */
    seabottom() :
      z(0)
    {}

    /**
     * Test match.
     *
     * \param  object               sea bottom
     * \return                      true if matches; else false
     */
    inline bool match(const seabottom& object) const
    {
      return (z == object.z);
    }

    double z;                  ///< Sea bottom [m]

    ClassDefNV(seabottom,1);
  };


  /**
   * Livetime of DAQ data.
   */
  struct DAQ {
  public:
    /**
     * Default constructor.
     */
    DAQ() :
      livetime_s(0.0)
    {}

    /**
     * Test match.
     *
     * \param  object               DAQ
     * \return                      true if matches; else false
     */
    inline bool match(const DAQ& object) const
    {
      return ((livetime_s == 0.0 && object.livetime_s == 0.0) ||
              (livetime_s >  0.0 && object.livetime_s >  0.0));
    }

    /**
     * Addition.
     *
     * \param  object               DAQ
     * \return                      this DAQ
     */
    inline DAQ& add(const DAQ& object)
    {
      livetime_s += object.livetime_s;

      return *this;
    }

    double livetime_s;         ///< Live time [s]

    ClassDefNV(DAQ,1);
  };


  /**
   * Scaling time of neutrino interaction generators.
   */
  struct tgen {
    /**
     * Default constructor.
     */
    tgen() : 
      numberOfSeconds(0)
    {}
    
    double numberOfSeconds;                ///< Time in seconds
    
    ClassDefNV(tgen,1);
  };


  /**
   * Primary particle.
   */
  struct primary {
    /**
     * Default constructor.
     */
    primary() : 
      type(0)
    {}

    /**
     * Comparison.
     *
     * \param  object               primary particle
     * \return                      true if this primary particle is less than given primary particle; else false
     */
    inline bool less(const primary& object) const
    {
      return type < object.type;
    }

    /**
     * Test match.
     *
     * \param  object               primary particle
     * \return                      true if matches; else false
     */
    inline bool match(const primary& object) const
    {
      return (type == object.type);
    }
   
    int type;                              ///< Particle type
    
    ClassDefNV(primary,1);
  };
  

  /**
   * End of event record.
   */
  struct end_event {
    /**
     * Default constructor.
     */
    end_event()
    {}

    ClassDefNV(end_event,1);
  };


  /**
   * Monte Carlo run header.
   *
   * This class extends the Head class so that the data from specific tags 
   * can be promoted to concrete data types.
   *
   * Note that for the copy of new JHead data (e.g. data not obtained via a previous JAANET::copy) to become effective, 
   * the key words in the corresponding map of the Head class should be set. \n
   * To this end, member method JHead::push can be used.
   */
  struct JHead : 
    public Head
  {
    /**
     * Generators.
     */
    static const std::string GENHEN;
    static const std::string GENIE;
    static const std::string GSEAGEN;
    static const std::string MUPAGE;
    static const std::string JSIRENE;
    static const std::string KM3;
    static const std::string KM3SIM;

    
    /**
     * Default constructor.
     */
    JHead()
    {}


    /**
     * Copy constructor.
     *
     * \param  header          header
     */
    JHead(const Head& header)
    {
      copy(header, *this);
    }


    /**
     * Virtual destructor.
     */
    virtual ~JHead()
    {}


    /**
     * Get header.
     *
     * \return                 header
     */
    const JHead& getHeader() const
    {
      return static_cast<const JHead&>(*this);
    }


    /**
     * Get header.
     *
     * \return                 header
     */
    JHead& getHeader()
    {
      return static_cast<JHead&>(*this);
    }


    /**
     * Set header.
     *
     * \param  header          header
     */
    void setHeader(const JHead& header)
    {
      static_cast<JHead&>(*this) = header;
    }


    /**
     * Check validity of given data member in Head.
     *
     * The validity is defined by the appearence of the name of the data member of the underlying map.
     *
     * \param  pd              pointer to data member
     * \return                 true if valid; else false
     */
    template<class T>
    inline bool is_valid(T JHead::*pd) const
    {
      return (this->pull(pd) != this->end());
    }


    /**
     * Pull given data member from Head.
     *
     * \param  pd              pointer to data member
     * \return                 iterator of Head
     */
    template<class T>
    inline const_iterator pull(T JHead::*pd) const
    {
      return this->find(JROOT::getDataMember(pd)->GetName());
    }


    /**
     * Pull given data member from Head.
     *
     * \param  pd              pointer to data member
     * \return                 iterator of Head
     */
    template<class T>
    inline iterator pull(T JHead::*pd)
    {
      return this->find(JROOT::getDataMember(pd)->GetName());
    }


    /**
     * Push given data member to Head.
     *
     * \param  pd              pointer to data member
     */
    template<class T>
    inline void push(T JHead::*pd) 
    {
      (*this)[JROOT::getDataMember(pd)->GetName()] = "";
    }


    /**
     * Remove given data member from Head.
     *
     * \param  pd              pointer to data member
     */
    template<class T>
    inline void erase(T JHead::*pd) 
    {
      iterator p = this->pull(pd);

      if (p != this->end()) {

        this->*pd = T();

        static_cast<Head*>(this)->erase(p);
      }
    }


    /**
     * Test match of headers.
     *
     * Note that if option is set to <tt>false</tt>,
     * the match applies only to data which have a corresponding entry in the underlying map of the given header.
     *
     * \param  header          second header
     * \param  option          option
     * \return                 true if matches; else false
     */
    inline bool match(const JHead& header, const bool option = true) const
    {
      return (match(*this, header, option, &JHead::cut_primary)       &&
              match(*this, header, option, &JHead::cut_seamuon)       &&
              match(*this, header, option, &JHead::cut_in)            &&
              match(*this, header, option, &JHead::cut_nu)            &&
              match(*this, header, option, &JHead::physics)           &&
              match(*this, header, option, &JHead::simul)             &&
              match(*this, header, option, &JHead::spectrum)          &&
              match(*this, header, option, &JHead::can)               &&
              match(*this, header, option, &JHead::fixedcan)          &&
              match(*this, header, option, &JHead::genvol)            &&
              match(*this, header, option, &JHead::coord_origin)      &&
              match(*this, header, option, &JHead::norma)             &&
              match(*this, header, option, &JHead::livetime)          &&
              match(*this, header, option, &JHead::seabottom)         &&
              match(*this, header, option, &JHead::primary)           &&
              match(*this, header, option, &JHead::DAQ));
    }


    /**
     * Comparison of headers.
     *
     * \param  header          header
     * \return                 true if this header less than given header; else false
     */
    inline bool less(const JHead& header) const
    {
      if      (primary.less(header.primary))
        return true;
      else if (header.primary.less(primary))
        return false;
      else
        return cut_primary.less(header.cut_primary);
    }


    /**
     * Addition of headers.
     *
     * \param  header          header
     * \return                 this header
     */
    inline JHead& add(const JHead& header)
    {
      if (match(header)) {

        genvol  .add(header.genvol);
        norma   .add(header.norma);
        livetime.add(header.livetime);
        DAQ     .add(header.DAQ);

      } else {

        THROW(JException, "JHead::add() headers do not match.");
      }

      return *this;
    }


    /**
     * Equal operator.
     *
     * Note that this operator uses the JHead::match method.
     *
     * \param  first        first  header
     * \param  second       second header
     * \return              true if two headers are equal; else false
     */
    friend inline bool operator==(const JHead& first,
                                  const JHead& second)
    {
      return first.match(second);
    }


    /**
     * Less than operator.
     *
     * \param  first        first  header
     * \param  second       second header
     * \return              true if first header is less than second header; else false
     */
    friend inline bool operator<(const JHead& first,
                                 const JHead& second)
    {
      return first.less(second);
    }


    JAANET::start_run            start_run;        // first data member
    JAANET::XSecFile             XSecFile;
    JAANET::drawing              drawing;
    JAANET::detector             detector;
    JAANET::muon_desc_file       muon_desc_file;
    JAANET::target               target;
    std::vector<JAANET::physics> physics;
    std::vector<JAANET::simul>   simul;
    JAANET::cut_primary          cut_primary;
    JAANET::cut_seamuon          cut_seamuon;
    JAANET::cut_in               cut_in;
    JAANET::cut_nu               cut_nu;
    JAANET::spectrum             spectrum;
    JAANET::can                  can;
    JAANET::fixedcan             fixedcan;
    JAANET::genvol               genvol;
    JAANET::coord_origin         coord_origin;
    JAANET::genhencut            genhencut;
    JAANET::norma                norma;
    JAANET::livetime             livetime;
    JAANET::seabottom            seabottom;
    JAANET::DAQ                  DAQ;
    JAANET::tgen                 tgen;
    JAANET::primary              primary;
    JAANET::end_event            end_event;        // last data member

    
    /**
     * Get equation parameters corresponding to Monte Carlo ASCII format, i.e:
     * <pre>
     *    <key>: <value> [<value>]*
     *    <key>: <value> [<value>]*
     * </pre>
     *
     * \return                         equation parameters
     */
    static inline JLANG::JEquationParameters& getEquationParameters()
    {
      static JLANG::JEquationParameters parameters(":", "\n", "", "");
      
      return parameters;
    }
    

    /**
     * Set equation parameters.
     *
     * \param  equation                equation parameters
     */
    static inline void setEquationParameters(const JLANG::JEquationParameters& equation)
    {
      getEquationParameters() = equation;
    }


    /**
     * Read header from input.
     *
     * \param  in         input stream
     * \return            input stream
     */
    std::istream& read(std::istream& in);


    /**
     * Write header to output.
     *
     * \param  out        output stream
     * \return            output stream
     */
    std::ostream& write(std::ostream& out) const;


    /**
     * Print header to output.
     *
     * \param  out        output stream
     * \return            output stream
     */
    std::ostream& print(std::ostream& out) const;


    ClassDef(JHead,3);

  private:
    /**
     * Test match.
     *
     * \param  first           first  object
     * \param  second          second object
     * \return                 true if matches; else false
     */
    template<class T>
    static inline bool match(const T& first,
                             const T& second)
    {
      return first.match(second);
    }


    /**
     * Test one container is subset of other container or vice versa.
     *
     * \param  first           first  object
     * \param  second          second object
     * \return                 true if first (second) is subset of second (first); else false
     */
    template<class T>
    static inline bool match(const std::vector<T>& first,
                             const std::vector<T>& second)
    {
      if (!first.empty() && !second.empty()) {

        size_t ns = 0;
      
        for (typename std::vector<T>::const_iterator ix = first.begin(); ix != first.end(); ++ix) {
          for (typename std::vector<T>::const_iterator iy = second.begin(); iy != second.end(); ++iy) {
            if (match(*ix, *iy)) {
              ++ns;
            }
          }
        }
      
        return ns == first.size() || ns == second.size();

      }

      return first.empty() && second.empty();
    }

    
    /**
     * Test match of given data member of headers.
     *
     * Note that if option is set to <tt>false</tt>,
     * the match applies only to data which have a corresponding tag in the underlying map of the second header.
     *
     * \param  first           first  header
     * \param  second          second header
     * \param  option          option
     * \param  pd              pointer to data member     
     * \return                 true if matches; else false
     */
    template<class T>
    static inline bool match(const JHead& first,
                             const JHead& second,
                             const bool   option,
                             T JHead::*pd)
    {
      if (option || second.pull(pd) != second.end())
        return match(first.*pd, second.*pd);
      else
        return true;
    }
  };


  /**
   * Equal operator.
   *
   * Note that this operator uses the JHead::match method.
   *
   * \param  first        first  header
   * \param  second       second header
   * \return              true if two headers are equal; else false
   */
  inline bool operator==(const Head& first,
                         const Head& second)
  {
    return JHead(first).match(JHead(second));
  }
  
  
  /**
   * Less than operator.
   *
   * \param  first        first  header
   * \param  second       second header
   * \return              true if first header is less than second header; else false
   */
  inline bool operator<(const Head& first,
                        const Head& second)
  {
    return JHead(first).less(JHead(second));
  }
}


/**
 * Read header from input.
 *
 * \param  in         input stream
 * \param  header     header
 * \return            input stream
 */
inline std::istream& operator>>(std::istream& in, JAANET::JHead& header)
{
  return header.read(in);
}


/**
 * Write header to output.
 *
 * \param  out        output stream
 * \param  header     header
 * \return            output stream
 */
inline std::ostream& operator<<(std::ostream& out, const JAANET::JHead& header)
{
  return header.write(out);
}
#endif