Support for the testbeam prototype geometry

[EinarElen]

Hi everyone,

I’ve been working on simulations for the testbeam prototype, which requires some modifications to LDMX-sw to support the prototype geometry, specifically in the HCal reconstruction/Detdescr geometry code. The main difference between the back Hcal and the prototype is that the number of bars (strips) depends on the layer number and that the width and height of a given scintillator layer are not going to be the same.

Currently, I am writing code with the assumption that it could be something that you would want to merge into mainline LDMX-sw, which means that I have to think about e.g. backwards compatibility with previous v3 config files. Alternatively, we could keep the prototype code on a dedicated branch indefinitely. I’m not sure what would be best.

Most of the changes that have to be done are in the DetDescr/HcalGeometry module. Functions like getZeroStrip, getNumStrips, and getHalfTotalWidth will need to return different values depending on the layer number, and functions like buildStripPositionMap need to be updated to make use of the new output of those functions.

We will also need to provide the C++ code with parameters describing the difference between layers. DetDescr/HcalGeometry.py already supports adding multiple geometry configurations that are selected at runtime based on the DetectorName variable in the detector GDML file, so adding a new geometry here in addition to the v12 geometry is straight-forward. However, the existing parameters don’t really work for this purpose.

Currently, most parameters are in the form of a vector where the indices represent the 5 sections of the LDMX Hcal. Since LDMX-sw doesn’t support parameters that are vectors of vectors, you can’t just repurpose these variables to be a vector of vectors where the first index represents the section and the second index represents the layer. If we could have that, we could in principle let the C++-layer figure out whether the parameter that was passed in was a regular vector or a vector of vectors, and thus handle both types detector geometries seamlessly.

I could have a go at adding such a feature myself, it does seem non-trivial but maybe it would be a feature that others would find useful.

Of course, if we were to keep the prototype branch separated indefinitely, we could just change the purpose of whatever variables we would see fit. Otherwise, the solutions I’ve been able to think of (and have used so far in my own work) are typically a lot uglier/verbose.

Finally, I would want to add some debugging support for the strip position map. Just enough to be useful for verifying that the strips are ending up at the positions you would expect.

[tomeichlersmith]

Hi Einar! This is definitely a big project, so thank you for starting this issue to help move it along!

You've already noticed the DetectorName variable in the configuration which allows the user to change the paramete values depending on the name of the geometry. Another feature I want to point out is that the set of parameters between different detector names don't need to be the same. This fact is because under-the-hood we are basically just keeping track of a Python dictionary which is converted to an instance of framework::config::Parameters at run-time. Nobody else has used this before, but it seems to definitely come in handy for your use case. This also avoids the issue of trying to translate a list of lists in Python to a vector of vectors in C++ (which absolutely begs for segmentation violations).

In order to change the parameters, you will need to maintain a certain set of consistent parameters that are required by the geometry-provider structure. I don't know this set off the top of my head, but detectors_valid is almost certainly one of them. This would involve modifying the python code to something like

class HcalReadoutGeometryBase :
    def __init__(self) :
        self.detectors_valid = []

    def __str__(self) :
        msg = f'HcalReadoutGeometry for {self.detectors_valid}'

class HcalReadoutGeometryFull(HcalReadoutGeometryBase) :
    def __init__(self) :
        super().__init__()
        # insert full geometry parameters here

    def __str__(self) :
        # old printout

class HcalReadoutGeometryTestBeam(HcalReadoutGeometryBase) :
    def __init__(self) :
        super().__init__()
        # insert test-beam geometry parameters here

    def __str__(self) :
        # printout parameters

class HcalGeometry :
    # class implementation ommitted for brevity

And then you will want to implement a make_testbeam similar to how make_v12 is implemented in HcalGeometry. I advocate that this is the preferred method so that we can have multiple geometries compiled together, making it much simpler to run comparisons.

This type of inheritance tree in the Python configuration could be mimic'd in the C++ by having the HcalGeometry class create a different sub-class implementation depending on the detectors_valid parameters passed. This would avoid lots of if...else trees elsewhere but may be more effort than it is worth. A currently-active PR #1003 uses a flag value and if...else trees if you want to see an example.

For some debugging support, I would suggest writing a simple "printer" that will print-out the strip/layer position map. We have done this in the ECal and it has proven to be extremely useful. This debugging is only helpful for developers since it can only be run manually and creates images that need to be viewed manually, but that is still a huge help when changing the geometry drastically.

This is already a lot of information, so I'm going to stop rambling now. I would be open to meeting with you if you want to have a more in-depth discussion.

[EinarElen]

Something I just realized I should have made a bit more clear, I have essentially already created a working version of this so that we could make use of the improvements that came with v3. So other than choosing strategies/design, most of the work of has already been done!

Hi Einar! This is definitely a big project, so thank you for starting this issue to help move it along!

