// HTrack.cpp

#include "HTrack.h"
#include "trfutil/trfstream.h"
#include "trfbase/PropStat.h"
#include "trfbase/Surface.h"

#include <set>

using std::multiset;

using trf::Surface;
using trf::ETrack;
using trf::PropStat;
using trf::Propagator;
using trf::McIdList;
using trf::TrackMcHitInfo;
using trf::HTrack;

//**********************************************************************
// member functions
//**********************************************************************

// Default constructor.

HTrack::HTrack()
: RefCounter(), _tre(), _chisq(0.0), _nfit(-1) { }

//**********************************************************************

// constructor

HTrack::HTrack(const ETrack& tre)
: RefCounter(), _tre(tre), _chisq(0.0), _nfit(0) { }

//**********************************************************************

// destructor

HTrack::~HTrack() { }

//**********************************************************************

// add a hit to the back of the list

void HTrack::add_hit(const HitPtr& phit) {
  _hits.push_back(phit);
}

//**********************************************************************

// Drop the last hit on the list.

void HTrack::drop_hit() {
  if ( _hits.size() == 0 ) {
    return;
  }
  _hits.pop_back();
  if ( _nfit > _hits.size() ) unset_fit();
}

//**********************************************************************

// Drop all the hits.

void HTrack::drop_hits() {
  _hits.erase( _hits.begin(), _hits.end() );
  if ( _nfit > _hits.size() ) unset_fit();
}

//**********************************************************************

// Set the fit.
// The track and chi-square are assigned the specified value and
// the number of fitted hits is set to the current number of hits.

void HTrack::set_fit(const ETrack& tre, double chisq) {
  _tre = tre;
  _chisq = chisq;
  _nfit = _hits.size();
}

//**********************************************************************

// Unset the fit.
// The number of hits is set to 0 and the chi-square is set to zero.
// The track is not modified (the last fit becomes the starting track).

void HTrack::unset_fit() {
  _nfit = -1;
  _chisq = 0.0;
}


// Return the number of measurements in the fit

int HTrack::get_num_meas() const {
  int nmeas = 0;
  HitList::const_iterator it;
  for (it = _hits.begin(); it!=_hits.end();++it){
    nmeas += (*it)->size();
  }
  return nmeas;
}

//Return the id of the "best" MC track for this track
TrackMcHitInfo HTrack::get_mcinfo() const {

  McIdList idlist;
  typedef multiset<McId> mset;
  mset mcidset;
  //Loop over the hits and accumulate the full list of McIds in the set
  HitList::const_iterator it;
  int nhits = _hits.size();
  int nshared=0;

  for (it = _hits.begin(); it!=_hits.end();++it){
    mcidset.insert((*it)->get_mc_ids().begin(),(*it)->get_mc_ids().end());
    McIdList::const_iterator mcl;
    int first=0;
    bool shared=false;
    for (mcl=(*it)->get_mc_ids().begin(); mcl!=(*it)->get_mc_ids().end(); ++mcl)
    {
      int id=(*mcl);
      if( first==0) first=abs(id);
      if( abs(id) != first ) shared=true;
    }
    if(shared) nshared+=1;
  }

  TrackMcHitInfo::McIdMap idmap;
  // Loop over the set and create the id map
  for(mset::const_iterator setit=mcidset.begin(); setit != mcidset.end(); ++setit){
    idmap[*setit] = mcidset.count(*setit);
  }

  // Return a TrackMcHitInfo object
  return TrackMcHitInfo(nhits,nshared,idmap);
}

//**********************************************************************

// Return the ID's of the intersection of the MC tracks contributing 
// to the clusters on this track.
McIdList HTrack::get_mc_ids() const {
  McIdList result;

  // Loop over hits.

  bool first = true;
  for(HitList::const_iterator ihit = _hits.begin(); ihit != _hits.end();
      ++ihit) {

    // On the first cluster, store all mcids in the result vector.

    if(first) {
      first = false;
      const McIdList& mcids1 = (*ihit)->get_mc_ids();
      for(McIdList::const_iterator i = mcids1.begin(); i != mcids1.end(); 
	  ++i)
	result.push_back(*i);
    }
    else {

      // On subsequent clusters, find the intersection between the previous
      // result and the current clsuter.

      McIdList new_result;
      const McIdList& mcids = (*ihit)->get_mc_ids();
      for(McIdList::const_iterator i = result.begin(); i != result.end();
	  ++i) {
	for(McIdList::const_iterator j = mcids.begin(); j != mcids.end();
	    ++j) {
	  if(*i == *j) {
	    new_result.push_back(*i);
	    break;
	  }
	}
      }
      result = new_result;
    }
  }
  return result;
}

//**********************************************************************

// Unset the fit and reset the starting track.
// The fit is unset and the track is set to the specified value.

void HTrack::unset_fit(const ETrack& tre) {
  unset_fit();
  _tre = tre;
}

//**********************************************************************

// output

void HTrack::ostr(ostream& stream) const {
  stream << begin_object;
  stream << "HTrack";
  stream << " is ";
  if ( ! is_fit() ) stream << "not ";
  stream << "fit:" << new_line;
  stream << _tre << new_line;
  stream << "Chi-square: " << _chisq << new_line;
  stream << "Track has " << _hits.size() << " hit";
  if ( _hits.size() != 1 ) stream << 's';
  if ( _hits.size() > 0 ) {
    stream << ':';
    for ( HitList::const_iterator ihit=_hits.begin();
          ihit!=_hits.end(); ++ihit ) {
      stream << new_line;
      const HitPtr& phit = *ihit;
      if ( phit ) stream << **ihit;
      else stream << "Invalid hit";
    }
  } else {
    stream << '.';
  }
  stream << end_object;
}

//**********************************************************************

// Propagate the fit track.
// Note this does not change the fit status because the propagation
// should be reversible.

PropStat HTrack::propagate(const Propagator& prop, const Surface& srf) {
  return prop.err_prop(_tre,srf);
}

//**********************************************************************

// Propagate the fit track in a specified direction.
// Note this does not change the fit status because the propagation
// should be reversible.

PropStat HTrack::propagate(const Propagator& prop, const Surface& srf,
                           Propagator::PropDir dir) {
  return prop.err_dir_prop(_tre,srf,dir);
}

//**********************************************************************
// free functions
//**********************************************************************

// output

ostream& trf::operator<<( ostream& stream, const HTrack& trh ) {
  trh.ostr( stream );
  return stream;
}

//**********************************************************************

// Equality.

bool operator==(const HTrack& lhs, const HTrack& rhs) {
  if ( lhs.get_hits()       != rhs.get_hits()       ) return false;
  if ( lhs.get_track()      != rhs.get_track()      ) return false;
  if ( lhs.get_chi_square() != rhs.get_chi_square() ) return false;
  if ( lhs.is_fit()         != rhs.is_fit()         ) return false;
  return true;
}

//**********************************************************************

// Inequality.

bool operator!=(const HTrack& lhs, const HTrack& rhs) {
  return ! ( lhs == rhs );
}

//**********************************************************************