CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/RecoPixelVertexing/PixelTrackFitting/src/PixelFitterByHelixProjections.cc

Go to the documentation of this file.
00001 #include "RecoPixelVertexing/PixelTrackFitting/interface/PixelFitterByHelixProjections.h"
00002 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00003 
00004 #include "FWCore/Framework/interface/EventSetup.h"
00005 #include "FWCore/Framework/interface/ESHandle.h"
00006 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00007 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00008 
00009 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
00010 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00011 #include "DataFormats/GeometryCommonDetAlgo/interface/GlobalError.h"
00012 
00013 #include "TrackingTools/DetLayers/interface/DetLayer.h"
00014 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
00015 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
00016 //#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
00017 #include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
00018 
00019 #include "MagneticField/Engine/interface/MagneticField.h"
00020 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00021 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
00022 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
00023 
00024 #include "CommonTools/Statistics/interface/LinearFit.h"
00025 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
00026 
00027 #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
00028 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
00029 #include "FWCore/Framework/interface/ESWatcher.h"
00030 
00031 
00032 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00033 
00034 #include "RZLine.h"
00035 #include "CircleFromThreePoints.h"
00036 #include "RecoPixelVertexing/PixelTrackFitting/interface/PixelTrackBuilder.h"
00037 #include "RecoPixelVertexing/PixelTrackFitting/interface/PixelTrackErrorParam.h"
00038 #include "DataFormats/GeometryVector/interface/Pi.h"
00039 
00040 using namespace std;
00041 
00042 namespace {
00043 
00044   int charge(const std::vector<GlobalPoint> & points) {
00045     // the cross product will tell me...
00046     float dir = (points[1].x()-points[0].x())*(points[2].y()-points[1].y())
00047       - (points[1].y()-points[0].y())*(points[2].x()-points[1].x());
00048     
00049     /*
00050       GlobalVector v21 = points[1]-points[0];
00051       GlobalVector v32 = points[2]-points[1];
00052       float dphi = v32.phi() - v21.phi();
00053       while (dphi >  Geom::fpi()) dphi -=  Geom::ftwoPi();
00054       while (dphi < -Geom::fpi()) dphi +=  Geom::ftwoPi();
00055       return (dphi > 0) ? -1 : 1;
00056     */
00057     return (dir>0) ? -1 : 1;
00058   }
00059 
00060   float cotTheta(const GlobalPoint& inner, const GlobalPoint& outer) {
00061     float dr = outer.perp()-inner.perp();
00062     float dz = outer.z()-inner.z();
00063     return (std::abs(dr) > 1.e-3f) ? dz/dr : 0;
00064   }
00065 
00066   inline float phi(float xC, float yC, int charge) {
00067     return  (charge>0) ? std::atan2(xC,-yC) :  std::atan2(-xC,yC);
00068   }
00069 
00070   float zip(float d0, float phi_p, float curv, 
00071             const GlobalPoint& pinner, const GlobalPoint& pouter) {
00072     //
00073     //phi = asin(r*rho/2) with asin(x) ~= x+x**3/(2*3) = x(1+x*x/6);
00074     //
00075     
00076     float phi0 = phi_p - Geom::fhalfPi();
00077     GlobalPoint pca(d0*std::cos(phi0), d0*std::sin(phi0),0.);
00078     
00079     float rho2 = curv*curv;
00080     float r1s = (pinner-pca).perp2();
00081     double phi1 = std::sqrt(r1s)*(curv*0.5f)*(1.f+r1s*(rho2/24.f));
00082     float r2s = (pouter-pca).perp2();
00083     double phi2 = std::sqrt(r2s)*(curv*0.5f)*(1.f+r2s*(rho2/24.f));
00084     double z1 = pinner.z();
00085     double z2 = pouter.z();
00086 
00087     if (fabs(curv)>1.e-5) 
00088       return z1 - phi1/(phi1-phi2)*(z1-z2);
00089     else {
00090       double dr = std::max(std::sqrt(r2s)-std::sqrt(r1s),1.e-5f);
00091       return z1-std::sqrt(r1s)*(z2-z1)/dr;
00092     }
00093   }
00094 }
00095   
00096 PixelFitterByHelixProjections::PixelFitterByHelixProjections(
00097    const edm::ParameterSet& cfg) 
00098  : theConfig(cfg), theTracker(0), theField(0), theTTRecHitBuilder(0) {}
00099 
00100 reco::Track* PixelFitterByHelixProjections::run(
00101     const edm::EventSetup& es,
00102     const std::vector<const TrackingRecHit * > & hits,
00103     const TrackingRegion & region) const
00104 {
00105   int nhits = hits.size();
00106   if (nhits <2) return 0;
00107 
00108   vector<GlobalPoint> points(nhits);
00109   vector<GlobalError> errors(nhits);
00110   vector<bool> isBarrel(nhits);
00111   
00112   static edm::ESWatcher<TrackerDigiGeometryRecord> watcherTrackerDigiGeometryRecord;
00113   if (!theTracker || watcherTrackerDigiGeometryRecord.check(es)) {
00114     edm::ESHandle<TrackerGeometry> trackerESH;
00115     es.get<TrackerDigiGeometryRecord>().get(trackerESH);
00116     theTracker = trackerESH.product();
00117   }
00118 
00119   static edm::ESWatcher<IdealMagneticFieldRecord>  watcherIdealMagneticFieldRecord;
00120   if (!theField || watcherIdealMagneticFieldRecord.check(es)) {
00121     edm::ESHandle<MagneticField> fieldESH;
00122     es.get<IdealMagneticFieldRecord>().get(fieldESH);
00123     theField = fieldESH.product();
00124   }
00125 
00126   static edm::ESWatcher<TransientRecHitRecord> watcherTransientRecHitRecord;
00127   if (!theTTRecHitBuilder || watcherTransientRecHitRecord.check(es)) {
00128     edm::ESHandle<TransientTrackingRecHitBuilder> ttrhbESH;
00129     std::string builderName = theConfig.getParameter<std::string>("TTRHBuilder");
00130     es.get<TransientRecHitRecord>().get(builderName,ttrhbESH);
00131     theTTRecHitBuilder = ttrhbESH.product();
00132   }
00133 
00134 
00135   for ( int i=0; i!=nhits; ++i) {
00136     TransientTrackingRecHit::RecHitPointer recHit = theTTRecHitBuilder->build(hits[i]);
00137     points[i]  = GlobalPoint( recHit->globalPosition().x()-region.origin().x(), 
00138                               recHit->globalPosition().y()-region.origin().y(),
00139                               recHit->globalPosition().z()-region.origin().z() 
00140                               );
00141     errors[i] = recHit->globalPositionError();
00142     isBarrel[i] = recHit->detUnit()->type().isBarrel();
00143   }
00144 
00145   CircleFromThreePoints circle = (nhits==2) ?
00146         CircleFromThreePoints( GlobalPoint(0.,0.,0.), points[0], points[1]) :
00147         CircleFromThreePoints(points[0],points[1],points[2]); 
00148 
00149   float valPhi, valTip, valPt;
00150 
00151   int iCharge = charge(points);
00152   float curvature = circle.curvature();
00153 
00154   if (curvature > 1.e-4) {
00155     float invPt = PixelRecoUtilities::inversePt( circle.curvature(), es);
00156     valPt = (invPt > 1.e-4f) ? 1.f/invPt : 1.e4f;
00157     CircleFromThreePoints::Vector2D center = circle.center();
00158     valTip = iCharge * (center.mag()-1.f/curvature);
00159     valPhi = phi(center.x(), center.y(), iCharge);
00160   } 
00161   else {
00162     valPt = 1.e4f; 
00163     GlobalVector direction(points[1]-points[0]);
00164     valPhi =  direction.phi(); 
00165     valTip = -points[0].x()*sin(valPhi) + points[0].y()*cos(valPhi); 
00166   }
00167 
00168   float valCotTheta = cotTheta(points[0],points[1]);
00169   float valEta = asinh(valCotTheta);
00170   float valZip = zip(valTip, valPhi, curvature, points[0],points[1]);
00171 
00172   PixelTrackErrorParam param(valEta, valPt);
00173   float errValPt  = param.errPt();
00174   float errValCot = param.errCot();
00175   float errValTip = param.errTip();
00176   float errValPhi = param.errPhi();
00177   float errValZip = param.errZip();
00178 
00179 
00180   float chi2 = 0;
00181   if (nhits > 2) {
00182     RZLine rzLine(points,errors,isBarrel);
00183     float cottheta, intercept, covss, covii, covsi; 
00184     rzLine.fit(cottheta, intercept, covss, covii, covsi);
00185     chi2 = rzLine.chi2(cottheta, intercept);         //FIXME: check which intercept to use!
00186   }
00187 
00188   PixelTrackBuilder builder;
00189   Measurement1D pt(valPt, errValPt);
00190   Measurement1D phi(valPhi, errValPhi);
00191   Measurement1D cotTheta(valCotTheta, errValCot);
00192   Measurement1D tip(valTip, errValTip);
00193   Measurement1D zip(valZip, errValZip);
00194 
00195   return builder.build(pt, phi, cotTheta, tip, zip, chi2, iCharge, hits, theField, region.origin() );
00196 }
00197 
00198 
00199 
00200