You've already noticed the DetectorName variable in the configuration which allows the user to change the paramete values depending on the name of the geometry. Another feature I want to point out is that the set of parameters between different detector names don't need to be the same. This fact is because under-the-hood we are basically just keeping track of a Python dictionary which is converted to an instance of framework::config::Parameters at run-time. Nobody else has used this before, but it seems to definitely come in handy for your use case. This also avoids the issue of trying to translate a list of lists in Python to a vector of vectors in C++ (which absolutely begs for segmentation violations).

I can imagine how this could be a bit dangerous, yes. What you are suggesting is essentially what i have been doing so far. It works quite well, the only downside I can think of is that it makes the initialization code on the C++ side a bit messier. It's nothing horrible and you can add some helper functions that reduce the repetition involved. Although, as you mention later you could make use of the detectors_valid parameter here. I'm not certain that would be shorter and/or more clear but it would be worth trying out.

In order to change the parameters, you will need to maintain a certain set of consistent parameters that are required by the geometry-provider structure. I don't know this set off the top of my head, but detectors_valid is almost certainly one of them. This would involve modifying the python code to something like

class HcalReadoutGeometryBase :
    def __init__(self) :
        self.detectors_valid = []

    def __str__(self) :
        msg = f'HcalReadoutGeometry for {self.detectors_valid}'

class HcalReadoutGeometryFull(HcalReadoutGeometryBase) :
    def __init__(self) :
        super().__init__()
        # insert full geometry parameters here

    def __str__(self) :
        # old printout

class HcalReadoutGeometryTestBeam(HcalReadoutGeometryBase) :
    def __init__(self) :
        super().__init__()
        # insert test-beam geometry parameters here

    def __str__(self) :
        # printout parameters

class HcalGeometry :
    # class implementation ommitted for brevity

And then you will want to implement a make_testbeam similar to how make_v12 is implemented in HcalGeometry. I advocate that this is the preferred method so that we can have multiple geometries compiled together, making it much simpler to run comparisons.

This is essentially what I have been doing, although I think that the version I have doesn't make use of inheritance.

This type of inheritance tree in the Python configuration could be mimic'd in the C++ by having the HcalGeometry class create a different sub-class implementation depending on the detectors_valid parameters passed. This would avoid lots of if...else trees elsewhere but may be more effort than it is worth. A currently-active PR #1003 uses a flag value and if...else trees if you want to see an example.

The method I've been using so far which I think could be a good alternative is to behind the scenes use the parameters for either geometry to set up a generic

Essentially, if you had a variable like NumStrips in the Hcal config representing the number of strips for the different sections and another in the prototype config representing the number of strips for each layer, looking something like this

# Somewhere in the Hcal config
 self.v12.NumStrips = [62,12,12,12,12]
 self.v12.NumLayers = [100,28,28,26,26]
 # Somewhere in the prototype config 
self.prototype.NumStrips = [8 for i in range(9)] + [12 for i in range(10)]

On the C++ side, you can represent both of these cases as a vector of vectors (or some other container of vectors/maybe even have different variables for the different sections) that you set up differently depending on which geometry you are provided. So in this case you would end up with something like this (although not initialized directly like here)

// Hcal array
NumStrips = [[62,62, 62, ..., 62],
                      [28, 28, 28, ..., 28],
                      ... 
                      ]
// Prototype array
NumStrips = [[ 8, 8, 8, 8, ..., 12, 12, 12],
                      [],[],[],[]
                      ]

One caveat that you might notice from the code above is that the vectors representing the other four Hcal sections will be empty vectors, so you could get into trouble if you were to do something like

int section = 2
int layer = 0
NumStrips[section][layer]

The rest of the code, primarily the code setting up the strip position map, can then have a single abstraction to work with. So you can loop over the number of layers the same way, the Hcal reconstruction code can happily ask how many strips are in a layer or what the midpoint of a bar with a given HcalID is, without the underlying code being filled with if/else statements that can be hard to read/end up becoming out of date.

I'm less worried about the segfault risk here than for allowing it in the parameters. In fact, this could allow you to use more safe looping mechanisms

// Without 
for (int i = 0; i < NumLayers; ++i) 
// maybe this access is out of bounds 
NumStrips[i]
// With
for (int i = 0; i < NumStrips[section].size(); ++i) 
// or 
for (auto strip : NumStrips[section]) 

You could of course still use a section number that isn't in bounds, but then you would be in trouble elsewhere in the code already.

For some debugging support, I would suggest writing a simple "printer" that will print-out the strip/layer position map. We have done this in the ECal and it has proven to be extremely useful. This debugging is only helpful for developers since it can only be run manually and creates images that need to be viewed manually, but that is still a huge help when changing the geometry drastically.

This is a lot fancier than what I have been having in my mind/my current repository, which has been more along the lines of a printer-style debugging facility. But at some point implementing something like what you have in the Ecal module would be very useful.

This is already a lot of information, so I'm going to stop rambling now. I would be open to meeting with you if you want to have a more in-depth discussion.

I've tried rambling a bit back to you here. Meeting sounds great, maybe we could find a time on slack?

[tomeichlersmith]

Merged this branch in manually.