// HitTestFit1.h

#ifndef HitTestFit1_H
#define HitTestFit1_H

// Hit class to test fitter.
// Prediction is the first track parameter.
// Each hit returns one prediction.
// Measurement and error are specified in hit constructor.

#include "trfbase/Hit.h"
#include <math.h>
#include <iostream.h>
#include "objstream/ObjData.hpp"
#include "trfutil/trfstream.h"
#include "trfbase/SurfTest.h"
#include "trfbase/ETrack.h"

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

// Hit.

ObjPtr create_HitTestFit1(const ObjData& data);

class HitTestFit1 : public trf::Hit {

public:  // static methods

  // Return the type name.
  static TypeName get_type_name() { return "HitTestFit1"; }

  // Return the creator.
  static ObjCreator get_creator() { return create_HitTestFit1; }

  // Return the type.
  static Type get_static_type() { return get_creator(); }

private:  // data
  double _msmt;               // measurement
  double _emsmt;              // square of measured error
  double _pred;               // prediction
  double _epred;              // square of prediction error

private:  // methods

  // Output stream.
  void ostr(ostream& stream) const {
    stream << begin_object;
    stream << "TestFit1 hit:" << new_line
           << "measure = " << measured_vector()(0) << " +/- "
           << sqrt(_emsmt) << new_line
           << "predict = " << predicted_vector()(0) << " +/- "
           << sqrt(_epred);
    stream << end_object;
  };

  // Equality.
  bool equal(const trf::Hit& hit) const {
     assert( hit.get_type() == get_type() );
     return get_cluster() == hit.get_cluster() && 
            _pred == ((const HitTestFit1&) hit)._pred;
   };

public:

  // Constructor.
  HitTestFit1(double msmt, double emsmt, const trf::ETrack& tre)
  : _msmt(msmt), _emsmt(emsmt) {
    update(tre);
  };

  // Return the type.
  Type get_type() const { return get_static_type(); }
  
  // Write the object data.
  ObjData write_data() const {
    ObjData data("TestFit1");
    data.add_double("msmt",_msmt);
    data.add_double("emsmt",_emsmt);
    data.add_double("pred",_pred);
    data.add_double("epred",_epred);
    data.add_bare_pointer("cluster",get_cluster());
    return data;
  }

  // Return hit chractersistics.
  int size() const { return 1; };
  trf::HitVector measured_vector() const { return trf::HitVector(_msmt); };
  trf::HitError measured_error() const { return trf::HitError(_emsmt); };
  trf::HitVector predicted_vector() const { return trf::HitVector(_pred); };
  trf::HitError predicted_error() const { return trf::HitError(_epred); };
  trf::HitDerivative dhit_dtrack() const {
    trf::HitDerivative hder(1);
    hder(0,0) = 1.0;
    return hder;
  }
  trf::HitVector difference_vector() const
    { return predicted_vector() - measured_vector(); };

  // Update.
  void update(const trf::ETrack& tre) {
    _pred = tre.get_vector()(0);
    _epred = tre.get_error()(0,0);
  };

};

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

// Cluster.

ObjPtr create_ClusTestFit1(const ObjData& data);

class ClusTestFit1 : public trf::Cluster {

public:  // static methods

  // Return the type name.
  static TypeName get_type_name() { return "ClusTestFit1"; }

  // Return the creator.
  static ObjCreator get_creator() { return create_ClusTestFit1; }

  // Return the type.
  static Type get_static_type() { return get_creator(); }

private:
  const SurfTest& _stst;     // surface
  double _msmt;               // measurement
  double _emsmt;              // square of measured error
  void ostr(ostream& stream) const {
    stream << "ClusTestFit1 with measurement " << _msmt << " +/- "
           << sqrt(_emsmt);
  };

  // Equality.
  bool equal(const trf::Cluster& rhs) const {
    assert( rhs.get_type() == get_type() );
    const ClusTestFit1& clus = (const ClusTestFit1&) rhs;
    return _msmt == clus._msmt && _emsmt == clus._emsmt;
  };

  // generate hit prediction ==> pass msmt and track
  // In this simple example, there is only one hit/cluster.
  trf::HitList _predict(const trf::ETrack& tre) const  {
    trf::HitList hits;
    hits.push_back( trf::HitPtr(new HitTestFit1(_msmt,_emsmt,tre)) );
    return hits;
  };

public:

  // constructor from msmt and its error
  ClusTestFit1(const SurfTest& srf,double msmt,double emsmt)
  : _stst(srf), _msmt(msmt), _emsmt(emsmt) { };

  // Return the type.
  Type get_type() const { return get_static_type(); }

  // Write the object data.
  ObjData write_data() const {
    ObjData data("TestFit1");
    data.add_bare_pointer("surface",&_stst);
    data.add_double("msmt",_msmt);
    data.add_double("emsmt",_emsmt);
    return data;
  }

  // Return the surface.
  const SurfTest& get_surface() const { return _stst; };

};

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

// Return a weighted average.
class WAvg {
private:
  double _sum0;
  double _sum1;
  double _sum2;
  int _count;
public:
  WAvg() : _sum0(0.0), _sum1(0.0), _sum2(0.0), _count(0) { };
  // input is value and square of the uncertainty
  void add_pair(double val, double err) {
    double weight = 1.0/err;
    _sum0 += weight;
    _sum1 += weight*val;
    _sum2 += weight*val*val;
    ++_count;
  };
  int get_count() const { return _count; };
  double get_average() const { return _sum1/_sum0; };
  double get_error() const { return 1.0/_sum0; };
  double get_chi_square() const { return _sum2 - _sum1*_sum1/_sum0; };
};

inline ostream& operator<<(ostream& stream, const WAvg& avg) {
  stream << begin_object;
  stream << "Weighted average = " << avg.get_average() 
         << " +/- " << sqrt(avg.get_error())
         << " (err = " << avg.get_error() << ")"
         << new_line
         << "Chi-square = " << avg.get_chi_square();
  stream << end_object;
  return stream;
}

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

// Creators.

ObjPtr create_HitTestFit1(const ObjData& data) {
  assert( data.get_object_type() == "HitTestFit1" );
  if ( data.get_object_type() != "HitTestFit1" ) return ObjPtr(0);
  double msmt = data.get_double("msmt");
  double emsmt = data.get_double("emsmt");
  double pred = data.get_double("pred");
  double epred = data.get_double("epred");
  trf::ClusterPtr pclu;
  data.get_bare_pointer("cluster",pclu);
  trf::SurfacePtr psrf( pclu->get_surface().new_surface() );
  trf::TrackVector vec;
  vec(0) = pred;
  trf::TrackError err;
  err(0,0) = epred;
  trf::ETrack tre(psrf,vec,err);
  trf::Hit* phit( new HitTestFit1(msmt,emsmt,tre) );
  phit->set_parent_pointer(pclu);
  return ObjPtr(phit);
}

ObjPtr create_ClusTestFit1(const ObjData& data) {
  assert( data.get_object_type() == "ClusTestFit1" );
  if ( data.get_object_type() != "ClusTestFit1" ) return ObjPtr(0);
  double msmt = data.get_double("msmt");
  double emsmt = data.get_double("emsmt");
  const SurfTest* psrf;
  data.get_bare_pointer("surface",psrf);
  return ObjPtr( new ClusTestFit1(*psrf,msmt,emsmt) );
}

#endif