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 <ctype.h>
0011 #include <edm4hep/Vector3f.h>
0012 #include <fmt/core.h>
0013 #include <algorithm>
0014 #include <cmath>
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_, std::shared_ptr<spdlog::logger> log_)
0021   : m_detName(detName_), m_det(det_), m_log(log_)
0022 {
0023   // capitalize m_detName
0024   std::transform(m_detName.begin(), m_detName.end(), m_detName.begin(), ::toupper);
0025 }
0026
0027 // generate list ACTS disc surfaces, for a given radiator
0028 std::vector<eicrecon::SurfaceConfig> richgeo::ActsGeo::TrackingPlanes(int radiator, int numPlanes) {
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, trackZmax;
0059     std::function<double(double)> trackRmin, 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) return rmax0 + snoutSlope * z0;
0072           else return rmax2;
0073         };
0074         break;
0075       default:
0076         m_log->error("unknown radiator number {}",numPlanes);
0077         return discs;
0078     }
0079     trackRmin = [&] (auto z) { return rmin0 + boreSlope * (z - zmin); };
0080
0081     // define discs: `numPlanes` z-equidistant planes *within* the radiator;
0082     /* NOTE: do not allow planes to be at radiator boundary
0083      * NOTE: alternative binning strategies do not seem to work well with the IRT algorithm;
0084      *       this one seems to work the best
0085      *
0086      * EXAMPLE: numPlanes=4
0087      *
0088      *    trackZmin         trackZmax
0089      *       :                   :
0090      *       :                   :
0091      *       +===================+....trackRmax
0092      *       [   |   |   |   |   ]
0093      *       [   |   |   |   |   ]
0094      *       [   <--planes--->   ]
0095      *       [   0   1   2   3   ]
0096      *       [   |   |   |   |   ]
0097      *       [   |   |   |   |   ]
0098      *       +===================+....trackRmin
0099      *       :   :       :   :
0100      *     ->:   :<-     :   :
0101      *     trackZstep    :   :
0102      *                 ->:   :<-
0103      *                 trackZStep
0104      */
0105     m_log->debug("Define ACTS disks for {} radiator: {} disks in z=[ {}, {} ]",
0106         RadiatorName(radiator), numPlanes, trackZmin, trackZmax);
0107     double trackZstep = std::abs(trackZmax-trackZmin) / (numPlanes+1);
0108     for(int i=0; i<numPlanes; i++) {
0109       auto z         = trackZmin + (i+1)*trackZstep;
0110       auto rmin      = trackRmin(z);
0111       auto rmax      = trackRmax(z);
0112       discs.push_back(eicrecon::DiscSurfaceConfig{"ForwardRICH_ID", z, rmin, rmax});
0113       m_log->debug("  disk {}: z={} r=[ {}, {} ]", i, z, rmin, rmax);
0114     }
0115   }
0116
0117   // pfRICH DD4hep-ACTS bindings --------------------------------------------------------------------
0118   else if(m_detName=="PFRICH") {
0119     m_log->error("TODO: pfRICH DD4hep-ACTS bindings have not yet been implemented");
0120   }
0121
0122   // ------------------------------------------------------------------------------------------------
0123   else m_log->error("ActsGeo is not defined for detector '{}'",m_detName);
0124   return discs;
0125 }
0126
0127 // generate a cut to remove any track points that should not be used
0128 std::function<bool(edm4eic::TrackPoint)> richgeo::ActsGeo::TrackPointCut(int radiator) {
0129
0130   // reject track points in dRICH gas that are beyond the dRICH mirrors
0131   // FIXME: assumes the full mirror spheres are much bigger than the dRICH
0132   // FIXME: needs to be generalized for dual or multi-mirror (per sector) design
0133   if(m_detName=="DRICH" && radiator==kGas) {
0134
0135     // get sphere centers
0136     std::vector<dd4hep::Position> mirror_centers;
0137     for(int isec = 0; isec < m_det->constant<int>("DRICH_num_sectors"); isec++)
0138       mirror_centers.emplace_back(
0139           m_det->constant<double>("DRICH_mirror_center_x_sec" + std::to_string(isec)) / dd4hep::mm,
0140           m_det->constant<double>("DRICH_mirror_center_y_sec" + std::to_string(isec)) / dd4hep::mm,
0141           m_det->constant<double>("DRICH_mirror_center_z_sec" + std::to_string(isec)) / dd4hep::mm
0142           );
0143     auto mirror_radius = m_det->constant<double>("DRICH_mirror_radius") / dd4hep::mm;
0144
0145     // beyond the mirror cut
0146     return [mirror_centers, mirror_radius] (edm4eic::TrackPoint p) {
0147       for(const auto& c : mirror_centers) {
0148         auto dist = std::hypot(
0149             c.x() - p.position.x,
0150             c.y() - p.position.y,
0151             c.z() - p.position.z
0152             );
0153         if(dist < mirror_radius)
0154           return true;
0155       }
0156       return false;
0157     };
0158   }
0159
0160   // otherwise return a cut which always passes
0161   return [] (edm4eic::TrackPoint p) { return true; };
0162
0163 }