#include <ThirdHitPrediction.h>

Public Types
typedef TkTrackingRegionsMargin< float >	Margin
typedef PixelRecoRange< float >	Range
Public Member Functions
void	getRanges (const DetLayer *layer, float phi[], float rz[])
void	getRanges (float rORz, float phi[], float rz[])
bool	isCompatibleWithMultipleScattering (GlobalPoint g3, std::vector< const TrackingRecHit * > h, std::vector< GlobalVector > &localDirs, const edm::EventSetup &es)
	ThirdHitPrediction (const TrackingRegion &region, GlobalPoint inner, GlobalPoint outer, const edm::EventSetup &es, double nSigMultipleScattering, double maxAngleRatio, std::string builderName)
	~ThirdHitPrediction ()
Private Member Functions
float	angleRatio (const Global2DVector &p3, const Global2DVector &c)
float	areaParallelogram (const Global2DVector &a, const Global2DVector &b)
void	calculateRanges (float rz3, float phi[2], float rz[2])
void	calculateRangesBarrel (float r3, float phi[2], float z[2], bool keep)
void	calculateRangesForward (float z3, float phi[2], float r[2], bool keep)
std::pair< float, float >	findArcIntersection (std::pair< float, float > a, std::pair< float, float > b, bool &keep)
std::pair< float, float >	findMinimalCircles (float r)
void	findRectangle (const float x[3], const float y[3], const float par[3], float phi[2], float z[2])
std::pair< float, float >	findTouchingCircles (float r)
void	fitParabola (const float x[3], const float y[3], float par[3])
void	initLayer (const DetLayer *layer)
void	invertCircle (Global2DVector &c, float &r)
void	invertPoint (Global2DVector &p)
void	printOut (char *text)
void	spinCloser (float phi[3])
Private Attributes
std::pair< float, float >	arc_0m
float	Bz
Global2DVector	c0
Global2DVector	dif
GlobalPoint	g1
GlobalPoint	g2
bool	keep
double	maxRatio
double	nSigma
Global2DVector	p1
Global2DVector	p2
float	r0
float	rm
bool	theBarrel
Range	theDetRange
bool	theForward
const DetLayer *	theLayer
PixelRecoLineRZ	theLine
Margin	theTolerance
const TransientTrackingRecHitBuilder *	theTTRecHitBuilder

Detailed Description

Definition at line 33 of file ThirdHitPrediction.h.

Member Typedef Documentation

typedef TkTrackingRegionsMargin<float> ThirdHitPrediction::Margin

Definition at line 36 of file ThirdHitPrediction.h.

typedef PixelRecoRange<float> ThirdHitPrediction::Range

Definition at line 35 of file ThirdHitPrediction.h.

Constructor & Destructor Documentation

ThirdHitPrediction::ThirdHitPrediction	(	const TrackingRegion &	region,
		GlobalPoint	inner,
		GlobalPoint	outer,
		const edm::EventSetup &	es,
		double	nSigMultipleScattering,
		double	maxAngleRatio,
		std::string	builderName
	)

Definition at line 24 of file ThirdHitPrediction.cc.

References edm::EventSetup::get(), SurfaceOrientation::inner, keep, autoMagneticFieldProducer_cfi::magfield, TrackingRegion::origin(), TrackingRegion::originRBound(), SurfaceOrientation::outer, p1, p2, TrackingRegion::ptMin(), submit::rm, PV3DBase< T, PVType, FrameType >::x(), and PV3DBase< T, PVType, FrameType >::y().

{
 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 ( )

Definition at line 63 of file ThirdHitPrediction.cc.

{
}

Member Function Documentation

float ThirdHitPrediction::angleRatio	(	const Global2DVector &	p3,
		const Global2DVector &	c
	)		`[private]`

Definition at line 192 of file ThirdHitPrediction.cc.

References trackerHits::c, mag2(), p1, and p2.

{
  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;
}

float ThirdHitPrediction::areaParallelogram	(	const Global2DVector &	a,
		const Global2DVector &	b
	)		`[private]`

Definition at line 185 of file ThirdHitPrediction.cc.

References PV2DBase< T, PVType, FrameType >::x(), and PV2DBase< T, PVType, FrameType >::y().

{
  return a.x() * b.y() - a.y() * b.x();
}

void ThirdHitPrediction::calculateRanges	(	float	rz3,
		float	phi[2],
		float	rz[2]
	)		`[private]`

Definition at line 323 of file ThirdHitPrediction.cc.

{
  // 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::calculateRangesBarrel	(	float	r3,
		float	phi[2],
		float	z[2],
		bool	keep
	)		`[private]`

Definition at line 214 of file ThirdHitPrediction.cc.

References angle(), trackerHits::c, funct::cos(), delta, i, mag(), p1, p2, p3, PV2DBase< T, PVType, FrameType >::phi(), funct::sin(), funct::sqr(), and mathSSE::sqrt().

{
  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
	)		`[private]`

Definition at line 270 of file ThirdHitPrediction.cc.

References angle(), trackerHits::c, funct::cos(), i, PV2DBase< T, PVType, FrameType >::mag(), mag2(), p1, p2, p3, PV2DBase< T, PVType, FrameType >::phi(), funct::sin(), x, and detailsBasic3DVector::y.

{
  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);
  }
}

