File indexing completed on 2024-09-28 07:03:04
0001
0002
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
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
0024 std::transform(m_detName.begin(), m_detName.end(), m_detName.begin(), ::toupper);
0025 }
0026
0027
0028 std::vector<eicrecon::SurfaceConfig> richgeo::ActsGeo::TrackingPlanes(int radiator, int numPlanes) {
0029
0030
0031 std::vector<eicrecon::SurfaceConfig> discs;
0032
0033
0034 if(m_detName=="DRICH") {
0035
0036
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
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
0054 auto boreSlope = (rmin1 - rmin0) / (zmax - zmin);
0055 auto snoutSlope = (rmax1 - rmax0) / snoutLength;
0056
0057
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
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
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
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
0128 std::function<bool(edm4eic::TrackPoint)> richgeo::ActsGeo::TrackPointCut(int radiator) {
0129
0130
0131
0132
0133 if(m_detName=="DRICH" && radiator==kGas) {
0134
0135
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
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
0161 return [] (edm4eic::TrackPoint p) { return true; };
0162
0163 }