00001 //----------------------------------------------------------------------- 00002 // File from RkFit module 00003 // 00004 // Filename : RkFitCylinder.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 "TrkReco/TRunge.h" 00018 #include "TrkReco/RkFitMaterial.h" 00019 #include "TrkReco/RkFitCylinder.h" 00020 00021 double RkFitCylinder::intersect( TRunge& track, 00022 HepPoint3D& x) const 00023 { 00024 double dPhi[4]; 00025 dPhi[0] = track.intersect_cylinder(ro_); 00026 if(dPhi[0] == 0) return -1; 00027 dPhi[1] = track.intersect_cylinder(ri_); 00028 if(dPhi[1] == 0) return -1; 00029 dPhi[2] = track.intersect_xy_plane(zf_); 00030 dPhi[3] = track.intersect_xy_plane(zb_); 00031 00032 int n[2]; 00033 int j = 0; 00034 for(int i = 0; i < 4 && j < 2; i++){ 00035 HepPoint3D xx = track.helix().x(dPhi[i]); 00036 if(isInside(xx)) n[j++] = i; 00037 } 00038 if(j < 2) return -1; 00039 00040 x = track.helix().x((dPhi[n[0]] + dPhi[n[1]]) * .5); 00041 00042 double tanl = track.helix().tanl(); 00043 //cout<<"RkFitCylinder: track radius"<<track.radius()<<" dphi0 " 00044 // <<dPhi[n[0]]<<" dphi1 "<<dPhi[n[1]]<<" tanl "<<tanl<<endl; 00045 return fabs(track.helix().radius() * (dPhi[n[0]] - dPhi[n[1]]) 00046 * sqrt(1 + tanl * tanl)); 00047 // return 0; 00048 } 00049 00050 00051 double RkFitCylinder::intersect( TRunge& track, 00052 HepPoint3D& x, const HepPoint3D& point) const 00053 { 00054 00055 const double ro = sqrt(point.x()*point.x()+point.y()*point.y()); 00056 00057 //std::cout<<" ro: "<<ro<<std::endl; 00058 00059 double dPhi[4]; 00060 dPhi[0] = track.intersect_cylinder(ro); 00061 if(dPhi[0] == 0) return -1; 00062 dPhi[1] = track.intersect_cylinder(ri_); 00063 if(dPhi[1] == 0) return -1; 00064 dPhi[2] = track.intersect_xy_plane(zf_); 00065 dPhi[3] = track.intersect_xy_plane(zb_); 00066 00067 //for(int ii=0; ii<4; ii++) 00068 //std::cout<<"dPhi["<<ii<<"]"<<dPhi[ii]<<std::endl; 00069 00070 int n[2]; 00071 int j = 0; 00072 for(int i = 0; i < 4 && j < 2; i++){ 00073 HepPoint3D xx = track.helix().x(dPhi[i]); 00074 if(isInside(xx)) n[j++] = i; 00075 } 00076 00077 if(j < 2) return -1; 00078 00079 x = track.helix().x((dPhi[n[0]] + dPhi[n[1]]) * .5); 00080 00081 double tanl = track.helix().tanl(); 00082 00083 return fabs(track.helix().radius() * (dPhi[n[0]] - dPhi[n[1]]) 00084 * sqrt(1 + tanl * tanl)); 00085 } 00086 00087 void RkFitCylinder::updateTrack(TRunge& track,double y[6]) const{ 00088 HepPoint3D x; 00089 double path = intersect(track, x); 00090 double mass=0.000511; 00091 if(path > 0){ 00092 // move pivot 00093 // multiple scattering and energy loss 00094 // if(muls_) track.ms(path, *material_, index_element); 00095 track.eloss(path, material_,mass,y,1); 00096 // track.Propagate(path,y); 00097 } 00098 00099 } 00100 00101 00102 bool RkFitCylinder::isInside(const HepPoint3D& x) const 00103 { 00104 double r = x.perp(); 00105 double z = x.z(); 00106 //std::cout<<"r: "<<r<<" z: "<<z<<" ri: "<<ri_<<" ro: "<<ro_<<" zb_: "<<zb_<<"zf: "<<zf_<<std::endl; 00107 00108 return (r >= ri_ - FLT_EPSILON && 00109 r <= ro_ + FLT_EPSILON && 00110 z >= zb_ - FLT_EPSILON && 00111 z <= zf_ + FLT_EPSILON); 00112 } 00113 00114 00115 bool RkFitCylinder::isInside2(const HepPoint3D& x) const 00116 { 00117 double r = x.perp(); 00118 double z = x.z(); 00119 //std::cout<<"r: "<<r<<" z: "<<z<<" ri: "<<ri_<<" ro: "<<ro_<<" zb_: "<<zb_<<"zf: "<<zf_<<std::endl; 00120 00121 return (r <= ro_ + FLT_EPSILON && 00122 z >= zb_ - FLT_EPSILON && 00123 z <= zf_ + FLT_EPSILON); 00124 } 00125
1.4.7