geometry/richgeo/IrtGeo.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 "IrtGeo.h"
0005
0006 #include <DD4hep/Volumes.h>
0007 #include <Evaluator/DD4hepUnits.h>
0008 #include <IRT/CherenkovRadiator.h>
0009 #include <Math/GenVector/DisplacementVector3D.h>
0010 #include <TGDMLMatrix.h>
0011 #include <TString.h>
0012 #include <TVector3.h>
0013 #include <fmt/core.h>
0014 #include <fmt/format.h>
0015 #include <cmath>
0016 #include <cstdint>
0017 #include <functional>
0018 #include <map>
0019 #include <numbers>
0020 #include <utility>
0021 #include <vector>
0022
0023 #include "services/geometry/richgeo/RichGeo.h"
0024
0025 // constructor: creates IRT-DD4hep bindings using main `Detector` handle `*det_`
0026 richgeo::IrtGeo::IrtGeo(std::string detName_, gsl::not_null<const dd4hep::Detector*> det_,
0027                         gsl::not_null<const dd4hep::rec::CellIDPositionConverter*> conv_,
0028                         std::shared_ptr<spdlog::logger> log_)
0029     : m_detName(std::move(detName_)), m_det(det_), m_converter(conv_), m_log(std::move(log_)) {
0030   Bind();
0031 }
0032
0033 // Bind() -----------------------------------------
0034 // set IRT and DD4hep geometry handles
0035 void richgeo::IrtGeo::Bind() {
0036   // DD4hep geometry handles
0037   m_detRich = m_det->detector(m_detName);
0038   m_posRich = m_detRich.placement().position();
0039
0040   // IRT geometry handles
0041   m_irtDetectorCollection = new CherenkovDetectorCollection();
0042   m_irtDetector           = m_irtDetectorCollection->AddNewDetector(m_detName.c_str());
0043 }
0044 // ------------------------------------------------
0045
0046 // define the `cell ID -> pixel position` converter, correcting to sensor surface
0047 void richgeo::IrtGeo::SetReadoutIDToPositionLambda() {
0048
0049   m_irtDetector->m_ReadoutIDToPosition =
0050       [&m_log = this->m_log, // capture logger by reference
0051        // capture instance members by value, so those owned by `this` are not mutable here
0052        cell_mask = this->m_irtDetector->GetReadoutCellMask(), converter = this->m_converter,
0053        sensor_info = this->m_sensor_info](auto cell_id) {
0054         // decode cell ID to get the sensor ID and pixel volume centroid
0055         auto sensor_id             = cell_id & cell_mask;
0056         auto pixel_volume_centroid = (1 / dd4hep::mm) * converter->position(cell_id);
0057         // get sensor info
0058         auto sensor_info_it = sensor_info.find(sensor_id);
0059         if (sensor_info_it == sensor_info.end()) {
0060           m_log->warn("cannot find sensor ID {} in IrtGeo; using pixel volume centroid instead",
0061                       sensor_id);
0062           return TVector3(pixel_volume_centroid.x(), pixel_volume_centroid.y(),
0063                           pixel_volume_centroid.z());
0064         }
0065         auto sensor_obj = sensor_info_it->second;
0066         // get pixel surface centroid, given sensor surface offset w.r.t centroid
0067         auto pixel_surface_centroid = pixel_volume_centroid + sensor_obj.surface_offset;
0068         // cross check: make sure pixel and sensor surface centroids are close enough
0069         auto dist = sqrt((pixel_surface_centroid - sensor_obj.surface_centroid).Mag2());
0070         if (dist > sensor_obj.size / std::numbers::sqrt2) {
0071           m_log->warn("dist(pixel,sensor) is too large: {} mm", dist);
0072         }
0073         return TVector3(pixel_surface_centroid.x(), pixel_surface_centroid.y(),
0074                         pixel_surface_centroid.z());
0075       };
0076 }
0077 // ------------------------------------------------
0078
0079 // fill table of refractive indices
0080 void richgeo::IrtGeo::SetRefractiveIndexTable() {
0081   m_log->debug("{:-^60}", " Refractive Index Tables ");
0082   for (auto rad_obj : m_irtDetector->Radiators()) {
0083     m_log->debug("{}:", rad_obj.first.Data());
0084     auto* const rad = rad_obj.second;
0085     const auto* rindex_matrix =
0086         m_det->material(rad->GetAlternativeMaterialName()).property("RINDEX");
0087     for (unsigned row = 0; row < rindex_matrix->GetRows(); row++) {
0088       auto energy = rindex_matrix->Get(row, 0) / dd4hep::eV;
0089       auto rindex = rindex_matrix->Get(row, 1);
0090       m_log->debug("  {:>5} eV   {:<}", energy, rindex);
0091       rad->m_ri_lookup_table.emplace_back(energy, rindex);
0092     }
0093   }
0094 }
0095
0096 // destructor
0097 richgeo::IrtGeo::~IrtGeo() {
0098   delete m_irtDetector;
0099   delete m_irtDetectorCollection;
0100 }