11 : thePosition(point), theRho(curvature), thePropDir(propDir) {
15 double px = direction.
x();
16 double py = direction.
y();
17 double pz = direction.
z();
18 double pt = px * px + py *
py;
19 double p =
sqrt(pt + pz * pz);
37 ceq[3] = 2 * helix2.
mag2();
63 ceq[3] = 2 * helix2.
dot(helix2p);
66 ceq[2] = 3 *
theDirection.
dot(lineDirection) * helix2.dot(lineDirection);
68 ceq[0] = deltaPos.
dot(helix1p);
86 for (
unsigned int i = 0;
i < nRaw;
i++) {
89 solutions[nDir++] = solutions[
i];
92 return std::make_pair(
false, 0.);
97 for (
unsigned int i = 0;
i < nDir;
i++) {
98 double st = solutions[
i];
99 double deri2 = (3 * ceq[3] * st + 2 * ceq[2]) * st + ceq[1];
101 solutions[nMin++] = st;
104 return std::make_pair(
false, 0.);
108 double dSt = solutions[0];
109 for (
unsigned int i = 1;
i <
nMin;
i++) {
110 if (fabs(solutions[
i]) < fabs(dSt))
121 if (fabs(ceq[3]) > FLT_MIN) {
123 double q = (ceq[2] * ceq[2] - 3 * ceq[3] * ceq[1]) / (ceq[3] * ceq[3]) / 9.;
124 double r = (2 * ceq[2] * ceq[2] * ceq[2] - 9 * ceq[3] * ceq[2] * ceq[1] + 27 * ceq[3] * ceq[3] * ceq[0]) /
125 (ceq[3] * ceq[3] * ceq[3]) / 54.;
126 double q3 = q * q *
q;
128 double phi = acos(r /
sqrt(q3)) / 3.;
129 double rootq =
sqrt(q);
130 for (
int i = 0;
i < 3;
i++) {
131 solutions[
i] = -2 * rootq *
cos(phi) - ceq[2] / ceq[3] / 3.;
132 phi += 2. / 3. *
M_PI;
136 double a =
pow(fabs(r) +
sqrt(r * r - q3), 1. / 3.);
139 double b = fabs(a) > FLT_MIN ? q / a : 0.;
140 solutions[0] = a + b - ceq[2] / ceq[3] / 3.;
148 else if (fabs(ceq[2]) > FLT_MIN) {
154 solutions[0] = -ceq[0] / ceq[1];
161 double deq1 = coeff[1] * coeff[1];
162 double deq2 = coeff[2] * coeff[0];
163 if (fabs(deq1) < FLT_MIN || fabs(deq2 / deq1) > 1.
e-6) {
167 double deq = deq1 + 2 * deq2;
170 double ceq = -0.5 * (coeff[1] + (coeff[1] > 0 ? 1 : -1) *
sqrt(deq));
171 solutions[0] = -2 * ceq / coeff[2];
172 solutions[1] = coeff[0] / ceq;
178 double ceq = coeff[1] / coeff[2];
179 double deq = deq2 / deq1;
180 deq *= (1 - deq / 2);
181 solutions[0] = -ceq * deq;
182 solutions[1] = ceq * (2 + deq);
DirectionType direction() const
T y() const
Cartesian y coordinate.
T x() const
Cartesian x coordinate.
PositionType position() const
T curvature(T InversePt, const edm::EventSetup &iSetup)
T z() const
Cartesian z coordinate.
Cos< T >::type cos(const T &t)
Power< A, B >::type pow(const A &a, const B &b)
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
T mag2() const
The vector magnitude squared. Equivalent to vec.dot(vec)
T dot(const Basic3DVector &rh) const
Scalar product, or "dot" product, with a vector of same type.