File indexing completed on 2024-09-28 07:03:04
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 <TGeoNode.h>
0017 #include <TRef.h>
0018 #include <TVector3.h>
0019 #include <fmt/core.h>
0020 #include <math.h>
0021 #include <stdint.h>
0022 #include <map>
0023 #include <string>
0024 #include <unordered_map>
0025 #include <utility>
0026
0027 #include "services/geometry/richgeo/RichGeo.h"
0028
0029 void richgeo::IrtGeoPFRICH::DD4hep_to_IRT() {
0030
0031
0032
0033
0034
0035
0036
0037
0038 auto vesselZmin = m_det->constant<double>("PFRICH_zmin") / dd4hep::mm;
0039 auto gasvolMaterial = m_det->constant<std::string>("PFRICH_gasvol_material");
0040 TVector3 normX(1, 0, 0);
0041 TVector3 normY(0, 1, 0);
0042 m_surfEntrance = new FlatSurface(TVector3(0, 0, vesselZmin), normX, normY);
0043 auto cv = m_irtDetectorCollection->SetContainerVolume(
0044 m_irtDetector,
0045 RadiatorName(kGas).c_str(),
0046 0,
0047 (G4LogicalVolume*)(0x0),
0048 nullptr,
0049 m_surfEntrance
0050 );
0051 cv->SetAlternativeMaterialName(gasvolMaterial.c_str());
0052
0053
0054
0055 auto cellMask = uint64_t(std::stoull(m_det->constant<std::string>("PFRICH_cell_mask")));
0056 m_irtPhotonDetector = new CherenkovPhotonDetector(nullptr, nullptr);
0057 m_irtDetector->SetReadoutCellMask(cellMask);
0058 m_irtDetectorCollection->AddPhotonDetector(
0059 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(0,0,1);
0116 auto surfaceOffset = sensorNorm.Unit() * (0.5*sensorThickness);
0117 auto posSensorSurface = posSensor + surfaceOffset;
0118
0119 richgeo::Sensor sensor_info;
0120 sensor_info.size = sensorSize;
0121 sensor_info.surface_centroid = posSensorSurface;
0122 sensor_info.surface_offset = surfaceOffset;
0123 m_sensor_info.insert({ imod, sensor_info });
0124
0125 double sensorLocalNormX[3] = {1.0, 0.0, 0.0};
0126 double sensorLocalNormY[3] = {0.0, 1.0, 0.0};
0127 double sensorGlobalNormX[3], sensorGlobalNormY[3];
0128 pvSensor.ptr()->LocalToMasterVect(sensorLocalNormX, sensorGlobalNormX);
0129 pvSensor.ptr()->LocalToMasterVect(sensorLocalNormY, sensorGlobalNormY);
0130
0131
0132
0133 dd4hep::Direction normXdir, normYdir;
0134 normXdir.SetCoordinates(sensorGlobalNormX);
0135 normYdir.SetCoordinates(sensorGlobalNormY);
0136 auto normZdir = normXdir.Cross(normYdir);
0137 auto testOrtho = normXdir.Dot(normYdir);
0138 auto testRadial = sensorNorm.Cross(normZdir).Mag2();
0139 if(abs(testOrtho)>1e-6 || abs(testRadial)>1e-6) {
0140 m_log->error(
0141 "sensor normal is wrong: normX.normY = {:f} |sensorNorm x normZdir|^2 = {:f}",
0142 testOrtho,
0143 testRadial
0144 );
0145 return;
0146 }
0147
0148
0149 m_sensorFlatSurface = new FlatSurface(
0150 TVector3(posSensorSurface.x(), posSensorSurface.y(), posSensorSurface.z()),
0151 TVector3(sensorGlobalNormX),
0152 TVector3(sensorGlobalNormY)
0153 );
0154 m_irtDetector->CreatePhotonDetectorInstance(
0155 0,
0156 m_irtPhotonDetector,
0157 imod,
0158 m_sensorFlatSurface
0159 );
0160 m_log->trace(
0161 "sensor: id={:#08X} pos=({:5.2f}, {:5.2f}, {:5.2f}) normX=({:5.2f}, {:5.2f}, {:5.2f}) normY=({:5.2f}, {:5.2f}, {:5.2f})",
0162 imod,
0163 posSensorSurface.x(), posSensorSurface.y(), posSensorSurface.z(),
0164 normXdir.x(), normXdir.y(), normXdir.z(),
0165 normYdir.x(), normYdir.y(), normYdir.z()
0166 );
0167
0168
0169
0170
0171 if(firstSensor) {
0172 m_irtDetector->GetRadiator(RadiatorName(kGas).c_str())->m_Borders[0].second = dynamic_cast<ParametricSurface*>(m_sensorFlatSurface);
0173 firstSensor = false;
0174 }
0175
0176 }
0177 }
0178
0179
0180 std::map<const char*, double> rIndices;
0181 rIndices.insert({RadiatorName(kGas).c_str(), 1.0013});
0182 rIndices.insert({RadiatorName(kAerogel).c_str(), 1.0190});
0183 rIndices.insert({"Filter", 1.5017});
0184 for (auto const& [rName, rIndex] : rIndices) {
0185 auto rad = m_irtDetector->GetRadiator(rName);
0186 if (rad)
0187 rad->SetReferenceRefractiveIndex(rIndex);
0188 }
0189
0190
0191 SetRefractiveIndexTable();
0192
0193
0194 SetReadoutIDToPositionLambda();
0195 }
0196
0197
0198 richgeo::IrtGeoPFRICH::~IrtGeoPFRICH() {
0199 delete m_surfEntrance;
0200 delete m_irtPhotonDetector;
0201 delete m_aerogelFlatSurface;
0202 delete m_filterFlatSurface;
0203 delete m_sensorFlatSurface;
0204 }