Updates to the hcal prototype geometry

[EinarElen]

The current ldmx-hcal-prototype-v1 geometry differs from what was actually used at the test beam so we need to add a couple of new GDML descriptions (corresponding to different configurations that were used) which also brings along some necessary updates to LDMX-sw (mostly for the Hcal submodule but also in DetDescr).

• The ordering of horizontal and vertical layers in the prototype differs from the regular Hcal
• Some measurements were carried out with the horizontal bars shifted which means that the center of the bars cannot be assumed to be at x/y = 0

[EinarElen]

@bryngemark I've pushed a version which has @PeterGy's latest work on https://github.com/LDMX-Software/ldmx-sw/tree/iss1063 in the Detectors/data/ldmx-hcal-prototype-v2.0 directory. There are some different versions of the TS available in the trigger_physvols file, but everything except plastic should be commented out by default.

from LDMX.Framework import ldmxcfg
from LDMX.SimCore import generators
from LDMX.SimCore import simulator


process = ldmxcfg.Process('test_TS')

process.maxEvents = 1000
process.run = 10
process.outputFiles = ['test_TS.root']

gun = generators.gun('particle_gun')
gun.particle = 'e-'
gun.direction = [0., 0., 1.]
gun.position = [0., 0., -1000.]
gun.energy = 4.

simulation = simulator.simulator('test_TS')
simulation.generators=[gun]
simulation.setDetector('ldmx-hcal-prototype-v2.0')

from LDMX.TrigScint.trigScint import TrigScintQIEDigiProducer
from LDMX.TrigScint.trigScint import TrigScintRecHitProducer

tsDigis = TrigScintQIEDigiProducer.up()
tsDigis.number_of_strips = 12
tsRecHits = TrigScintRecHitProducer.up()


process.sequence=[simulation,
                  tsDigis,
                  tsRecHits
                  ]

I stole this from what @PeterGy has been using and it does something although you and Peter are probably better than me at knowing if it is sensible :)

[bryngemark]

very helpful, thanks @EinarElen! this looks reasonable.

[EinarElen]

I've uploaded the basic configuration (plastic TS, no shifted bars) of the April testbeam to https://github.com/LDMX-Software/ldmx-sw/tree/iss1063 together with an updated DetDescr readout geometry. See updated test config below.

However, if you try to run this by default you will crash when running the HcalRecProducer. This is because the ordering of the layers has changed. In the April testbeam, odd layers were vertical while in the standard geometries odd layers in the Hcal are horizontal. To fix this, there are two places (at least) that need to be able to take this into account.

In
https://github.com/LDMX-Software/SimCore/blob/1fd615f297b6daad5a80ba41e29596f1c64cf633/src/SimCore/HcalSD.cxx#L136-L146

  int stripID = -1;
  // Odd layers have bars horizontal
  // Even layers have bars vertical
  // 5cm wide bars are HARD-CODED
  if (section == ldmx::HcalID::BACK && layer % 2 == 1)
    stripID = int((localPosition.y() + scint->GetYHalfLength()) / 50.0);
  else if (section == ldmx::HcalID::BACK && layer % 2 == 0)
    stripID = int((localPosition.x() + scint->GetXHalfLength()) / 50.0);
  else
    stripID = int((localPosition.z() + scint->GetZHalfLength()) / 50.0);

in https://github.com/LDMX-Software/Hcal/blob/49ee21d7634d353ee538470aaa4d89d1d2cebccf/src/Hcal/HcalRecProducer.cxx#L254-L259

      // set position along the bar
      if ((id_posend.layer() % 2) == 1) {
        position.SetX(position_bar);
      } else {
        position.SetY(position_bar);
      }

in

ldmx-sw/DetDescr/src/DetDescr/HcalGeometry.cxx