pair< float, float > ThirdHitPrediction::findArcIntersection	(	std::pair< float, float >	a,
		std::pair< float, float >	b,
		bool &	keep
	)		`[private]`

Definition at line 111 of file ThirdHitPrediction.cc.

References trackerHits::c, M_PI, max(), and min.

{ 
  // 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;
}

pair< float, float > ThirdHitPrediction::findMinimalCircles ( float r ) [private]

Definition at line 85 of file ThirdHitPrediction.cc.

References a, and funct::sqr().

{ 
  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;
}

void ThirdHitPrediction::findRectangle	(	const float	x[3],
		const float	y[3],
		const float	par[3],
		float	phi[2],
		float	z[2]
	)		`[private]`

Definition at line 164 of file ThirdHitPrediction.cc.

References max(), min, and funct::sqr().

{ 
  // 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;
  }
}

pair< float, float > ThirdHitPrediction::findTouchingCircles ( float r ) [private]

Definition at line 97 of file ThirdHitPrediction.cc.

References a, trackerHits::c, mag2(), p2, phi, and funct::sqr().

{
  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;
}

void ThirdHitPrediction::fitParabola	(	const float	x[3],
		const float	y[3],
		float	par[3]
	)		`[private]`

Definition at line 142 of file ThirdHitPrediction.cc.

References indexGen::s2, and funct::sqr().

{ 
  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::getRanges	(	const DetLayer *	layer,
		float	phi[],
		float	rz[]
	)

Definition at line 336 of file ThirdHitPrediction.cc.

References M_PI, max(), and min.

Referenced by PixelTripletLowPtGenerator::hitTriplets().

{
  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	rORz,
		float	phi[],
		float	rz[]
	)

Definition at line 375 of file ThirdHitPrediction.cc.

{
  calculateRanges(rz3, phi,rz);
}

void ThirdHitPrediction::initLayer ( const DetLayer * layer ) [private]

Definition at line 470 of file ThirdHitPrediction.cc.

References Reference_intrackfit_cff::barrel, BoundSurface::bounds(), Reference_intrackfit_cff::endcap, DetLayer::location(), GeometricSearchDet::position(), CosmicsPD_Skims::radius, Cylinder::radius(), BarrelDetLayer::specificSurface(), BarrelDetLayer::surface(), ForwardDetLayer::surface(), Bounds::thickness(), and PV3DBase< T, PVType, FrameType >::z().

{
  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);
  }
}

void ThirdHitPrediction::invertCircle	(	Global2DVector &	c,
		float &	r
	)		`[private]`

Definition at line 68 of file ThirdHitPrediction.cc.

References mag2(), p1, asciidump::s, and funct::sqr().

{
  float s = dif.mag2() / ((c - p1).mag2() - sqr(r));

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

void ThirdHitPrediction::invertPoint ( Global2DVector & p ) [private]

Definition at line 77 of file ThirdHitPrediction.cc.

References mag2(), p1, and asciidump::s.

{
  float s = dif.mag2() / (p - p1).mag2();

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

bool ThirdHitPrediction::isCompatibleWithMultipleScattering	(	GlobalPoint	g3,
		std::vector< const TrackingRecHit * >	h,
		std::vector< GlobalVector > &	localDirs,
		const edm::EventSetup &	es
	)

Definition at line 382 of file ThirdHitPrediction.cc.

References beta, trackerHits::c, PixelRecoUtilities::curvature(), dir, m_pi, mag2(), max(), L1TEmulatorMonitor_cff::p, p1, p2, p3, slope, funct::sqr(), mathSSE::sqrt(), v, PV3DBase< T, PVType, FrameType >::x(), PV2DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), PV2DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by PixelTripletLowPtGenerator::hitTriplets().

{
  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*>::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::printOut ( char * text ) [private]

void ThirdHitPrediction::spinCloser ( float phi[3] ) [private]

Definition at line 203 of file ThirdHitPrediction.cc.

References M_PI.

{
  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;
}