A revision of the HelixClass

[dudarboh]

Brief description

After the discussion with @yradkhorrami we discovered some issues in the HelixClass which I will describe here below. All problems are connected to the fact that HelixClass initializers all rely on the assumption that reference point = (0, 0, 0) even in the situations when it is really not supposed to.
This might lead to ambiguity and misunderstanding all over the place.

HelixClass::Initialize_Canonical()

Issue

A helix in the canonical representation is defined by the reference point and 5 parameters: phi0, omega, tanL, d0, z0.

This initializer is useful to define the helix from the trackState which has all parameters defined.
However this initializer doesn't take reference point as an input argument and assumes reference point always to be (0,0,0)
which can lead to a confusion:

auto ts = track->getTrackState(AtLastHit);
HelixClass helix;
helix.Initialize_Canonical( ts->getPhi0(), ts->getD0(), ts->getZ0(), ts->getOmega(), ts->geTanLambda(), 3.5);

Here the helixClass treats passed d0, z0 arguments as with respect to the (0, 0, 0). However in the trackState they are written with respect to the arbitrary reference_point=ts->getReferencePoint() which is not necessarily (0, 0, 0).
This will result with a helix of the same curvature and dip, but completely different positioning as one would expect from the trackState.

By simply copying d0 and z0 parameters we completely transform their geometrical meaning:

helix.getD0() == ts.getD0() // true
helix.getZ0() == ts.getZ0() // true

As the ts.getZ0() is the distance along z between helix PCA and the ts->getReferencePoint()[2]
While the helix.getZ0() is the distance along z between helix PCA and the z=0.
So they shouldn't be identical if we talk about "the same" helix.

In the current state to get the desired helix one would have to be aware to recalculate d0 and z0 from the track state assuming reference point is (0, 0, 0), e.g.:

auto ts = track->getTrackState(AtLastHit);
HelixClass helix;
double newD0 = getD0AtIP(ts->getReferencePoint(), ts->Phi0(), ts->getD0(), ts->getZ0), ...);
double newZ0 = getZ0AtIP(ts->getReferencePoint(), ts->Phi0(), ts->getD0(), ts->getZ0), ...);

helix.Initialize_Canonical( ts->getPhi0(), newD0, newZ0, ts->getOmega(), ts->geTanLambda(), bField);

which is not necessarily trivial...

Potential way to improve

Provide reference point as an input argument to avoid ambiguity.

HelixClass::Initialize_VP()

Issue

This initializer is useful if one has:

• (x, y, z) coordinates of the point on the helix
• momentum of the helix at the given point

It uses the given point as a reference point to construct the helix:

