File indexing completed on 2024-04-06 12:28:30
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
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
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;
0135 par[1] = -s2 / s3;
0136 par[0] = -s4 / s3;
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
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
0149 float xe = -par[1] / (2 * par[2]);
0150 float ye = par[0] - sqr(par[1]) / (4 * par[2]);
0151
0152
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.);
0195 invertCircle(c3, r3);
0196
0197 float angle[3];
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]));
0207 float lambda = delta * vec - sqrt(sqr(delta * vec) - delta * delta + sqr(r3));
0208
0209 Global2DVector p3 = p2 + lambda * vec;
0210 invertPoint(p3);
0211 phi3[i] = p3.phi();
0212
0213 float ratio;
0214
0215 if (keep && i == 1) {
0216 ratio = (p2 - p3).mag() / (p1 - p2).mag();
0217 } else {
0218 Global2DVector c = p2 - vec * (vec * (p2 - p1));
0219 invertPoint(c);
0220 c = 0.5 * (p1 + c);
0221
0222 ratio = angleRatio(p3, c);
0223 }
0224
0225 z3[i] = g2.z() + (g2.z() - g1.z()) * ratio;
0226 }
0227
0228 spinCloser(phi3);
0229
0230
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];
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) {
0253 p3 = p2 + ratio * (p2 - p1);
0254 } else {
0255 Global2DVector vec(cos(angle[i]), sin(angle[i]));
0256
0257 Global2DVector c = p2 - vec * (vec * (p2 - p1));
0258 invertPoint(c);
0259 c = 0.5 * (p1 + c);
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
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
0286 phi[0] = 0.;
0287 rz[0] = 0.;
0288 phi[1] = 0.;
0289 rz[1] = 0.;
0290
0291
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
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();
0373 float coshEta = sqrt(1 + sqr(cotTheta));
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
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
0414 float sigma_ms = sigma_z * coshEta;
0415
0416
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 }