File indexing completed on 2025-07-06 07:55:52
0001
0002
0003
0004 #include "IrtGeoDRICH.h"
0005
0006 #include <DD4hep/DetElement.h>
0007 #include <DD4hep/Fields.h>
0008 #include <DD4hep/Objects.h>
0009 #include <DDRec/DetectorData.h>
0010 #include <Evaluator/DD4hepUnits.h>
0011 #include <IRT/CherenkovDetector.h>
0012 #include <IRT/CherenkovDetectorCollection.h>
0013 #include <IRT/CherenkovRadiator.h>
0014 #include <IRT/G4Object.h>
0015 #include <Math/GenVector/Cartesian3D.h>
0016 #include <Math/GenVector/DisplacementVector3D.h>
0017 #include <TRef.h>
0018 #include <cstdint>
0019 #include <fmt/core.h>
0020 #include <map>
0021 #include <stdexcept>
0022 #include <string>
0023 #include <unordered_map>
0024 #include <utility>
0025
0026 void richgeo::IrtGeoDRICH::DD4hep_to_IRT() {
0027
0028
0029
0030
0031
0032
0033
0034
0035 auto nSectors = m_det->constant<int>("DRICH_num_sectors");
0036 auto vesselZmin = m_det->constant<double>("DRICH_zmin") / dd4hep::mm;
0037 auto vesselWindowThickness = m_det->constant<double>("DRICH_window_thickness") / dd4hep::mm;
0038 auto gasvolMaterial = m_det->constant<std::string>("DRICH_gasvol_material");
0039 TVector3 normX(1, 0, 0);
0040 TVector3 normY(0, -1, 0);
0041 m_surfEntrance =
0042 new FlatSurface(TVector3(0, 0, vesselZmin + vesselWindowThickness), normX, normY);
0043 for (int isec = 0; isec < nSectors; isec++) {
0044 auto* cv = m_irtDetectorCollection->SetContainerVolume(
0045 m_irtDetector,
0046 RadiatorName(kGas).c_str(),
0047 isec,
0048 (G4LogicalVolume*)nullptr,
0049 nullptr,
0050 m_surfEntrance
0051 );
0052 cv->SetAlternativeMaterialName(gasvolMaterial.c_str());
0053 }
0054
0055
0056
0057 auto cellMask = uint64_t(std::stoull(m_det->constant<std::string>("DRICH_cell_mask")));
0058 m_irtPhotonDetector = new CherenkovPhotonDetector(nullptr, nullptr);
0059 m_irtDetector->SetReadoutCellMask(cellMask);
0060 m_irtDetectorCollection->AddPhotonDetector(m_irtDetector,
0061 nullptr,
0062 m_irtPhotonDetector
0063 );
0064 m_log->debug("cellMask = {:#X}", cellMask);
0065
0066
0067
0068
0069
0070
0071
0072 auto aerogelZpos = m_det->constant<double>("DRICH_aerogel_zpos") / dd4hep::mm;
0073 auto aerogelThickness = m_det->constant<double>("DRICH_aerogel_thickness") / dd4hep::mm;
0074 auto aerogelMaterial = m_det->constant<std::string>("DRICH_aerogel_material");
0075 auto filterZpos = m_det->constant<double>("DRICH_filter_zpos") / dd4hep::mm;
0076 auto filterThickness = m_det->constant<double>("DRICH_filter_thickness") / dd4hep::mm;
0077 auto filterMaterial = m_det->constant<std::string>("DRICH_filter_material");
0078 m_aerogelFlatSurface = new FlatSurface(TVector3(0, 0, aerogelZpos), normX, normY);
0079 m_filterFlatSurface = new FlatSurface(TVector3(0, 0, filterZpos), normX, normY);
0080 for (int isec = 0; isec < nSectors; isec++) {
0081 auto* aerogelFlatRadiator = m_irtDetectorCollection->AddFlatRadiator(
0082 m_irtDetector,
0083 RadiatorName(kAerogel).c_str(),
0084 isec,
0085 (G4LogicalVolume*)(0x1),
0086 nullptr,
0087 m_aerogelFlatSurface,
0088 aerogelThickness
0089 );
0090 auto* filterFlatRadiator = m_irtDetectorCollection->AddFlatRadiator(
0091 m_irtDetector,
0092 "Filter",
0093 isec,
0094 (G4LogicalVolume*)(0x2),
0095 nullptr,
0096 m_filterFlatSurface,
0097 filterThickness
0098 );
0099 aerogelFlatRadiator->SetAlternativeMaterialName(aerogelMaterial.c_str());
0100 filterFlatRadiator->SetAlternativeMaterialName(filterMaterial.c_str());
0101 }
0102 m_log->debug("aerogelZpos = {:f} mm", aerogelZpos);
0103 m_log->debug("filterZpos = {:f} mm", filterZpos);
0104 m_log->debug("aerogel thickness = {:f} mm", aerogelThickness);
0105 m_log->debug("filter thickness = {:f} mm", filterThickness);
0106
0107
0108 for (int isec = 0; isec < nSectors; isec++) {
0109 std::string secName = "sec" + std::to_string(isec);
0110
0111 auto mirrorRadius = m_det->constant<double>("DRICH_mirror_radius") / dd4hep::mm;
0112 dd4hep::Position mirrorCenter(
0113 m_det->constant<double>("DRICH_mirror_center_x_" + secName) / dd4hep::mm,
0114 m_det->constant<double>("DRICH_mirror_center_y_" + secName) / dd4hep::mm,
0115 m_det->constant<double>("DRICH_mirror_center_z_" + secName) / dd4hep::mm);
0116 m_mirrorSphericalSurface = new SphericalSurface(
0117 TVector3(mirrorCenter.x(), mirrorCenter.y(), mirrorCenter.z()), mirrorRadius);
0118 m_mirrorOpticalBoundary =
0119 new OpticalBoundary(m_irtDetector->GetContainerVolume(),
0120 m_mirrorSphericalSurface,
0121 false
0122 );
0123 m_irtDetector->AddOpticalBoundary(isec, m_mirrorOpticalBoundary);
0124 m_log->debug("");
0125 m_log->debug(" SECTOR {:d} MIRROR:", isec);
0126 m_log->debug(" mirror x = {:f} mm", mirrorCenter.x());
0127 m_log->debug(" mirror y = {:f} mm", mirrorCenter.y());
0128 m_log->debug(" mirror z = {:f} mm", mirrorCenter.z());
0129 m_log->debug(" mirror R = {:f} mm", mirrorRadius);
0130
0131
0132 auto* rad = m_irtDetector->GetRadiator(RadiatorName(kGas).c_str());
0133 if (rad != nullptr) {
0134 rad->m_Borders[isec].second = m_mirrorSphericalSurface;
0135 } else {
0136 throw std::runtime_error("Gas radiator not built in IrtGeo");
0137 }
0138
0139
0140 m_log->trace(" SENSORS:");
0141 m_log->trace(
0142 "--------------------------------------------------------------------------------------");
0143 m_log->trace(
0144 "name ID sector pos_x pos_y pos_z normX_x normX_y normX_z normY_x normY_y normY_z");
0145 m_log->trace(
0146 "--------------------------------------------------------------------------------------");
0147 auto sensorThickness = m_det->constant<double>("DRICH_sensor_thickness") / dd4hep::mm;
0148 auto sensorSize = m_det->constant<double>("DRICH_sensor_size") / dd4hep::mm;
0149 for (auto const& [de_name, detSensor] : m_detRich.children()) {
0150 if (de_name.find("sensor_de_" + secName) != std::string::npos) {
0151
0152
0153 const auto sensorID = detSensor.id();
0154 auto* const detSensorPars = detSensor.extension<dd4hep::rec::VariantParameters>(true);
0155 if (detSensorPars == nullptr) {
0156 throw std::runtime_error(
0157 fmt::format("sensor '{}' does not have VariantParameters", de_name));
0158 }
0159
0160 auto posSensor =
0161 GetVectorFromVariantParameters<dd4hep::Position>(detSensorPars, "pos") / dd4hep::mm;
0162
0163 auto normXdir = GetVectorFromVariantParameters<dd4hep::Direction>(detSensorPars, "normX");
0164 auto normYdir = GetVectorFromVariantParameters<dd4hep::Direction>(detSensorPars, "normY");
0165 auto normZdir = normXdir.Cross(normYdir);
0166
0167 auto surfaceOffset = normZdir.Unit() * (0.5 * sensorThickness);
0168
0169
0170 richgeo::Sensor sensor_info;
0171 sensor_info.size = sensorSize;
0172 sensor_info.surface_centroid = posSensor;
0173 sensor_info.surface_offset = surfaceOffset;
0174 m_sensor_info.insert({sensorID, sensor_info});
0175
0176 m_sensorFlatSurface = new FlatSurface(TVector3(posSensor.x(), posSensor.y(), posSensor.z()),
0177 TVector3(normXdir.x(), normXdir.y(), normXdir.z()),
0178 TVector3(normYdir.x(), normYdir.y(), normYdir.z()));
0179 m_irtDetector->CreatePhotonDetectorInstance(isec,
0180 m_irtPhotonDetector,
0181 sensorID,
0182 m_sensorFlatSurface
0183 );
0184 m_log->trace("{} {:#X} {} {:5.2f} {:5.2f} {:5.2f} {:5.2f} {:5.2f} {:5.2f} {:5.2f} "
0185 "{:5.2f} {:5.2f}",
0186 de_name, sensorID, isec, posSensor.x(), posSensor.y(), posSensor.z(),
0187 normXdir.x(), normXdir.y(), normXdir.z(), normYdir.x(), normYdir.y(),
0188 normYdir.z());
0189 }
0190 }
0191
0192 }
0193
0194
0195 std::map<const std::string, double> rIndices;
0196 rIndices.insert({RadiatorName(kGas), 1.00076});
0197 rIndices.insert({RadiatorName(kAerogel), 1.0190});
0198 rIndices.insert({"Filter", 1.5017});
0199 for (auto const& [rName, rIndex] : rIndices) {
0200 auto* rad = m_irtDetector->GetRadiator(rName.c_str());
0201 if (rad != nullptr) {
0202 rad->SetReferenceRefractiveIndex(rIndex);
0203 }
0204 }
0205
0206
0207 SetRefractiveIndexTable();
0208
0209
0210 SetReadoutIDToPositionLambda();
0211 }
0212 TVector3 richgeo::IrtGeoDRICH::GetSensorSurfaceNorm(CellIDType id) {
0213 TVector3 sensorNorm;
0214 auto cellMask = uint64_t(std::stoull(m_det->constant<std::string>("DRICH_cell_mask")));
0215 auto sensor_info = this->m_sensor_info;
0216 auto sID = id & cellMask;
0217 auto sensor_info_it = sensor_info.find(sID);
0218 if (sensor_info_it != sensor_info.end()) {
0219 auto sensor_obj = sensor_info_it->second;
0220 auto normZdir = sensor_obj.surface_offset.Unit();
0221 sensorNorm.SetX(static_cast<double>(normZdir.x()));
0222 sensorNorm.SetY(static_cast<double>(normZdir.y()));
0223 sensorNorm.SetZ(static_cast<double>(normZdir.z()));
0224 } else {
0225 m_log->error("Cannot find sensor {} in IrtGeoDRICH::GetSensorSurface", id);
0226 throw std::runtime_error("sensor not found in IrtGeoDRIC::GetSensorSurfaceNormal");
0227 }
0228 return sensorNorm;
0229 }
0230
0231 richgeo::IrtGeoDRICH::~IrtGeoDRICH() {
0232 delete m_surfEntrance;
0233 delete m_irtPhotonDetector;
0234 delete m_aerogelFlatSurface;
0235 delete m_filterFlatSurface;
0236 delete m_mirrorSphericalSurface;
0237 delete m_mirrorOpticalBoundary;
0238 delete m_sensorFlatSurface;
0239 }