#ifndef DET_HH_INCLUDED
#define DET_HH_INCLUDED

#include "Evt.hh"
#include "AAObject.hh"

#include "constants.hh"

class EventFile;

#include <map>
#include <string>
#include "foreach.hh"
#include "stringutil.hh"
//using namespace stringutil;

using stringutil::str;
using stringutil::Filedesc;
using stringutil::form;
using stringutil::file_exists;

using namespace std;

/*! a PMT */

struct Pmt : public AAObject
{
  int id;                    //dox: global PMT identifyer
  unsigned int channel_id;   //dox: local PMT identifyer 0-31
  int dom_id;                //dox: ID of the dom this PMT sits in
  Vec pos;                   //dox: position (global)
  Vec dir;                   //dox: direction
  double t0;                 //dox: PMT t0 calibration
  double rate;               //dox: the singles rate on this pmt
  int status;                //dox: 0 means PMT OK for analysis; 1 = PMT NOT OK 

  int flag;     //dox: can use this for flagging purposes
  Pmt() : id(0), channel_id(0), dom_id(0), t0(0), rate(0), status(0), flag(0) {}

  /*! set the position of the hit and compute the hit-time from the TDC and the T0
      (time slewing correction not yet done) */

  void dress( Hit& h , bool set_t = true , bool use_t = false, bool v = false ) const
  {
    h.pos = pos;
    h.dir = dir;
    h.channel_id = channel_id; // set all ids since we
    h.dom_id     = dom_id;     // dont' know how dress was called.
    h.pmt_id     = id;

    if ( set_t )
    {
      double t = use_t ? h.t : h.tdc; 

      h.t =  t + t0;
      if ( t == 0.0 && v ) stringutil::print("warning: hit with time == 0 ", h );
    }

  }

  void print( std::ostream& out = std::cout ) const
  {
    out << "Pmt: id=" + str(id) + ", dom_id=" + str(dom_id) + ", channel_id=" + str(channel_id) + ", pos=";
    pos.print( out );
  }

  // string __str__() const
  // {
  //   ostringstream s;
  //   print(s);
  //   return s.str();
  // }


  ClassDef(Pmt, 4);
};

struct Dom : public AAObject
{
  int id;                //dox: global Dom identifyer
  Vec pos;               //dox: position

  int line_id;           //dox: line (a.k.a. string) id
  int floor_id;          //dox: floor id
  vector<Pmt> pmts;      //dox: The index in the vector is the channel-id.
  int flag;
  int status;            //dox: see https://git.km3net.de/common/km3net-dataformat/-/blob/master/definitions/module_status.csv 

  Dom(): id(0), line_id(0), floor_id(0), flag(0), status(0) {}

  // string __str__() const
  // {
  //   ostringstream s;
  //   print(s);
  //   return s.str();
  // }

  void print( std::ostream& out = std::cout ) const
  {
    out << "Dom: " + str(id) + ", line=" + str(line_id) + ", floor=" + str(floor_id) + " pos=";
    pos.print( out );
  }

  ClassDef(Dom, 5);
};


struct Det;
bool read_detfile_detx( string filename, Det& det );
bool read_detfile_det( string filename, Det& det );

/*! Class representing the Detector geometry and configuration data.
 */

struct Det : public TObject
{
  int id;
  int version;
  TTimeStamp validity_range_start;
  TTimeStamp validity_range_stop;

  /*! geographical location of the detector - see KM3NeT_SOFT_2016_003 */
  string utm_reference_elipsoid;
  int utm_zone;
  double utm_ref_easting;
  double utm_ref_northing;
  double utm_ref_z;

  double longitude, latitude;
  double meridian_convergence_angle;


  // the doms map owns the doms, which own the pmts
  std::map<int, Dom>  doms; //dox: DOMs by dom-id.

  // the pmts map only serves to provide fast lookup
  std::map<int, Pmt*> pmts; //dox: PMT pointers by global id

  void init_pmts() // call this when the doms map will not change anymore
  {
    foreach_map(i, d, doms)
    {
      foreach( p, d.pmts ) pmts[p.id] = &p;
      (void)i; // unused varible warning
    }
  }



  Dom& get_dom( Pmt& pmt) { return doms[ pmt.dom_id ]; }

