// PropDCACyl.cpp

#include "PropDCACyl.h"
#include <iostream>
#include <cassert>
#include "objstream/ObjData.hpp"
#include "objstream/ObjTable.hpp"
#include "trfutil/trfstream.h"
#include "trfutil/TRFMath.h"
#include "trfbase/TrackVector.h"
#include "trfbase/VTrack.h"
#include "trfbase/ETrack.h"
#include "trfbase/PropStat.h"
#include "trfcyl/SurfCylinder.h"
#include "SurfDCA.h"

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

#ifndef DEFECT_CMATH_NOT_STD
using std::sqrt;
using std::cos;
using std::sin;
using std::atan2;
using std::asin;
#endif

using namespace trf;

//**********************************************************************
// Local definitions.
//**********************************************************************

namespace {

// Creator.

ObjPtr create(const ObjData& data) {
  assert( data.get_object_type() == "PropDCACyl" );
  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 << "PropDCACyl: can't get bfield parameter." << endl;
    abort();
  }
  return ObjPtr( new PropDCACyl(bfield) );
}

}

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

// Assign track parameter indices

enum {
  IRSIGNED = SurfDCA::IRSIGNED,
  IZ_DCA   = SurfDCA::IZ,
  IPHID    = SurfDCA::IPHID,
  ITLM_DCA = SurfDCA::ITLM,
  IQPT_DCA = SurfDCA::IQPT,

  IPHI     = SurfCylinder::IPHI,
  IZ_CYL   = SurfCylinder::IZ,
  IALF     = SurfCylinder::IALF,
  ITLM_CYL = SurfCylinder::ITLM,
  IQPT_CYL = SurfCylinder::IQPT
};

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

// Private class STCalc.
//
// An STCalc_ object calculates sT (the signed transverse path length)
// and its derivatives w.r.t. alf1 and crv1.  It is constructed from
// the starting (r1, phi1, alf1, crv1) and final track parameters
// (r2, phi2, alf2) assuming these are consistent.  Methods are
// provided to retrieve sT and the two derivatives.

class ST_DCACyl {
private:
  bool _big_crv;
  double _st;
  //  double _dst_dalf2;
  double _dst_dr1;
  double _dst_dcrv1;
  double _dst_dphi2;
  double _dst_dphi2_or;
public:
  // constructor
  ST_DCACyl();
  ST_DCACyl(double r1, double phi1, double alf1, double crv,
         double r2, double phi2, double alf2);
  double st() const { return _st; };
  double d_st_dr1_sign(double sign,   double d_phi2_dr1,   double d_alf2_dr1   ) const;
  double d_st_dr1(double d_phi2_dr1,   double d_alf2_dr1   ) const;
  double d_st_dcrv1( double d_phi2_dcrv1, double d_alf2_dcrv1 ) const;
  double d_st_dalf1_or( double r1,double dphi2_dalf1 , double dalf2_dalf1  ) const;
  double _crv1;
};

ST_DCACyl::ST_DCACyl() { }

