HcalGeometry.cxx

#include "DetDescr/HcalGeometry.h"
 
#include <assert.h>
 
#include <algorithm>
#include <iostream>
 
#include "Framework/Exception/Exception.h"
 
namespace ldmx {
 
HcalGeometry::HcalGeometry(const framework::config::Parameters& ps)
    : framework::ConditionsObject(HcalGeometry::CONDITIONS_OBJECT_NAME) {
  scint_thickness_ = ps.get<double>("scint_thickness");
  scint_width_ = ps.get<double>("scint_width");
  zero_layer_ = ps.get<std::vector<double>>("zero_layer");
  layer_thickness_ = ps.get<std::vector<double>>("layer_thickness");
  num_layers_ = ps.get<std::vector<int>>("num_layers");
  num_sections_ = ps.get<int>("num_sections");
  ecal_dx_ = ps.get<double>("ecal_dx");
  ecal_dy_ = ps.get<double>("ecal_dy");
  verbose_ = ps.get<int>("verbose");
  back_horizontal_parity_ = ps.get<int>("back_horizontal_parity");
  side_3d_readout_ = ps.get<int>("side_3d_readout");
  y_offset_ = ps.get<double>("y_offset");
 
  auto detectors_valid = ps.get<std::vector<std::string>>("detectors_valid");
  // If one of the strings in detectors_valid is "ldmx-hcal-prototype", we
  // will use prototype geometry initialization
  is_prototype_ = std::find_if(detectors_valid.cbegin(), detectors_valid.cend(),
                               [](const auto detector) {
                                 return detector.find("ldmx-hcal-prototype") !=
                                        std::string::npos;
                               }) != detectors_valid.cend();
 
  num_strips_ = ps.get<std::vector<std::vector<int>>>("num_strips");
  half_total_width_ =
      ps.get<std::vector<std::vector<double>>>("half_total_width");
  zero_strip_ = ps.get<std::vector<std::vector<double>>>("zero_strip");
  scint_length_ = ps.get<std::vector<std::vector<double>>>("scint_length");
 
  // optional explicit layout, empty for the LDMX and prototype geometries
  layer_z_ = ps.get<std::vector<std::vector<double>>>("layer_z", {});
  layer_long_axis_ =
      ps.get<std::vector<std::vector<int>>>("layer_long_axis", {});
  layer_bar_center_ =
      ps.get<std::vector<std::vector<double>>>("layer_bar_center", {});
  layer_strip_zero_ =
      ps.get<std::vector<std::vector<double>>>("layer_strip_zero", {});
  double_ended_ = ps.get<std::vector<int>>("double_ended", {});
  if (hasExplicitLayout()) {
    for (int section = 0; section < num_sections_; ++section) {
      const auto n{static_cast<std::size_t>(num_layers_.at(section))};
      if (layer_z_.at(section).size() != n ||
          layer_long_axis_.at(section).size() != n ||
          layer_bar_center_.at(section).size() != n ||
          layer_strip_zero_.at(section).size() != n ||
          double_ended_.size() != static_cast<std::size_t>(num_sections_)) {
        EXCEPTION_RAISE(
            "HcalGeometry",
            "Explicit layout does not match num_layers in section " +
                std::to_string(section));
      }
    }
  }
 
  buildStripPositionMap();
 
  if (verbose_ > 0) {
    printPositionMap();
  }
}
double HcalGeometry::getDistanceAlongBar(ldmx::HcalID id, double x, double y,
                                         double z) const {
  const auto center{getStripCenterPosition(id)};
  switch (getScintillatorOrientation(id)) {
    case ScintillatorOrientation::horizontal:
      return x - center.X();
    case ScintillatorOrientation::vertical:
      return y - center.Y();
    default:
      return z - center.Z();
  }
}
 
std::vector<double> HcalGeometry::rotateGlobalToLocalBarPosition(
    const std::vector<double>& globalPosition, const ldmx::HcalID& id) const {
  const auto orientation{getScintillatorOrientation(id)};
  // explicit layout bars lie flat in z like the back HCal
  const unsigned int section{
      hasExplicitLayout()
          ? static_cast<unsigned int>(ldmx::HcalID::HcalSection::BACK)
          : id.section()};
  switch (section) {
    case ldmx::HcalID::HcalSection::BACK:
      switch (orientation) {
        case ScintillatorOrientation::horizontal:
          return {globalPosition[2], globalPosition[1], globalPosition[0]};
        case ScintillatorOrientation::vertical:
          return {globalPosition[2], globalPosition[0], globalPosition[1]};
        default:  // Should not be possible with current geometries
          EXCEPTION_RAISE("InvalidRotation",
                          "Attempted to rotate into an invalid "
                          "orientation for a scintillator bar!");
      }
    case ldmx::HcalID::HcalSection::TOP:
      [[fallthrough]];
    case ldmx::HcalID::HcalSection::BOTTOM:
      switch (orientation) {
        case ScintillatorOrientation::horizontal:
          return {globalPosition[1], globalPosition[2], globalPosition[0]};
        case ScintillatorOrientation::depth:
          return {globalPosition[1], globalPosition[0], globalPosition[2]};
        default:  // Should not be possible with current geometries
          EXCEPTION_RAISE("InvalidRotation",
                          "Attempted to rotate into an invalid "
                          "orientation for a scintillator bar!");
      }
    case ldmx::HcalID::HcalSection::LEFT:
      [[fallthrough]];
    case ldmx::HcalID::HcalSection::RIGHT:
      switch (orientation) {
        case ScintillatorOrientation::vertical:
          return {globalPosition[0], globalPosition[2], globalPosition[1]};
        case ScintillatorOrientation::depth:
          return globalPosition;
        default:  // Should not be possible with current geometries
          EXCEPTION_RAISE("InvalidRotation",
                          "Attempted to rotate into an invalid "
                          "orientation for a scintillator bar!");
      }
    default:
      // Can only reach this part if we somehow didn't match any of the options
      // above. This could happen if someone introduces a new geometry but
      // doesn't patch this part.
      EXCEPTION_RAISE("InvalidRotation",
                      "Attempted to rotate into an invalid "
                      "orientation for a scintillator bar!");
  }
}
 
HcalGeometry::ScintillatorOrientation HcalGeometry::getScintillatorOrientation(
    const ldmx::HcalID id) const {
  if (hasExplicitLayout()) {
    return layer_long_axis_.at(id.section()).at(id.layer() - 1) == 0
               ? ScintillatorOrientation::horizontal
               : ScintillatorOrientation::vertical;
  }
  if (hasSide3DReadout()) {
    // v14 or later detector
    switch (id.section()) {
      case ldmx::HcalID::HcalSection::TOP:
      case ldmx::HcalID::HcalSection::BOTTOM:
        // Odd layers are in z/depth direction, even are in the x/horizontal
        // direction
        return id.layer() % 2 == 0 ? ScintillatorOrientation::horizontal
                                   : ScintillatorOrientation::depth;
 
      case ldmx::HcalID::HcalSection::LEFT:
      case ldmx::HcalID::HcalSection::RIGHT:
        // Odd layers are in the z/depth direction, even are in the y/vertical
        // direction
        return id.layer() % 2 == 0 ? ScintillatorOrientation::vertical
                                   : ScintillatorOrientation::depth;
      case ldmx::HcalID::HcalSection::BACK:
        // Configurable
        return id.layer() % 2 == back_horizontal_parity_
                   ? ScintillatorOrientation::horizontal
                   : ScintillatorOrientation::vertical;
    }  // V14 or later detector
  }
  if (isPrototype()) {
    // The prototype only has the back section. However, the orientation
    // depends on the configuration so we delegate to the
    // back_horizontal_parity parameter
    return id.layer() % 2 == back_horizontal_parity_
               ? ScintillatorOrientation::horizontal
               : ScintillatorOrientation::vertical;
  }  // Prototype detector
  // v13/v12
  switch (id.section()) {
    // For the v13 side hcal, the bars in each section have the same
    // orientation
    case ldmx::HcalID::HcalSection::TOP:
    case ldmx::HcalID::HcalSection::BOTTOM:
      return ScintillatorOrientation::horizontal;
    case ldmx::HcalID::HcalSection::LEFT:
    case ldmx::HcalID::HcalSection::RIGHT:
      return ScintillatorOrientation::vertical;
    case ldmx::HcalID::HcalSection::BACK:
      // Configurable
      return id.layer() % 2 == back_horizontal_parity_
                 ? ScintillatorOrientation::horizontal
                 : ScintillatorOrientation::vertical;
  }  // v13/v12 detector
  // Can only reach this part if we somehow didn't match any of the options
  // above. This could happen if someone introduces a new geometry but doesn't
  // patch this part.
  EXCEPTION_RAISE("InvalidRotation",
                  "Attempted to rotate into an invalid "
                  "orientation for a scintillator bar!");
}
void HcalGeometry::printPositionMap(int section) const {
  // Note that layer numbering starts at 1 rather than 0
  for (int layer = 1; layer <= num_layers_[section]; ++layer) {
    for (int strip = 0; strip < getNumStrips(section, layer); ++strip) {
      HcalID id(section, layer, strip);
      auto center_position = getStripCenterPosition(id);
      auto x = center_position.X();
      auto y = center_position.Y();
      auto z = center_position.Z();
      std::cout << id << ": Center position: (" << x << ", " << y << ", " << z
                << ")\n";
    }
  }
}
 
void HcalGeometry::buildStripPositionMap() {
  // We hard-code the number of sections as seen in HcalID
  for (unsigned int section = 0; section < num_sections_; section++) {
    for (unsigned int layer = 1; layer <= num_layers_[section]; layer++) {
      for (unsigned int strip = 0; strip < getNumStrips(section, layer);
           strip++) {
        // initialize values
        double x{-99999}, y{-99999}, z{-99999};
 
        // get hcal section
        ldmx::HcalID::HcalSection hcalsection =
            (ldmx::HcalID::HcalSection)section;
 
        const ldmx::HcalID id{section, layer, strip};
        const auto orientation{getScintillatorOrientation(id)};
        // the center of a layer: (layer-1) * (layer_thickness) +
        // scint_thickness/2
        double layercenter =
            (layer - 1) * layer_thickness_.at(section) + 0.5 * scint_thickness_;
 
        // the center of a strip: (strip + 0.5) * (strip_dx)
        double stripcenter = (strip + 0.5) * scint_width_;
 
        if (hasExplicitLayout()) {
          const auto i_layer{layer - 1};
          const double along{layer_bar_center_.at(section).at(i_layer)};
          const double across{layer_strip_zero_.at(section).at(i_layer) +
                              strip * scint_width_};
          if (orientation == ScintillatorOrientation::horizontal) {
            x = along;
            y = across;
          } else {
            x = across;
            y = along;
          }
          z = layer_z_.at(section).at(i_layer);
        } else if (hcalsection == ldmx::HcalID::HcalSection::BACK) {
          // z position: zero-layer(z) + layer_z + scint_thickness / 2
          z = zero_layer_.at(section) + layercenter;
 
          if (orientation == ScintillatorOrientation::horizontal) {
            y = stripcenter - getZeroStrip(section, layer);
            x = 0;
          } else {
            x = stripcenter - getZeroStrip(section, layer);
            y = 0;
          }
        } else {
          if (side_3d_readout_) {
            /*
             *
             * For 3D readout:
             * - odd layers have strips in z
             * - even layers have strips in x(y) for top-bottom (left-right)
             * sections
             * - odd layers have strips occupying width of scintillator in x(y)
             * - even layers have strips occupying width of scintillator in z
             *
             */
            switch (hcalsection) {
              case ldmx::HcalID::HcalSection::BACK:
                // Handled earlier in the code!
              case ldmx::HcalID::HcalSection::LEFT:
              case ldmx::HcalID::HcalSection::RIGHT:
                if (orientation == ScintillatorOrientation::vertical) {
                  x = zero_layer_[section] + 0.5 * scint_thickness_ +
                      (layer - 1) * layer_thickness_[section];
                  y = ecal_dy_ -
                      (getScintillatorLength({id.section(), 2, id.strip()}) -
                       getScintillatorLength(id)) /
                          2;
                  z = getZeroStrip(section, layer) +
                      (strip + 0.5) * getScintillatorWidth();
                } else if (orientation == ScintillatorOrientation::depth) {
                  x = zero_layer_[section] + 0.5 * scint_thickness_ +
                      layer_thickness_[section] * (layer - 1);
                  y = -ecal_dy_ / 2 + (strip + 0.5) * getScintillatorWidth();
                  z = getZeroStrip(section, layer + 1) +
                      getScintillatorLength(id) / 2;
                }
                if (section == ldmx::HcalID::HcalSection::LEFT) {
                  y *= -1;
                  x *= -1;
                }
                break;
 
              case ldmx::HcalID::HcalSection::BOTTOM:
              case ldmx::HcalID::HcalSection::TOP:
                if (orientation == ScintillatorOrientation::horizontal) {
                  //
                  // Second half of the expression is the difference between the
                  // longest strips (first module) and the current module.
                  //
                  // 22 mm extra for space for 1 absorber and one air box
                  x = -ecal_dx_ / 2 - 2 - 20 +
                      (getScintillatorLength({id.section(), 2, id.strip()}) -
                       getScintillatorLength(id)) /
                          2;
                  y = zero_layer_[section] + 0.5 * scint_thickness_ +
                      (layer - 1) * layer_thickness_[section];
                  z = getZeroStrip(section, layer) +
                      (strip + 0.5) * getScintillatorWidth();
                }
                if (orientation == ScintillatorOrientation::depth) {
                  x = (ecal_dx_ / 2) - (strip + 0.5) * getScintillatorWidth();
                  y = zero_layer_[section] + 0.5 * scint_thickness_ +
                      layer_thickness_[section] * (layer - 1);
                  z = getZeroStrip(section, layer + 1) +
                      getScintillatorLength(id) / 2;
                }
                if (section == ldmx::HcalID::HcalSection::BOTTOM) {
                  y *= -1;
                  x *= -1;
                }
                break;
            }
 
          } else {
            // z position: zero-strip(z) + strip_center(z)
            z = getZeroStrip(section, layer) + stripcenter;
            if (hcalsection == ldmx::HcalID::HcalSection::TOP or
                hcalsection == ldmx::HcalID::HcalSection::BOTTOM) {
              y = zero_layer_.at(section) + layercenter;
              x = getHalfTotalWidth(section, layer);
              if (hcalsection == ldmx::HcalID::HcalSection::BOTTOM) {
                y *= -1;
                x *= -1;
              }
 
            } else {
              x = zero_layer_.at(section) + layercenter;
              y = getHalfTotalWidth(section, layer);
              if (hcalsection == ldmx::HcalID::HcalSection::RIGHT) {
                x *= -1;
                y *= -1;
              }
            }
          }
        }
 
        y += y_offset_;
        ROOT::Math::XYZVector pos;
        pos.SetXYZ(x, y, z);
        strip_position_map_[ldmx::HcalID(section, layer, strip)] = pos;
      }  // loop over strips
    }  // loop over layers
  }  // loop over sections
}  // strip position map
 
}  // namespace ldmx