// HitCylPhiGenerator_t.cpp

// Test HitCylPhiGenerator.

#include "HitCylPhiGenerator.h"
#include <iostream>
#include <string>
#include <cassert>
#include "objstream/ObjData.hpp"
#include "trfutil/trfstream.h"
#include "SurfCylinder.h"
#include "trfbase/ETrack.h"
#include "trfbase/Hit.h"
#include "trfutil/TRFMath.h"
#include "trfbase/CrossStat.h"

using std::cout;
using std::cerr;
using std::endl;
using std::string;
#ifndef DEFECT_CMATH_NOT_STD
using std::sqrt;
using std::fabs;
#endif

using namespace trf;

ObjPtr create(const ObjData& data);

// SurfCylinder subclass for testing.
// z > 100 is considered out of bounds.
class TestSurfCylinder : public SurfCylinder {
public:  // static methods
  // Return the type name.
  static TypeName get_type_name() { return "TestSurfCylinder"; }
  // Return the creator.
  static ObjCreator get_creator() { return create; };
  // Return the type.
  static Type get_static_type() { return get_creator(); }
  double _zmax;
public:
  TestSurfCylinder(double radius) : SurfCylinder(radius) { };
  CrossStat status(const ETrack& tre) const {
    return status( (const VTrack&) tre );
  };
  // Return the type.
  Type get_type() const { return get_static_type(); }
  // Write the object data.
  ObjData write_data() const { return ObjData(get_type_name()); }
  CrossStat status(const VTrack& trv) const {
    CrossStat xstat = pure_status(trv);
    if ( ! xstat.at() ) return xstat;
    double z = trv.get_vector(SurfCylinder::IZ);
    if ( z > 100.0 ) return CrossStat(CrossStat::OUT_OF_BOUNDS);
    return CrossStat(CrossStat::IN_BOUNDS);
  };
  Surface* new_surface() const { 
    return new TestSurfCylinder(get_radius());
  };
};

ObjPtr create(const ObjData& data) {
  return ObjPtr(0);
}

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

int main( ) {

  string component = "HitCylPhiGenerator";
  string ok_prefix = component + " (I): ";
  string error_prefix = component + " test (E): ";

  cout << ok_prefix
       << "---------- Testing component " + component
          + ". ----------" << endl;

  // Make sure assert is enabled.
  bool assert_flag = false;
  assert ( ( assert_flag = true, assert_flag ) );
  if ( ! assert_flag ) {
    cerr << "Assert is disabled" << endl;
    return 1;
  }

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

  cout << trf_format;
  cout << ok_prefix << "Test constructor." << endl;
  double radius = 25.0;
  double dphi = 0.01;
  TestSurfCylinder tcy(radius);
  HitCylPhiGenerator gen(tcy,dphi);
  cout << gen << endl;

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

  cout << ok_prefix << "Generate a list of hits." << endl;
  double phi0 = 1.2345;
  SurfacePtr psrf( new SurfCylinder(radius) );
  TrackVector vec;
  vec(SurfCylinder::IPHI) = phi0;
  TrackError err;
  VTrack trv(psrf,vec);
  ETrack tre(trv,err);
  ClusterList clusters;
  int nclus = 20;
  for ( int i=0; i<nclus; ++i ) {
    Cluster* pclu = gen.new_cluster(trv);
    assert( pclu != 0 );
    clusters.push_back( pclu );
  }
  // Verify and display the list.
  HitList hits;
  ClusterList::iterator iclu;
  for ( iclu=clusters.begin(); iclu!=clusters.end(); ++iclu ) {
    ClusterPtr pclu = *iclu;
    assert( pclu != 0 );
    HitList newhits = pclu->predict(tre);
    assert( newhits.size() == 1 );
    hits.push_back( newhits.front() );
  }
  assert( hits.size() == nclus );
  
  //********************************************************************

  cout << ok_prefix << "Generate list of measured values." << endl;
  double avg = 0.0;
  double sdev = 0.0;
  HitList::iterator ihit;
  for ( ihit=hits.begin(); ihit!=hits.end(); ++ihit ) {
    HitPtr phit = *ihit;
    assert( phit.pointer() != 0 );
    double mval = phit->measured_vector()(0);
    double merr = phit->measured_error()(0,0);
    double pval = phit->predicted_vector()(0);
    double perr = phit->predicted_error()(0,0);
    assert( is_equal(pval,phi0) );
    assert( is_equal(merr,dphi*dphi) );
    cout << mval << endl;
    avg += mval;
    sdev += (mval-pval)*(mval-pval);
  }
  avg /= nclus;
  sdev = sqrt(sdev/nclus);
  cout << " Avg: " << avg << endl;
  cout << "Sdev: " << sdev << endl;
  assert( fabs( avg - phi0 ) < dphi );
  assert( fabs( sdev - dphi ) < dphi );

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

  cout << ok_prefix << "Generate a hit out-of-bounds." << endl;
  vec(SurfCylinder::IZ) = 200.0;
  trv.set_vector(vec);
  assert( ! gen.new_cluster(trv) );

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

  cout << ok_prefix << "Print final generator state." << endl;
  cout << gen << endl;

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

  cout << ok_prefix
       << "------------- All tests passed. -------------" << endl;
  return 0;

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

}