ST_DCACyl::ST_DCACyl(double r1, double phi1, double alf1, double crv1,
               double r2, double phi2, double alf2) {

  _crv1 = crv1;
  assert( r1 >= 0.0 );
  assert( r2 >= 0.0 );
  double rmax = r1+r2;

  // Calculate the change in xy direction
  double phi_dir_diff = fmod2(phi2+alf2-phi1-alf1,TWOPI);
  assert( fabs(phi_dir_diff) <= PI );

  // Evaluate whether |C| is "big"
  _big_crv = rmax*fabs(crv1) > 0.001;

  // If the curvature is big we can use
  // sT = (phi_dir2 - phi_dir1)/crv1
  if ( _big_crv ) {
    assert( crv1 != 0.0 );
    _st = phi_dir_diff/crv1;
  }

  // Otherwise, we calculate the straight-line distance
  // between the points and use an approximate correction
  // for the (small) curvature.
  else {

    // evaluate the distance
    double d = sqrt( r1*r1 + r2*r2 - 2.0*r1*r2*cos(phi2-phi1) );
    double arg = 0.5*d*crv1;
    double arg2 = arg*arg;
    float st_minus_d = d*arg2*( 1.0/6.0 + 3.0/40.0*arg2 );
    _st = d + st_minus_d;

    // evaluate the sign
    // We define a metric xsign = abs( (dphid-d*C)/(d*C) ).
    // Because sT*C = dphid and d = abs(sT):
    // xsign = 0 for sT > 0
    // xsign = 2 for sT < 0
    // Numerical roundoff will smear these predictions.
    double sign = 0.0;
    if ( crv1*_st ) {
      double xsign = fabs( (phi_dir_diff - _st*crv1) / (_st*crv1) );
      if ( xsign < 0.5 ) sign = 1.0;
      if ( xsign > 1.5  &&  xsign < 3.0 ) sign = -1.0;
    }
    // If the above is indeterminate, assume zero curvature.
    // In this case abs(alpha) decreases monotonically
    // with sT.  Track passing through origin has alpha = 0 on one
    // side and alpha = +/-pi on the other.  If both points are on
    // the same side, we use dr/ds > 0 for |alpha|<pi/2.
    if ( ! sign ) {
      sign = 1.0;
      if ( fabs(alf2) > fabs(alf1) ) sign = -1.0;
      if ( fabs(alf2) == fabs(alf1) ) {
        if ( fabs(alf2) < PI2 ) {
          if ( r2 < r1 ) sign = -1.0;
        } else {
          if ( r2 > r1 ) sign = -1.0;
        }
      }
    }

    // Correct _st using the above sign.
    assert( fabs(sign) == 1.0 );
    _st = sign*_st;

    // save derivatives
    //    _dst_dalf2 = (_st/d)*2.0*(r1-r2*cos(phi2-phi1));
    //_dst_dphi2 = (_st/d)*2.0*r1*r2*sin(phi2-phi1);
    if ( is_zero(d) ) {
      _dst_dcrv1 = 0.0;
      _dst_dphi2 = sign*r1;
      _dst_dphi2_or = sign;
      _dst_dr1 =  0.0;
    } else {
      _dst_dcrv1 = sign*d*d*arg*( 1.0/6.0 + 3.0/20.0*arg2);
      double root = (1.0 + 0.5*arg*arg + 3.0/8.0*arg*arg*arg*arg );
      _dst_dphi2_or = sign*(r2*sin(phi2-phi1))*root/d;
      _dst_dphi2 = _dst_dphi2_or*r1;
      _dst_dr1 =   sign*(r1-r2*cos(phi2-phi1))*root/d;
    }
  }

}

double ST_DCACyl::d_st_dr1_sign(double sign_alf1,double dphi2_dr1_sign,double dalf2_dr1_sign) const {
  if ( _big_crv ) return ( dphi2_dr1_sign + dalf2_dr1_sign ) / _crv1;
  else return   _dst_dphi2 * dphi2_dr1_sign + _dst_dr1*sign_alf1;
}
double ST_DCACyl::d_st_dr1(double dphi2_dr1,double dalf2_dr1) const {
  if ( _big_crv ) return ( dphi2_dr1 + dalf2_dr1 ) / _crv1;
  else return   _dst_dphi2 * dphi2_dr1 + _dst_dr1;
}

double ST_DCACyl::d_st_dcrv1(double dphi2_dcrv1,double dalf2_dcrv1) const {
  if ( _big_crv ) return ( dphi2_dcrv1 + dalf2_dcrv1 - _st ) / _crv1;
  else return   _dst_dphi2 * dphi2_dcrv1 + _dst_dcrv1;
}

double ST_DCACyl::d_st_dalf1_or(double r1,double dphi2_dalf1_m1_or, double dalf2_dalf1_or) const {
  if ( _big_crv ) return ( dphi2_dalf1_m1_or + dalf2_dalf1_or) / _crv1;
  else return _dst_dphi2_or * (dphi2_dalf1_m1_or*r1+1);
}


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

