ThirdHitPrediction.cc

Go to the documentation of this file.

#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "MagneticField/Engine/interface/MagneticField.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "RecoPixelVertexing/PixelLowPtUtilities/interface/ThirdHitPrediction.h"
#include "RecoPixelVertexing/PixelTrackFitting/interface/CircleFromThreePoints.h"
#include "RecoTracker/TkMSParametrization/interface/MultipleScatteringParametrisation.h"
#include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
#include "TrackingTools/DetLayers/interface/DetLayer.h"
#include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
#include "TrackingTools/Records/interface/TransientRecHitRecord.h"
#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"

inline float sqr(float x) { return x*x; }

using namespace std;

/*****************************************************************************/
ThirdHitPrediction::ThirdHitPrediction
  (const TrackingRegion & region, GlobalPoint inner, GlobalPoint outer,
   const edm::EventSetup& es,
   double nSigMultipleScattering, double maxAngleRatio,
   string builderName)
{
 using namespace edm;
 ESHandle<MagneticField> magfield;
 es.get<IdealMagneticFieldRecord>().get(magfield);

 edm::ESHandle<TransientTrackingRecHitBuilder> ttrhbESH;
 es.get<TransientRecHitRecord>().get(builderName,ttrhbESH);
 theTTRecHitBuilder = ttrhbESH.product();

 Bz = fabs(magfield->inInverseGeV(GlobalPoint(0,0,0)).z());

 c0 = Global2DVector(region.origin().x(),
                     region.origin().y());

 r0 = region.originRBound();
 rm = region.ptMin() / Bz;

 g1 = inner;
 g2 = outer;

 p1 = Global2DVector(g1.x(), g1.y());
 p2 = Global2DVector(g2.x(), g2.y());

 dif = p1 - p2;

 // Prepare circles of minimal pt (rm) and cylinder of origin (r0)
 keep = true;
 arc_0m = findArcIntersection(findMinimalCircles (rm),
                              findTouchingCircles(r0), keep);

 nSigma   = nSigMultipleScattering;
 maxRatio = maxAngleRatio;
}

/*****************************************************************************/
ThirdHitPrediction::~ThirdHitPrediction()
{
}

/*****************************************************************************/
void ThirdHitPrediction::invertCircle(Global2DVector& c,float& r)
{
  float s = dif.mag2() / ((c - p1).mag2() - sqr(r));

  c = p1 + (c - p1)*s;
  r *= fabsf(s);
}

/*****************************************************************************/
void ThirdHitPrediction::invertPoint(Global2DVector& p)
{
  float s = dif.mag2() / (p - p1).mag2();

  p = p1 + (p - p1)*s;
}

/*****************************************************************************/
pair<float,float> ThirdHitPrediction::findMinimalCircles(float r)
{ 
  pair<float,float> a(0.,0.);

  if(dif.mag2() <  2 * sqr(r))
    a = pair<float,float>( dif.phi(),
                           0.5*acos(1 - 0.5 * dif.mag2()/sqr(r)) );

  return a;
}

/*****************************************************************************/
pair<float,float> ThirdHitPrediction::findTouchingCircles(float r)
{
  Global2DVector c = c0;
  invertCircle(c,r);

  pair<float,float> a(0.,0.);
  a = pair<float,float>( (c - p2).phi(),
                          0.5*acos(1 - 2*sqr(r)/(c - p2).mag2()) );

  return a;
}

/*****************************************************************************/
pair<float,float> ThirdHitPrediction::findArcIntersection
  (pair<float,float> a, pair<float,float> b, bool& keep)
{ 
  // spin closer
  while(b.first < a.first - M_PI) b.first += 2*M_PI;
  while(b.first > a.first + M_PI) b.first -= 2*M_PI;

  float min,max;

  if(a.first - a.second > b.first - b.second)
    min = a.first - a.second;
  else
  { min = b.first - b.second; keep = false; }

  if(a.first + a.second < b.first + b.second)
    max = a.first + a.second;
  else
  { max = b.first + b.second; keep = false; }
  
  pair<float,float> c(0.,0.);

  if(min < max)
  {
    c.first  = 0.5*(max + min);
    c.second = 0.5*(max - min);
  }

  return c;
}

