geometry/richgeo/ActsGeo.cc

0001 // Copyright (C) 2022, 2023, Christopher Dilks
0002 // Subject to the terms in the LICENSE file found in the top-level directory.
0003
0004 #include "ActsGeo.h"
0005
0006 #include <DD4hep/Objects.h>
0007 #include <Evaluator/DD4hepUnits.h>
0008 #include <Math/GenVector/Cartesian3D.h>
0009 #include <Math/GenVector/DisplacementVector3D.h>
0010 #include <algorithm>
0011 #include <cmath>
0012 #include <edm4hep/Vector3f.h>
0013 #include <fmt/core.h>
0014 #include <utility>
0015
0016 #include "algorithms/tracking/TrackPropagationConfig.h"
0017 #include "services/geometry/richgeo/RichGeo.h"
0018
0019 // constructor
0020 richgeo::ActsGeo::ActsGeo(std::string detName_, gsl::not_null<const dd4hep::Detector*> det_,
0021                           std::shared_ptr<spdlog::logger> log_)
0022     : m_detName(std::move(detName_)), m_det(det_), m_log(std::move(log_)) {}
0023
0024 // generate list ACTS disc surfaces, for a given radiator
0025 std::vector<eicrecon::SurfaceConfig> richgeo::ActsGeo::TrackingPlanes(int radiator,
0026                                                                       int numPlanes) const {
0027
0028   // output list of surfaces
0029   std::vector<eicrecon::SurfaceConfig> discs;
0030
0031   // dRICH DD4hep-ACTS bindings --------------------------------------------------------------------
0032   if (m_detName == "DRICH") {
0033
0034     // vessel constants
0035     auto zmin  = m_det->constant<double>("DRICH_zmin");
0036     auto zmax  = m_det->constant<double>("DRICH_zmax");
0037     auto rmin0 = m_det->constant<double>("DRICH_rmin0");
0038     auto rmin1 = m_det->constant<double>("DRICH_rmin1");
0039     auto rmax0 = m_det->constant<double>("DRICH_rmax0");
0040     auto rmax1 = m_det->constant<double>("DRICH_rmax1");
0041     auto rmax2 = m_det->constant<double>("DRICH_rmax2");
0042
0043     // radiator constants
0044     auto snoutLength      = m_det->constant<double>("DRICH_snout_length");
0045     auto aerogelZpos      = m_det->constant<double>("DRICH_aerogel_zpos");
0046     auto aerogelThickness = m_det->constant<double>("DRICH_aerogel_thickness");
0047     auto filterZpos       = m_det->constant<double>("DRICH_filter_zpos");
0048     auto filterThickness  = m_det->constant<double>("DRICH_filter_thickness");
0049     auto window_thickness = m_det->constant<double>("DRICH_window_thickness");
0050
0051     // radial wall slopes
0052     auto boreSlope  = (rmin1 - rmin0) / (zmax - zmin);
0053     auto snoutSlope = (rmax1 - rmax0) / snoutLength;
0054
0055     // get z and radial limits where we will expect charged particles in the RICH
0056     double trackZmin = NAN;
0057     double trackZmax = NAN;
0058     std::function<double(double)> trackRmin;
0059     std::function<double(double)> trackRmax;
0060     switch (radiator) {
0061     case kAerogel:
0062       trackZmin = aerogelZpos - aerogelThickness / 2;
0063       trackZmax = aerogelZpos + aerogelThickness / 2;
0064       trackRmax = [&](auto z) { return rmax0 + snoutSlope * (z - zmin); };
0065       break;
0066     case kGas:
0067       trackZmin = filterZpos + filterThickness / 2;
0068       trackZmax = zmax - window_thickness;
0069       trackRmax = [&](auto z) {
0070         auto z0 = z - zmin;
0071         if (z0 < snoutLength) {
0072           return rmax0 + snoutSlope * z0;
0073         }
0074         return rmax2;
0075       };
0076       break;
0077     default:
0078       m_log->error("unknown radiator number {}", numPlanes);
0079       return discs;
0080     }
0081     trackRmin = [&](auto z) { return rmin0 + boreSlope * (z - zmin); };
0082
0083     // define discs: `numPlanes` z-equidistant planes *within* the radiator;
0084     /* NOTE: do not allow planes to be at radiator boundary
0085      * NOTE: alternative binning strategies do not seem to work well with the IRT algorithm;
0086      *       this one seems to work the best
0087      *
0088      * EXAMPLE: numPlanes=4
0089      *
0090      *    trackZmin         trackZmax
0091      *       :                   :
0092      *       :                   :
0093      *       +===================+....trackRmax
0094      *       [   |   |   |   |   ]
0095      *       [   |   |   |   |   ]
0096      *       [   <--planes--->   ]
0097      *       [   0   1   2   3   ]
0098      *       [   |   |   |   |   ]
0099      *       [   |   |   |   |   ]
0100      *       +===================+....trackRmin
0101      *       :   :       :   :
0102      *     ->:   :<-     :   :
0103      *     trackZstep    :   :
0104      *                 ->:   :<-
0105      *                 trackZStep
0106      */
0107     m_log->debug("Define ACTS disks for {} radiator: {} disks in z=[ {}, {} ]",
0108                  RadiatorName(radiator), numPlanes, trackZmin, trackZmax);
0109     double trackZstep = std::abs(trackZmax - trackZmin) / (numPlanes + 1);
0110     for (int i = 0; i < numPlanes; i++) {
0111       auto z    = trackZmin + (i + 1) * trackZstep;
0112       auto rmin = trackRmin(z);
0113       auto rmax = trackRmax(z);
0114       discs.emplace_back(eicrecon::DiscSurfaceConfig{"ForwardRICH_ID", z, rmin, rmax});
0115       m_log->debug("  disk {}: z={} r=[ {}, {} ]", i, z, rmin, rmax);
0116     }
0117   }
0118
0119   // pfRICH DD4hep-ACTS bindings --------------------------------------------------------------------
0120   else if (m_detName == "RICHEndcapN") {
0121     m_log->error("TODO: pfRICH DD4hep-ACTS bindings have not yet been implemented");
0122   }
0123
0124   // ------------------------------------------------------------------------------------------------
0125   else {
0126     m_log->error("ActsGeo is not defined for detector '{}'", m_detName);
0127   }
0128   return discs;
0129 }
0130
0131 // generate a cut to remove any track points that should not be used
0132 std::function<bool(edm4eic::TrackPoint)> richgeo::ActsGeo::TrackPointCut(int radiator) const {
0133
0134   // reject track points in dRICH gas that are beyond the dRICH mirrors
0135   // FIXME: assumes the full mirror spheres are much bigger than the dRICH
0136   // FIXME: needs to be generalized for dual or multi-mirror (per sector) design
0137   if (m_detName == "DRICH" && radiator == kGas) {
0138
0139     // get sphere centers
0140     std::vector<dd4hep::Position> mirror_centers;
0141     for (int isec = 0; isec < m_det->constant<int>("DRICH_num_sectors"); isec++) {
0142       mirror_centers.emplace_back(
0143           m_det->constant<double>("DRICH_mirror_center_x_sec" + std::to_string(isec)) / dd4hep::mm,
0144           m_det->constant<double>("DRICH_mirror_center_y_sec" + std::to_string(isec)) / dd4hep::mm,
0145           m_det->constant<double>("DRICH_mirror_center_z_sec" + std::to_string(isec)) / dd4hep::mm);
0146     }
0147     auto mirror_radius = m_det->constant<double>("DRICH_mirror_radius") / dd4hep::mm;
0148
0149     // beyond the mirror cut
0150     return [mirror_centers, mirror_radius](edm4eic::TrackPoint p) {
0151       return std::ranges::any_of(mirror_centers, [&p, &mirror_radius](const auto& c) {
0152         auto dist = std::hypot(c.x() - p.position.x, c.y() - p.position.y, c.z() - p.position.z);
0153         return dist < mirror_radius;
0154       });
0155     };
0156   }
0157
0158   // otherwise return a cut which always passes
0159   return [](edm4eic::TrackPoint /* p */) { return true; };
0160 }