PixelLowPtUtilities/src/ThirdHitPrediction.cc

0001 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0002 #include "MagneticField/Engine/interface/MagneticField.h"
0003 #include "RecoTracker/PixelLowPtUtilities/interface/ThirdHitPrediction.h"
0004 #include "RecoTracker/PixelTrackFitting/interface/CircleFromThreePoints.h"
0005 #include "RecoTracker/TkMSParametrization/interface/MultipleScatteringParametrisationMaker.h"
0006 #include "RecoTracker/TkMSParametrization/interface/MultipleScatteringParametrisation.h"
0007 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0008 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0009 #include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
0010 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0011 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0012
0013 inline float sqr(float x) { return x * x; }
0014
0015 using namespace std;
0016
0017 /*****************************************************************************/
0018 ThirdHitPrediction::ThirdHitPrediction(const TrackingRegion& region,
0019                                        GlobalPoint inner,
0020                                        GlobalPoint outer,
0021                                        const MagneticField& magfield,
0022                                        const TransientTrackingRecHitBuilder& ttrhBuilder,
0023                                        double nSigMultipleScattering,
0024                                        double maxAngleRatio) {
0025   using namespace edm;
0026   theTTRecHitBuilder = &ttrhBuilder;
0027
0028   Bz = fabs(magfield.inInverseGeV(GlobalPoint(0, 0, 0)).z());
0029
0030   c0 = Global2DVector(region.origin().x(), region.origin().y());
0031
0032   r0 = region.originRBound();
0033   rm = region.ptMin() / Bz;
0034
0035   g1 = inner;
0036   g2 = outer;
0037
0038   p1 = Global2DVector(g1.x(), g1.y());
0039   p2 = Global2DVector(g2.x(), g2.y());
0040
0041   dif = p1 - p2;
0042
0043   // Prepare circles of minimal pt (rm) and cylinder of origin (r0)
0044   keep = true;
0045   arc_0m = findArcIntersection(findMinimalCircles(rm), findTouchingCircles(r0), keep);
0046
0047   nSigma = nSigMultipleScattering;
0048   maxRatio = maxAngleRatio;
0049 }
0050
0051 /*****************************************************************************/
0052 ThirdHitPrediction::~ThirdHitPrediction() {}
0053
0054 /*****************************************************************************/
0055 void ThirdHitPrediction::invertCircle(Global2DVector& c, float& r) {
0056   float s = dif.mag2() / ((c - p1).mag2() - sqr(r));
0057
0058   c = p1 + (c - p1) * s;
0059   r *= fabsf(s);
0060 }
0061
0062 /*****************************************************************************/
0063 void ThirdHitPrediction::invertPoint(Global2DVector& p) {
0064   float s = dif.mag2() / (p - p1).mag2();
0065
0066   p = p1 + (p - p1) * s;
0067 }
0068
0069 /*****************************************************************************/
0070 pair<float, float> ThirdHitPrediction::findMinimalCircles(float r) {
0071   pair<float, float> a(0., 0.);
0072
0073   if (dif.mag2() < 2 * sqr(r))
0074     a = pair<float, float>(dif.phi(), 0.5 * acos(1 - 0.5 * dif.mag2() / sqr(r)));
0075
0076   return a;
0077 }
0078
0079 /*****************************************************************************/
0080 pair<float, float> ThirdHitPrediction::findTouchingCircles(float r) {
0081   Global2DVector c = c0;
0082   invertCircle(c, r);
0083
0084   pair<float, float> a(0., 0.);
0085   a = pair<float, float>((c - p2).phi(), 0.5 * acos(1 - 2 * sqr(r) / (c - p2).mag2()));
0086
0087   return a;
0088 }
0089
0090 /*****************************************************************************/
0091 pair<float, float> ThirdHitPrediction::findArcIntersection(pair<float, float> a, pair<float, float> b, bool& keep) {
0092   // spin closer
0093   while (b.first < a.first - M_PI)
0094     b.first += 2 * M_PI;
0095   while (b.first > a.first + M_PI)
0096     b.first -= 2 * M_PI;
0097
0098   float min, max;
0099
0100   if (a.first - a.second > b.first - b.second)
0101     min = a.first - a.second;
0102   else {
0103     min = b.first - b.second;
0104     keep = false;
0105   }
0106
0107   if (a.first + a.second < b.first + b.second)
0108     max = a.first + a.second;
0109   else {
0110     max = b.first + b.second;
0111     keep = false;
0112   }
0113
0114   pair<float, float> c(0., 0.);
0115
0116   if (min < max) {
0117     c.first = 0.5 * (max + min);
0118     c.second = 0.5 * (max - min);
0119   }
0120
0121   return c;
0122 }
0123
0124 /*****************************************************************************/
0125 void ThirdHitPrediction::fitParabola(const float x[3], const float y[3], float par[3]) {
0126   float s2 = sqr(x[0]) * (y[1] - y[2]) + sqr(x[1]) * (y[2] - y[0]) + sqr(x[2]) * (y[0] - y[1]);
0127
0128   float s1 = x[0] * (y[1] - y[2]) + x[1] * (y[2] - y[0]) + x[2] * (y[0] - y[1]);
0129
0130   float s3 = (x[0] - x[1]) * (x[1] - x[2]) * (x[2] - x[0]);
0131   float s4 =
0132       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]);
0133
0134   par[2] = s1 / s3;   // a2
0135   par[1] = -s2 / s3;  // a1
0136   par[0] = -s4 / s3;  // a0
0137 }
0138
0139 /*****************************************************************************/
0140 void ThirdHitPrediction::findRectangle(
0141     const float x[3], const float y[3], const float par[3], float phi[2], float z[2]) {
0142   // Initial guess
0143   phi[0] = min(x[0], x[2]);
0144   z[0] = min(y[0], y[2]);
0145   phi[1] = max(x[0], x[2]);
0146   z[1] = max(y[0], y[2]);
0147
0148   // Extremum: position and value
0149   float xe = -par[1] / (2 * par[2]);
0150   float ye = par[0] - sqr(par[1]) / (4 * par[2]);
0151
0152   // Check if extremum is inside the phi range
0153   if (phi[0] < xe && xe < phi[1]) {
0154     if (ye < z[0])
0155       z[0] = ye;
0156     if (ye > z[1])
0157       z[1] = ye;
0158   }
0159 }
0160
0161 /*****************************************************************************/
0162 float ThirdHitPrediction::areaParallelogram(const Global2DVector& a, const Global2DVector& b) {
0163   return a.x() * b.y() - a.y() * b.x();
0164 }
0165
0166 /*****************************************************************************/
0167 float ThirdHitPrediction::angleRatio(const Global2DVector& p3, const Global2DVector& c) {
0168   float rad2 = (p1 - c).mag2();
0169
0170   float a12 = asin(fabsf(areaParallelogram(p1 - c, p2 - c)) / rad2);
0171   float a23 = asin(fabsf(areaParallelogram(p2 - c, p3 - c)) / rad2);
0172
0173   return a23 / a12;
0174 }
0175
0176 /*****************************************************************************/
0177 void ThirdHitPrediction::spinCloser(float phi[3]) {
0178   while (phi[1] < phi[0] - M_PI)
0179     phi[1] += 2 * M_PI;
0180   while (phi[1] > phi[0] + M_PI)
0181     phi[1] -= 2 * M_PI;
0182
0183   while (phi[2] < phi[1] - M_PI)
0184     phi[2] += 2 * M_PI;
0185   while (phi[2] > phi[1] + M_PI)
0186     phi[2] -= 2 * M_PI;
0187 }
0188
0189 /*****************************************************************************/
0190 void ThirdHitPrediction::calculateRangesBarrel(float r3, float phi[2], float z[2], bool keep) {
0191   pair<float, float> arc_all = findArcIntersection(arc_0m, findTouchingCircles(r3), keep);
0192
0193   if (arc_all.second != 0.) {
0194     Global2DVector c3(0., 0.);  // barrel at r3
0195     invertCircle(c3, r3);       // inverted
0196
0197     float angle[3];  // prepare angles
0198     angle[0] = arc_all.first - arc_all.second;
0199     angle[1] = arc_all.first;
0200     angle[2] = arc_all.first + arc_all.second;
0201
0202     float phi3[3], z3[3];
0203     Global2DVector delta = c3 - p2;
0204
0205     for (int i = 0; i < 3; i++) {
0206       Global2DVector vec(cos(angle[i]), sin(angle[i]));  // unit vector
0207       float lambda = delta * vec - sqrt(sqr(delta * vec) - delta * delta + sqr(r3));
0208
0209       Global2DVector p3 = p2 + lambda * vec;  // inverted third hit
0210       invertPoint(p3);                        // third hit
0211       phi3[i] = p3.phi();                     // phi of third hit
0212
0213       float ratio;
0214
0215       if (keep && i == 1) {  // Straight line
0216         ratio = (p2 - p3).mag() / (p1 - p2).mag();
0217       } else {                                            // Circle
0218         Global2DVector c = p2 - vec * (vec * (p2 - p1));  // inverted antipodal
0219         invertPoint(c);                                   // antipodal
0220         c = 0.5 * (p1 + c);                               // center
0221
0222         ratio = angleRatio(p3, c);
0223       }
0224
0225       z3[i] = g2.z() + (g2.z() - g1.z()) * ratio;  // z of third hit
0226     }
0227
0228     spinCloser(phi3);
0229
0230     // Parabola on phi - z
0231     float par[3];
0232     fitParabola(phi3, z3, par);
0233     findRectangle(phi3, z3, par, phi, z);
0234   }
0235 }
0236
0237 /*****************************************************************************/
0238 void ThirdHitPrediction::calculateRangesForward(float z3, float phi[2], float r[2], bool keep) {
0239   float angle[3];  // prepare angles
0240   angle[0] = arc_0m.first - arc_0m.second;
0241   angle[1] = arc_0m.first;
0242   angle[2] = arc_0m.first + arc_0m.second;
0243
0244   float ratio = (z3 - g2.z()) / (g2.z() - g1.z());
0245
0246   if (0 < ratio && ratio < maxRatio) {
0247     float phi3[3], r3[3];
0248
0249     for (int i = 0; i < 3; i++) {
0250       Global2DVector p3;
0251
0252       if (keep && i == 1) {  // Straight line
0253         p3 = p2 + ratio * (p2 - p1);
0254       } else {                                             // Circle
0255         Global2DVector vec(cos(angle[i]), sin(angle[i]));  // unit vector
0256
0257         Global2DVector c = p2 - vec * (vec * (p2 - p1));  // inverted antipodal
0258         invertPoint(c);                                   // antipodal
0259         c = 0.5 * (p1 + c);                               // center
0260
0261         float rad2 = (p1 - c).mag2();
0262
0263         float a12 = asin(areaParallelogram(p1 - c, p2 - c) / rad2);
0264         float a23 = ratio * a12;
0265
0266         p3 = c + Global2DVector((p2 - c).x() * cos(a23) - (p2 - c).y() * sin(a23),
0267                                 (p2 - c).x() * sin(a23) + (p2 - c).y() * cos(a23));
0268       }
0269
0270       phi3[i] = p3.phi();
0271       r3[i] = p3.mag();
0272     }
0273
0274     spinCloser(phi3);
0275
0276     // Parabola on phi - z
0277     float par[3];
0278     fitParabola(phi3, r3, par);
0279     findRectangle(phi3, r3, par, phi, r);
0280   }
0281 }
0282
0283 /*****************************************************************************/
0284 void ThirdHitPrediction::calculateRanges(float rz3, float phi[2], float rz[2]) {
0285   // Clear
0286   phi[0] = 0.;
0287   rz[0] = 0.;
0288   phi[1] = 0.;
0289   rz[1] = 0.;
0290
0291   // Calculate
0292   if (theBarrel)
0293     calculateRangesBarrel(rz3, phi, rz, keep);
0294   else
0295     calculateRangesForward(rz3, phi, rz, keep);
0296 }
0297
0298 /*****************************************************************************/
0299 void ThirdHitPrediction::getRanges(const DetLayer* layer, float phi[], float rz[]) {
0300   theLayer = layer;
0301
0302   if (layer)
0303     initLayer(layer);
0304
0305   float phi_inner[2], rz_inner[2];
0306   calculateRanges(theDetRange.min(), phi_inner, rz_inner);
0307
0308   float phi_outer[2], rz_outer[2];
0309   calculateRanges(theDetRange.max(), phi_outer, rz_outer);
0310
0311   if ((phi_inner[0] == 0. && phi_inner[1] == 0.) || (phi_outer[0] == 0. && phi_outer[1] == 0.)) {
0312     phi[0] = 0.;
0313     phi[1] = 0.;
0314
0315     rz[0] = 0.;
0316     rz[1] = 0.;
0317   } else {
0318     while (phi_outer[0] > phi_inner[0] + M_PI) {
0319       phi_outer[0] -= 2 * M_PI;
0320       phi_outer[1] -= 2 * M_PI;
0321     }
0322
0323     while (phi_outer[0] < phi_inner[0] - M_PI) {
0324       phi_outer[0] += 2 * M_PI;
0325       phi_outer[1] += 2 * M_PI;
0326     }
0327
0328     phi[0] = min(phi_inner[0], phi_outer[0]);
0329     phi[1] = max(phi_inner[1], phi_outer[1]);
0330
0331     rz[0] = min(rz_inner[0], rz_outer[0]);
0332     rz[1] = max(rz_inner[1], rz_outer[1]);
0333   }
0334 }
0335
0336 /*****************************************************************************/
0337 void ThirdHitPrediction::getRanges(float rz3, float phi[], float rz[]) { calculateRanges(rz3, phi, rz); }
0338
0339 /*****************************************************************************/
0340 bool ThirdHitPrediction::isCompatibleWithMultipleScattering(GlobalPoint g3,
0341                                                             const vector<const TrackingRecHit*>& h,
0342                                                             vector<GlobalVector>& globalDirs,
0343                                                             const MultipleScatteringParametrisationMaker& msmaker) {
0344   Global2DVector p1(g1.x(), g1.y());
0345   Global2DVector p2(g2.x(), g2.y());
0346   Global2DVector p3(g3.x(), g3.y());
0347
0348   CircleFromThreePoints circle(g1, g2, g3);
0349
0350   if (circle.curvature() != 0.) {
0351     Global2DVector c(circle.center().x(), circle.center().y());
0352
0353     float rad2 = (p1 - c).mag2();
0354     float a12 = asin(fabsf(areaParallelogram(p1 - c, p2 - c)) / rad2);
0355     float a23 = asin(fabsf(areaParallelogram(p2 - c, p3 - c)) / rad2);
0356
0357     float slope = (g2.z() - g1.z()) / a12;
0358
0359     float rz3 = g2.z() + slope * a23;
0360     float delta_z = g3.z() - rz3;
0361
0362     // Transform to tt
0363     vector<TransientTrackingRecHit::RecHitPointer> th;
0364     for (vector<const TrackingRecHit*>::const_iterator ih = h.begin(); ih != h.end(); ih++)
0365       th.push_back(theTTRecHitBuilder->build(*ih));
0366
0367     float sigma1_le2 = max(th[0]->parametersError()[0][0], th[0]->parametersError()[1][1]);
0368     float sigma2_le2 = max(th[1]->parametersError()[0][0], th[1]->parametersError()[1][1]);
0369
0370     float sigma_z2 = (1 + a23 / a12) * (1 + a23 / a12) * sigma2_le2 + (a23 / a12) * (a23 / a12) * sigma1_le2;
0371
0372     float cotTheta = slope * circle.curvature();  // == sinhEta
0373     float coshEta = sqrt(1 + sqr(cotTheta));      // == 1/sinTheta
0374
0375     float pt = Bz / circle.curvature();
0376     float p = pt * coshEta;
0377
0378     float m_pi = 0.13957018;
0379     float beta = p / sqrt(sqr(p) + sqr(m_pi));
0380
0381     MultipleScatteringParametrisation msp = msmaker.parametrisation(theLayer);
0382     PixelRecoPointRZ rz2(g2.perp(), g2.z());
0383
0384     float sigma_z = msp(pt, cotTheta, rz2) / beta;
0385
0386     // Calculate globalDirs
0387     float sinTheta = 1. / coshEta;
0388     float cosTheta = cotTheta * sinTheta;
0389
0390     int dir;
0391     if (areaParallelogram(p1 - c, p2 - c) > 0)
0392       dir = 1;
0393     else
0394       dir = -1;
0395
0396     float curvature = circle.curvature();
0397
0398     {
0399       Global2DVector v = (p1 - c) * curvature * dir;
0400       globalDirs.push_back(GlobalVector(-v.y() * sinTheta, v.x() * sinTheta, cosTheta));
0401     }
0402
0403     {
0404       Global2DVector v = (p2 - c) * curvature * dir;
0405       globalDirs.push_back(GlobalVector(-v.y() * sinTheta, v.x() * sinTheta, cosTheta));
0406     }
0407
0408     {
0409       Global2DVector v = (p3 - c) * curvature * dir;
0410       globalDirs.push_back(GlobalVector(-v.y() * sinTheta, v.x() * sinTheta, cosTheta));
0411     }
0412
0413     // Multiple scattering
0414     float sigma_ms = sigma_z * coshEta;
0415
0416     // Local error squared
0417     float sigma_le2 = max(th[2]->parametersError()[0][0], th[2]->parametersError()[1][1]);
0418
0419     return (delta_z * delta_z / (sigma_ms * sigma_ms + sigma_le2 + sigma_z2) < nSigma * nSigma);
0420   }
0421
0422   return false;
0423 }
0424
0425 /*****************************************************************************/
0426 void ThirdHitPrediction::initLayer(const DetLayer* layer) {
0427   if (layer->location() == GeomDetEnumerators::barrel) {
0428     theBarrel = true;
0429     theForward = false;
0430     const BarrelDetLayer& bl = dynamic_cast<const BarrelDetLayer&>(*layer);
0431     float halfThickness = bl.surface().bounds().thickness() / 2;
0432     float radius = bl.specificSurface().radius();
0433     theDetRange = Range(radius - halfThickness, radius + halfThickness);
0434   } else if (layer->location() == GeomDetEnumerators::endcap) {
0435     theBarrel = false;
0436     theForward = true;
0437     const ForwardDetLayer& fl = dynamic_cast<const ForwardDetLayer&>(*layer);
0438     float halfThickness = fl.surface().bounds().thickness() / 2;
0439     float zLayer = fl.position().z();
0440     theDetRange = Range(zLayer - halfThickness, zLayer + halfThickness);
0441   }
0442 }