File indexing completed on 2024-09-28 07:03:04
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 <fmt/core.h>
0019 #include <stdint.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 = new FlatSurface(TVector3(0, 0, vesselZmin + vesselWindowThickness), normX, normY);
0042 for (int isec=0; isec<nSectors; isec++) {
0043 auto *cv = m_irtDetectorCollection->SetContainerVolume(
0044 m_irtDetector,
0045 RadiatorName(kGas).c_str(),
0046 isec,
0047 (G4LogicalVolume*)(0x0),
0048 nullptr,
0049 m_surfEntrance
0050 );
0051 cv->SetAlternativeMaterialName(gasvolMaterial.c_str());
0052 }
0053
0054
0055
0056 auto cellMask = uint64_t(std::stoull(m_det->constant<std::string>("DRICH_cell_mask")));
0057 m_irtPhotonDetector = new CherenkovPhotonDetector(nullptr, nullptr);
0058 m_irtDetector->SetReadoutCellMask(cellMask);
0059 m_irtDetectorCollection->AddPhotonDetector(
0060 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 );
0117 m_mirrorSphericalSurface = new SphericalSurface(TVector3(mirrorCenter.x(), mirrorCenter.y(), mirrorCenter.z()), mirrorRadius);
0118 m_mirrorOpticalBoundary = new OpticalBoundary(
0119 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) rad->m_Borders[isec].second = m_mirrorSphericalSurface;
0134 else throw std::runtime_error("Gas radiator not built in IrtGeo");
0135
0136
0137 m_log->trace(" SENSORS:");
0138 m_log->trace("--------------------------------------------------------------------------------------");
0139 m_log->trace("name ID sector pos_x pos_y pos_z normX_x normX_y normX_z normY_x normY_y normY_z");
0140 m_log->trace("--------------------------------------------------------------------------------------");
0141 auto sensorThickness = m_det->constant<double>("DRICH_sensor_thickness") / dd4hep::mm;
0142 auto sensorSize = m_det->constant<double>("DRICH_sensor_size") / dd4hep::mm;
0143 for(auto const& [de_name, detSensor] : m_detRich.children()) {
0144 if(de_name.find("sensor_de_"+secName)!=std::string::npos) {
0145
0146
0147 const auto sensorID = detSensor.id();
0148 const auto detSensorPars = detSensor.extension<dd4hep::rec::VariantParameters>(true);
0149 if(detSensorPars==nullptr)
0150 throw std::runtime_error(fmt::format("sensor '{}' does not have VariantParameters", de_name));
0151
0152 auto posSensor = GetVectorFromVariantParameters<dd4hep::Position>(detSensorPars, "pos") / dd4hep::mm;
0153
0154 auto normXdir = GetVectorFromVariantParameters<dd4hep::Direction>(detSensorPars, "normX");
0155 auto normYdir = GetVectorFromVariantParameters<dd4hep::Direction>(detSensorPars, "normY");
0156 auto normZdir = normXdir.Cross(normYdir);
0157
0158 auto surfaceOffset = normZdir.Unit() * (0.5*sensorThickness);
0159
0160
0161 richgeo::Sensor sensor_info;
0162 sensor_info.size = sensorSize;
0163 sensor_info.surface_centroid = posSensor;
0164 sensor_info.surface_offset = surfaceOffset;
0165 m_sensor_info.insert({ sensorID, sensor_info });
0166
0167 m_sensorFlatSurface = new FlatSurface(
0168 TVector3(posSensor.x(), posSensor.y(), posSensor.z()),
0169 TVector3(normXdir.x(), normXdir.y(), normXdir.z()),
0170 TVector3(normYdir.x(), normYdir.y(), normYdir.z())
0171 );
0172 m_irtDetector->CreatePhotonDetectorInstance(
0173 isec,
0174 m_irtPhotonDetector,
0175 sensorID,
0176 m_sensorFlatSurface
0177 );
0178 m_log->trace(
0179 "{} {:#X} {} {:5.2f} {:5.2f} {:5.2f} {:5.2f} {:5.2f} {:5.2f} {:5.2f} {:5.2f} {:5.2f}",
0180 de_name, sensorID, isec,
0181 posSensor.x(), posSensor.y(), posSensor.z(),
0182 normXdir.x(), normXdir.y(), normXdir.z(),
0183 normYdir.x(), normYdir.y(), normYdir.z()
0184 );
0185 }
0186 }
0187
0188 }
0189
0190
0191 std::map<const std::string, double> rIndices;
0192 rIndices.insert({RadiatorName(kGas), 1.00076});
0193 rIndices.insert({RadiatorName(kAerogel), 1.0190});
0194 rIndices.insert({"Filter", 1.5017});
0195 for (auto const& [rName, rIndex] : rIndices) {
0196 auto rad = m_irtDetector->GetRadiator(rName.c_str());
0197 if (rad)
0198 rad->SetReferenceRefractiveIndex(rIndex);
0199 }
0200
0201
0202 SetRefractiveIndexTable();
0203
0204
0205 SetReadoutIDToPositionLambda();
0206 }
0207 TVector3 richgeo::IrtGeoDRICH::GetSensorSurfaceNorm(CellIDType id){
0208 TVector3 sensorNorm;
0209 auto cellMask = uint64_t(std::stoull(m_det->constant<std::string>("DRICH_cell_mask")));
0210 auto sensor_info = this->m_sensor_info;
0211 auto sID = id & cellMask;
0212 auto sensor_info_it = sensor_info.find(sID);
0213 if(sensor_info_it!=sensor_info.end()){
0214 auto sensor_obj = sensor_info_it->second;
0215 auto normZdir = sensor_obj.surface_offset.Unit();
0216 sensorNorm.SetX(static_cast<double>(normZdir.x()));
0217 sensorNorm.SetY(static_cast<double>(normZdir.y()));
0218 sensorNorm.SetZ(static_cast<double>(normZdir.z()));
0219 }
0220 else{
0221 m_log->error("Cannot find sensor {} in IrtGeoDRICH::GetSensorSurface", id);
0222 throw std::runtime_error("sensor not found in IrtGeoDRIC::GetSensorSurfaceNormal");
0223 }
0224 return sensorNorm;
0225 }
0226
0227 richgeo::IrtGeoDRICH::~IrtGeoDRICH() {
0228 delete m_surfEntrance;
0229 delete m_irtPhotonDetector;
0230 delete m_aerogelFlatSurface;
0231 delete m_filterFlatSurface;
0232 delete m_mirrorSphericalSurface;
0233 delete m_mirrorOpticalBoundary;
0234 delete m_sensorFlatSurface;
0235 }