PixelFitterByConformalMappingAndLine.cc

Go to the documentation of this file.

#include "RecoPixelVertexing/PixelTrackFitting/interface/PixelFitterByConformalMappingAndLine.h"

#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"

#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "DataFormats/GeometryCommonDetAlgo/interface/GlobalError.h"


#include "MagneticField/Engine/interface/MagneticField.h"
#include "CommonTools/Statistics/interface/LinearFit.h"
#include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"


#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "TrackingTools/DetLayers/interface/DetLayer.h"
#include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"

#include "Geometry/CommonDetUnit/interface/GeomDet.h"
#include "Geometry/CommonDetUnit/interface/GeomDetType.h"

#include "ConformalMappingFit.h"
#include "RecoPixelVertexing/PixelTrackFitting/interface/PixelTrackBuilder.h"
#include "RecoPixelVertexing/PixelTrackFitting/interface/RZLine.h"

#include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
#include "RecoTracker/TkMSParametrization/interface/LongitudinalBendingCorrection.h"

using namespace std;

template <class T> T sqr( T t) {return t*t;}


PixelFitterByConformalMappingAndLine::PixelFitterByConformalMappingAndLine(const edm::EventSetup *es,
                                                                           const TransientTrackingRecHitBuilder *ttrhBuilder,
                                                                           const TrackerGeometry *tracker,
                                                                           const MagneticField *field,
                                                                           double fixImpactParameter,
                                                                           bool useFixImpactParameter):
  theES(es),
  theTTRHBuilder(ttrhBuilder),
  theTracker(tracker),
  theField(field),
  theFixImpactParameter(fixImpactParameter),
  theUseFixImpactParameter(useFixImpactParameter)
{ }

std::unique_ptr<reco::Track> PixelFitterByConformalMappingAndLine::run(
    const std::vector<const TrackingRecHit * > & hits,
    const TrackingRegion & region) const
{

  int nhits = hits.size();

  vector<GlobalPoint> points;
  vector<GlobalError> errors;
  vector<bool> isBarrel;
  

  for (vector<const TrackingRecHit*>::const_iterator ih=hits.begin();  ih!=hits.end(); ih++) {
    TransientTrackingRecHit::RecHitPointer recHit = theTTRHBuilder->build(*ih);
    points.push_back( recHit->globalPosition() );
    errors.push_back( recHit->globalPositionError() );
    isBarrel.push_back( recHit->detUnit()->type().isBarrel() );
  }

//    if (useMultScatt) {
//      MultipleScatteringParametrisation ms(hits[i].layer());
//      float cotTheta = (p.z()-zVtx)/p.perp();
//      err += sqr( ms( pt, cotTheta, PixelRecoPointRZ(0.,zVtx) ) );
//    }

  //
  // simple fit to get pt, phi0 used for precise calcul.
  //
  typedef ConformalMappingFit::PointXY PointXY;
  vector<PointXY> xy; vector<float> errRPhi2;
  for (int i=0; i < nhits; ++i) {
    const GlobalPoint & point = points[i];
    xy.push_back(PointXY( point.x()-region.origin().x(), point.y()-region.origin().y()));
    float phiErr2 = errors[i].phierr(point);
    errRPhi2.push_back( point.perp2()*phiErr2);
  }
  ConformalMappingFit parabola(xy, errRPhi2);
  if (theUseFixImpactParameter) 
    parabola.fixImpactParmaeter(theFixImpactParameter);
  else if (nhits < 3) parabola.fixImpactParmaeter(0.);
 

  Measurement1D curv = parabola.curvature();
  float invPt = PixelRecoUtilities::inversePt( curv.value(), *theES);
  float valPt =  (invPt > 1.e-4) ? 1./invPt : 1.e4;
  float errPt =PixelRecoUtilities::inversePt(curv.error(), *theES) * sqr(valPt);
  Measurement1D pt (valPt,errPt);
  Measurement1D phi = parabola.directionPhi();
  Measurement1D tip = parabola.impactParameter();

  //
  // precalculate theta to correct errors:
  //
  vector<float> r(nhits),z(nhits),errZ(nhits);
  float simpleCot = ( points.back().z()-points.front().z() )/ (points.back().perp() - points.front().perp() );
  for (int i=0; i< nhits; ++i) {
    const GlobalPoint & point = points[i]; 
    const GlobalError & error = errors[i];
    r[i] = sqrt( sqr(point.x()-region.origin().x()) + sqr(point.y()-region.origin().y()) );
    r[i] += pixelrecoutilities::LongitudinalBendingCorrection(pt.value(), *theES)(r[i]);
    z[i] = point.z()-region.origin().z();  
    errZ[i] =  (isBarrel[i]) ? sqrt(error.czz()) : sqrt( error.rerr(point) )*simpleCot;
  }

  //
  // line fit (R-Z plane)
  //
  RZLine rzLine(r,z,errZ);
  
//
// parameters for track builder 
//
  Measurement1D zip(rzLine.intercept(), sqrt(rzLine.covii()));
  Measurement1D cotTheta(rzLine.cotTheta(), sqrt(rzLine.covss()));
  float chi2 = parabola.chi2() +  rzLine.chi2();
  int charge = parabola.charge();


  PixelTrackBuilder builder;
  return std::unique_ptr<reco::Track>(builder.build(pt, phi, cotTheta, tip, zip, chi2, charge, hits,  theField, region.origin()));
}