  int get_dom_id (Pmt& pmt)    { return  get_dom(pmt).id;      } // for symmetry
  int get_line_id( Pmt& pmt )  { return  get_dom(pmt).line_id; }
  int get_floor_id( Pmt& pmt ) { return  get_dom(pmt).floor_id;}
  int get_status( Pmt& pmt )   { return  get_dom(pmt).status;}

  Pmt& get_pmt( int dom_id, int channel_id )
  {
    return doms[dom_id].pmts[channel_id];
  }

  Det(): id(0) {}
  
  Det (string filename )
  {
    read(filename);
  }

  Det(EventFile& evtfile);

  void remove_dom_if( bool (*pred)(const Dom&) )
  {
    remove_element_if( doms, pred );
  }


  bool read( string filename )
  {
    *this = Det(); // clean slate

    if (!file_exists( filename) )
    {
      stringutil::print ("Error : detector file", filename, " does not exist!" );
      return false;
    }

    string ext = Filedesc(filename).ext;
    if ( ext == "det"  ) return read_detfile_det( filename , *this );

    if ( ext == "detx" )
    {
      bool r = read_detfile_detx( filename , *this );
      set_lonlat_from_utm();
      return r;
    }

    return false;
  }

  /*! fills a map, lines_floor_map()[line-nuber][floornumber] gives a pointer to the dom */

  map<int, map< int, Dom* > > lines_floor_map() 
  {
    map<int, map< int, Dom* > > r;
    foreach_map( dom_id, dom, doms ) r[dom.line_id][dom.floor_id] = &dom;
    return r;
  }


  // Table domtable()
  // {
  //   Table T;
  //   auto m = lines_floor_map()

  //   auto lines  = keys( m );
  //   auto floors = keys( m.begin().second );

  //   T << name;
  //   for( int line : keys( m ) ) T << line;

  //   foreach_map_map( line, floor, pdom, m ) {


  //   }



  // }

  
  void print( std::ostream& out = std::cout ) const
  {
    out << "Det with " << pmts.size() << " pmts in " << doms.size() << "doms at lon/lat " << longitude * aa::deg << " " << latitude * aa::deg << "(deg)";

  }

  void prnt()
  {
    stringutil::print("Detector_______________________________");

    stringutil::print ("size of containers:");
    stringutil::print ("- pmts : ", pmts.size() );
    stringutil::print ("- doms : ", doms.size() );

    stringutil::print("--------first doms--------------");

    int i = 0;
    foreach_map(id, dom, doms)
    {
      stringutil::print ("id, npmts", id, dom.pmts.size() );
      foreach( p, dom.pmts ) cout << p.id << " ";
      cout << endl;

      if (i++ > 10) break;
    }
  }


  /*! The center of gravity of the doms */

  Vec get_cg()
  {
    Vec r;
    foreach_map( id, dom, doms )
    {
      r += dom.pos;
      (void) id; // usused variable warning
    }
    return r / doms.size();
  }

  /*! Shift every dom by v */

  void shift(const Vec& v)
  {
    foreach_map( id, dom, doms )
    {
      (void)id;
      dom.pos += v;
      foreach( pmt, dom.pmts ) pmt.pos += v;
    }
  }


  /*! return a vector with the maximal postional, rotational, and time difference
      between the position,directionsa and times of the hits and those in the detector */

  vector<double> compare_hits ( const vector<Hit>& hits )
  {
    double posdif(0), dirdif(0), tdif(0);
   
    vector<Hit> hits2 = hits;
    apply( hits2 );

    for ( unsigned i = 0 ; i < hits.size() ; i++ )
    {
      const Hit& h1 = hits[i];
      const Hit& h2 = hits2[i];

      posdif = max( posdif, (h1.pos - h2.pos ).len() );
      dirdif = max( dirdif, angle_between( h1.dir, h2.dir ) );
      tdif   = max( tdif  , abs( h1.t - h2.t ));
    }

    return {posdif, dirdif, tdif };
  }


  /*! check if the set of hits agrees with the current detector for what concerns
      position and time. For the hit-times, we expect some difference due to time-slewing */

