00001 //----------------------------------------------------------------------- 00002 // File from KalFit module 00003 // 00004 // Filename : KalFitCylinder.cc 00005 //------------------------------------------------------------------------ 00006 // Description : 00007 // Cylinder is an Element whose shape is a cylinder. 00008 //------------------------------------------------------------------------ 00009 // Modif : 00010 //------------------------------------------------------------------------ 00011 #include <float.h> 00012 #include "CLHEP/Geometry/Point3D.h" 00013 #ifndef ENABLE_BACKWARDS_COMPATIBILITY 00014 typedef HepGeom::Point3D<double> HepPoint3D; 00015 #endif 00016 //#include "TrackUtil/Helix.h" 00017 #include "KalFitAlg/helix/Helix.h" 00018 00019 #include "KalFitAlg/KalFitTrack.h" 00020 #include "KalFitAlg/KalFitMaterial.h" 00021 #include "KalFitAlg/KalFitCylinder.h" 00022 00023 double KalFitCylinder::intersect(const KalFitTrack& track, 00024 HepPoint3D& x) const 00025 { 00026 double dPhi[4]; 00027 dPhi[0] = track.intersect_cylinder(ro_); 00028 if(dPhi[0] == 0) return -1; 00029 dPhi[1] = track.intersect_cylinder(ri_); 00030 if(dPhi[1] == 0) return -1; 00031 dPhi[2] = track.intersect_xy_plane(zf_); 00032 dPhi[3] = track.intersect_xy_plane(zb_); 00033 00034 int n[2]; 00035 int j = 0; 00036 for(int i = 0; i < 4 && j < 2; i++){ 00037 HepPoint3D xx = track.x(dPhi[i]); 00038 if(isInside(xx)) n[j++] = i; 00039 } 00040 if(j < 2) return -1; 00041 00042 x = track.x((dPhi[n[0]] + dPhi[n[1]]) * .5); 00043 00044 double tanl = track.tanl(); 00045 //cout<<"KalFitCylinder: track radius"<<track.radius()<<" dphi0 " 00046 // <<dPhi[n[0]]<<" dphi1 "<<dPhi[n[1]]<<" tanl "<<tanl<<endl; 00047 return fabs(track.radius() * (dPhi[n[0]] - dPhi[n[1]]) 00048 * sqrt(1 + tanl * tanl)); 00049 } 00050 00051 00052 double KalFitCylinder::intersect(const KalFitTrack& track, 00053 HepPoint3D& x, const HepPoint3D& point) const 00054 { 00055 00056 const double ro = sqrt(point.x()*point.x()+point.y()*point.y()); 00057 00058 //std::cout<<" ro: "<<ro<<std::endl; 00059 00060 double dPhi[4]; 00061 dPhi[0] = track.intersect_cylinder(ro); 00062 if(dPhi[0] == 0) return -1; 00063 dPhi[1] = track.intersect_cylinder(ri_); 00064 if(dPhi[1] == 0) return -1; 00065 dPhi[2] = track.intersect_xy_plane(zf_); 00066 dPhi[3] = track.intersect_xy_plane(zb_); 00067 00068 //for(int ii=0; ii<4; ii++) 00069 //std::cout<<"dPhi["<<ii<<"]"<<dPhi[ii]<<std::endl; 00070 00071 int n[2]; 00072 int j = 0; 00073 for(int i = 0; i < 4 && j < 2; i++){ 00074 HepPoint3D xx = track.x(dPhi[i]); 00075 if(isInside(xx)) n[j++] = i; 00076 } 00077 00078 if(j < 2) return -1; 00079 00080 x = track.x((dPhi[n[0]] + dPhi[n[1]]) * .5); 00081 00082 double tanl = track.tanl(); 00083 00084 return fabs(track.radius() * (dPhi[n[0]] - dPhi[n[1]]) 00085 * sqrt(1 + tanl * tanl)); 00086 } 00087 00088 00089 00090 00091 bool KalFitCylinder::isInside(const HepPoint3D& x) const 00092 { 00093 double r = x.perp(); 00094 double z = x.z(); 00095 //std::cout<<"r: "<<r<<" z: "<<z<<" ri: "<<ri_<<" ro: "<<ro_<<" zb_: "<<zb_<<"zf: "<<zf_<<std::endl; 00096 00097 return (r >= ri_ - FLT_EPSILON && 00098 r <= ro_ + FLT_EPSILON && 00099 z >= zb_ - FLT_EPSILON && 00100 z <= zf_ + FLT_EPSILON); 00101 } 00102 00103 00104 bool KalFitCylinder::isInside2(const HepPoint3D& x) const 00105 { 00106 double r = x.perp(); 00107 double z = x.z(); 00108 //std::cout<<"r: "<<r<<" z: "<<z<<" ri: "<<ri_<<" ro: "<<ro_<<" zb_: "<<zb_<<"zf: "<<zf_<<std::endl; 00109 00110 return (r <= ro_ + FLT_EPSILON && 00111 z >= zb_ - FLT_EPSILON && 00112 z <= zf_ + FLT_EPSILON); 00113 } 00114
1.4.7