PropStat
dca_cyl_propagate_(       double             _bfac,
                          VTrack&             trv,
                    const Surface&            srf,
                    const Propagator::PropDir dir,
                          TrackDerivative* pderiv =0 ) {

  // construct return status
  PropStat pstat;

  // fetch the originating Surface and check it is a DCA surface 
  const Surface& srf0 = *trv.get_surface();
  assert( srf0.get_pure_type() == SurfDCA::get_static_type() );
  if ( srf0.get_pure_type( ) != SurfDCA::get_static_type() )
    return pstat;
  const SurfDCA& srf_dca = (const SurfDCA&)srf0;

  // check that the destination surface is a Cylinder
  assert( srf.get_pure_type() == SurfCylinder::get_static_type() );
  if ( srf.get_pure_type( ) != SurfCylinder::get_static_type() )
    return pstat;
  const SurfCylinder& srf_cyl = (const SurfCylinder&) srf;

  // fetch the originating TrackVector
  TrackVector vec_dca = trv.get_vector();
  double r1p    = fabs(vec_dca(IRSIGNED));             // r           
  double z1    = vec_dca(IZ_DCA);                // z
  double phid1p = vec_dca(IPHID);                 // phi_direction
  double tlam1 = vec_dca(ITLM_DCA);              // tan(lambda)
  double qpt1  = vec_dca(IQPT_DCA);              // q/pT

  double xv = srf_dca.get_x();
  double yv = srf_dca.get_y();

  // calculate alf1 and phi1
  double sign_alf1p;
  double sign_r1p=0.;
  if ( !is_zero(vec_dca(IRSIGNED)) ) {
    sign_alf1p  = vec_dca(IRSIGNED)/fabs(vec_dca(IRSIGNED));
    sign_r1p = sign_alf1p;
  }
  else
    sign_alf1p  = 0.0;

  if( (!is_zero(xv) || !is_zero(yv)) &&  is_zero(sign_alf1p) )
    sign_alf1p=1;

  double alf1p = sign_alf1p * PI2;                 // alpha
  double phi1p = phid1p - alf1p;                    // phi_position
  phi1p = fmod2( phi1p, TWOPI );


  // calculate salf1, calf1 and crv1
  double salf1p=0.;
  salf1p = alf1p>0. ? 1.:( alf1p < 0. ? -1. : 0.) ;
    
  double cosphi1p=cos(phi1p);
  double sinphi1p=sin(phi1p);

  double x1=r1p*cosphi1p+xv;
  double y1=r1p*sinphi1p+yv;

  double phi1 = atan2(y1,x1);       // phi position in (xv,yv) coord system
  double r1 = sqrt(x1*x1+y1*y1);    // r1 in (xv,yv) coord system
  double alf1 = alf1p + phi1p - phi1; // alf1 in (xv,yv) coord system
  if( is_zero(x1) && is_zero(y1)  ) {
    phi1 = phid1p;
    alf1 = 0.;
  }
  if( sign_r1p == 0 && !is_zero(r1) )
    sign_r1p=1;

  // calculate salf1, calf1 and crv1
  double salf1 = sin( alf1 );
  double calf1 = cos( alf1 );  

  if( is_zero(alf1) ) {
    salf1=0.;
    calf1=1.;
  }
  if( is_equal(fabs(alf1),PI2) ) {
    salf1= alf1 > 0 ? 1. : -1. ;
    calf1=0.;
  }
 
  double  crv1 = _bfac * qpt1;
  double sign_crv = 1.;

  // fetch r2 of the destination Cylinder
  double r2 = srf_cyl.get_parameter(SurfCylinder::RADIUS);

  // calculate tlam2, qpt2 and crv2 of the destination Cylinder
  double tlam2 = tlam1;
  double qpt2  = qpt1;
  double crv2  = crv1;

  // calculate salf2 of the destination Cylinder
  double salf2 = r1/r2*salf1 + 0.5*crv1/r2*(r2*r2-r1*r1);
   
  // If salf2 is close to 1 or -1, set it to that value.
  double diff = fabs( fabs(salf2) - 1.0 );
  if ( diff < 1.e-10 ) salf2 = salf2>0 ? 1.0 : -1.0;
  // if salf2 > 1, track does not cross cylinder
  if ( fabs(salf2) > 1.0 ) return pstat;

  // there are two possibilities for alf2
  double alf21 = asin( salf2 );
  double alf22 = alf21>0 ? PI-alf21 : -PI-alf21;
  double calf21 = cos( alf21 );
  double calf22 = cos( alf22 );
  //double phi21 = phi1 + atan2( -sign_crv*calf21, r2*crv2-salf2 );
  //double phi22 = phi1 + atan2( -sign_crv*calf22, r2*crv2-salf2 );
  double cnst = atan2(salf1-r1*crv1,calf1);
  if( is_equal(fabs(alf1),PI2) ) {
    cnst = salf1-r1*crv1 > 0 ? PI2 : -PI2;
    cnst = (r1==0.) ? 0. : cnst;
  }
  double phi21 = phi1 + cnst 
                      - atan2( salf2-r2*crv2, calf21 );
  double phi22 = phi1 + cnst
                      - atan2( salf2-r2*crv2, calf22 );

  if( is_zero(crv1) ) {
    phi21 = phi1 + alf1  - alf21;
    phi22 = phi1 + alf1  - alf22;
  }

  if ( is_zero(calf21) ) {
    phi21 = phi1;
    phi22 = phi1;
  }

  // construct an sT object for each solution
  ST_DCACyl sto1(r1,phi1,alf1,crv1,r2,phi21,alf21);
  ST_DCACyl sto2(r1,phi1,alf1,crv1,r2,phi22,alf22);
  // check the two solutions are nonzero and have opposite sign
  // or at least one is nonzero

  // choose the correct solution
  bool use_first_solution;
  switch (dir) {
  case Propagator::NEAREST:
    use_first_solution = abs(sto2.st()) > abs(sto1.st());
    break;
  case Propagator::NEAREST_MOVE:
    use_first_solution = abs(sto2.st()) > abs(sto1.st());
    break;
  case Propagator::FORWARD:
    use_first_solution = sto1.st() > 0.0;
    break;
  case Propagator::BACKWARD:
    use_first_solution = sto1.st() < 0.0;
    break;
  default:
    cerr << "PropCyl::_vec_propagate: Unknown direction." << endl;
    abort();
  }

  // assign phi2, alf2 and sto2 for the chosen solution
  double phi2, alf2;
  ST_DCACyl sto;
  double calf2;
  if ( use_first_solution ) {
    sto = sto1;
    phi2 = phi21;
    alf2 = alf21;
    calf2 = calf21;
   } else {
    sto = sto2;
    phi2 = phi22;
    alf2 = alf22;
    calf2 = calf22;
   }

  // check alpha range
  assert( fabs(alf2) <= PI );

  // fetch sT
  double st = sto.st();
  double s = st*sqrt(1+tlam1*tlam1);

  // calculate z2 of the destination Cylinder
  double z2 = z1 + st * tlam1;

  // construct the destination TrackVector
  TrackVector vec_cyl;
  vec_cyl(IPHI)     = phi2;
  vec_cyl(IZ_CYL)   = z2;
  vec_cyl(IALF)     = alf2;
  vec_cyl(ITLM_CYL) = tlam2;
  vec_cyl(IQPT_CYL) = qpt2;

  // set the surface of trv to the destination Cylinder
  trv.set_surface( SurfacePtr(srf_cyl.new_pure_surface()) );

  // set the vector of trv to the destination TrackVector (Cyl. coord.)
  trv.set_vector(vec_cyl);

  // set the direction of trv
  if ( fabs(alf2) <= PI2 ) trv.set_forward();
  else trv.set_backward();

  // set the return status
  pstat.set_path_distance(s);

  // exit now if user did not ask for error matrix.
  if ( ! pderiv ) return pstat;


  // calculate derivatives
 
  //double dalf2_dr1   = (1.0-r1*crv1*salf1)/(r2*calf2);
  // commented out by SK ( DCA to DCA(xv,yv) change)
  double salf12 = sign_r1p*salf1;
  if( sign_r1p == 0 && salf1 == 0. )
    salf12 = 1;
  double dalf2_dr1_sign = (salf12-crv1*r1*sign_r1p)/r2/calf2;
  double dalf2_dr1 = (salf1-crv1*r1)/r2/calf2;
  double dalf2_dalf1_or = 1/r2*calf1/calf2; // over r1
  double dalf2_dalf1 = r1*dalf2_dalf1_or;
  double dalf2_dcrv1 = (r2*r2-r1*r1)/(2.0*r2*calf2);

  //double dphi2_dr1 = -sign_crv*
  //               (r2*crv2*salf2-1.0)*(r1*crv1-salf1)/
  //               (r2*calf2*(1.0 + r2*r2*crv2*crv2 - 2.0*r2*crv2*salf2));
  // commented out on 10.16.2001 by SK ( DCA to DCA(xv,yv) change)
  //double dphi2_dr1 = -sign_crv *
  //               (r2*crv2*salf2-1.0)*(sign_alf1*r1*crv1-1.)/
  //               (r2*calf2*(1.0 + r2*r2*crv2*crv2 - 2.0*r2*crv2*salf2));
  // 
  //double dphi2_dcrv1 = -sign_crv*
  //           ((1.0-r2*crv2*salf2)*(r2*r2-r1*r1)-2.0*r2*r2*calf2*calf2)/
  //           (2.0*r2*calf2*(1.0 + r2*r2*crv2*crv2 - 2.0*r2*crv2*salf2));

  double dphi2_dphi1=1.;
  double dphi2_dr1_sign = -crv1*calf1*sign_r1p/(1.-2.*r1*crv1*salf1+r1*r1*crv1*crv1)
          - dalf2_dr1_sign*(1-salf2*crv2*r2)/(1.-2.*r2*crv2*salf2+r2*r2*crv2*crv2);

  double dphi2_dr1 = -crv1*calf1/(1.-2.*r1*crv1*salf1+r1*r1*crv1*crv1)
          - dalf2_dr1*(1-salf2*crv2*r2)/(1.-2.*r2*crv2*salf2+r2*r2*crv2*crv2);
  double dphi2_dalf1_m1_or = (-salf1*crv1-r1*crv1*crv1+2.*crv1*salf1)/(1.-2.*r1*crv1*salf1+r1*r1*crv1*crv1)
    - dalf2_dalf1_or*(1-salf2*crv2*r2)/(1.-2.*r2*crv2*salf2+r2*r2*crv2*crv2); // minus 1
  double dphi2_dalf1 = dphi2_dalf1_m1_or*r1 + 1;
  double dphi2_dcrv1 = -r1*calf1/(1.-2.*r1*crv1*salf1+r1*r1*crv1*crv1)
         -(dalf2_dcrv1*(1.-r2*crv2*salf2) - r2*calf2)/(1.-2.*r2*crv2*salf2+r2*r2*crv2*crv2);
  

  //  double dst_dr1   = sto.d_st_dr1(sign_alf1*dphi2_dr1,sign_alf1*dalf2_dr1  );
  double dst_dr1_sign   = sto.d_st_dr1_sign(sign_r1p,dphi2_dr1_sign,dalf2_dr1_sign);  
  double dst_dr1   = sto.d_st_dr1(dphi2_dr1,dalf2_dr1);  
  double dst_dcrv1 = sto.d_st_dcrv1(dphi2_dcrv1, dalf2_dcrv1);
  double dst_dalf1_or = sto.d_st_dalf1_or(r1,dphi2_dalf1_m1_or, dalf2_dalf1_or);

  // build the derivative matrix
  // derivatives to prime coordinates
  
  double dphi1_dphi1p_tr = r1p*cos(phi1p-phi1); // times r1
  double dphi1_dr1p_tr = sign_r1p*sin(phi1p-phi1); // times r1

  double dalf1_dphi1p_tr = r1 - dphi1_dphi1p_tr;
  double dalf1_dr1p_tr = -dphi1_dr1p_tr;

  double dr1_dr1p_sign = cos(phi1p-phi1);
  double dr1_dphi1p = r1p*sin(phi1-phi1p);


  double dphi2_dphi1p = dphi2_dr1*dr1_dphi1p +dphi2_dalf1 - dphi1_dphi1p_tr*dphi2_dalf1_m1_or;
  double dphi2_dr1p = dphi2_dr1_sign*dr1_dr1p_sign - dphi2_dalf1_m1_or*dphi1_dr1p_tr;
  
  double dalf2_dphi1p = dalf2_dr1*dr1_dphi1p + dalf2_dalf1_or*dalf1_dphi1p_tr;
  double dalf2_dr1p = dalf2_dr1_sign*dr1_dr1p_sign + dalf2_dalf1_or*dalf1_dr1p_tr;

  double dst_dphi1p = dst_dr1*dr1_dphi1p + dst_dalf1_or*dalf1_dphi1p_tr;
  double dst_dr1p = dst_dr1_sign*dr1_dr1p_sign + dst_dalf1_or*dalf1_dr1p_tr;

  TrackDerivative& deriv = *pderiv;
  deriv.fill( 0.0 );

  //deriv(IPHI,     IRSIGNED) = sign_alf1 * dphi2_dr1;
  deriv(IPHI,     IRSIGNED) = dphi2_dr1p;
  deriv(IPHI,     IZ_DCA)   = 0.0;
  deriv(IPHI,     IPHID)    = dphi2_dphi1p;
  deriv(IPHI,     ITLM_DCA) = 0.0;
  deriv(IPHI,     IQPT_DCA) = _bfac * dphi2_dcrv1;

  deriv(IZ_CYL,   IRSIGNED) = tlam1 * dst_dr1p;
  deriv(IZ_CYL,   IZ_DCA)   = 1.0;
  deriv(IZ_CYL,   IPHID)    = tlam1 * dst_dphi1p;
  deriv(IZ_CYL,   ITLM_DCA) = st;
  deriv(IZ_CYL,   IQPT_DCA) = tlam1 * _bfac * dst_dcrv1;

  //deriv(IALF,     IRSIGNED) = sign_alf1 * dalf2_dr1;
  deriv(IALF,     IRSIGNED) = dalf2_dr1p;
  deriv(IALF,     IZ_DCA)   = 0.0;
  deriv(IALF,     IPHID)    = dalf2_dphi1p;
  deriv(IALF,     ITLM_DCA) = 0.0;
  deriv(IALF,     IQPT_DCA) = _bfac * dalf2_dcrv1;

  deriv(ITLM_CYL, IRSIGNED) = 0.0;
  deriv(ITLM_CYL, IZ_DCA)   = 0.0;
  deriv(ITLM_CYL, IPHID)    = 0.0;
  deriv(ITLM_CYL, ITLM_DCA) = 1.0;
  deriv(ITLM_CYL, IQPT_DCA) = 0.0;

  deriv(IQPT_CYL, IRSIGNED) = 0.0;
  deriv(IQPT_CYL, IZ_DCA)   = 0.0;
  deriv(IQPT_CYL, IPHID)    = 0.0;
  deriv(IQPT_CYL, ITLM_DCA) = 0.0;
  deriv(IQPT_CYL, IQPT_DCA) = 1.0;

  return pstat;

}


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

