CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/RecoPixelVertexing/PixelLowPtUtilities/src/TrackFitter.cc

Go to the documentation of this file.
00001 #include "RecoPixelVertexing/PixelLowPtUtilities/interface/TrackFitter.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 "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
00017 
00018 #include "MagneticField/Engine/interface/MagneticField.h"
00019 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00020 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
00021 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
00022 
00023 #include "CommonTools/Statistics/interface/LinearFit.h"
00024 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
00025 
00026 #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
00027 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
00028 
00029 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00030 
00031 #include "RecoPixelVertexing/PixelTrackFitting/src/RZLine.h"
00032 #include "RecoPixelVertexing/PixelTrackFitting/src/CircleFromThreePoints.h"
00033 #include "RecoPixelVertexing/PixelTrackFitting/interface/PixelTrackBuilder.h"
00034 
00035 using namespace std;
00036 
00037 /*****************************************************************************/
00038 TrackFitter::TrackFitter
00039   (const edm::ParameterSet& cfg) :
00040    theConfig(cfg), theTracker(0), theField(0), theTTRecHitBuilder(0)
00041 {
00042 }
00043 
00044 /*****************************************************************************/
00045 reco::Track* TrackFitter::run
00046   (const edm::EventSetup& es,
00047    const std::vector<const TrackingRecHit *> & hits,
00048    const TrackingRegion & region) const
00049 {
00050   int nhits = hits.size();
00051   if(nhits <2) return 0;
00052 
00053   vector<GlobalPoint> points;
00054   vector<GlobalError> errors;
00055   vector<bool> isBarrel;
00056   
00057   if (!theField || !theTracker || !theTTRecHitBuilder)
00058   {
00059     edm::ESHandle<TrackerGeometry> trackerESH;
00060     es.get<TrackerDigiGeometryRecord>().get(trackerESH);
00061     theTracker = trackerESH.product();
00062 
00063     edm::ESHandle<MagneticField> fieldESH;
00064     es.get<IdealMagneticFieldRecord>().get(fieldESH);
00065     theField = fieldESH.product();
00066 
00067     edm::ESHandle<TransientTrackingRecHitBuilder> ttrhbESH;
00068     std::string builderName = theConfig.getParameter<std::string>("TTRHBuilder");
00069     es.get<TransientRecHitRecord>().get(builderName,ttrhbESH);
00070     theTTRecHitBuilder = ttrhbESH.product();
00071   }
00072 
00073   for (vector<const TrackingRecHit*>::const_iterator ih = hits.begin();
00074                                                      ih!= hits.end(); ih++)
00075   {
00076     TransientTrackingRecHit::RecHitPointer recHit =
00077       theTTRecHitBuilder->build(*ih);
00078 
00079     points.push_back(recHit->globalPosition());
00080     errors.push_back(recHit->globalPositionError());
00081     isBarrel.push_back(recHit->detUnit()->type().isBarrel() );
00082   }
00083 
00084   CircleFromThreePoints circle = (nhits==2) ?
00085         CircleFromThreePoints(GlobalPoint(0.,0.,0.), points[0], points[1]) :
00086         CircleFromThreePoints(points[0],points[1],points[2]); 
00087 
00088   int charge = getCharge(points);
00089   float curvature = circle.curvature();
00090 
00091   // pt
00092   float invPt = PixelRecoUtilities::inversePt(curvature, es);
00093   float valPt = (invPt > 1.e-4) ? 1./invPt : 1.e4;
00094   float errPt = 0.055*valPt + 0.017*valPt*valPt;
00095 
00096   CircleFromThreePoints::Vector2D center = circle.center();
00097 
00098   // tip
00099   float valTip = charge * (center.mag()-1/curvature);
00100   // zip
00101   float valZip = getZip(valTip, curvature, points[0],points[1]);
00102   // phi
00103   float valPhi = getPhi(center.x(), center.y(), charge);
00104   // cot(theta), update zip
00105   float valCotTheta =
00106      getCotThetaAndUpdateZip(points[0],points[1], 1/curvature,
00107                              valPhi,valTip,valZip);
00108 
00109   // errors
00110   float errTip, errZip;
00111   getErrTipAndErrZip(valPt, points.back().eta(), errTip,errZip);
00112   float errPhi      = 0.002;
00113   float errCotTheta = 0.002;
00114 
00115   float chi2 = 0;
00116   if(nhits > 2)
00117   {
00118     RZLine rzLine(points,errors,isBarrel);
00119     float      cotTheta, intercept, covss, covii, covsi; 
00120     rzLine.fit(cotTheta, intercept, covss, covii, covsi);
00121     chi2 = rzLine.chi2(cotTheta, intercept);
00122     //FIXME: check which intercept to use!
00123   }
00124 
00125   // build pixel track
00126   PixelTrackBuilder builder;
00127 
00128   Measurement1D pt      (valPt,       errPt);
00129   Measurement1D phi     (valPhi,      errPhi);
00130   Measurement1D cotTheta(valCotTheta, errCotTheta);
00131   Measurement1D tip     (valTip,      errTip);
00132   Measurement1D zip     (valZip,      errZip);
00133 
00134   return builder.build(pt, phi, cotTheta, tip, zip, chi2,
00135                        charge, hits, theField);
00136 }
00137 
00138 /*****************************************************************************/
00139 int TrackFitter::getCharge
00140   (const vector<GlobalPoint> & points) const
00141 {
00142    GlobalVector v21 = points[1]-points[0];
00143    GlobalVector v32 = points[2]-points[1];
00144    float dphi = v32.phi() - v21.phi();
00145    while (dphi >  M_PI) dphi -= 2*M_PI;
00146    while (dphi < -M_PI) dphi += 2*M_PI;
00147    return (dphi > 0) ? -1 : 1;
00148 }
00149 
00150 /*****************************************************************************/
00151 float TrackFitter::getCotThetaAndUpdateZip
00152   (const GlobalPoint& inner, const GlobalPoint& outer,
00153    float radius, float phi, float d0, float& zip) const
00154 {
00155    float chi = phi - M_PI_2;
00156    GlobalPoint IP(d0*cos(chi), d0*sin(chi),zip);
00157 
00158    float phi1 = 2*asin(0.5*(inner - IP).perp()/radius);
00159    float phi2 = 2*asin(0.5*(outer - IP).perp()/radius);
00160 
00161    float dr = radius*(phi2 - phi1);
00162    float dz = outer.z()-inner.z();
00163 
00164    // Recalculate ZIP
00165    zip = (inner.z()*phi2 - outer.z()*phi1)/(phi2 - phi1);
00166 
00167    return (fabs(dr) > 1.e-3) ? dz/dr : 0;
00168 }
00169 
00170 /*****************************************************************************/
00171 float TrackFitter::getPhi
00172   (float xC, float yC, int charge) const
00173 {
00174   float phiC;
00175 
00176   if (charge>0) phiC = atan2(xC,-yC);
00177            else phiC = atan2(-xC,yC);
00178 
00179   return phiC;
00180 }
00181 
00182 /*****************************************************************************/
00183 float TrackFitter::getZip
00184   (float d0, float curv,
00185    const GlobalPoint& inner, const GlobalPoint& outer) const
00186 {
00187   // phi = asin(r*rho/2) with asin(x) ~= x+x**3/(2*3)
00188   float rho3 = curv*curv*curv;
00189 
00190   float r1 = inner.perp();
00191   double phi1 = r1*curv/2 + inner.perp2()*r1*rho3/48.;
00192 
00193   float r2 = outer.perp();
00194   double phi2 = r2*curv/2 + outer.perp2()*r2*rho3/48.;
00195 
00196   double z1 = inner.z();
00197   double z2 = outer.z();
00198 
00199   return z1 - phi1/(phi1-phi2)*(z1-z2);
00200 }
00201 
00202 /*****************************************************************************/
00203 void TrackFitter::getErrTipAndErrZip
00204   (float pt, float eta, float & errTip, float & errZip) const 
00205 {
00206   float coshEta = cosh(eta);
00207 
00208   { // transverse
00209     float c_ms = 0.0115; //0.0115;
00210     float s_le = 0.0095; //0.0123;
00211     float s_ms2 = c_ms*c_ms / (pt*pt) * coshEta;
00212 
00213     errTip = sqrt(s_le*s_le + s_ms2                  );
00214   }
00215 
00216   { // z
00217     float c_ms = 0.0070;
00218     float s_le = 0.0135;
00219 
00220     errZip = sqrt( (s_le*s_le + c_ms*c_ms/(pt*pt)) * coshEta*coshEta*coshEta);
00221   }
00222 }
00223