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 <float.h>
#include "CLHEP/Geometry/Point3D.h"
#ifndef ENABLE_BACKWARDS_COMPATIBILITY
typedef HepGeom::Point3D<double> HepPoint3D;
#endif
#include "TrackUtil/Helix.h"
#include "TrkReco/TRunge.h"
#include "TrkReco/RkFitMaterial.h"
#include "TrkReco/RkFitCylinder.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);
}