void HelixClass::Initialize_VP(float * pos, float * mom, float q, float B) {
_referencePoint[0] = pos[0];
_referencePoint[1] = pos[1];
_referencePoint[2] = pos[2];

One would expect that output d0 and z0 paremeters must be 0 as the reference point is on the helix by definition.
However, while calculating d0 and z0 initializer comes backs to assuming reference point as (0,0,0) without mentioning it anywhere:

_d0 = double(q)*radius - q_sign * double(sqrt(xCentre*xCentre+yCentre*yCentre)); // not 0
_z0 = pos[2] + _radius*_tanLambda*q*(deltaPhi + _const_2pi*nCircles); // not 0

This would fail the user, if he plans to use d0, z0 helix parameters further in the program, e.g. to draw the helix,
which @yradkhorrami has experienced in the following example:

float pos[3] = getVertexPosition();
float mom[3] = getMomentum();

HelixClass helix;
helix.Initialize_VP(pos, mom, charge, bField);

double Px = pt * std::cos(  helix->getPhi() );
double Py = pt * std::sin(  helix->getPhi() );
double Pz = pt * helix->getTanLambda() ;
double Xs = helix->getReferencePoint()[0] -  helix->getD0() * std::sin( helix->getPhi() ) ; // d0 should be 0
double Ys = helix->getReferencePoint()[1] +  helix->getD0() * std::cos( helix->getPhi() ) ; // d0 should be 0
double Zs = helix->getReferencePoint()[2] +  helix->getZ0() ; // z0 should be 0

DDMarlinCED::drawHelix(bField, charge, Xs, Ys, Zs, Px, Py, Pz, ...);
// results in the different helix from:
// DDMarlinCED::drawHelix(bField, charge, pos[0], pos[1], pos[2], mom[0], mom[1], mom[2], ...);

Potential way to improve

Cross check all the equations and remove any assumptions of reference point as 0,0,0 where it is definitely shouldn't be and set d0 and z0 as 0.

HelixClass::Initialize_BZ()

Issue

This initializer is good when we know exact geometrical properties of the helix:

• Coordinates of the center: xCentre, yCentre
• Radius radius
• Dip bZ
• starting phase of the helix phi0

One should immediately note very inconvenient naming of the phi0 parameter.
Inside the Initialize_BZ() the phi0 parameter is a phase shift of the initial helix along z and is NOT a track parameter as in previous two initializers, which is a great point to become confused in this context.

We also see inconsistency between helix.getReferencePoint() which seems to be ON the helix:

  _referencePoint[2] = zBegin;
  _referencePoint[0] = xCentre + radius*cos(bZ*zBegin+phi0);
  _referencePoint[1] = yCentre + radius*sin(bZ*zBegin+phi0);

and output d0 and z0 parameters that have nothing to do with the helix reference point. It is hard to understand for me from the equations with respect to what point d0 and z0 are actually calculated in this case...

_d0 = -_xAtPCA*sin(_phi0) + _yAtPCA*cos(_phi0);
_z0 = _referencePoint[2] - (deltaPhi + _const_2pi*nCircles)/bZ;  

I am very confused about the most of the equations in the Initialize_BZ(). So some discussion with experts would be nice. E.g.:
_phiAtPCA = atan2(-_yCentre,-_xCentre); does it make sense?

Potential way to improve

Cross check all of the equations, rename phi0 to something different, e.g. phase0

HelixClass_double ?

This is just a comment to bring a discussion on why we have separate class for the HelixClass with double. Is it bad to transit HelixClass from float to double? Overload/template? Isn't there any better way to implement this?

List of packages using HelixClass

1. Garlic/GarlicExtendedTrack
   providing a comment:
   (231) // get helix-plane intersection by hand (the one from HelixClass seems strange/buggy...)
2. LCFIVertex/ConversionTagger
3. FPCCDSiliconTracking_MarlinTrk
4. V0Finder
5. (Its quite a lot, I will finish it later at some point...)

Any thoughts, comments @tmadlener, @gaede?

[jfklama]

A similar issue appears in the function HelixClass::getDistanceToHelix(HelixClass * helix, float * pos, float * mom), not only because it obviously will not give correct results for a reference point other than (0,0,0). Also, the following set of conditions

if (deltaPhi21 < 0 && charge2 < 0) {
  deltaPhi21 += _const_2pi;
}
else if (deltaPhi21 > 0 && charge2 > 0) { 
  deltaPhi21 -= _const_2pi;
}

if (deltaPhi22 < 0 && charge2 < 0) {
  deltaPhi22 += _const_2pi;
}
else if (deltaPhi22 > 0 && charge2 > 0) { 
  deltaPhi22 -= _const_2pi;
}

and the corresponding one for the other helix again assume that the ref. point is (0,0,0). In a different case (e.g. for a ref. point in the first/last hit), deltaPhi should be shifted if its sign is opposite to the charge sign. This can be easily fixed by simply adding the alternative set of conditions for a case when ref != (0,0,0). I also think that the shift should be done only if deltaPhi > M_PI.

[dudarboh]

As this topic is popped up in my feed I will also add

that I think Yasser also encountered that getDistanceToHelix gave very unreasonable results, so he wrote his own version iirc.

Also, while, HelixClass is not properly working people will create their own versions of getDistanceToHelix, getDistanceToPlane, getDistanceToPoint all over the code, like here:

https://github.com/iLCSoft/Garlic/blob/32edf27735ac1a58bbb986535b82ce0c18c3dc9e/src/GarlicExtendedTrack.cc#L231
(I found only this one, because it explicitly documented in the comment, but there might be more)

which is, of course, not ideal in the long term, as it hard to catch them all.

[jfklama]

For me it works fine just after these two fixes (Initialize_Canonical with a ref. point from input + deltaPhi), so it is fairly easy, at least in this case. But on the other hand it will require initialisation changes in the packages which use HelixClass.

Do you know if Yasser found some problems also for a referencePoint = (0, 0, 0)? Because I assumed they only occur for other ref. points; at least if you go through the getDistanceToHelix code it seems reasonable for (0, 0, 0).

[dudarboh]

I think he used it for some analysis involving tracks from 2nd and/or 3rd vertices. So I think so the main issue was indeed in reference point not being at 0.
I would believe that it should work for (0,0,0), but I wasn't investigating code in depth myself.

[tmadlener]

Hi @dudarboh, @jfklama,

@gaede and me just had a discussion about this issue and we agree that the current status is far from ideal and should be addressed. However, there is overall no easy solution to this, given the number of places where HelixClass actually shows up in iLCSoft see here. Touching all these places and making sure that nothing gets broken in the process is quite a large effort and we don't see who could do it at the moment.

One thing that we would like to point out is that in aidaTT, there is similar and in principle more modern Helix utility functionality. We are not yet sure how large the overlap between the two implementations actually is and whether the aidaTT version addresses all your needs at this point. If you have time you could have a look into that to see if you can find the necessary functionality and use that instead, to at least not introduce more usage of the HelixClass into existing code. Obviously, this mainly applies for new things, rather than already existing pieces of code that use the HelixClass.

In the long run the goal would be to converge on only one implementation of this functionality (and the current goal for that would be aidaTT). From the list of usages above this will be a major effort and we will most likely have to do it piece by piece. The overall strategy would be somewhat along the following lines:

• Cover the existing HelixClass with unit tests to capture the current functionality. The main purpose is to essentially programmatically document the current functionality, rather than making sure that this is the most sensible thing.
• See if the current HelixClass can be implemented via the aidaTT functionality under the hood. With the unit tests it should at least be possible to find out whether there are some changes in the results.
• Mark the HelixClass as deprecated and start replacing it.

For now starting with the unit tests would be the first step. Here I would start with providing the necessary "harness", but we would also definitely need your help to actually put together meaningful unit tests.