/*****************************************************************************/
void ThirdHitPrediction::fitParabola
  (const float x[3], const float y[3], float par[3])
{ 
  float s2 = sqr(x[0]) * (y[1] - y[2]) +
             sqr(x[1]) * (y[2] - y[0]) +
             sqr(x[2]) * (y[0] - y[1]);
  
  float s1 =     x[0]  * (y[1] - y[2]) +
                 x[1]  * (y[2] - y[0]) +
                 x[2]  * (y[0] - y[1]);

  float s3 = (x[0] - x[1]) * (x[1] - x[2]) * (x[2] - x[0]);
  float s4 = x[0]*x[1]*y[2] * (x[0] - x[1]) +
             x[0]*y[1]*x[2] * (x[2] - x[0]) +
             y[0]*x[1]*x[2] * (x[1] - x[2]);

  par[2] =  s1 / s3; // a2
  par[1] = -s2 / s3; // a1
  par[0] = -s4 / s3; // a0
}

/*****************************************************************************/
void ThirdHitPrediction::findRectangle
  (const float x[3], const float y[3], const float par[3],
   float phi[2],float z[2])
{ 
  // Initial guess
  phi[0] = min(x[0],x[2]); z[0] = min(y[0],y[2]);
  phi[1] = max(x[0],x[2]); z[1] = max(y[0],y[2]);

  // Extremum: position and value
  float xe = -par[1]/(2*par[2]);
  float ye = par[0] - sqr(par[1])/(4*par[2]);

  // Check if extremum is inside the phi range
  if(phi[0] < xe  && xe < phi[1])
  {
    if(ye < z[0]) z[0] = ye;
    if(ye > z[1]) z[1] = ye;
  }
}

/*****************************************************************************/
float ThirdHitPrediction::areaParallelogram
  (const Global2DVector& a, const Global2DVector& b)
{
  return a.x() * b.y() - a.y() * b.x();
}

/*****************************************************************************/
float ThirdHitPrediction::angleRatio
  (const Global2DVector& p3, const Global2DVector& c)
{
  float rad2 = (p1 - c).mag2();

  float a12 = asin(fabsf(areaParallelogram(p1 - c, p2 - c)) / rad2);
  float a23 = asin(fabsf(areaParallelogram(p2 - c, p3 - c)) / rad2);

  return a23/a12;
}

/*****************************************************************************/
void ThirdHitPrediction::spinCloser(float phi[3])
{
  while(phi[1] < phi[0] - M_PI) phi[1] += 2*M_PI;
  while(phi[1] > phi[0] + M_PI) phi[1] -= 2*M_PI;

  while(phi[2] < phi[1] - M_PI) phi[2] += 2*M_PI;
  while(phi[2] > phi[1] + M_PI) phi[2] -= 2*M_PI;
}

/*****************************************************************************/
void ThirdHitPrediction::calculateRangesBarrel
  (float r3, float phi[2],float z[2], bool keep)
{
  pair<float,float> arc_all =
    findArcIntersection(arc_0m, findTouchingCircles(r3), keep);

  if(arc_all.second != 0.)
  {
    Global2DVector c3(0.,0.); // barrel at r3
    invertCircle(c3,r3);      // inverted

    float angle[3];           // prepare angles
    angle[0] = arc_all.first - arc_all.second;
    angle[1] = arc_all.first;
    angle[2] = arc_all.first + arc_all.second;

    float phi3[3], z3[3];
    Global2DVector delta = c3 - p2;

    for(int i=0; i<3; i++)
    {
      Global2DVector vec(cos(angle[i]), sin(angle[i])); // unit vector
      float lambda = delta*vec - sqrt(sqr(delta*vec) - delta*delta + sqr(r3));

      Global2DVector p3 = p2 + lambda * vec;  // inverted third hit 
      invertPoint(p3);                        // third hit
      phi3[i] = p3.phi();                     // phi of third hit

      float ratio;

      if(keep && i==1)
      { // Straight line
        ratio = (p2 - p3).mag() / (p1 - p2).mag();
      }
      else
      { // Circle
        Global2DVector c = p2 - vec * (vec * (p2 - p1)); // inverted antipodal
        invertPoint(c);                                  // antipodal
        c = 0.5*(p1 + c);                                // center

        ratio = angleRatio(p3,c);
      }

      z3[i] = g2.z() + (g2.z() - g1.z()) * ratio;        // z of third hit
    }

    spinCloser(phi3);

    // Parabola on phi - z
    float par[3];
    fitParabola  (phi3,z3, par);
    findRectangle(phi3,z3, par, phi,z);
  }
}