  bool check( const vector<Hit>& hits , bool set_t = true, 
              double maxdist = 1e-5, double maxang = 1e-5, double maxt = 20  )
  {
    auto diff = compare_hits( hits );

    stringutil::print( "detector check postion-diff, direction-diff, time-diff :", diff[0], diff[1] , diff[2] );

    if ( diff[0] > maxdist || diff[1] > maxang || diff[2] > maxt ) return false;
    return true;
  }

  /*! Apply the detector to the hits. If set_t is true, the hit-times
      will be corrected for the offset if the PMTs */

  void apply( vector<Hit>& hits , bool set_t = true )
  {
    if ( hits.size() == 0 ) return;

    int have_glob_pmt_ids = 0;
    int have_dom_chan_ids = 0;
    int have_tdc_values = 0;
    int have_t_values = 0;
    bool use_t = false;

    foreach ( h, hits )
    {
      have_glob_pmt_ids += h.pmt_id;
      have_dom_chan_ids += (int)h.dom_id;
      have_tdc_values += h.tdc != 0.0;
      have_t_values += h.t != 0.0;
    }

    if ( have_dom_chan_ids != 0 && have_glob_pmt_ids != 0 )
    {
      static bool once = true;

      if ( once )
      {
        stringutil::print("Both global and local pmts ids are set");
        once = false;
      }
    }

    if ( !(have_glob_pmt_ids || have_dom_chan_ids) )
    {
      stringutil::print("cannot figure out how to apply Det object to hits");
      stringutil::print("because neither global or local pmts ids are set");
      for ( auto& h : hits ) stringutil::print(h);
      return;
    }

    if ( set_t && have_t_values && !have_tdc_values ) 
    {
      static bool once = true;
      if (once) stringutil::print("will use legacy Hit::t as hit time as uncalibrated hit time, rather than Hit::tdc");
      use_t = true;
      once = false;
    }


    if ( set_t && have_tdc_values && !have_t_values )
    {  
      static bool once = true;
      if (once) stringutil::print("using Hit::tdc as uncalibrated hit time.");
      once = false;
    }


    foreach (h, hits )
    {
      Pmt* tube = nullptr;

      if ( have_glob_pmt_ids )
      {
        tube = pmts[ h.pmt_id ];

        if ( have_dom_chan_ids ) // if we have both, check
        {
          if ( tube -> dom_id     != h.dom_id ||
               tube -> channel_id != h.channel_id )
          {
            stringutil::print ("Error: incompatible dom/channel/pmt id combination incompatible with detector.");
            h.print(cout);
          }
        }
      }

      if ( have_dom_chan_ids )
      {
        tube = & ( doms[h.dom_id].pmts[ (int) h.channel_id ] );
      }

      if (tube) tube -> dress( h , set_t , use_t );
    }
  }

  void apply( Evt& evt )
  {
    apply( evt.hits    , true  );
    apply( evt.mc_hits , false );
  }

  /*! set all rates of all pmts to zero (or some other nubmer) */
  void zero_rates(double value = 0 )
  {
    for (auto& d : doms ) {
      for ( auto& pmt : d.second.pmts )
      {
        pmt.rate = value;
      }
    }
  }


  //=============== geographical coordinates functions ==========================

  /*! compute the longitude and latitude from the UTM data */
  static vector<double> compute_lonlat_from_utm( double utmX, double utmY, int zone );

  /*! Return latitude of the central meridian in the UTM zone */
  static double longitude_of_central_meridian(int utmzone);

  /*! Compute UTM zone from the latitude */
  static int get_utm_zone(double lat);

  /*! Compute the convergence angle. */
  static double compute_meridian_convergence_angle( double longitude, double latitude);


  /*! Compute and set the MCA */
  void set_meridian_convergence_angle();

  /*! set longitude and lattitude, and utm zone and convergence-angle */
  void set_lonlat( double lon, double lat)
  {
    longitude = lon;
    latitude  = lat;
    utm_zone  = get_utm_zone( lon );
    meridian_convergence_angle = compute_meridian_convergence_angle( lon, lat );
  }

