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