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