#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/JUUID.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::JUUID;
  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);
  };


  /**
   * Calibration.
   */
  struct calibration :
    public String
  {
    static const std::string statical()  { return "statical"; }
    static const std::string dynamical() { return "dynamical"; }

    ClassDefNV(calibration,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    &&
	      volume == object.volume);
    }

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


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

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

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

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

      return *this;
    }

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

    ClassDefNV(K40,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.
   */
  class JHead : 
    public Head
  {
  public:
    /**
     * Default constructor.
     */
    JHead()
    {
      createUUID();
    }


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


    /**
     * Create UUID if not already set.
     */
    void createUUID()
    {
      if (!is_valid(&JHead::UUID)) {
	this->UUID = JUUID::rndm();
	this->push(&JHead::UUID);
      }
    }

    
    /**
     * 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()] = "";
    }


    /**
     * Push all data members to Head.
     */
    void push();


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

      if (p != this->end()) {
	static_cast<Head*>(this)->erase(p);
      }
    }


    /**
     * Get matching fields.
     *
     * \param  header          header
     * \return                 header with matching fields
     */
    inline JHead getMatch(const JHead& header) const
    {
#define IF_MATCH(A, B, C, D)	        \
      if (match(B,C,D)) { A.push(D);  }	\
      else              { A.erase(D); }
      
      JHead buffer(*this);

      buffer.clear();

      buffer.createUUID();

      IF_MATCH(buffer, *this, header, &JHead::cut_primary);
      IF_MATCH(buffer, *this, header, &JHead::cut_seamuon);
      IF_MATCH(buffer, *this, header, &JHead::cut_in);
      IF_MATCH(buffer, *this, header, &JHead::cut_nu);
      IF_MATCH(buffer, *this, header, &JHead::simul);
      IF_MATCH(buffer, *this, header, &JHead::physics);
      IF_MATCH(buffer, *this, header, &JHead::spectrum);
      IF_MATCH(buffer, *this, header, &JHead::can);
      IF_MATCH(buffer, *this, header, &JHead::fixedcan);
      IF_MATCH(buffer, *this, header, &JHead::genvol);
      IF_MATCH(buffer, *this, header, &JHead::coord_origin);

      IF_MATCH(buffer, *this, header, &JHead::norma);
      IF_MATCH(buffer, *this, header, &JHead::livetime);

      IF_MATCH(buffer, *this, header, &JHead::seabottom);
      IF_MATCH(buffer, *this, header, &JHead::depth);
      IF_MATCH(buffer, *this, header, &JHead::tgen);
      //IF_MATCH(buffer, *this, header, &JHead::time_interval);

      IF_MATCH(buffer, *this, header, &JHead::primary);
      IF_MATCH(buffer, *this, header, &JHead::flux);
      IF_MATCH(buffer, *this, header, &JHead::DAQ);      
      IF_MATCH(buffer, *this, header, &JHead::K40);
      IF_MATCH(buffer, *this, header, &JHead::target);

      return buffer;
      
#undef IF_MATCH
    }


    /**
     * Get number of matching fields.
     *
     * \param  header          header
     * \return                 number of matching header fields
     */
    inline size_t getNumberOfMatches(const JHead& header) const
    {
      const JHead head = getMatch(header);

      return head.size();
    }


    /**
     * 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) == getNumberOfMatches(*this);
    }


    /**
     * 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);
      RETURN_IF_DIFFERENT(this->target,      header.target);

      // compare compatibility

      if (is_valid(this->livetime) == is_valid(header.livetime)  &&
	  is_valid(this->DAQ)      == is_valid(header.DAQ)       &&
	  is_valid(this->K40)      == is_valid(header.K40))      { 
	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)) {

	this->createUUID();
	this->UUID = JUUID::rndm();
	
	genvol  .add(header.genvol);
	norma   .add(header.norma);
	livetime.add(header.livetime);
	DAQ     .add(header.DAQ);
	K40     .add(header.K40);

      } 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
    JUUID                         UUID;             // header unique identifier
    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::calibration           calibration;
    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::K40                   K40;
    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 size_t getMaximumNumberOfMatches()
    {
      static JHead header;

      header.push();
      
      static const size_t 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,5);

  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 (first .is_valid(pd) &&
	      second.is_valid(pd) &&
	      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.
   *
   * Note that this operator uses the JHead::less method.
   *
   * \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