RkFitCylinder.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------
// File from RkFit module
//
// Filename : RkFitCylinder.cc
//------------------------------------------------------------------------
// Description :
// Cylinder is an Element whose shape is a cylinder.
//------------------------------------------------------------------------
// Modif :
//------------------------------------------------------------------------
#include "CLHEP/Geometry/Point3D.h"
#include <float.h>
#ifndef ENABLE_BACKWARDS_COMPATIBILITY
typedef HepGeom::Point3D<double> HepPoint3D;
#endif
#include "TrackUtil/Helix.h"
#include "TrkReco/RkFitCylinder.h"
#include "TrkReco/RkFitMaterial.h"
#include "TrkReco/TRunge.h"
 
double RkFitCylinder::intersect( TRunge& track, HepPoint3D& x ) const {
  double dPhi[4];
  dPhi[0] = track.intersect_cylinder( ro_ );
  if ( dPhi[0] == 0 ) return -1;
  dPhi[1] = track.intersect_cylinder( ri_ );
  if ( dPhi[1] == 0 ) return -1;
  dPhi[2] = track.intersect_xy_plane( zf_ );
  dPhi[3] = track.intersect_xy_plane( zb_ );
 
  int n[2];
  int j = 0;
  for ( int i = 0; i < 4 && j < 2; i++ )
  {
    HepPoint3D xx = track.helix().x( dPhi[i] );
    if ( isInside( xx ) ) n[j++] = i;
  }
  if ( j < 2 ) return -1;
 
  x = track.helix().x( ( dPhi[n[0]] + dPhi[n[1]] ) * .5 );
 
  double tanl = track.helix().tanl();
  // cout<<"RkFitCylinder: track radius"<<track.radius()<<" dphi0 "
  //     <<dPhi[n[0]]<<" dphi1 "<<dPhi[n[1]]<<" tanl "<<tanl<<endl;
  return fabs( track.helix().radius() * ( dPhi[n[0]] - dPhi[n[1]] ) *
               sqrt( 1 + tanl * tanl ) );
  //   return 0;
}
 
double RkFitCylinder::intersect( TRunge& track, HepPoint3D& x,
                                 const HepPoint3D& point ) const {
 
  const double ro = sqrt( point.x() * point.x() + point.y() * point.y() );
 
  // std::cout<<" ro: "<<ro<<std::endl;
 
  double dPhi[4];
  dPhi[0] = track.intersect_cylinder( ro );
  if ( dPhi[0] == 0 ) return -1;
  dPhi[1] = track.intersect_cylinder( ri_ );
  if ( dPhi[1] == 0 ) return -1;
  dPhi[2] = track.intersect_xy_plane( zf_ );
  dPhi[3] = track.intersect_xy_plane( zb_ );
 
  // for(int ii=0; ii<4; ii++)
  // std::cout<<"dPhi["<<ii<<"]"<<dPhi[ii]<<std::endl;
 
  int n[2];
  int j = 0;
  for ( int i = 0; i < 4 && j < 2; i++ )
  {
    HepPoint3D xx = track.helix().x( dPhi[i] );
    if ( isInside( xx ) ) n[j++] = i;
  }
 
  if ( j < 2 ) return -1;
 
  x = track.helix().x( ( dPhi[n[0]] + dPhi[n[1]] ) * .5 );
 
  double tanl = track.helix().tanl();
 
  return fabs( track.helix().radius() * ( dPhi[n[0]] - dPhi[n[1]] ) *
               sqrt( 1 + tanl * tanl ) );
}
 
void RkFitCylinder::updateTrack( TRunge& track, double y[6] ) const {
  HepPoint3D x;
  double     path = intersect( track, x );
  double     mass = 0.000511;
  if ( path > 0 )
  {
    // move pivot
    // multiple scattering and energy loss
    // if(muls_) track.ms(path, *material_, index_element);
    track.eloss( path, material_, mass, y, 1 );
    //    track.Propagate(path,y);
  }
}
 
bool RkFitCylinder::isInside( const HepPoint3D& x ) const {
  double r = x.perp();
  double z = x.z();
  // std::cout<<"r: "<<r<<" z: "<<z<<" ri: "<<ri_<<" ro: "<<ro_<<" zb_: "<<zb_<<"zf:
  // "<<zf_<<std::endl;
 
  return ( r >= ri_ - FLT_EPSILON && r <= ro_ + FLT_EPSILON && z >= zb_ - FLT_EPSILON &&
           z <= zf_ + FLT_EPSILON );
}
 
bool RkFitCylinder::isInside2( const HepPoint3D& x ) const {
  double r = x.perp();
  double z = x.z();
  // std::cout<<"r: "<<r<<" z: "<<z<<" ri: "<<ri_<<" ro: "<<ro_<<" zb_: "<<zb_<<"zf:
  // "<<zf_<<std::endl;
 
  return ( r <= ro_ + FLT_EPSILON && z >= zb_ - FLT_EPSILON && z <= zf_ + FLT_EPSILON );
}