// CheckCylDca.cpp

#include "CheckCylDca.h"
#include <iostream>
#include "objstream/ObjData.hpp"
#include "objstream/ObjTable.hpp"
#include "trfutil/trfstream.h"
#include "trffit/MTrack.h"
#include "trfcyl/SurfCylinder.h"
#include "spacegeom/SpacePoint.h"
#include "trfcut/CutRecorder.hpp"

using std::cerr;
using std::endl;
using std::string;

using trf::TrackVector;
using trf::ETrack;
using trf::CutRecord;
using trf::CutRecorder;
using trf::CheckCylDca;

//**********************************************************************
// Free functions.
//**********************************************************************

namespace {

// Creator.

ObjPtr create(const ObjData& data) {
  assert( data.get_object_type() == "CheckCylDca" );
  double max_dca = data.get_double("max_dca");
  ObjDoublePtr bfield;
  ObjData::Type dtype = data.get_type("bfield");
  if(dtype == ObjData::DOUBLE)
    bfield = new ObjDouble(data.get_double("bfield"));
  else if(dtype == ObjData::SHARE_PTR)
    bfield.assign_with_dynamic_cast(data.get_share_pointer("bfield"));
  else {
    cerr << "PropCylDCA: can't get bfield parameter." << endl;
    abort();
  }
  string cut_name;
  if(data.has("cut_name"))
    cut_name = data.get_string("cut_name");
  return ObjPtr( new CheckCylDca(max_dca,bfield,cut_name) );
}

}

//**********************************************************************
// Member functions.
//**********************************************************************

// Output stream.

void CheckCylDca::ostr(ostream& stream) const {
  stream << begin_object;
  stream << "Check distance of closest approach in r-phi <= " << _max_dca;
  stream << "in magnetic field of strength " << (*_bfield)() << endl;
  stream << end_object;
}

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

// Constructor.

CheckCylDca::CheckCylDca(double max_dca, ObjDoublePtr bfield, 
			 const string& cut_name)
  : _max_dca(max_dca), _bfield(bfield) {
  CutRecorder* p = CutRecorder::instance();
  assert(p != 0);
  _cut_id = p->get_cutid(cut_name);
}

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

// Return the creator.

ObjCreator CheckCylDca::get_creator() {
  return create;
}

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

// Write the object data.

ObjData CheckCylDca::write_data() const {
  ObjData data( get_type_name() );
  data.add_double( "max_dca", get_max_dca() );
  data.add_double( "bfield", get_bfield() );
  return data;
}

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

// return track status.

bool CheckCylDca::status_with_record(const MTrack& trm, 
				     CutRecord* prec) const {
  //Calculate distance of closest approach in r-phi at beam line...

  // Check that we really have a cylinder
  const ETrack& tre=trm.get_track();
  assert( tre.get_surface()->get_pure_type()  ==
	  SurfCylinder::get_static_type()     );
  if( tre.get_surface()->get_pure_type()  !=
      SurfCylinder::get_static_type()     )
  {
    cerr << "CheckCylDca: Surface Type is NOT Cylinder" << endl;
    abort();
  } 
  
  // Fetch the originating TrackVector
  TrackVector vec_scy = tre.get_vector();
  
  // Get the track parameters

  double phipos  = vec_scy(0);        // phi_position
  double alfa    = vec_scy(2);        // phi_dir - phi_pos

  if(vec_scy(4)==0.) return false;
  double qpt     = 1./vec_scy(4);     // signed pT of the track
  double radius  = trm.get_track().space_point().rxy(); //radius of cylinder

  double phidir = alfa+phipos;

  double x0 = radius*cos(phipos); // x position of the track on the cylinder
  double y0 = radius*sin(phipos); // y position of the track on the cylinder
  double dca;
  double bfield=-(*_bfield)();
  if(bfield!=0.) { 
    double rtrk = qpt/(0.003*bfield); // Radius of curvature for the track
    
    double xc = x0 - rtrk*sin(phidir); // x of center of track circle
    double yc = y0 + rtrk*cos(phidir); // y of center of track circle
    double sign;
    rtrk>0 ? sign=1. : sign = -1.;
    dca = rtrk - sign*sqrt(xc*xc+yc*yc);
  }
  else {
    dca = radius*sin(alfa);  // No field, straight-line extrapolation
  }

  // Cut

  bool passed = abs(dca) <= _max_dca;

  // Fill cut record, if defined.

  if(prec && _cut_id.is_valid()) {
    prec->add_cut(_cut_id, dca);
    prec->set_passed(passed && prec->passed());
  }

  return passed;
}

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