/*****************************************************************************/
void ThirdHitPrediction::calculateRangesForward
  (float z3, float phi[2],float r[2], bool keep)
{
  float angle[3];           // prepare angles
  angle[0] = arc_0m.first - arc_0m.second;
  angle[1] = arc_0m.first;
  angle[2] = arc_0m.first + arc_0m.second;

  float ratio = (z3 - g2.z()) / (g2.z() - g1.z());

  if(0 < ratio  && ratio < maxRatio)
  {
    float phi3[3], r3[3];

    for(int i=0; i<3; i++)
    {
      Global2DVector p3;

      if(keep && i==1)
      { // Straight line
        p3 = p2 + ratio * (p2 - p1);
      }
      else
      { // Circle
        Global2DVector vec(cos(angle[i]), sin(angle[i])); // unit vector
  
        Global2DVector c = p2 - vec * (vec * (p2 - p1));  // inverted antipodal
        invertPoint(c);                                   // antipodal
        c = 0.5*(p1 + c);                                 // center

        float rad2 = (p1 - c).mag2();

        float a12 = asin(areaParallelogram(p1 - c, p2 - c) / rad2);
        float a23 = ratio * a12;

        p3 = c + Global2DVector((p2-c).x()*cos(a23) - (p2-c).y()*sin(a23),
                                (p2-c).x()*sin(a23) + (p2-c).y()*cos(a23));
      }

      phi3[i] = p3.phi();
      r3[i]   = p3.mag();
    }

    spinCloser(phi3);

    // Parabola on phi - z
    float par[3];
    fitParabola  (phi3,r3, par);
    findRectangle(phi3,r3, par, phi,r);
  }
}

/*****************************************************************************/
void ThirdHitPrediction::calculateRanges
  (float rz3, float phi[2],float rz[2])
{
  // Clear
  phi[0] = 0.; rz[0] = 0.;
  phi[1] = 0.; rz[1] = 0.;

  // Calculate
  if(theBarrel) calculateRangesBarrel (rz3, phi,rz, keep);
           else calculateRangesForward(rz3, phi,rz, keep);
}

/*****************************************************************************/
void ThirdHitPrediction::getRanges
  (const DetLayer *layer, float phi[],float rz[])
{
  theLayer = layer;
  
  if (layer) initLayer(layer);

  float phi_inner[2],rz_inner[2];
  calculateRanges(theDetRange.min(), phi_inner,rz_inner);

  float phi_outer[2],rz_outer[2];
  calculateRanges(theDetRange.max(), phi_outer,rz_outer);

  if( (phi_inner[0] == 0. && phi_inner[1] == 0.) ||
      (phi_outer[0] == 0. && phi_outer[1] == 0.) )
  {
    phi[0] = 0.;
    phi[1] = 0.;

     rz[0] = 0.;
     rz[1] = 0.;
  }
  else
  {
    while(phi_outer[0] > phi_inner[0] + M_PI)
    { phi_outer[0] -= 2*M_PI; phi_outer[1] -= 2*M_PI; }

    while(phi_outer[0] < phi_inner[0] - M_PI)
    { phi_outer[0] += 2*M_PI; phi_outer[1] += 2*M_PI; }

    phi[0] = min(phi_inner[0],phi_outer[0]);
    phi[1] = max(phi_inner[1],phi_outer[1]);
  
     rz[0] = min( rz_inner[0], rz_outer[0]);
     rz[1] = max( rz_inner[1], rz_outer[1]);
  }
}

/*****************************************************************************/
void ThirdHitPrediction::getRanges
  (float rz3, float phi[],float rz[])
{
  calculateRanges(rz3, phi,rz);
}

