Initialize_Canonical() and getDistanceToHelix() updates for reference points other than (0,0,0)

[jfklama]

BEGINRELEASENOTES

Fixes of some of the issues mentioned in #24

• Added overloaded HelixClass::Initialize_Canonical() with a reference point passed as an input, not assuming that it is (0,0,0)
• In HelixClass::getDistanceToHelix(), for the arc length calculation, added the cases when the ref. point is passed and when it is nonzero

ENDRELEASENOTES

[tmadlener]

Hi @jfklama,

Thanks for this. One request I would have from my side is the following: Could you either remove the commented std::cout debug prints, or make them "proper" debug outputs with streamlog_out(DEBUG) (potentially choosing an appropriate debug level here to make sure it only appears with rather detailed output).

@dudarboh could you also have a look at the implementation to see if this indeed solves some of the issues you raised?

[dudarboh]

Hi @jfklama,
see some comments inline, and note that I am also not an expert, so these are my guesses which are also should be thoroughly checked.

I can't really comment on your changes for float HelixClass::getDistanceToHelix(HelixClass * helix, float * pos, float * mom) because I have no idea what it does algorithm-wise. Maybe you can write this more clear in doc-string inside the header, so one could make a rough understanding how does it check for the distance.

Also I see you re-calculate x(y)AtPCA many times, which should be already stored inside each helix.
Maybe it makes sense to write "getters" for each private member and use them? Or make members public at all...

Also, do these fixes work for tracks which are created with other broken constructors?
HelixClass::Initialize_BZ() and HelixClass::Initialize_VP() which put a reference point somewhere on the helix and have non-zero d0, z0 parameters?

@tmadlener, we also should not forget that this ...
https://github.com/iLCSoft/MarlinUtil/blob/master/source/src/HelixClass_double.cc
exists..

We should either keep them both consistent with each other, or remove one of them. Otherwise it is even a bigger mess..

Also, I think HelixClass is a perfect place to have tests that will check all possible combinations of different reference point positions, charges, coordinate quarters, even magnetic field, etc.

Debugging all the trigonometry 100 times is not fun. Better to do it once and put in some tests.

[tmadlener]

@tmadlener, we also should not forget that this ...
https://github.com/iLCSoft/MarlinUtil/blob/master/source/src/HelixClass_double.cc
exists..

Very good point. I can have a look at merging those into one templated class (in a separate PR).

Also, I think HelixClass is a perfect place to have tests that will check all possible combinations of different reference point positions, charges, coordinate quarters, even magnetic field, etc.

Debugging all the trigonometry 100 times is not fun. Better to do it once and put in some tests.

I fully agree with this. Covering everything from the Helix class now is probably a lot of work, but at least adding tests for newly added functionality or bug fixes should be rather straight forward. In principle there is a very simple test program already present here, but we can also try and set up "proper unit testing" rather quickly here.

[jfklama]

Dear @tmadlener, @dudarboh

I've made the additional changes - removed / changed outputs, added general description in the header and "getters" for x(y)AtPCA.

Also, do these fixes work for tracks which are created with other broken constructors?
HelixClass::Initialize_BZ() and HelixClass::Initialize_VP() which put a reference point somewhere on the helix and have non-zero d0, z0 parameters?

I have never used other initializers and didn't touch them. getDistanceToHelix() apparently doesn't use d0 and z0, so as long as other parameters provided by these initializers are correct (most importantly radius, x(y)Centre, tanLambda), the results should also be fine. However, I must admit I was considering only ref. points at the IP, first/last hits and maybe calorimeter, so I'm not sure; I can try to check this.

Regarding "proper" unit testing, unfortunately I have no experience with that - I can try to help somehow, but probably I would need to see some examples first...

[tmadlener]

@jfklama I have rebased this onto the latest version. Is this pull request still relevant? Could you also briefly check that I didn't introduce any obvious mistakes during rebasing?

[tmadlener]

@jfklama I rebased this again to pick up the unit test setup, so it should now be a lot easier to add a few tests that make sure that things work as expected.

