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