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 <sstream>

#include "TObject.h"

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

#include "JLang/JException.hh"
#include "JLang/JEquationParameters.hh"
#include "JLang/JLangToolkit.hh"

#include "JTools/JRange.hh"

#include "JROOT/JRootClass.hh"


/**
 * \author mdejong
 */

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

namespace JAANET {

  using JLANG::JException;
  using JLANG::JEquationParameters;
  
  /**
   * Type definition of range.
   */
  struct JRange_t :
    public JTOOLS::JRange<double>
  {
    /**
     * Default constructor.
     */
    JRange_t() :
      JTOOLS::JRange<double>()
    {}


    /**
     * Constructor.
     *
     * \param  x                lower limit
     * \param  y                upper limit
     */
    JRange_t(double x,
             double y) :
      JTOOLS::JRange<double>(x, y)
    {}
  };

  
  static const char AANET_TAG_SEPARATOR  =  '_';     //!< Separator for tag extension of multiple tags in Head ("_<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;
  }


  /**
   * Get tag with aanet extension "_<counter>" for identical tags.
   *
   * \param  tag             tag
   * \param  counter         counter
   * \return                 tag
   */
  inline std::string getTag(const std::string& tag, const int counter)
  {
    std::ostringstream os;

    os << tag << AANET_TAG_SEPARATOR << counter;

    return os.str();
  }


  /**
   * 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 !(*this).less(object) && !(object).less(*this);
    }

    /**
     * Comparison.
     *
     * \param  object               string
     * \return                      true if this string less than given string; else false
     */
    inline bool less(const String& object) const
    {
      using namespace std;

      istringstream i0(this ->buffer);
      istringstream i1(object.buffer);

      vector<string> v0;
      vector<string> v1;

      copy(istream_iterator<string>(i0), istream_iterator<string>(), back_inserter(v0));
      copy(istream_iterator<string>(i1), istream_iterator<string>(), back_inserter(v1));

      return v0 < v1;
    }

    /**
     * Read string from input stream.
     *
     * \param  in                   input stream
     * \param  object               string
     * \return                      input stream
     */
    friend inline std::istream& operator>>(std::istream& in, String& object)
    {
      return std::getline(in, object.buffer);
    }

    /**
     * Write string to output stream.
     *
     * \param  out                  output stream
     * \param  object               string
     * \return                      output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const String& object)
    {
      return out << object.buffer;
    }

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

    ClassDefNV(String,1);
  };


  /**
   * Detector file.
   */
  struct detector
  {
    /**
     * Default constructor.
     */
    detector()
    {}

    /**
     * Test match.
     *
     * \param  object               detector
     * \return                      true if matches; else false
     */
    inline bool match(const detector& object) const
    {
      return this->program == object.program && this->filename == object.filename; 
    }

    /**
     * Read detector from input stream.
     *
     * \param  in                   input stream
     * \param  object               detector
     * \return                      input stream
     */
    friend inline std::istream& operator>>(std::istream& in, detector& object)
    {
      using namespace JPP;

      in >> object.program >> object.filename;

      if (!in || is_integer(object.filename)) { // legacy
        object.filename = object.program;
        object.program  = "?";
      }

      return in;
    }

    /**
     * Write detector to output stream.
     *
     * \param  out                  output stream
     * \param  object               detector
     * \return                      output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const detector& object)
    {
      return out << object.program << ' ' << object.filename;
    }

    std::string program;
    std::string filename;
    
    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() :
      E   (0.0, 0.0),
      cosT(0.0, 0.0)
    {}

    /**
     * Constructor.
     * 
     * \param  _E            energy range
     * \param  _cosT         cosine zenith angle range
     */
    cut(const JRange_t& _E,
        const JRange_t& _cosT) :
      E   (_E),
      cosT(_cosT)
    {}

    /**
     * Constructor.
     *
     * \param  Emin                energy range lower bound
     * \param  Emax                energy range upper bound
     * \param  cosTmin             cosine zenith angle lower bound
     * \param  cosTmax             cosine zenith angle upper bound
     */
    cut(const double Emin,
        const double Emax,
        const double cosTmin,
        const double cosTmax) :
      E   (Emin,    Emax),
      cosT(cosTmin, cosTmax)
    {}

    /**
     * Comparison.
     *
     * \param  object               cut
     * \return                      true if this cut less than given cut; else false
     */
    inline bool less(const cut& object) const
    {
      if (E.getLowerLimit() == object.E.getLowerLimit()) {

        if (E.getUpperLimit() == object.E.getUpperLimit()) {

          if (cosT.getLowerLimit() == object.cosT.getLowerLimit()) {
            return cosT.getUpperLimit() < object.cosT.getUpperLimit();
          } else {
            return cosT.getLowerLimit() < object.cosT.getLowerLimit();
          }

        } else {

          return E.getUpperLimit() < object.E.getUpperLimit();
        }

      } else {

        return E.getLowerLimit() < object.E.getLowerLimit();
      }
    }

    /**
     * Test match.
     *
     * \param  object               cut
     * \return                      true if matches; else false
     */
    inline bool match(const cut& object) const
    {
      return (E   .equals(object.E) &&
              cosT.equals(object.cosT));
    }
    
    JRange_t E;     ///< Energy range                [GeV]
    JRange_t cosT;  ///< Cosine zenith angle range   
    
    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()
    {}

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

    /**
     * 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);
  };


  /**
   * Physics information.
   */
  struct physics : 
    public String
  { 
    ClassDefNV(physics,1);
  };


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

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

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

    /**
     * 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 fixed 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 &&
              ycenter == 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)
    {}

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

    /**
     * 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)
    {}

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

    /**
     * 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()
    {}

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

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

    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);
  };


  /**
   * Depth.
   */
  struct depth {
    /**
     * Default constructor.
     */
    depth() :
      z(0)
    {}

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

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

    ClassDefNV(depth,1);
  };


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

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

    /**
     * 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);
  };


  /**
   * Time duration of event generation.
   */
  struct tgen {
    /**
     * Default constructor.
     */
    tgen() : 
      numberOfSeconds(0)
    {}

    /**
     * Test match.
     *
     * \param  object               time duration
     * \return                      true if matches; else false
     */
    inline bool match(const tgen& object) const
    {
      return this->numberOfSeconds == object.numberOfSeconds;
    }
    
    double numberOfSeconds;                ///< Time in seconds
    
    ClassDefNV(tgen,1);
  };


  /**
   * UTC time interval for event generation.
   */
  struct time_interval {
    /**
     * Default constructor.
     */
    time_interval() : 
      t1(0),
      t2(0)
    {}

    /**
     * Test match.
     *
     * \param  object               time interval
     * \return                      true if matches; else false
     */
    inline bool match(const time_interval& object) const
    {
      return (this->t1 == object.t1 && this->t2 == object.t2);
    }
    
    double t1;                             ///< Start time in seconds
    double t2;                             ///< Stop  time in seconds
    
    ClassDefNV(time_interval,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
  {
    /**
     * 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 JHead.
     *
     * The validity is defined by the presence of the name of the data member in 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());
    }


    /**
     * Check validity of given data member in JHead.
     *
     * The validity is defined by difference between actual and default value.
     *
     * \param  object          object 
     * \return                 true if valid; else false
     */
    template<class T>
    static bool is_valid(const T& object)
    {
      static const T value;

      return (object.less(value) || value.less(object));
    }


    /**
     * 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);
      }
    }


    /**
     * Get number of matching fields.
     *
     * \param  header          second header
     * \return                 number of matching header fields
     */
    inline int getNumberOfMatches(const JHead& header) const
    {
      return ((int) match(*this, header, &JHead::cut_primary)       +
              (int) match(*this, header, &JHead::cut_seamuon)       +
              (int) match(*this, header, &JHead::cut_in)            +
              (int) match(*this, header, &JHead::cut_nu)            +
              (int) match(*this, header, &JHead::simul)             +
              (int) match(*this, header, &JHead::spectrum)          +
              (int) match(*this, header, &JHead::can)               +
              (int) match(*this, header, &JHead::fixedcan)          +
              (int) match(*this, header, &JHead::genvol)            +
              (int) match(*this, header, &JHead::coord_origin)      +
              (int) match(*this, header, &JHead::norma)             +
              (int) match(*this, header, &JHead::livetime)          +   
              (int) match(*this, header, &JHead::seabottom)         +
              (int) match(*this, header, &JHead::depth)             +
              (int) match(*this, header, &JHead::tgen)              +
              (int) match(*this, header, &JHead::time_interval)     +
              (int) match(*this, header, &JHead::primary)           +
              (int) match(*this, header, &JHead::flux)              +   
              (int) match(*this, header, &JHead::DAQ));
    }


    /**
     * Test match of headers.
     *
     * \param  header          second header
     * \return                 true if all header fields match; else false
     */
    inline bool match(const JHead& header) const
    {
      return getNumberOfMatches(header) == getMaximumNumberOfMatches();
    }


    /**
     * Comparison of headers.
     *
     * \param  header          header
     * \return                 true if this header less than given header; else false
     */
    inline bool less(const JHead& header) const
    {
#define RETURN_IF_DIFFERENT(A, B)       \
      if (less(A,B)) { return true;  }  \
      if (less(B,A)) { return false; } 

      // compare physics

      RETURN_IF_DIFFERENT(this->physics,     header.physics);
      
      // compare simulation
      
      RETURN_IF_DIFFERENT(this->simul,       header.simul);
      
      // compare generation data
      
      RETURN_IF_DIFFERENT(this->primary,     header.primary);
      RETURN_IF_DIFFERENT(this->flux,        header.flux);
      RETURN_IF_DIFFERENT(this->spectrum,    header.spectrum);
      RETURN_IF_DIFFERENT(this->cut_primary, header.cut_primary);
      RETURN_IF_DIFFERENT(this->cut_seamuon, header.cut_seamuon);
      RETURN_IF_DIFFERENT(this->cut_in,      header.cut_in);
      RETURN_IF_DIFFERENT(this->cut_nu,      header.cut_nu);
      RETURN_IF_DIFFERENT(this->genvol,      header.genvol);

      // compare compatibility

      if (is_valid(this->livetime) == is_valid(header.livetime)  &&
          is_valid(this->DAQ)      == is_valid(header.DAQ))      { 
        return false;
      }

      THROW(JException, "JHead::less() headers do not compare.");

#undef RETURN_IF_DIFFERENT
    }


    /**
     * 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;
    std::vector<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::depth                 depth;
    JAANET::DAQ                   DAQ;
    JAANET::tgen                  tgen;
    JAANET::time_interval         time_interval;
    JAANET::primary               primary;
    std::vector<JAANET::flux>     flux;
    JAANET::end_event             end_event;        // last data member


    /**
     * Get maximum number of matching header fields.
     *
     * \return                 maximum number of matching header fields
     */
    static inline const int getMaximumNumberOfMatches()
    {
      static const JHead header;
      static const int   value = header.getNumberOfMatches(header);
      
      return value;
    }

    
    /**
     * 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:
    /**
     * Comparison.
     *
     * \param  first           first  object
     * \param  second          second object
     * \return                 true if first less than second; else false
     */
    template<class T>
    static inline bool less(const T& first,
                            const T& second)
    {
      return first.less(second);
    }


    /**
     * 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);
    }


    /**
     * Comparison of containers.
     * It is assumed that the containers are ordered in the same way.
     *
     * \param  first           first  object
     * \param  second          second object
     * \return                 true if first is less than second; else false
     */
    template<class T>
    static inline bool less(const std::vector<T>& first,
                            const std::vector<T>& second)
    {
      if (first.size() == second.size()) {

        for (size_t i = 0; i != first.size(); ++i) {
          if (less(first[i], second[i])) {
            return true;
          }
        }

        return false;

      } else {

        return first.size() < second.size();
      }
    }


    /**
     * Test is containers match.
     * It is assumed that the containers are ordered in the same way.
     *
     * \param  first           first  object
     * \param  second          second object
     * \return                 true if matches; else false
     */
    template<class T>
    static inline bool match(const std::vector<T>& first,
                             const std::vector<T>& second)
    {
      for (size_t i = 0; i != first.size() && i != second.size(); ++i) {
        if (!match(first[i], second[i])) {
          return false;
        }
      }

      return first.size() == second.size();
    }

    
    /**
     * Test match of given data member of headers.
     *
     * \param  first           first  header
     * \param  second          second header
     * \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,
                             T JHead::*pd)
    {
      return match(first.*pd, second.*pd);
    }
  };
  
  
  /**
   * 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