[jfklama]

@tmadlener I have found some small inconsistencies, it should fine now.

[tmadlener]

Thanks for checking (and fixing my mistakes). Would it be possible to add some (or at least one) test for this? There is a basic example of how to do this in source/tests/unittests/TestHelixClass.cpp

Essentially what you need to do is:

• Create a helix and initialize it with this new constructor (choosing some appropriate values for the parameters. If necessary simply invent values)
• Check that the properties of the initialized helix are as expected using the REQUIRE macro.
  • Since this constructor is about reference points different from (0, 0, 0), properties that are affected by the reference point should get the most attention in the tests.

[jfklama]

I've added some tests of calculating distance between helices. It seems that my changes worked at least in my case (highly displaced, non-pointing tracks) but in general the behavior can still be unexpected. However, the backward compatibility looks fine. I assume that if someone used HelixClass so far, they tested the performance (like e.g. in the case of V0Finder) so we should keep it.

It's quite surprising that some tests are passed by one HelixClass version and failed by the other...

[tmadlener]

Hi @jfklama, thanks for adding some tests. It seems that at least in the CI that we are running here quite a few of them are failing at the moment: https://github.com/iLCSoft/MarlinUtil/actions/runs/3432147649/jobs/5734311976#step:4:628 Are these working locally for you?

I have an additional style comment below.

[jfklama]

Are these working locally for you?

No, these tests correspond to the '(0,0,0) reference point' case with large d0 and z0 parameters which I wanted to add for completeness, but of course I didn't take into account it will make the checks here fail. Also, calculating distance with Initialize_VP and zero ref. point give wrong results. These tests also fail if I use the HelixClass version without my changes, but I don't have any idea yet on how to fix this without breaking the backward compatibility. Should I remove these tests?

Honestly, with the form of this algorithm, I am not sure whether it even fundamentally makes sense to use (0,0,0) as a reference point if tracks are highly displaced. Would it make sense to put some warnings in the code for those case scenarios?

[tmadlener]

OK. So if I understand correctly, you added tests for your newly added changes, and these are working as expected (at least for the use cases you intend to use them)? But then you also added some more tests, trying to cover existing functionality, and these are currently broken. If that is the case, I would leave the unittests for your expected behavior and remove the ones from that try to use the new developments together with existing functionality.

We will deal with the non-trivial things from the existing functionality in separate PRs. Those should be kept rather small, so that they are not too much work to create and review.

[dudarboh]

@jfklama
Maybe my inline comments got missed in the active discussion, but to highlight, as I mentioned above:

Two critical bugs

There are two critical bugs:

1. In the new Initialize_Canonical()

// this is currently implemented and it is simply wrong
//  _xCentre = _referencePoint[0] + _radius*cos(_phi0-_const_pi2*_charge);
// _yCentre = _referencePoint[1] + _radius*sin(_phi0-_const_pi2*_charge);
// this is how it should be
_xCentre = _xAtPCA + _radius*cos(_phi0-_const_pi2*_charge);
_yCentre = _yAtPCA + _radius*sin(_phi0-_const_pi2*_charge);

2. As I also just now have seen: in the old Initialize_Canonical()

// this is currently implemented and it is simply wrong
// _referencePoint[0] = _xAtPCA;
// _referencePoint[1] = _yAtPCA;
// _referencePoint[2] = _z0;
// this is how it should be
_referencePoint[0] = 0.;
_referencePoint[1] = 0.;
_referencePoint[2] = 0.;

but changing latter might have an impact on currently running code, so one has to thoroughly check if this part is used anywhere @tmadlener...

