File indexing completed on 2025-07-11 07:53:43
0001
0002
0003
0004 #include "IrtGeoPFRICH.h"
0005
0006 #include <DD4hep/DetElement.h>
0007 #include <DD4hep/Fields.h>
0008 #include <DD4hep/Volumes.h>
0009 #include <Evaluator/DD4hepUnits.h>
0010 #include <IRT/CherenkovDetector.h>
0011 #include <IRT/CherenkovDetectorCollection.h>
0012 #include <IRT/CherenkovRadiator.h>
0013 #include <IRT/G4Object.h>
0014 #include <Math/GenVector/Cartesian3D.h>
0015 #include <Math/GenVector/DisplacementVector3D.h>
0016 #include <RtypesCore.h>
0017 #include <TGeoNode.h>
0018 #include <TRef.h>
0019 #include <TVector3.h>
0020 #include <fmt/core.h>
0021 #include <cmath>
0022 #include <cstdint>
0023 #include <map>
0024 #include <string>
0025 #include <unordered_map>
0026 #include <utility>
0027
0028 #include "services/geometry/richgeo/RichGeo.h"
0029
0030 void richgeo::IrtGeoPFRICH::DD4hep_to_IRT() {
0031
0032
0033
0034
0035
0036
0037
0038
0039 auto vesselZmin = m_det->constant<double>("PFRICH_zmin") / dd4hep::mm;
0040 auto gasvolMaterial = m_det->constant<std::string>("PFRICH_gasvol_material");
0041 TVector3 normX(1, 0, 0);
0042 TVector3 normY(0, 1, 0);
0043 m_surfEntrance = new FlatSurface(TVector3(0, 0, vesselZmin), normX, normY);
0044 auto* cv = m_irtDetectorCollection->SetContainerVolume(
0045 m_irtDetector,
0046 RadiatorName(kGas).c_str(),
0047 0,
0048 (G4LogicalVolume*)nullptr,
0049 nullptr,
0050 m_surfEntrance
0051 );
0052 cv->SetAlternativeMaterialName(gasvolMaterial.c_str());
0053
0054
0055
0056 auto cellMask = uint64_t(std::stoull(m_det->constant<std::string>("PFRICH_cell_mask")));
0057 m_irtPhotonDetector = new CherenkovPhotonDetector(nullptr, nullptr);
0058 m_irtDetector->SetReadoutCellMask(cellMask);
0059 m_irtDetectorCollection->AddPhotonDetector(m_irtDetector,
0060 nullptr,
0061 m_irtPhotonDetector
0062 );
0063 m_log->debug("cellMask = {:#X}", cellMask);
0064
0065
0066
0067
0068
0069 auto aerogelZpos = m_det->constant<double>("PFRICH_aerogel_zpos") / dd4hep::mm;
0070 auto aerogelThickness = m_det->constant<double>("PFRICH_aerogel_thickness") / dd4hep::mm;
0071 auto aerogelMaterial = m_det->constant<std::string>("PFRICH_aerogel_material");
0072 auto filterZpos = m_det->constant<double>("PFRICH_filter_zpos") / dd4hep::mm;
0073 auto filterThickness = m_det->constant<double>("PFRICH_filter_thickness") / dd4hep::mm;
0074 auto filterMaterial = m_det->constant<std::string>("PFRICH_filter_material");
0075 m_aerogelFlatSurface = new FlatSurface(TVector3(0, 0, aerogelZpos), normX, normY);
0076 m_filterFlatSurface = new FlatSurface(TVector3(0, 0, filterZpos), normX, normY);
0077 auto* aerogelFlatRadiator = m_irtDetectorCollection->AddFlatRadiator(
0078 m_irtDetector,
0079 RadiatorName(kAerogel).c_str(),
0080 0,
0081 (G4LogicalVolume*)(0x1),
0082 nullptr,
0083 m_aerogelFlatSurface,
0084 aerogelThickness
0085 );
0086 auto* filterFlatRadiator = m_irtDetectorCollection->AddFlatRadiator(
0087 m_irtDetector,
0088 "Filter",
0089 0,
0090 (G4LogicalVolume*)(0x2),
0091 nullptr,
0092 m_filterFlatSurface,
0093 filterThickness
0094 );
0095 aerogelFlatRadiator->SetAlternativeMaterialName(aerogelMaterial.c_str());
0096 filterFlatRadiator->SetAlternativeMaterialName(filterMaterial.c_str());
0097 m_log->debug("aerogelZpos = {:f} mm", aerogelZpos);
0098 m_log->debug("filterZpos = {:f} mm", filterZpos);
0099 m_log->debug("aerogel thickness = {:f} mm", aerogelThickness);
0100 m_log->debug("filter thickness = {:f} mm", filterThickness);
0101
0102
0103 auto sensorThickness = m_det->constant<double>("PFRICH_sensor_thickness") / dd4hep::mm;
0104 auto sensorSize = m_det->constant<double>("PFRICH_sensor_size") / dd4hep::mm;
0105 bool firstSensor = true;
0106 for (auto const& [de_name, detSensor] : m_detRich.children()) {
0107 if (de_name.find("sensor_de") != std::string::npos) {
0108
0109
0110 auto imod = detSensor.id();
0111
0112 auto pvSensor = detSensor.placement();
0113 auto posSensor = (1 / dd4hep::mm) * (m_posRich + pvSensor.position());
0114
0115 dd4hep::Direction sensorNorm(
0116 0, 0,
0117 1);
0118 auto surfaceOffset = sensorNorm.Unit() * (0.5 * sensorThickness);
0119 auto posSensorSurface = posSensor + surfaceOffset;
0120
0121 richgeo::Sensor sensor_info;
0122 sensor_info.size = sensorSize;
0123 sensor_info.surface_centroid = posSensorSurface;
0124 sensor_info.surface_offset = surfaceOffset;
0125 m_sensor_info.insert({imod, sensor_info});
0126
0127 double sensorLocalNormX[3] = {1.0, 0.0, 0.0};
0128 double sensorLocalNormY[3] = {0.0, 1.0, 0.0};
0129 double sensorGlobalNormX[3];
0130 double sensorGlobalNormY[3];
0131 pvSensor.ptr()->LocalToMasterVect(
0132 static_cast<const Double_t*>(sensorLocalNormX),
0133 static_cast<Double_t*>(
0134 sensorGlobalNormX));
0135 pvSensor.ptr()->LocalToMasterVect(static_cast<const Double_t*>(sensorLocalNormY),
0136 static_cast<Double_t*>(sensorGlobalNormY));
0137
0138
0139
0140 dd4hep::Direction normXdir;
0141 dd4hep::Direction normYdir;
0142 normXdir.SetCoordinates(static_cast<const Double_t*>(sensorGlobalNormX));
0143 normYdir.SetCoordinates(static_cast<const Double_t*>(sensorGlobalNormY));
0144 auto normZdir =
0145 normXdir.Cross(normYdir);
0146 auto testOrtho = normXdir.Dot(normYdir);
0147 auto testRadial = sensorNorm.Cross(normZdir)
0148 .Mag2();
0149 if (std::abs(testOrtho) > 1e-6 || std::abs(testRadial) > 1e-6) {
0150 m_log->error(
0151 "sensor normal is wrong: normX.normY = {:f} |sensorNorm x normZdir|^2 = {:f}",
0152 testOrtho, testRadial);
0153 return;
0154 }
0155
0156
0157 m_sensorFlatSurface = new FlatSurface(
0158 TVector3(posSensorSurface.x(), posSensorSurface.y(), posSensorSurface.z()),
0159 TVector3(sensorGlobalNormX), TVector3(sensorGlobalNormY));
0160 m_irtDetector->CreatePhotonDetectorInstance(0,
0161 m_irtPhotonDetector,
0162 imod,
0163 m_sensorFlatSurface
0164 );
0165 m_log->trace("sensor: id={:#08X} pos=({:5.2f}, {:5.2f}, {:5.2f}) normX=({:5.2f}, {:5.2f}, "
0166 "{:5.2f}) normY=({:5.2f}, {:5.2f}, {:5.2f})",
0167 imod, posSensorSurface.x(), posSensorSurface.y(), posSensorSurface.z(),
0168 normXdir.x(), normXdir.y(), normXdir.z(), normYdir.x(), normYdir.y(),
0169 normYdir.z());
0170
0171
0172
0173
0174 if (firstSensor) {
0175 m_irtDetector->GetRadiator(RadiatorName(kGas).c_str())->m_Borders[0].second =
0176 dynamic_cast<ParametricSurface*>(m_sensorFlatSurface);
0177 firstSensor = false;
0178 }
0179
0180 }
0181 }
0182
0183
0184 std::map<const char*, double> rIndices;
0185 rIndices.insert({RadiatorName(kGas).c_str(), 1.0013});
0186 rIndices.insert({RadiatorName(kAerogel).c_str(), 1.0190});
0187 rIndices.insert({"Filter", 1.5017});
0188 for (auto const& [rName, rIndex] : rIndices) {
0189 auto* rad = m_irtDetector->GetRadiator(rName);
0190 if (rad != nullptr) {
0191 rad->SetReferenceRefractiveIndex(rIndex);
0192 }
0193 }
0194
0195
0196 SetRefractiveIndexTable();
0197
0198
0199 SetReadoutIDToPositionLambda();
0200 }
0201
0202
0203 richgeo::IrtGeoPFRICH::~IrtGeoPFRICH() {
0204 delete m_surfEntrance;
0205 delete m_irtPhotonDetector;
0206 delete m_aerogelFlatSurface;
0207 delete m_filterFlatSurface;
0208 delete m_sensorFlatSurface;
0209 }