// Constructor

PropDCACyl::PropDCACyl(ObjDoublePtr bfield)
  : PropDirected(Propagator::FORWARD),
    _bfield(bfield) { }

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

// Destructor

PropDCACyl::~PropDCACyl() { }

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

// Clone

Propagator* PropDCACyl::new_propagator() const {
  return new PropDCACyl( _bfield );
}

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

// Return the creator.

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

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

// Write the object data.

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

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

// propagate a track without error in the specified direction

PropStat PropDCACyl::
vec_dir_prop(VTrack& trv, const Surface& srf, PropDir dir,
             TrackDerivative* pderiv ) const {
  double bfac = -BFAC*(*_bfield)();
  return dca_cyl_propagate_( bfac, trv, srf, dir, pderiv );
}

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

// propagate a track with error in the specified direction

// PropStat PropDCACyl::
// err_dir_prop(ETrack& tre, const Surface& srf, PropDir dir ) const {
//   TrackDerivative deriv;
//   PropStat pstat =  dca_cyl_propagate_( _bfac, tre, srf, dir, &deriv );
//   if ( ! pstat.success() ) return pstat;
//   TrackError err_dca = tre.get_error();
//   TrackError err_cyl = err_dca.Xform(deriv);
//   tre.set_error(err_cyl);
//   return pstat;
// }

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

// return the magnetic field

double PropDCACyl::get_bfield() const {
  return (*_bfield)();
}

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

// output stream

void PropDCACyl::ostr( ostream& stream ) const {
  stream << begin_object;
  stream << "DCA propagation to a Cylinder with constant "
         << get_bfield() << " Tesla field";
  stream << end_object;
}

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