CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HelixBarrelCylinderCrossing.cc
Go to the documentation of this file.
4 
6 
7 #include <iostream>
8 #include <cmath>
9 
10 template <typename T>
11 inline T sqr(const T& t) {return t*t;}
12 
15  const GlobalVector& startingDir,
16  double rho, PropagationDirection propDir,
17  const Cylinder& cyl)
18 {
19  // assumes the cylinder is centered at 0,0
20  double R = cyl.radius();
21 
22  // protect for zero curvature case
23  const double sraightLineCutoff = 1.e-7;
24  if (fabs(rho)*R < sraightLineCutoff &&
25  fabs(rho)*startingPos.perp() < sraightLineCutoff) {
26  // switch to straight line case
27  StraightLineCylinderCrossing slc( cyl.toLocal(startingPos),
28  cyl.toLocal(startingDir), propDir);
29  std::pair<bool,double> pl = slc.pathLength( cyl);
30  if (pl.first) {
31  theSolExists = true;
32  theS = pl.second;
33  thePos = cyl.toGlobal(slc.position(theS));
34  theDir = startingDir;
35  }
36  else theSolExists = false;
37  return; // all needed data members have been set
38  }
39 
40  double R2cyl = R*R;
41  double pt = startingDir.perp();
42  Point center( startingPos.x()-startingDir.y()/(pt*rho),
43  startingPos.y()+startingDir.x()/(pt*rho));
44  double p2 = startingPos.perp2();
45  bool solveForX;
46  double B, C, E, F;
47  if (fabs(center.x()) > fabs(center.y())) {
48  solveForX = false;
49  E = (R2cyl - p2) / (2.*center.x());
50  F = center.y()/center.x();
51  B = 2.*( startingPos.y() - F*startingPos.x() - E*F);
52  C = 2.*E*startingPos.x() + E*E + p2 - R2cyl;
53  }
54  else {
55  solveForX = true;
56  E = (R2cyl - p2) / (2.*center.y());
57  F = center.x()/center.y();
58  B = 2.*( startingPos.x() - F*startingPos.y() - E*F);
59  C = 2.*E*startingPos.y() + E*E + p2 - R2cyl;
60  }
61 
62  RealQuadEquation eq( 1+F*F, B, C);
63  if (!eq.hasSolution) {
64  theSolExists = false;
65  return;
66  }
67 
68  Vector d1, d2;;
69  if (solveForX) {
70  d1 = Point(eq.first, E-F*eq.first);
71  d2 = Point(eq.second, E-F*eq.second);
72  }
73  else {
74  d1 = Point( E-F*eq.first, eq.first);
75  d2 = Point( E-F*eq.second, eq.second);
76  }
77 
78  chooseSolution(d1, d2, startingPos, startingDir, propDir);
79  if (!theSolExists) return;
80 
81  double ipabs = 1./startingDir.mag();
82  double sinTheta = pt * ipabs;
83  double cosTheta = startingDir.z() * ipabs;
84 
85  double dMag = theD.mag();
86  double tmp = 0.5 * dMag * rho;
87  if (std::abs(tmp)>1.) tmp = ::copysign(1.,tmp);
88  theS = theActualDir * 2.* asin( tmp ) / (rho*sinTheta);
89  thePos = GlobalPoint( startingPos.x() + theD.x(),
90  startingPos.y() + theD.y(),
91  startingPos.z() + theS*cosTheta);
92 
93  if (theS < 0) tmp = -tmp;
94  double sinPhi = 2.*tmp*sqrt(1.-tmp*tmp);
95  double cosPhi = 1.-2.*tmp*tmp;
96  theDir = DirectionType(startingDir.x()*cosPhi-startingDir.y()*sinPhi,
97  startingDir.x()*sinPhi+startingDir.y()*cosPhi,
98  startingDir.z());
99 }
100 
102  const PositionType& startingPos,
103  const DirectionType& startingDir,
104  PropagationDirection propDir)
105 {
106  double momProj1 = startingDir.x()*d1.x() + startingDir.y()*d1.y();
107  double momProj2 = startingDir.x()*d2.x() + startingDir.y()*d2.y();
108 
109  if ( propDir == anyDirection ) {
110  theSolExists = true;
111  if (d1.mag2()<d2.mag2()) {
112  theD = d1;
113  theActualDir = (momProj1 > 0) ? 1 : -1;
114  }
115  else {
116  theD = d2;
117  theActualDir = (momProj2 > 0) ? 1 : -1;
118  }
119  }
120  else {
121  int propSign = propDir==alongMomentum ? 1 : -1;
122  if (momProj1*momProj2 < 0) {
123  // if different signs return the positive one
124  theSolExists = true;
125  theD = (momProj1*propSign > 0) ? d1 : d2;
126  theActualDir = propSign;
127  }
128  else if (momProj1*propSign > 0) {
129  // if both positive, return the shortest
130  theSolExists = true;
131  theD = (d1.mag2()<d2.mag2()) ? d1 : d2;
132  theActualDir = propSign;
133  }
134  else theSolExists = false;
135  }
136 }
137 
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:78
void chooseSolution(const Point &p1, const Point &p2, const PositionType &startingPos, const DirectionType &startingDir, PropagationDirection propDir)
T perp() const
Definition: PV3DBase.h:66
double_binary B
Definition: DDStreamer.cc:235
Definition: DDAxes.h:10
std::pair< bool, double > pathLength(const Cylinder &cyl) const
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
T y() const
Definition: PV3DBase.h:57
#define abs(x)
Definition: mlp_lapack.h:159
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
Definition: Basic2DVector.h:71
T perp2() const
Definition: PV3DBase.h:65
PropagationDirection
tuple d1
Definition: debug_cff.py:7
T mag() const
Definition: PV3DBase.h:61
Scalar radius() const
Radius of the cylinder.
Definition: Cylinder.h:55
T sqrt(T t)
Definition: SSEVec.h:28
LocalPoint toLocal(const GlobalPoint &gp) const
T z() const
Definition: PV3DBase.h:58
HelixBarrelCylinderCrossing(const GlobalPoint &startingPos, const GlobalVector &startingDir, double rho, PropagationDirection propDir, const Cylinder &cyl)
double p2[4]
Definition: TauolaWrapper.h:90
T y() const
Cartesian y coordinate.
Definition: Basic2DVector.h:65
T mag2() const
The vector magnitude squared. Equivalent to vec.dot(vec)
Definition: Basic2DVector.h:68
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
Square< F >::type sqr(const F &f)
Definition: Square.h:13
T x() const
Definition: PV3DBase.h:56
T x() const
Cartesian x coordinate.
Definition: Basic2DVector.h:62