Lines 110 to 116 in f49a7f8

	if (layer % 2 == 1) {
	y = stripcenter - getZeroStrip(section, layer);
	x = 0;
	} else {
	x = stripcenter - getZeroStrip(section, layer);
	y = 0;
	}

          /**
            Now compute, y(x) position for even(odd) layers, relative to the
            center of detector. Strips enumeration starts from -y(-x)
            stripcenter will be large for +y(+x) and the halfwidth of the strip
            needs to be subtracted The halfwidth of the scintillator is given by
            ZeroStrip_. The x(y) position is set to the center of the strip (0).
          */
          if (layer % 2 == 1) {
            y = stripcenter - getZeroStrip(section, layer);
            x = 0;
          } else {
            x = stripcenter - getZeroStrip(section, layer);
            y = 0;
          }

and in https://github.com/LDMX-Software/Hcal/blob/49ee21d7634d353ee538470aaa4d89d1d2cebccf/src/Hcal/HcalDigiProducer.cxx#L171

        distance_along_bar = (layer % 2) ? position[0] : position[1];

If you want to start testing running simulations with what's currently in this branch, you can just manually patch these three places and you should (probably) be fine to start testing things.

The easiest way to support these geometry differences that I can think of would be to rely on the HcalGeometry condition as the single source of truth. Then you could just have a member function in https://github.com/LDMX-Software/ldmx-sw/blob/trunk/DetDescr/include/DetDescr/HcalGeometry.h along the lines of (although maybe with some more elegant but equivalent expression that I couldn't think of right now)

bool layerIsHorizontal(const int layer) const { 
    if (someRuntimeConfig) { 
     return layer % 2 == 0; 
     } else { 
     return layer % 2 == 1;
     }
}

then the places that currently has the layer ordering hard coded could rely on the geometry condition and do something along the lines of

      // set position along the bar
      if (hcalGeometry.layerIsHorizontal(id_posend.layer())) {
        position.SetX(position_bar);
      } else {
        position.SetY(position_bar);
      }

For the Hcal parts, this trivial to do. However, this would add a new dependency for SimCore on import LDMX.Hcal.HcalGeometry. Maybe that's ok? There's already a dependence on the corresponding the Ecal geometry in EcalSD.

from LDMX.Framework import ldmxcfg
from LDMX.SimCore import generators
from LDMX.SimCore import simulator


process = ldmxcfg.Process('test')

process.maxEvents = 1000
process.run = 10
process.outputFiles = ['test.root']

gun = generators.gun('particle_gun')
gun.particle = 'e-'
gun.direction = [0., 0., 1.]
gun.position = [0., 0., -1000.]
gun.energy = 4.

simulation = simulator.simulator('test_TS')
simulation.generators=[gun]
simulation.setDetector('ldmx-hcal-prototype-v2.0')

from LDMX.TrigScint.trigScint import TrigScintQIEDigiProducer
from LDMX.TrigScint.trigScint import TrigScintRecHitProducer

tsDigis = TrigScintQIEDigiProducer.up()
tsDigis.number_of_strips = 12
tsRecHits = TrigScintRecHitProducer.up()


import LDMX.Hcal.HcalGeometry
from LDMX.Hcal import digi
from LDMX.Hcal import hcal
from LDMX.Hcal import hcal_hardcoded_conditions

digis = digi.HcalDigiProducer()
reco = digi.HcalRecProducer()
process.sequence=[simulation,
                  digis,
                  reco,
                  tsDigis,
                  tsRecHits
                  ]

[EinarElen]

I went ahead and tested doing things this way and it wasn't particularly complicated to do. So if it makes sense, I can go ahead and implement it.

[bryngemark]

Another option is to have the parity of the layers configurable, as in horizontalParity_=1, and then check,

if layer % 2 == horizontalParity_
...

Could still be set in the python script or the conditions but you wouldn't have to check for which configuration you're using.

[EinarElen]

Oh yeah that's much better. The question of putting it as part of the condition or as a part of the configuration is hard, i would lean towards putting it in the condition since, in principle at least, once you choose your geometry it would be a bug to select the wrong parameter. On the other hand though, having the partity change invisibly when you change detector geometry might be the kind of thing that makes someone pull their hair out if it isn't obvious