/*****************************************************************************/
bool ThirdHitPrediction::isCompatibleWithMultipleScattering
  (GlobalPoint g3, const vector<const TrackingRecHit*>& h,
   vector<GlobalVector>& globalDirs, const edm::EventSetup& es)
{
  Global2DVector p1(g1.x(),g1.y());
  Global2DVector p2(g2.x(),g2.y());
  Global2DVector p3(g3.x(),g3.y());

  CircleFromThreePoints circle(g1,g2,g3);

  if(circle.curvature() != 0.)
  {
  Global2DVector c (circle.center().x(), circle.center().y());

  float rad2 = (p1 - c).mag2();
  float a12 = asin(fabsf(areaParallelogram(p1 - c, p2 - c)) / rad2);
  float a23 = asin(fabsf(areaParallelogram(p2 - c, p3 - c)) / rad2);

  float slope = (g2.z() - g1.z()) / a12;

  float rz3 = g2.z() + slope * a23;
  float delta_z = g3.z() - rz3;

  // Transform to tt
  vector<TransientTrackingRecHit::RecHitPointer> th;
  for(vector<const TrackingRecHit*>::const_iterator ih = h.begin(); ih!= h.end(); ih++)
    th.push_back(theTTRecHitBuilder->build(*ih));

  float sigma1_le2 = max(th[0]->parametersError()[0][0],
                         th[0]->parametersError()[1][1]);
  float sigma2_le2 = max(th[1]->parametersError()[0][0],
                         th[1]->parametersError()[1][1]);

  float sigma_z2 = (1 + a23/a12)*(1 + a23/a12) * sigma2_le2 +
                   (    a23/a12)*(    a23/a12) * sigma1_le2;

  float cotTheta = slope * circle.curvature(); // == sinhEta
  float coshEta  = sqrt(1 + sqr(cotTheta));    // == 1/sinTheta

  float pt = Bz / circle.curvature();
  float p  = pt * coshEta;

  float m_pi = 0.13957018;
  float beta = p / sqrt(sqr(p) + sqr(m_pi));

  MultipleScatteringParametrisation msp(theLayer,es);
  PixelRecoPointRZ rz2(g2.perp(), g2.z());

  float sigma_z = msp(pt, cotTheta, rz2) / beta;

  // Calculate globalDirs
  float sinTheta =       1. / coshEta;
  float cosTheta = cotTheta * sinTheta;

  int dir;
  if(areaParallelogram(p1 - c, p2 - c) > 0) dir = 1; else dir = -1;

  float curvature = circle.curvature();

  {
   Global2DVector v = (p1 - c)*curvature*dir;
   globalDirs.push_back(GlobalVector(-v.y()*sinTheta,v.x()*sinTheta,cosTheta));
  }

  {
   Global2DVector v = (p2 - c)*curvature*dir;
   globalDirs.push_back(GlobalVector(-v.y()*sinTheta,v.x()*sinTheta,cosTheta));
  }

  {
   Global2DVector v = (p3 - c)*curvature*dir;
   globalDirs.push_back(GlobalVector(-v.y()*sinTheta,v.x()*sinTheta,cosTheta));
  }

  // Multiple scattering
  float sigma_ms  = sigma_z * coshEta;

  // Local error squared
  float sigma_le2 = max(th[2]->parametersError()[0][0],
                        th[2]->parametersError()[1][1]);

  return (delta_z*delta_z / (sigma_ms*sigma_ms + sigma_le2 + sigma_z2)
          < nSigma * nSigma);
  }

  return false;
}

/*****************************************************************************/
void ThirdHitPrediction::initLayer(const DetLayer *layer)
{
  if ( layer->location() == GeomDetEnumerators::barrel) {
    theBarrel = true;
    theForward = false;
    const BarrelDetLayer& bl = dynamic_cast<const BarrelDetLayer&>(*layer);
    float halfThickness  = bl.surface().bounds().thickness()/2;
    float radius = bl.specificSurface().radius();
    theDetRange = Range(radius-halfThickness, radius+halfThickness);
  } else if ( layer->location() == GeomDetEnumerators::endcap) {
    theBarrel= false;
    theForward = true;
    const ForwardDetLayer& fl = dynamic_cast<const ForwardDetLayer&>(*layer);
    float halfThickness  = fl.surface().bounds().thickness()/2;
    float zLayer = fl.position().z() ;
    theDetRange = Range(zLayer-halfThickness, zLayer+halfThickness);
  }
}