  /*! use own UTM data to set longitude and lattitude, and utm zone and convergence-angle */
  void set_lonlat_from_utm()
  {
    auto c = compute_lonlat_from_utm( utm_ref_easting, utm_ref_northing, utm_zone );
    int oldzone = utm_zone;

    cout << " zone " << utm_zone << endl;

    set_lonlat( c[0], c[1] );
    if (utm_zone != oldzone ) dbg ("internal error computing utm coordinates");
  }

  ClassDef( Det, 7 );
};


/*! Return Summary information on the geometry */

struct DetDim
{
  Vec origin;
  Vec coord_min;
  Vec coord_max;
  Det& detector;
  uint dom_count;

  DetDim( Det& det ): detector( det ) {

    dom_count = 0;
    coord_max = det.doms[1].pos;
    coord_min = det.doms[1].pos;
    for ( auto& dom_pair : det.doms ) {

      Dom dom = dom_pair.second;

      if ( dom.pos.z < 50 ) continue; // Get rid of annoying bottom DOM.
      coord_max.x = coord_max.x > dom.pos.x ? coord_max.x : dom.pos.x;
      coord_min.x = coord_min.x < dom.pos.x ? coord_min.x : dom.pos.x;

      coord_max.y = coord_max.y > dom.pos.y ? coord_max.y : dom.pos.y;
      coord_min.y = coord_min.y < dom.pos.y ? coord_min.y : dom.pos.y;

      coord_max.z = coord_max.z > dom.pos.z ? coord_max.z : dom.pos.z;
      coord_min.z = coord_min.z < dom.pos.z ? coord_min.z : dom.pos.z;

      dom_count++;

    }

    origin = Vec( center_x(), center_y(), coord_min.z );
    cout << "Origin set at " << origin.x << "x, " << origin.y << "y, " << origin.z << "z" << endl;

  }

  uint get_dom_count() const {
    return dom_count;
  }



  // Vec get_dom_pos( uint dom_id ) const {
  //   try {
  //     if ( dom_id <= 20 ) {
  //       cout << "dom id: " << dom_id << ", pos: " << detector.doms[dom_id].pos.x << ", " << detector.doms[dom_id].pos.y << ", " << detector.doms[dom_id].pos.z << endl;
  //     }
  //     return detector.doms[dom_id].pos;
  //   } catch ( const std::exception& e ) {
  //     return Vec(-1, -1, -1);
  //   }
  // }

  double height() const {
    return coord_max.z - coord_min.z;
  }

  double radius() const {
    return ( coord_max.x - coord_min.x + coord_max.y - coord_min.y ) / 4.0;
  }

  double center_x() const {
    return ( coord_max.x - coord_min.x ) / 2 + coord_min.x;
  }

  double center_y() const {
    return ( coord_max.y - coord_min.y ) / 2 + coord_min.y;
  }

  double center_z() const {
    return ( coord_max.z - coord_min.z ) / 2 + coord_min.z;
  }

  Vec get_origin() const {
    return origin;
  }
};

inline double get_dist_to_edge( const Trk trk, const DetDim detdim) {
  /*
    Get the smallest distance to the surface of the detector. Negative if inside, positive if outside.
  */

  Vec origin = detdim.get_origin();

  const double r = detdim.radius();

  const double D_wall = ( trk.pos - Vec( origin.x, origin.y, trk.pos.z ) ).len() - r;
  const double D_bot  = detdim.coord_min.z - trk.pos.z;
  const double D_top  = trk.pos.z - detdim.coord_max.z;

  if ( D_top > 0 ) {
    if ( D_wall > 0 ) { //above outside
      return ( sqrt( D_wall * D_wall + D_top * D_top ) );
    } else { //straight above
      return ( D_top );
    }
  } else if ( D_bot > 0 ) {
    if ( D_wall > 0 ) { //below outside
      return ( sqrt( D_wall * D_wall + D_bot * D_bot ) );
    } else { //straight below
      return ( D_bot );
    }
  } else {
    if ( D_wall > 0 || ( D_wall > D_bot && D_wall > D_top ) ) {
      return ( D_wall ); //nearest wall
    } else if ( D_top > D_bot ) {
      return ( D_top ); //nearest top inside
    } else {
      return ( D_bot ); //nearest bot inside
    }
  };
};


#include "JDAQSummaryslice.hh"

void read(Det& det, const KM3NETDAQ::JDAQSummaryslice& ss );


#endif