28 declareInterface<VertexPointEstimator>(
this);
38 return StatusCode::SUCCESS;
43 return StatusCode::SUCCESS;
69 decors[
"deltaPhiTracks"],
76 return {
"deltaPhiTracks",
"DR1R2"};
102 double minArcLength(0.);
103 double maxArcLength(0.);
143 for (
int I=0;
I<2; ++
I) {
158 DZSVI[
I] = 2.0*DZMAX;
166 const AmgVector(5)& perParam1 = per1->parameters();
168 const AmgVector(5)& perParam2 = per2->parameters();
174 double theta1 = perParam1[
Trk::theta]; cotTheta[0] = 1./tan(theta1);
175 double theta2 = perParam2[
Trk::theta]; cotTheta[1] = 1./tan(theta2);
178 RC[0] = 10.*0.5/(qOverPt[0]*aconst); RC[1] = 10.*0.5/(qOverPt[1]*aconst);
180 XC[0] = -1*U*sin(
phi[0]);
181 YC[0] = U*cos(
phi[0]);
184 XC[1] = -1*U*sin(
phi[1]);
185 YC[1] = U*cos(
phi[1]);
191 double theta2 = perParam2[
Trk::theta]; cotTheta[0] = 1./tan(theta2);
192 double theta1 = perParam1[
Trk::theta]; cotTheta[1] = 1./tan(theta1);
195 RC[0] = 10.*0.5/(qOverPt[0]*aconst); RC[1] = 10.*0.5/(qOverPt[1]*aconst);
196 double U = RC[0] + d0[0];
197 XC[0] = -1*U*sin(
phi[0]);
198 YC[0] = U*cos(
phi[0]);
201 XC[1] = -1*U*sin(
phi[1]);
202 YC[1] = U*cos(
phi[1]);
207 double DX = XC[1] - XC[0];
208 double DY = YC[1] - YC[0];
209 double D = DX*DX + DY*DY;
212 ATH_MSG_DEBUG(
"Concentric circles, should not happen return (0,0,0)");
216 U = D - RA[0] - RA[1];
220 deltaR = U;
double PHI = -99999.;
221 double hl =
areaVar(XC[0], YC[0], RA[0], XC[1], YC[1], RA[1], PHI);
225 double Y0 = (-1*XC[0]*YC[1] + XC[1]*YC[0])/D;
226 for (
int I=0;
I<2; ++
I) {
227 AB[
I] = COST*XC[
I] + SINT*YC[
I];
229 U = (XC[1] + XC[0])*(XC[1] - XC[0]) + (YC[1] + YC[0])*(YC[1] - YC[0]);
230 double V = (RA[1] + RA[0])*(RA[1] - RA[0]);
231 double XX = 0.5*(U - V)/D;
232 if ((XX - AB[0])*(XX - AB[0]) > 0.) {
233 U = sqrt((XX - AB[0])*(XX - AB[0]));
235 double YY2 = (RA[0] + U)*(RA[0] - U);
243 for (
int I=0;
I<2; ++
I) {
245 XSVI[
I] = COST*XX - SINT*U;
246 YSVI[
I] = SINT*XX + COST*U;
255 U = D - RA[0] - RA[1];
261 if (AB[0] < AB[1]) XX = AB[0] + RA[0];
266 V = AB[0] - RA[0] - XX;
270 V = AB[1] + RA[1] - XX;
274 XSVI[0] = COST*XX - SINT*Y0;
275 YSVI[0] = SINT*XX + COST*Y0;
279 ATH_MSG_DEBUG(
"XY distance of minimum approach is too large, return (0,0,0)");
285 for (
int J=0; J<nsol; ++J) {
286 U = (XSVI[J] - XC[0])/RC[0];
287 V = -1*(YSVI[J] - YC[0])/RC[0];
288 U = atan2(U,V) -
phi[0];
292 ZZ1[J] = z0[0] + SS1[J]*cotTheta[0];
293 U = (XSVI[J] - XC[1])/RC[1];
294 V = -1*(YSVI[J] - YC[1])/RC[1];
295 U = atan2(U,V) -
phi[1];
299 ZZ2[J] = z0[1] + SS1[J]*cotTheta[1];
300 RVI[J] = sqrt(XSVI[J]*XSVI[J] + YSVI[J]*YSVI[J]);
301 DZSVI[J] = ZZ2[J] - ZZ1[J];
303 ZSVI[0] = 0.5*(ZZ1[0] + ZZ2[0]);
304 ZSVI[1] = 0.5*(ZZ1[1] + ZZ2[1]);
306 if (std::min(RVI[0],RVI[1]) > RVMAX) {
311 if (std::min(fabs(DZSVI[0]),fabs(DZSVI[1])) > DZMAX) {
316 double A = std::min(SS1[0],SS2[0]);
317 double B = std::min(SS1[1],SS2[1]);
318 if (std::max(
A,B) < minArcLength || std::max(
A,B) > maxArcLength) {
325 if ( nsol == 2 && (fabs(DZSVI[1]) < fabs(DZSVI[0])) ) J = 1;
330 intPoint(0)=X; intPoint(1)=Y; intPoint(2)=Z;
357 double ret = -999999;
364 double h = 0.5*(sqrt(
pow(xi1-xi2,2) +
pow(yi1-yi2,2) ));
365 double l = sqrt(
pow(xc1-xc2,2) +
pow(yc1-yc2,2) );
368 double norm1 = sqrt((xi1-xc1)*(xi1-xc1)+(yi1-yc1)*(yi1-yc1));
369 double norm2 = sqrt((xi1-xc2)*(xi1-xc2)+(yi1-yc2)*(yi1-yc2));
370 if((norm1!=0.) && (norm2!=0.)){
374 double xa = (xi1-xc1)*norm1;
375 double ya = (yi1-yc1)*norm1;
376 double xb = (xi1-xc2)*norm2;
377 double yb = (yi1-yc2)*norm2;
378 double costheta = xa*xb + ya*yb;
379 phi =
M_PI-std::acos(costheta);
388 if (
x < -1.)
return M_PI;
389 if (
x > 1.)
return 0;
397 double ret = -999999;
410 h = 0.5*(sqrt(
pow(xi1-xi2,2) +
pow(yi1-yi2,2) ));
412 double l = sqrt(
pow(xc1-xc2,2) +
pow(yc1-yc2,2) );
420 double norm1 = sqrt(
pow(xi1-xc1,2)+
pow(yi1-yc1,2));
421 double norm2 = sqrt(
pow(xi1-xc2,2)+
pow(yi1-yc2,2));
422 if ((norm1 != 0) && (norm2 != 0)){
424 double xa = (xi1 - xc1)/norm1;
425 double ya = (yi1 - yc1)/norm1;
426 double xb = (xi1 - xc2)/norm2;
427 double yb = (yi1 - yc2)/norm2;
428 double costheta = xa*xb +ya*yb;
439 double xc2,
double yc2,
double r2,
440 double& xi1,
double& yi1,
441 double& xi2,
double& yi2)
454 double A = (xc1 - xc2) / (yc2- yc1);
455 double B = (r1*r1 - r2*r2 - xc1*xc1 + xc2*xc2 - yc1*yc1 + yc2*yc2) / 2. / ( yc2 -yc1);
457 double b = 2*
A*B - 2*xc1 -2*
A*yc1;
458 double c = B*B - 2*B*yc1 + xc1*xc1 + yc1*yc1 - r1*r1;
467 double A = (yc1 - yc2) / (xc2- xc1);
468 double B = (r1*r1 - r2*r2 - xc1*xc1 + xc2*xc2 - yc1*yc1 + yc2*yc2) / 2. / ( xc2 -xc1);
470 double b = 2*
A*B - 2*yc1 -2*
A*xc1;
471 double c = B*B - 2*B*xc1 + xc1*xc1 + yc1*yc1 - r1*r1;
491 double discr = b*b - 4*
a*c;
492 if (discr < 0)
return false;
493 y1 = (-b+sqrt(discr))/2./
a;
494 y2 = (-b-sqrt(discr))/2./
a;
502 double g,
double h) {
506 return fabs(0.5* ( (
a*e*i + d*
h*c + b*f*g) - (g*e*c + d*b*i +
a*f*
h) ) );
Scalar deltaPhi(const MatrixBase< Derived > &vec) const
Scalar deltaR(const MatrixBase< Derived > &vec) const
Scalar phi() const
phi method
constexpr int pow(int base, int exp) noexcept
Header file for AthHistogramAlgorithm.
static double areaVar(double, double, double, double, double, double, double &)
static std::vector< std::string > decorKeys()
Return list of keys used for decorations.
BooleanProperty m_returnOnError
static bool secondDegree(double, double, double, double &, double &)
DoubleArrayProperty m_maxDZ
Amg::Vector3D intersectionImpl(const Trk::Perigee *per1, const Trk::Perigee *per2, unsigned int flag, int &errorcode, float &deltaPhi, float &deltaR) const
internal implementation
static bool circleIntersection(double, double, double, double, double, double, double &, double &, double &, double &)
static const double s_bmagnt
DoubleArrayProperty m_maxHl
virtual StatusCode initialize() override
VertexPointEstimator(const std::string &type, const std::string &name, const IInterface *parent)
DoubleArrayProperty m_minDr
DoubleArrayProperty m_minArcLength
static const InterfaceID & interfaceID()
DoubleArrayProperty m_maxDr
DoubleArrayProperty m_maxR
static double areaTriangle(double, double, double, double, double, double)
Amg::Vector3D getCirclesIntersectionPoint(const Trk::Perigee *per1, const Trk::Perigee *per2, unsigned int flag, int &errorcode) const
Get intersection point of two track helices.
virtual StatusCode finalize() override
DoubleArrayProperty m_maxPhi
DoubleArrayProperty m_maxArcLength
std::map< std::string, float > Values_t
DoubleArrayProperty m_maxDR
Eigen::Matrix< double, 3, 1 > Vector3D
static const InterfaceID IID_IVertexPointEstimator("InDet::VertexPointEstimator", 1, 0)
double internal_acos(double x)
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
hold the test vectors and ease the comparison