CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
ThirdHitPrediction.cc
Go to the documentation of this file.
12 
13 inline float sqr(float x) { return x * x; }
14 
15 using namespace std;
16 
17 /*****************************************************************************/
21  const MagneticField& magfield,
22  const TransientTrackingRecHitBuilder& ttrhBuilder,
24  double maxAngleRatio) {
25  using namespace edm;
26  theTTRecHitBuilder = &ttrhBuilder;
27 
28  Bz = fabs(magfield.inInverseGeV(GlobalPoint(0, 0, 0)).z());
29 
30  c0 = Global2DVector(region.origin().x(), region.origin().y());
31 
32  r0 = region.originRBound();
33  rm = region.ptMin() / Bz;
34 
35  g1 = inner;
36  g2 = outer;
37 
38  p1 = Global2DVector(g1.x(), g1.y());
39  p2 = Global2DVector(g2.x(), g2.y());
40 
41  dif = p1 - p2;
42 
43  // Prepare circles of minimal pt (rm) and cylinder of origin (r0)
44  keep = true;
45  arc_0m = findArcIntersection(findMinimalCircles(rm), findTouchingCircles(r0), keep);
46 
48  maxRatio = maxAngleRatio;
49 }
50 
51 /*****************************************************************************/
53 
54 /*****************************************************************************/
56  float s = dif.mag2() / ((c - p1).mag2() - sqr(r));
57 
58  c = p1 + (c - p1) * s;
59  r *= fabsf(s);
60 }
61 
62 /*****************************************************************************/
64  float s = dif.mag2() / (p - p1).mag2();
65 
66  p = p1 + (p - p1) * s;
67 }
68 
69 /*****************************************************************************/
70 pair<float, float> ThirdHitPrediction::findMinimalCircles(float r) {
71  pair<float, float> a(0., 0.);
72 
73  if (dif.mag2() < 2 * sqr(r))
74  a = pair<float, float>(dif.phi(), 0.5 * acos(1 - 0.5 * dif.mag2() / sqr(r)));
75 
76  return a;
77 }
78 
79 /*****************************************************************************/
80 pair<float, float> ThirdHitPrediction::findTouchingCircles(float r) {
82  invertCircle(c, r);
83 
84  pair<float, float> a(0., 0.);
85  a = pair<float, float>((c - p2).phi(), 0.5 * acos(1 - 2 * sqr(r) / (c - p2).mag2()));
86 
87  return a;
88 }
89 
90 /*****************************************************************************/
91 pair<float, float> ThirdHitPrediction::findArcIntersection(pair<float, float> a, pair<float, float> b, bool& keep) {
92  // spin closer
93  while (b.first < a.first - M_PI)
94  b.first += 2 * M_PI;
95  while (b.first > a.first + M_PI)
96  b.first -= 2 * M_PI;
97 
98  float min, max;
99 
100  if (a.first - a.second > b.first - b.second)
101  min = a.first - a.second;
102  else {
103  min = b.first - b.second;
104  keep = false;
105  }
106 
107  if (a.first + a.second < b.first + b.second)
108  max = a.first + a.second;
109  else {
110  max = b.first + b.second;
111  keep = false;
112  }
113 
114  pair<float, float> c(0., 0.);
115 
116  if (min < max) {
117  c.first = 0.5 * (max + min);
118  c.second = 0.5 * (max - min);
119  }
120 
121  return c;
122 }
123 
124 /*****************************************************************************/
125 void ThirdHitPrediction::fitParabola(const float x[3], const float y[3], float par[3]) {
126  float s2 = sqr(x[0]) * (y[1] - y[2]) + sqr(x[1]) * (y[2] - y[0]) + sqr(x[2]) * (y[0] - y[1]);
127 
128  float s1 = x[0] * (y[1] - y[2]) + x[1] * (y[2] - y[0]) + x[2] * (y[0] - y[1]);
129 
130  float s3 = (x[0] - x[1]) * (x[1] - x[2]) * (x[2] - x[0]);
131  float s4 =
132  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]);
133 
134  par[2] = s1 / s3; // a2
135  par[1] = -s2 / s3; // a1
136  par[0] = -s4 / s3; // a0
137 }
138 
139 /*****************************************************************************/
141  const float x[3], const float y[3], const float par[3], float phi[2], float z[2]) {
142  // Initial guess
143  phi[0] = min(x[0], x[2]);
144  z[0] = min(y[0], y[2]);
145  phi[1] = max(x[0], x[2]);
146  z[1] = max(y[0], y[2]);
147 
148  // Extremum: position and value
149  float xe = -par[1] / (2 * par[2]);
150  float ye = par[0] - sqr(par[1]) / (4 * par[2]);
151 
152  // Check if extremum is inside the phi range
153  if (phi[0] < xe && xe < phi[1]) {
154  if (ye < z[0])
155  z[0] = ye;
156  if (ye > z[1])
157  z[1] = ye;
158  }
159 }
160 
161 /*****************************************************************************/
163  return a.x() * b.y() - a.y() * b.x();
164 }
165 
166 /*****************************************************************************/
168  float rad2 = (p1 - c).mag2();
169 
170  float a12 = asin(fabsf(areaParallelogram(p1 - c, p2 - c)) / rad2);
171  float a23 = asin(fabsf(areaParallelogram(p2 - c, p3 - c)) / rad2);
172 
173  return a23 / a12;
174 }
175 
176 /*****************************************************************************/
177 void ThirdHitPrediction::spinCloser(float phi[3]) {
178  while (phi[1] < phi[0] - M_PI)
179  phi[1] += 2 * M_PI;
180  while (phi[1] > phi[0] + M_PI)
181  phi[1] -= 2 * M_PI;
182 
183  while (phi[2] < phi[1] - M_PI)
184  phi[2] += 2 * M_PI;
185  while (phi[2] > phi[1] + M_PI)
186  phi[2] -= 2 * M_PI;
187 }
188 
189 /*****************************************************************************/
190 void ThirdHitPrediction::calculateRangesBarrel(float r3, float phi[2], float z[2], bool keep) {
191  pair<float, float> arc_all = findArcIntersection(arc_0m, findTouchingCircles(r3), keep);
192 
193  if (arc_all.second != 0.) {
194  Global2DVector c3(0., 0.); // barrel at r3
195  invertCircle(c3, r3); // inverted
196 
197  float angle[3]; // prepare angles
198  angle[0] = arc_all.first - arc_all.second;
199  angle[1] = arc_all.first;
200  angle[2] = arc_all.first + arc_all.second;
201 
202  float phi3[3], z3[3];
203  Global2DVector delta = c3 - p2;
204 
205  for (int i = 0; i < 3; i++) {
206  Global2DVector vec(cos(angle[i]), sin(angle[i])); // unit vector
207  float lambda = delta * vec - sqrt(sqr(delta * vec) - delta * delta + sqr(r3));
208 
209  Global2DVector p3 = p2 + lambda * vec; // inverted third hit
210  invertPoint(p3); // third hit
211  phi3[i] = p3.phi(); // phi of third hit
212 
213  float ratio;
214 
215  if (keep && i == 1) { // Straight line
216  ratio = (p2 - p3).mag() / (p1 - p2).mag();
217  } else { // Circle
218  Global2DVector c = p2 - vec * (vec * (p2 - p1)); // inverted antipodal
219  invertPoint(c); // antipodal
220  c = 0.5 * (p1 + c); // center
221 
222  ratio = angleRatio(p3, c);
223  }
224 
225  z3[i] = g2.z() + (g2.z() - g1.z()) * ratio; // z of third hit
226  }
227 
228  spinCloser(phi3);
229 
230  // Parabola on phi - z
231  float par[3];
232  fitParabola(phi3, z3, par);
233  findRectangle(phi3, z3, par, phi, z);
234  }
235 }
236 
237 /*****************************************************************************/
238 void ThirdHitPrediction::calculateRangesForward(float z3, float phi[2], float r[2], bool keep) {
239  float angle[3]; // prepare angles
240  angle[0] = arc_0m.first - arc_0m.second;
241  angle[1] = arc_0m.first;
242  angle[2] = arc_0m.first + arc_0m.second;
243 
244  float ratio = (z3 - g2.z()) / (g2.z() - g1.z());
245 
246  if (0 < ratio && ratio < maxRatio) {
247  float phi3[3], r3[3];
248 
249  for (int i = 0; i < 3; i++) {
251 
252  if (keep && i == 1) { // Straight line
253  p3 = p2 + ratio * (p2 - p1);
254  } else { // Circle
255  Global2DVector vec(cos(angle[i]), sin(angle[i])); // unit vector
256 
257  Global2DVector c = p2 - vec * (vec * (p2 - p1)); // inverted antipodal
258  invertPoint(c); // antipodal
259  c = 0.5 * (p1 + c); // center
260 
261  float rad2 = (p1 - c).mag2();
262 
263  float a12 = asin(areaParallelogram(p1 - c, p2 - c) / rad2);
264  float a23 = ratio * a12;
265 
266  p3 = c + Global2DVector((p2 - c).x() * cos(a23) - (p2 - c).y() * sin(a23),
267  (p2 - c).x() * sin(a23) + (p2 - c).y() * cos(a23));
268  }
269 
270  phi3[i] = p3.phi();
271  r3[i] = p3.mag();
272  }
273 
274  spinCloser(phi3);
275 
276  // Parabola on phi - z
277  float par[3];
278  fitParabola(phi3, r3, par);
279  findRectangle(phi3, r3, par, phi, r);
280  }
281 }
282 
283 /*****************************************************************************/
284 void ThirdHitPrediction::calculateRanges(float rz3, float phi[2], float rz[2]) {
285  // Clear
286  phi[0] = 0.;
287  rz[0] = 0.;
288  phi[1] = 0.;
289  rz[1] = 0.;
290 
291  // Calculate
292  if (theBarrel)
293  calculateRangesBarrel(rz3, phi, rz, keep);
294  else
295  calculateRangesForward(rz3, phi, rz, keep);
296 }
297 
298 /*****************************************************************************/
299 void ThirdHitPrediction::getRanges(const DetLayer* layer, float phi[], float rz[]) {
300  theLayer = layer;
301 
302  if (layer)
303  initLayer(layer);
304 
305  float phi_inner[2], rz_inner[2];
306  calculateRanges(theDetRange.min(), phi_inner, rz_inner);
307 
308  float phi_outer[2], rz_outer[2];
309  calculateRanges(theDetRange.max(), phi_outer, rz_outer);
310 
311  if ((phi_inner[0] == 0. && phi_inner[1] == 0.) || (phi_outer[0] == 0. && phi_outer[1] == 0.)) {
312  phi[0] = 0.;
313  phi[1] = 0.;
314 
315  rz[0] = 0.;
316  rz[1] = 0.;
317  } else {
318  while (phi_outer[0] > phi_inner[0] + M_PI) {
319  phi_outer[0] -= 2 * M_PI;
320  phi_outer[1] -= 2 * M_PI;
321  }
322 
323  while (phi_outer[0] < phi_inner[0] - M_PI) {
324  phi_outer[0] += 2 * M_PI;
325  phi_outer[1] += 2 * M_PI;
326  }
327 
328  phi[0] = min(phi_inner[0], phi_outer[0]);
329  phi[1] = max(phi_inner[1], phi_outer[1]);
330 
331  rz[0] = min(rz_inner[0], rz_outer[0]);
332  rz[1] = max(rz_inner[1], rz_outer[1]);
333  }
334 }
335 
336 /*****************************************************************************/
337 void ThirdHitPrediction::getRanges(float rz3, float phi[], float rz[]) { calculateRanges(rz3, phi, rz); }
338 
339 /*****************************************************************************/
341  const vector<const TrackingRecHit*>& h,
342  vector<GlobalVector>& globalDirs,
344  Global2DVector p1(g1.x(), g1.y());
345  Global2DVector p2(g2.x(), g2.y());
346  Global2DVector p3(g3.x(), g3.y());
347 
348  CircleFromThreePoints circle(g1, g2, g3);
349 
350  if (circle.curvature() != 0.) {
351  Global2DVector c(circle.center().x(), circle.center().y());
352 
353  float rad2 = (p1 - c).mag2();
354  float a12 = asin(fabsf(areaParallelogram(p1 - c, p2 - c)) / rad2);
355  float a23 = asin(fabsf(areaParallelogram(p2 - c, p3 - c)) / rad2);
356 
357  float slope = (g2.z() - g1.z()) / a12;
358 
359  float rz3 = g2.z() + slope * a23;
360  float delta_z = g3.z() - rz3;
361 
362  // Transform to tt
363  vector<TransientTrackingRecHit::RecHitPointer> th;
364  for (vector<const TrackingRecHit*>::const_iterator ih = h.begin(); ih != h.end(); ih++)
365  th.push_back(theTTRecHitBuilder->build(*ih));
366 
367  float sigma1_le2 = max(th[0]->parametersError()[0][0], th[0]->parametersError()[1][1]);
368  float sigma2_le2 = max(th[1]->parametersError()[0][0], th[1]->parametersError()[1][1]);
369 
370  float sigma_z2 = (1 + a23 / a12) * (1 + a23 / a12) * sigma2_le2 + (a23 / a12) * (a23 / a12) * sigma1_le2;
371 
372  float cotTheta = slope * circle.curvature(); // == sinhEta
373  float coshEta = sqrt(1 + sqr(cotTheta)); // == 1/sinTheta
374 
375  float pt = Bz / circle.curvature();
376  float p = pt * coshEta;
377 
378  float m_pi = 0.13957018;
379  float beta = p / sqrt(sqr(p) + sqr(m_pi));
380 
381  MultipleScatteringParametrisation msp = msmaker.parametrisation(theLayer);
382  PixelRecoPointRZ rz2(g2.perp(), g2.z());
383 
384  float sigma_z = msp(pt, cotTheta, rz2) / beta;
385 
386  // Calculate globalDirs
387  float sinTheta = 1. / coshEta;
388  float cosTheta = cotTheta * sinTheta;
389 
390  int dir;
391  if (areaParallelogram(p1 - c, p2 - c) > 0)
392  dir = 1;
393  else
394  dir = -1;
395 
396  float curvature = circle.curvature();
397 
398  {
399  Global2DVector v = (p1 - c) * curvature * dir;
400  globalDirs.push_back(GlobalVector(-v.y() * sinTheta, v.x() * sinTheta, cosTheta));
401  }
402 
403  {
404  Global2DVector v = (p2 - c) * curvature * dir;
405  globalDirs.push_back(GlobalVector(-v.y() * sinTheta, v.x() * sinTheta, cosTheta));
406  }
407 
408  {
409  Global2DVector v = (p3 - c) * curvature * dir;
410  globalDirs.push_back(GlobalVector(-v.y() * sinTheta, v.x() * sinTheta, cosTheta));
411  }
412 
413  // Multiple scattering
414  float sigma_ms = sigma_z * coshEta;
415 
416  // Local error squared
417  float sigma_le2 = max(th[2]->parametersError()[0][0], th[2]->parametersError()[1][1]);
418 
419  return (delta_z * delta_z / (sigma_ms * sigma_ms + sigma_le2 + sigma_z2) < nSigma * nSigma);
420  }
421 
422  return false;
423 }
424 
425 /*****************************************************************************/
427  if (layer->location() == GeomDetEnumerators::barrel) {
428  theBarrel = true;
429  theForward = false;
430  const BarrelDetLayer& bl = dynamic_cast<const BarrelDetLayer&>(*layer);
431  float halfThickness = bl.surface().bounds().thickness() / 2;
432  float radius = bl.specificSurface().radius();
433  theDetRange = Range(radius - halfThickness, radius + halfThickness);
434  } else if (layer->location() == GeomDetEnumerators::endcap) {
435  theBarrel = false;
436  theForward = true;
437  const ForwardDetLayer& fl = dynamic_cast<const ForwardDetLayer&>(*layer);
438  float halfThickness = fl.surface().bounds().thickness() / 2;
439  float zLayer = fl.position().z();
440  theDetRange = Range(zLayer - halfThickness, zLayer + halfThickness);
441  }
442 }
void findRectangle(const float x[3], const float y[3], const float par[3], float phi[2], float z[2])
def rm(path, rec=False)
Definition: eostools.py:363
virtual const Surface::PositionType & position() const
Returns position of the surface.
void invertCircle(Global2DVector &c, float &r)
void initLayer(const DetLayer *layer)
Vector2DBase< float, GlobalTag > Global2DVector
const BoundSurface & surface() const final
GeometricSearchDet interface.
const BoundSurface & surface() const final
The surface of the GeometricSearchDet.
T z() const
Definition: PV3DBase.h:61
PixelRecoRange< float > Range
static const double slope[3]
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
void calculateRangesForward(float z3, float phi[2], float r[2], bool keep)
void calculateRangesBarrel(float r3, float phi[2], float z[2], bool keep)
float angleRatio(const Global2DVector &p3, const Global2DVector &c)
T curvature(T InversePt, const MagneticField &field)
constexpr std::array< uint8_t, layerIndexSize > layer
virtual float thickness() const =0
T x() const
Definition: PV3DBase.h:59
T y() const
Definition: PV3DBase.h:60
MultipleScatteringParametrisation parametrisation(const DetLayer *layer, X0Source x0Source=X0Source::useX0AtEta) const
T sqrt(T t)
Definition: SSEVec.h:19
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
float areaParallelogram(const Global2DVector &a, const Global2DVector &b)
std::pair< float, float > findMinimalCircles(float r)
#define M_PI
ThirdHitPrediction(const TrackingRegion &region, GlobalPoint inner, GlobalPoint outer, const MagneticField &magfield, const TransientTrackingRecHitBuilder &ttrhBuilder, double nSigMultipleScattering, double maxAngleRatio)
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
T mag2() const
The vector magnitude squared. Equivalent to vec.dot(vec)
virtual const BoundCylinder & specificSurface() const final
Extension of the interface.
double b
Definition: hdecay.h:118
float sigma_z
void spinCloser(float phi[3])
HLT enums.
std::pair< float, float > findArcIntersection(std::pair< float, float > a, std::pair< float, float > b, bool &keep)
double a
Definition: hdecay.h:119
float x
void invertPoint(Global2DVector &p)
void fitParabola(const float x[3], const float y[3], float par[3])
void calculateRanges(float rz3, float phi[2], float rz[2])
void getRanges(const DetLayer *layer, float phi[], float rz[])
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
std::pair< float, float > findTouchingCircles(float r)
bool isCompatibleWithMultipleScattering(GlobalPoint g3, const std::vector< const TrackingRecHit *> &h, std::vector< GlobalVector > &localDirs, const MultipleScatteringParametrisationMaker &msmaker)
Global3DVector GlobalVector
Definition: GlobalVector.h:10
#define m_pi
Definition: RPCConst.cc:8
T angle(T x1, T y1, T z1, T x2, T y2, T z2)
Definition: angle.h:11
const Bounds & bounds() const
Definition: Surface.h:87
float sqr(float x)