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