`getDistanceToHelix() requires both the center of the helix and reference point:

  FloatT x01 = getXC();
  FloatT y01 = getYC();
  FloatT x02 = helix->getXC();
  FloatT y02 = helix->getYC();
 ...
  FloatT ref1[3];
  FloatT ref2[3];
  for (int i=0;i<3;++i) {
    ref1[i]=_referencePoint[i];
    ref2[i]=helix->getReferencePoint()[i];
  }

So, without fixing two bugs above, even if the algorithm inside getDistanceToHelix() is correct, there is no chance it will produce meaningful results with these constructors...

[jfklama]

Hi @dudarboh

There are two critical bugs:

1. In the new `Initialize_Canonical()`

// this is currently implemented and it is simply wrong
//  _xCentre = _referencePoint[0] + _radius*cos(_phi0-_const_pi2*_charge);
// _yCentre = _referencePoint[1] + _radius*sin(_phi0-_const_pi2*_charge);
// this is how it should be
_xCentre = _xAtPCA + _radius*cos(_phi0-_const_pi2*_charge);
_yCentre = _yAtPCA + _radius*sin(_phi0-_const_pi2*_charge);

Thanks, I missed it - most likely, because I used reference points close to the helix, so the difference was not noticeable. This is also a consequence of the "old" implementation (which you point out later), where x_yPCA were assigned to the reference point. I can quickly correct this.

2. As I also just now have seen: in the old `Initialize_Canonical()`

// this is currently implemented and it is simply wrong
// _referencePoint[0] = _xAtPCA;
// _referencePoint[1] = _yAtPCA;
// _referencePoint[2] = _z0;
// this is how it should be
_referencePoint[0] = 0.;
_referencePoint[1] = 0.;
_referencePoint[2] = 0.;

but changing latter might have an impact on currently running code, so one has to thoroughly check if this part is used anywhere @tmadlener...

This I have left on purpose, exactly because these variables are used in the HelixClass::getDistanceToHelix(), which, in turn, is used (at least) in the V0Finder. Its performance has been studied and it seems to give reasonable results: https://agenda.linearcollider.org/event/7839/contributions/40920/attachments/32702/49744/2018_02_18_SWAna_Premeeting.pdf

So I didn't want to break it, even though I would guess your suggestion should improve it (but this would require tests). I don't know what effect this fix would have on other packages.

[dudarboh]

Ah, ok, I see.

I should admit, I completely lost the understanding how does it work.
But if it works for your case, then good :)

[jfklama]

But if it works for your case, then good :)

Hopefully not only for my case...

I think that the old implementation was "approximately" good, as long as one used ref. point = (0,0,0) and rather small d0, z0.

On the other hand, maybe we should consider fixing also the current Initialize_Canonical, as it fundamentally doesn't make sense, instead of building on top of the bugged code? If other packages used it and gave reasonable results, either they should also work (or improve) after the corrections, or there is something wrong with themselves as well...

[dudarboh]

@jfklama I think we should consider fixing Initilize_Canonical indeed :)
But unfortunately this is a mess. Looking also at the initial code of the getDistanceToHelix() I don't understand some of its parts.. So I am not sure if fixing bugged code in the constructor will not break some underlying assumptions inside getDistanceToHelix() and completely break everything..

As we all agree, this whole business needs a complete ovehaul. Which we will not do for the next patch.

So, we either:

1. make few small positive changes which are backwards compatible for the next patch.
2. leave this mess as it is for the next patch and maybe try to assess how much effort it is to fix everything completely long-term, e.g. for the next-next patch.

@tmadlener, your take on this?

[tmadlener]

Apologies for not coming back to this earlier. This seems to be becoming a larger and larger project the more we (well you) look into it. Maybe we should try to get something into the next release that is useful for you and then try to see how to actually tackle this in a more systematic way.

I think currently not breaking existing code is the most important concern. So anything that could change that should be avoided (at least for the upcoming tag and patch release). As far as I can tell, the current version does that, right? I.e. it only adds a new constructor, but essentially leaves the rest of the functionality as it is. Or could it potentially also change outcomes for other constructors?

I have left a few comments inline. One thing that is not yet entirely clear to me (and here you are probably the experts), is how tight the floating point comparisons can be made. Ideally they would be very tight, and really essentially testing numerical equality, rather than "equality in a physics sense" (i.e. probably doesn't change the outcome of an analysis).

Another option to at least make it possible for you to continue with your work would be the following: As the long term plan is to go to the aidaTT helix utilities in any case, you could also directly start using that instead of trying to tack on the necessary functionality on to here, where we know, that we will have to change it in any case in the future. I am not entirely sure how much work that entails, but maybe you could have a quick look to see if that will take days or if it is more likely going to be weeks.

[tmadlener]

Just for slighly better readability of the tests, could some of these magic numbers be put into variables with a meaningful way? That would make the comparisons below a bit easier to interpret as well.

[tmadlener]

Suggested change

	FloatT vtx_momentum1[3];
	FloatT vtx_momentum2[3];
	FloatT vertex1[3];
	FloatT vertex2[3];
	std::array<FloatT, 3> vtx_momentum1;
	std::array<FloatT, 3> vtx_momentum2;
	std::array<FloatT, 3> vertex1;
	std::array<FloatT, 3> vertex2;

Just for consistency throughout the tests. You probably have to use a .data() below then to pass things.

[tmadlener]

same as above.

[tmadlener]

Is 2 mm(?) a reasonable difference for it still to be considered non-zero? Or is this mainly to make these tests pass? Is it possible to calculate the actual expected distance in this case, and check that with a much smaller epsilon?

[jfklama]

It's partly to make the tests pass, but also, honestly, I don't see any way to calculate this, other than as currently implemented in getDistanceToHelix(). If the helices do not cross (which of course is the case in real life), the calculation becomes highly non-trivial, as one can see from complexity of the algorithm. So either we accept the current algorithm prediction as the "correct" and use it in the test, or assume the distance should be small and accept everything close to zero within some margin.

I chose the latter option because: 1) I think 2 mm is a good precision for such a high track displacement 2) from my experience, when there is a bug in the algorithm, the calculated distance explodes, 3) even though I went through this code and I believe I understand how it works, it is still not mine and it might happen that I am missing something (as with the x_yCentre bug pointed out by @dudarboh)

[jfklama]

Hi @tmadlener, sorry for a late response.

I think currently not breaking existing code is the most important concern. So anything that could change that should be avoided (at least for the upcoming tag and patch release). As far as I can tell, the current version does that, right? I.e. it only adds a new constructor, but essentially leaves the rest of the functionality as it is. Or could it potentially also change outcomes for other constructors?

In principle yes, for the Initialize_Canonical(). If someone uses other constructors + getDistanceToHelix(), it also should work, but unfortunately I cannot guarantee that, as I haven't performed any proper tests - especially for the Initialize_BZ(). I have no experience with the other constructors and I am not sure I correctly understand them. Unless this combination "other constructor" + getDistanceToHelix() appears anywhere, the current functionality will remain intact.

I have left a few comments inline. One thing that is not yet entirely clear to me (and here you are probably the experts), is how tight the floating point comparisons can be made. Ideally they would be very tight, and really essentially testing numerical equality, rather than "equality in a physics sense" (i.e. probably doesn't change the outcome of an analysis).

Thanks for the comments. Unfortunately, regarding the floating point numbers, that is not clear for me as well. I have seen some strange behaviour depending on precision of the input variables...

Another option to at least make it possible for you to continue with your work would be the following: As the long term plan is to go to the aidaTT helix utilities in any case, you could also directly start using that instead of trying to tack on the necessary functionality on to here, where we know, that we will have to change it in any case in the future. I am not entirely sure how much work that entails, but maybe you could have a quick look to see if that will take days or if it is more likely going to be weeks.

For me the main problem with the aidaTT helix utilities is that it doesn't have any tool to calculate distance/PCA between helices ;) I can see some functions for finding intersections which could be potentially used for that, but at first glance I have no idea how they work. So it could take some time...

To be very honest, for my work I can use my implementation of the HelixClass as well, even if it's sub-optimal. I've created that pull request so that in the future no one runs into the same problems as me, but it seems that we still have more problems than solutions. Since the long-term goal is to switch to aidaTT anyway, maybe it is a waste of time to further tinker with this, especially if maintaining the current functionality is the main priority?