PMT-only Vertex Reconstruction

[S81D]

Describe your changes

VertexLeastSquares is a PMT-based vertex reconstruction tool that can be used for low energy events (NC-like, neutrons, etc...) without an MRD track to assist in the reconstruction. More info can be found on the included README, and in Andrew's slides.

It can be used on MC hits, MC waveform-reconstructed hits, and data hits. It is designed to run over clusters to yield a reconstructed vertex (x,y,z) and stdev to the timing residuals.

Checklist before submitting your PR

• This PR implements a single change (one new/modified Tool, or a set of changes to implement one new/modified feature)
• This PR alters the minimum number of files to affect this change
• If this PR includes a new Tool, a README and minimal demonstration ToolChain is provided
• If a new Tool/ToolChain requires model or configuration files, their paths are not hard-coded, and means of generating those files is described in the readme, with examples provided on /pnfs/annie/persistent
• For every new usage, there is a reason the data must be on the heap
• For every new there is a delete, unless I explicitly know why (e.g. ROOT or a BoostStore takes ownership)
  ^ store should take care of the new used in the tool (line 134, 135) but it would be nice to double check.

[jminock]

Did not spot any memory leaks or logical issues. Everything is well-documented (thank you!) and should not impede any other Tools or ToolChains

[marc1uk]

Code in general looks good, just a few places where we can reduce a lot of duplication with some templating or lambdas.

[marc1uk]

it's like reading code by @brichards64 😆

[marc1uk]

There's an if/else branch here of ~50 odd lines on each branch, but as far as I can see the only difference between them is whether to add the buffers to yMin, yMax, and max_radius. It seems you could just do

if (fExternalSeeding) {
   yMin -= fYBuffer;
   yMax x= fYBuffer;
   max_radius += fRadialBuffer;
}

and totally eliminate any code duplication.

[S81D]

Pushed this change. Will work on the others soon.

[marc1uk]

can you guess what I'm gonna say...? 🙃

This is wholly identical to EvalAtGuessVertex, other than our constant headache of taking vector<Hit> instead of vector<MCHit>. There are two simple ways to deal with this here:

1. make EvalAtGaussVertex a templated method. Only downside of this is the implementation needs to move to the header file.
2. use a lambda. Downside is that lambdas can't be members so wrapping is a bit awkward. One way could be as follows:

void VertexLeastSquares::EvalAtGuessVertex(const std::vector<Hit>& hits, const std::vector<MCHit>& mchits, ... ){
    auto process = [](auto& hits){
        // all current code in EvalAtGaussVertexMC can go here unchanged
    }
   
    if(hits) process(hits);
    else     process(mchits);
}

This keeps the EvalAtGuessVertex code within a member function without moving the implementation to the header, while taking advantage of lambdas' ability to accept auto to work with both types.
I'm happy to accept either (or something else), but please do one or the other.

[S81D]

I'll try these out, thank you for the recommendation. My apologies for much of the duplicated code - it's been built up from simply an MC implementation and it was easier at the time to add a separate block.

[marc1uk]

is a vertex necessary? can the map simply not have an entry?
with c++17 we could use std::optional.... but for now, perhaps as a minimum using an outside the tank position would be helpful.... 😕

[S81D]

I was thinking it follows the "blank" entries of -9999 established by the other MC tools (MCReco etc...). -99 in x,y,z is sufficiently far enough away from the tank.

[marc1uk]

while(true){
   ++loop_num;
   if(loop_num > max_loop_num) break;
}

.... so a for loop then? 😅

[marc1uk]

should this comment be 'reject this guess'

[S81D]

you're right - done.

[marc1uk]

another duplicate function

[marc1uk]

i'm assuming this can also be merged with FilterHits