CMS 3D CMS Logo

ThirdHitPrediction.cc
Go to the documentation of this file.
2 
4 
8 
11 
14 
17 
18 inline float sqr(float x) { return x*x; }
19 
20 using namespace std;
21 
22 /*****************************************************************************/
25  const edm::EventSetup& es,
27  string builderName)
28 {
29  using namespace edm;
30  ESHandle<MagneticField> magfield;
31  es.get<IdealMagneticFieldRecord>().get(magfield);
32 
34  es.get<TransientRecHitRecord>().get(builderName,ttrhbESH);
35  theTTRecHitBuilder = ttrhbESH.product();
36 
37  Bz = fabs(magfield->inInverseGeV(GlobalPoint(0,0,0)).z());
38 
39  c0 = Global2DVector(region.origin().x(),
40  region.origin().y());
41 
42  r0 = region.originRBound();
43  rm = region.ptMin() / Bz;
44 
45  g1 = inner;
46  g2 = outer;
47 
48  p1 = Global2DVector(g1.x(), g1.y());
49  p2 = Global2DVector(g2.x(), g2.y());
50 
51  dif = p1 - p2;
52 
53  // Prepare circles of minimal pt (rm) and cylinder of origin (r0)
54  keep = true;
55  arc_0m = findArcIntersection(findMinimalCircles (rm),
56  findTouchingCircles(r0), keep);
57 
59  maxRatio = maxAngleRatio;
60 }
61 
62 /*****************************************************************************/
64 {
65 }
66 
67 /*****************************************************************************/
69 {
70  float s = dif.mag2() / ((c - p1).mag2() - sqr(r));
71 
72  c = p1 + (c - p1)*s;
73  r *= fabsf(s);
74 }
75 
76 /*****************************************************************************/
78 {
79  float s = dif.mag2() / (p - p1).mag2();
80 
81  p = p1 + (p - p1)*s;
82 }
83 
84 /*****************************************************************************/
85 pair<float,float> ThirdHitPrediction::findMinimalCircles(float r)
86 {
87  pair<float,float> a(0.,0.);
88 
89  if(dif.mag2() < 2 * sqr(r))
90  a = pair<float,float>( dif.phi(),
91  0.5*acos(1 - 0.5 * dif.mag2()/sqr(r)) );
92 
93  return a;
94 }
95 
96 /*****************************************************************************/
98 {
100  invertCircle(c,r);
101 
102  pair<float,float> a(0.,0.);
103  a = pair<float,float>( (c - p2).phi(),
104  0.5*acos(1 - 2*sqr(r)/(c - p2).mag2()) );
105 
106  return a;
107 }
108 
109 /*****************************************************************************/
111  (pair<float,float> a, pair<float,float> b, bool& keep)
112 {
113  // spin closer
114  while(b.first < a.first - M_PI) b.first += 2*M_PI;
115  while(b.first > a.first + M_PI) b.first -= 2*M_PI;
116 
117  float min,max;
118 
119  if(a.first - a.second > b.first - b.second)
120  min = a.first - a.second;
121  else
122  { min = b.first - b.second; keep = false; }
123 
124  if(a.first + a.second < b.first + b.second)
125  max = a.first + a.second;
126  else
127  { max = b.first + b.second; keep = false; }
128 
129  pair<float,float> c(0.,0.);
130 
131  if(min < max)
132  {
133  c.first = 0.5*(max + min);
134  c.second = 0.5*(max - min);
135  }
136 
137  return c;
138 }
139 
140 /*****************************************************************************/
142  (const float x[3], const float y[3], float par[3])
143 {
144  float s2 = sqr(x[0]) * (y[1] - y[2]) +
145  sqr(x[1]) * (y[2] - y[0]) +
146  sqr(x[2]) * (y[0] - y[1]);
147 
148  float s1 = x[0] * (y[1] - y[2]) +
149  x[1] * (y[2] - y[0]) +
150  x[2] * (y[0] - y[1]);
151 
152  float s3 = (x[0] - x[1]) * (x[1] - x[2]) * (x[2] - x[0]);
153  float s4 = x[0]*x[1]*y[2] * (x[0] - x[1]) +
154  x[0]*y[1]*x[2] * (x[2] - x[0]) +
155  y[0]*x[1]*x[2] * (x[1] - x[2]);
156 
157  par[2] = s1 / s3; // a2
158  par[1] = -s2 / s3; // a1
159  par[0] = -s4 / s3; // a0
160 }
161 
162 /*****************************************************************************/
164  (const float x[3], const float y[3], const float par[3],
165  float phi[2],float z[2])
166 {
167  // Initial guess
168  phi[0] = min(x[0],x[2]); z[0] = min(y[0],y[2]);
169  phi[1] = max(x[0],x[2]); z[1] = max(y[0],y[2]);
170 
171  // Extremum: position and value
172  float xe = -par[1]/(2*par[2]);
173  float ye = par[0] - sqr(par[1])/(4*par[2]);
174 
175  // Check if extremum is inside the phi range
176  if(phi[0] < xe && xe < phi[1])
177  {
178  if(ye < z[0]) z[0] = ye;
179  if(ye > z[1]) z[1] = ye;
180  }
181 }
182 
183 /*****************************************************************************/
185  (const Global2DVector& a, const Global2DVector& b)
186 {
187  return a.x() * b.y() - a.y() * b.x();
188 }
189 
190 /*****************************************************************************/
193 {
194  float rad2 = (p1 - c).mag2();
195 
196  float a12 = asin(fabsf(areaParallelogram(p1 - c, p2 - c)) / rad2);
197  float a23 = asin(fabsf(areaParallelogram(p2 - c, p3 - c)) / rad2);
198 
199  return a23/a12;
200 }
201 
202 /*****************************************************************************/
204 {
205  while(phi[1] < phi[0] - M_PI) phi[1] += 2*M_PI;
206  while(phi[1] > phi[0] + M_PI) phi[1] -= 2*M_PI;
207 
208  while(phi[2] < phi[1] - M_PI) phi[2] += 2*M_PI;
209  while(phi[2] > phi[1] + M_PI) phi[2] -= 2*M_PI;
210 }
211 
212 /*****************************************************************************/
214  (float r3, float phi[2],float z[2], bool keep)
215 {
216  pair<float,float> arc_all =
217  findArcIntersection(arc_0m, findTouchingCircles(r3), keep);
218 
219  if(arc_all.second != 0.)
220  {
221  Global2DVector c3(0.,0.); // barrel at r3
222  invertCircle(c3,r3); // inverted
223 
224  float angle[3]; // prepare angles
225  angle[0] = arc_all.first - arc_all.second;
226  angle[1] = arc_all.first;
227  angle[2] = arc_all.first + arc_all.second;
228 
229  float phi3[3], z3[3];
230  Global2DVector delta = c3 - p2;
231 
232  for(int i=0; i<3; i++)
233  {
234  Global2DVector vec(cos(angle[i]), sin(angle[i])); // unit vector
235  float lambda = delta*vec - sqrt(sqr(delta*vec) - delta*delta + sqr(r3));
236 
237  Global2DVector p3 = p2 + lambda * vec; // inverted third hit
238  invertPoint(p3); // third hit
239  phi3[i] = p3.phi(); // phi of third hit
240 
241  float ratio;
242 
243  if(keep && i==1)
244  { // Straight line
245  ratio = (p2 - p3).mag() / (p1 - p2).mag();
246  }
247  else
248  { // Circle
249  Global2DVector c = p2 - vec * (vec * (p2 - p1)); // inverted antipodal
250  invertPoint(c); // antipodal
251  c = 0.5*(p1 + c); // center
252 
253  ratio = angleRatio(p3,c);
254  }
255 
256  z3[i] = g2.z() + (g2.z() - g1.z()) * ratio; // z of third hit
257  }
258 
259  spinCloser(phi3);
260 
261  // Parabola on phi - z
262  float par[3];
263  fitParabola (phi3,z3, par);
264  findRectangle(phi3,z3, par, phi,z);
265  }
266 }
267 
268 /*****************************************************************************/
270  (float z3, float phi[2],float r[2], bool keep)
271 {
272  float angle[3]; // prepare angles
273  angle[0] = arc_0m.first - arc_0m.second;
274  angle[1] = arc_0m.first;
275  angle[2] = arc_0m.first + arc_0m.second;
276 
277  float ratio = (z3 - g2.z()) / (g2.z() - g1.z());
278 
279  if(0 < ratio && ratio < maxRatio)
280  {
281  float phi3[3], r3[3];
282 
283  for(int i=0; i<3; i++)
284  {
286 
287  if(keep && i==1)
288  { // Straight line
289  p3 = p2 + ratio * (p2 - p1);
290  }
291  else
292  { // Circle
293  Global2DVector vec(cos(angle[i]), sin(angle[i])); // unit vector
294 
295  Global2DVector c = p2 - vec * (vec * (p2 - p1)); // inverted antipodal
296  invertPoint(c); // antipodal
297  c = 0.5*(p1 + c); // center
298 
299  float rad2 = (p1 - c).mag2();
300 
301  float a12 = asin(areaParallelogram(p1 - c, p2 - c) / rad2);
302  float a23 = ratio * a12;
303 
304  p3 = c + Global2DVector((p2-c).x()*cos(a23) - (p2-c).y()*sin(a23),
305  (p2-c).x()*sin(a23) + (p2-c).y()*cos(a23));
306  }
307 
308  phi3[i] = p3.phi();
309  r3[i] = p3.mag();
310  }
311 
312  spinCloser(phi3);
313 
314  // Parabola on phi - z
315  float par[3];
316  fitParabola (phi3,r3, par);
317  findRectangle(phi3,r3, par, phi,r);
318  }
319 }
320 
321 /*****************************************************************************/
323  (float rz3, float phi[2],float rz[2])
324 {
325  // Clear
326  phi[0] = 0.; rz[0] = 0.;
327  phi[1] = 0.; rz[1] = 0.;
328 
329  // Calculate
330  if(theBarrel) calculateRangesBarrel (rz3, phi,rz, keep);
331  else calculateRangesForward(rz3, phi,rz, keep);
332 }
333 
334 /*****************************************************************************/
336  (const DetLayer *layer, float phi[],float rz[])
337 {
338  theLayer = layer;
339 
340  if (layer) initLayer(layer);
341 
342  float phi_inner[2],rz_inner[2];
343  calculateRanges(theDetRange.min(), phi_inner,rz_inner);
344 
345  float phi_outer[2],rz_outer[2];
346  calculateRanges(theDetRange.max(), phi_outer,rz_outer);
347 
348  if( (phi_inner[0] == 0. && phi_inner[1] == 0.) ||
349  (phi_outer[0] == 0. && phi_outer[1] == 0.) )
350  {
351  phi[0] = 0.;
352  phi[1] = 0.;
353 
354  rz[0] = 0.;
355  rz[1] = 0.;
356  }
357  else
358  {
359  while(phi_outer[0] > phi_inner[0] + M_PI)
360  { phi_outer[0] -= 2*M_PI; phi_outer[1] -= 2*M_PI; }
361 
362  while(phi_outer[0] < phi_inner[0] - M_PI)
363  { phi_outer[0] += 2*M_PI; phi_outer[1] += 2*M_PI; }
364 
365  phi[0] = min(phi_inner[0],phi_outer[0]);
366  phi[1] = max(phi_inner[1],phi_outer[1]);
367 
368  rz[0] = min( rz_inner[0], rz_outer[0]);
369  rz[1] = max( rz_inner[1], rz_outer[1]);
370  }
371 }
372 
373 /*****************************************************************************/
375  (float rz3, float phi[],float rz[])
376 {
377  calculateRanges(rz3, phi,rz);
378 }
379 
380 /*****************************************************************************/
382  (GlobalPoint g3, const vector<const TrackingRecHit*>& h,
383  vector<GlobalVector>& globalDirs, const edm::EventSetup& es)
384 {
385  Global2DVector p1(g1.x(),g1.y());
386  Global2DVector p2(g2.x(),g2.y());
387  Global2DVector p3(g3.x(),g3.y());
388 
389  CircleFromThreePoints circle(g1,g2,g3);
390 
391  if(circle.curvature() != 0.)
392  {
393  Global2DVector c (circle.center().x(), circle.center().y());
394 
395  float rad2 = (p1 - c).mag2();
396  float a12 = asin(fabsf(areaParallelogram(p1 - c, p2 - c)) / rad2);
397  float a23 = asin(fabsf(areaParallelogram(p2 - c, p3 - c)) / rad2);
398 
399  float slope = (g2.z() - g1.z()) / a12;
400 
401  float rz3 = g2.z() + slope * a23;
402  float delta_z = g3.z() - rz3;
403 
404  // Transform to tt
405  vector<TransientTrackingRecHit::RecHitPointer> th;
406  for(vector<const TrackingRecHit*>::const_iterator ih = h.begin(); ih!= h.end(); ih++)
407  th.push_back(theTTRecHitBuilder->build(*ih));
408 
409  float sigma1_le2 = max(th[0]->parametersError()[0][0],
410  th[0]->parametersError()[1][1]);
411  float sigma2_le2 = max(th[1]->parametersError()[0][0],
412  th[1]->parametersError()[1][1]);
413 
414  float sigma_z2 = (1 + a23/a12)*(1 + a23/a12) * sigma2_le2 +
415  ( a23/a12)*( a23/a12) * sigma1_le2;
416 
417  float cotTheta = slope * circle.curvature(); // == sinhEta
418  float coshEta = sqrt(1 + sqr(cotTheta)); // == 1/sinTheta
419 
420  float pt = Bz / circle.curvature();
421  float p = pt * coshEta;
422 
423  float m_pi = 0.13957018;
424  float beta = p / sqrt(sqr(p) + sqr(m_pi));
425 
427  PixelRecoPointRZ rz2(g2.perp(), g2.z());
428 
429  float sigma_z = msp(pt, cotTheta, rz2) / beta;
430 
431  // Calculate globalDirs
432  float sinTheta = 1. / coshEta;
433  float cosTheta = cotTheta * sinTheta;
434 
435  int dir;
436  if(areaParallelogram(p1 - c, p2 - c) > 0) dir = 1; else dir = -1;
437 
438  float curvature = circle.curvature();
439 
440  {
441  Global2DVector v = (p1 - c)*curvature*dir;
442  globalDirs.push_back(GlobalVector(-v.y()*sinTheta,v.x()*sinTheta,cosTheta));
443  }
444 
445  {
446  Global2DVector v = (p2 - c)*curvature*dir;
447  globalDirs.push_back(GlobalVector(-v.y()*sinTheta,v.x()*sinTheta,cosTheta));
448  }
449 
450  {
451  Global2DVector v = (p3 - c)*curvature*dir;
452  globalDirs.push_back(GlobalVector(-v.y()*sinTheta,v.x()*sinTheta,cosTheta));
453  }
454 
455  // Multiple scattering
456  float sigma_ms = sigma_z * coshEta;
457 
458  // Local error squared
459  float sigma_le2 = max(th[2]->parametersError()[0][0],
460  th[2]->parametersError()[1][1]);
461 
462  return (delta_z*delta_z / (sigma_ms*sigma_ms + sigma_le2 + sigma_z2)
463  < nSigma * nSigma);
464  }
465 
466  return false;
467 }
468 
469 /*****************************************************************************/
471 {
472  if ( layer->location() == GeomDetEnumerators::barrel) {
473  theBarrel = true;
474  theForward = false;
475  const BarrelDetLayer& bl = dynamic_cast<const BarrelDetLayer&>(*layer);
476  float halfThickness = bl.surface().bounds().thickness()/2;
477  float radius = bl.specificSurface().radius();
478  theDetRange = Range(radius-halfThickness, radius+halfThickness);
479  } else if ( layer->location() == GeomDetEnumerators::endcap) {
480  theBarrel= false;
481  theForward = true;
482  const ForwardDetLayer& fl = dynamic_cast<const ForwardDetLayer&>(*layer);
483  float halfThickness = fl.surface().bounds().thickness()/2;
484  float zLayer = fl.position().z() ;
485  theDetRange = Range(zLayer-halfThickness, zLayer+halfThickness);
486  }
487 }
const double beta
dbl * delta
Definition: mlp_gen.cc:36
float originRBound() const
bounds the particle vertex in the transverse plane
void findRectangle(const float x[3], const float y[3], const float par[3], float phi[2], float z[2])
void invertCircle(Global2DVector &c, float &r)
T y() const
Definition: PV2DBase.h:46
void initLayer(const DetLayer *layer)
GlobalPoint const & origin() const
PixelRecoRange< float > Range
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
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)
virtual Location location() const =0
Which part of the detector (barrel, endcap)
T mag() const
Definition: PV2DBase.h:49
void calculateRangesBarrel(float r3, float phi[2], float z[2], bool keep)
T y() const
Definition: PV3DBase.h:63
const Bounds & bounds() const
Definition: Surface.h:120
virtual const BoundCylinder & specificSurface() const final
Extension of the interface.
float angleRatio(const Global2DVector &p3, const Global2DVector &c)
const int keep
GlobalVector inInverseGeV(const GlobalPoint &gp) const
Field value ad specified global point, in 1/Gev.
Definition: MagneticField.h:41
Geom::Phi< T > phi() const
Definition: PV2DBase.h:51
T curvature(T InversePt, const edm::EventSetup &iSetup)
susybsm::MuonSegmentRefProd msp
Definition: classes.h:34
T sqrt(T t)
Definition: SSEVec.h:18
T z() const
Definition: PV3DBase.h:64
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
T mag2() const
The vector magnitude squared. Equivalent to vec.dot(vec)
float areaParallelogram(const Global2DVector &a, const Global2DVector &b)
std::pair< float, float > findMinimalCircles(float r)
T min(T a, T b)
Definition: MathUtil.h:58
double p2[4]
Definition: TauolaWrapper.h:90
#define M_PI
virtual float thickness() const =0
virtual const Surface::PositionType & position() const
Returns position of the surface.
const T & get() const
Definition: EventSetup.h:58
float ptMin() const
minimal pt of interest
double b
Definition: hdecay.h:120
ThirdHitPrediction(const TrackingRegion &region, GlobalPoint inner, GlobalPoint outer, const edm::EventSetup &es, double nSigMultipleScattering, double maxAngleRatio, std::string builderName)
void spinCloser(float phi[3])
HLT enums.
const BoundSurface & surface() const final
The surface of the GeometricSearchDet.
double p1[4]
Definition: TauolaWrapper.h:89
std::pair< float, float > findArcIntersection(std::pair< float, float > a, std::pair< float, float > b, bool &keep)
double a
Definition: hdecay.h:121
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])
dbl *** dir
Definition: mlp_gen.cc:35
bool isCompatibleWithMultipleScattering(GlobalPoint g3, const std::vector< const TrackingRecHit * > &h, std::vector< GlobalVector > &localDirs, const edm::EventSetup &es)
void getRanges(const DetLayer *layer, float phi[], float rz[])
std::pair< float, float > findTouchingCircles(float r)
rm
Definition: submit.py:76
T x() const
Definition: PV2DBase.h:45
T x() const
Definition: PV3DBase.h:62
T const * product() const
Definition: ESHandle.h:86
const BoundSurface & surface() const final
GeometricSearchDet interface.
Global3DVector GlobalVector
Definition: GlobalVector.h:10
Vector2DBase< float, GlobalTag > Global2DVector
double p3[4]
Definition: TauolaWrapper.h:91
#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
float sqr(float x)