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/format.h>
0012 #include <algorithm>
0013 #include <cmath>
0014 #include <utility>
0015 #include <variant>
0016
0017 #include "algorithms/tracking/TrackPropagationConfig.h"
0018 #include "services/geometry/richgeo/RichGeo.h"
0019
0020 // constructor
0021 richgeo::ActsGeo::ActsGeo(std::string detName_, gsl::not_null<const dd4hep::Detector*> det_,
0022                           std::shared_ptr<spdlog::logger> log_)
0023     : m_detName(std::move(detName_)), m_det(det_), m_log(std::move(log_)) {}
0024
0025 // generate list ACTS disc surfaces, for a given radiator
0026 std::vector<eicrecon::SurfaceConfig> richgeo::ActsGeo::TrackingPlanes(int radiator,
0027                                                                       int numPlanes) const {
0028
0029   // output list of surfaces
0030   std::vector<eicrecon::SurfaceConfig> discs;
0031
0032   // dRICH DD4hep-ACTS bindings --------------------------------------------------------------------
0033   if (m_detName == "DRICH") {
0034
0035     // vessel constants
0036     auto zmin  = m_det->constant<double>("DRICH_zmin");
0037     auto zmax  = m_det->constant<double>("DRICH_zmax");
0038     auto rmin0 = m_det->constant<double>("DRICH_rmin0");
0039     auto rmin1 = m_det->constant<double>("DRICH_rmin1");
0040     auto rmax0 = m_det->constant<double>("DRICH_rmax0");
0041     auto rmax1 = m_det->constant<double>("DRICH_rmax1");
0042     auto rmax2 = m_det->constant<double>("DRICH_rmax2");
0043
0044     // radiator constants
0045     auto snoutLength      = m_det->constant<double>("DRICH_snout_length");
0046     auto aerogelZpos      = m_det->constant<double>("DRICH_aerogel_zpos");
0047     auto aerogelThickness = m_det->constant<double>("DRICH_aerogel_thickness");
0048     auto filterZpos       = m_det->constant<double>("DRICH_filter_zpos");
0049     auto filterThickness  = m_det->constant<double>("DRICH_filter_thickness");
0050     auto window_thickness = m_det->constant<double>("DRICH_window_thickness");
0051
0052     // radial wall slopes
0053     auto boreSlope  = (rmin1 - rmin0) / (zmax - zmin);
0054     auto snoutSlope = (rmax1 - rmax0) / snoutLength;
0055
0056     // get z and radial limits where we will expect charged particles in the RICH
0057     double trackZmin = NAN;
0058     double trackZmax = NAN;
0059     std::function<double(double)> trackRmin;
0060     std::function<double(double)> trackRmax;
0061     switch (radiator) {
0062     case kAerogel:
0063       trackZmin = aerogelZpos - aerogelThickness / 2;
0064       trackZmax = aerogelZpos + aerogelThickness / 2;
0065       trackRmax = [&](auto z) { return rmax0 + snoutSlope * (z - zmin); };
0066       break;
0067     case kGas:
0068       trackZmin = filterZpos + filterThickness / 2;
0069       trackZmax = zmax - window_thickness;
0070       trackRmax = [&](auto z) {
0071         auto z0 = z - zmin;
0072         if (z0 < snoutLength) {
0073           return rmax0 + snoutSlope * z0;
0074         }
0075         return rmax2;
0076       };
0077       break;
0078     default:
0079       m_log->error("unknown radiator number {}", numPlanes);
0080       return discs;
0081     }
0082     trackRmin = [&](auto z) { return rmin0 + boreSlope * (z - zmin); };
0083
0084     // define discs: `numPlanes` z-equidistant planes *within* the radiator;
0085     /* NOTE: do not allow planes to be at radiator boundary
0086      * NOTE: alternative binning strategies do not seem to work well with the IRT algorithm;
0087      *       this one seems to work the best
0088      *
0089      * EXAMPLE: numPlanes=4
0090      *
0091      *    trackZmin         trackZmax
0092      *       :                   :
0093      *       :                   :
0094      *       +===================+....trackRmax
0095      *       [   |   |   |   |   ]
0096      *       [   |   |   |   |   ]
0097      *       [   <--planes--->   ]
0098      *       [   0   1   2   3   ]
0099      *       [   |   |   |   |   ]
0100      *       [   |   |   |   |   ]
0101      *       +===================+....trackRmin
0102      *       :   :       :   :
0103      *     ->:   :<-     :   :
0104      *     trackZstep    :   :
0105      *                 ->:   :<-
0106      *                 trackZStep
0107      */
0108     m_log->debug("Define ACTS disks for {} radiator: {} disks in z=[ {}, {} ]",
0109                  RadiatorName(radiator), numPlanes, trackZmin, trackZmax);
0110     double trackZstep = std::abs(trackZmax - trackZmin) / (numPlanes + 1);
0111     for (int i = 0; i < numPlanes; i++) {
0112       auto z    = trackZmin + (i + 1) * trackZstep;
0113       auto rmin = trackRmin(z);
0114       auto rmax = trackRmax(z);
0115       discs.emplace_back(eicrecon::DiscSurfaceConfig{
0116           .id = "ForwardRICH_ID", .zmin = z, .rmin = rmin, .rmax = 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 }