// HitCylPhiZGenerator_t.cpp

// Test HitCylPhiZGenerator.

#include "HitCylPhiZGenerator.h"
#include <iostream>
#include <string>
#include <cassert>
#include "trfutil/trfstream.h"
#include "SurfCylinder.h"
#include "trfbase/ETrack.h"
#include "trfbase/Hit.h"
#include "trfutil/TRFMath.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;

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

int main( ) {

  string component = "HitCylPhiZGenerator";
  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 stereo = 0.1;
  double dphiz = 0.01;
  HitCylPhiZGenerator gen(radius,stereo,dphiz);
  cout << gen << endl;

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

  cout << ok_prefix << "Generate a list of hits." << endl;
  double phi0 = 1.2345;
  SurfacePtr psrf( new SurfCylinder(radius) );
  TrackVector vec;
  vec(0) = phi0;
  TrackError err;
  VTrack trv(psrf,vec);
  ETrack tre(trv,err);
  ClusterList clusters;
  int nclus = 20;
  for ( int i=0; i<nclus; ++i ) 
    clusters.push_back( gen.new_cluster(trv) );
  // 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,dphiz*dphiz) );
    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 ) < dphiz );
  assert( fabs( sdev - dphiz ) < dphiz );

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

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

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

}