geometry/richgeo/IrtGeoDRICH.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 "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   // begin envelope
0029   /* FIXME: have no connection to GEANT G4LogicalVolume pointers; however all is needed
0030    * is to make them unique so that std::map work internally; resort to using integers,
0031    * who cares; material pointer can seemingly be '0', and effective refractive index
0032    * for all radiators will be assigned at the end by hand; FIXME: should assign it on
0033    * per-photon basis, at birth, like standalone GEANT code does;
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); // normal vectors
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,              // Cherenkov detector
0046         RadiatorName(kGas).c_str(), // name
0047         isec,                       // path
0048         (G4LogicalVolume*)nullptr,  // G4LogicalVolume (inaccessible? use an integer instead)
0049         nullptr,                    // G4RadiatorMaterial (inaccessible?)
0050         m_surfEntrance              // surface
0051     );
0052     cv->SetAlternativeMaterialName(gasvolMaterial.c_str());
0053   }
0054
0055   // photon detector
0056   // - FIXME: args (G4Solid,G4Material) inaccessible?
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,      // Cherenkov detector
0061                                              nullptr,            // G4LogicalVolume (inaccessible?)
0062                                              m_irtPhotonDetector // photon detector
0063   );
0064   m_log->debug("cellMask = {:#X}", cellMask);
0065
0066   // aerogel + filter
0067   /* AddFlatRadiator will create a pair of flat refractive surfaces internally;
0068    * FIXME: should make a small gas gap at the upstream end of the gas volume;
0069    * FIXME: do we need a sector loop?
0070    * FIXME: airgap radiator?
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,                  // Cherenkov detector
0083         RadiatorName(kAerogel).c_str(), // name
0084         isec,                           // path
0085         (G4LogicalVolume*)(0x1),        // G4LogicalVolume (inaccessible? use an integer instead)
0086         nullptr,                        // G4RadiatorMaterial
0087         m_aerogelFlatSurface,           // surface
0088         aerogelThickness                // surface thickness
0089     );
0090     auto* filterFlatRadiator = m_irtDetectorCollection->AddFlatRadiator(
0091         m_irtDetector,           // Cherenkov detector
0092         "Filter",                // name
0093         isec,                    // path
0094         (G4LogicalVolume*)(0x2), // G4LogicalVolume (inaccessible? use an integer instead)
0095         nullptr,                 // G4RadiatorMaterial
0096         m_filterFlatSurface,     // surface
0097         filterThickness          // surface thickness
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   // sector loop
0108   for (int isec = 0; isec < nSectors; isec++) {
0109     std::string secName = "sec" + std::to_string(isec);
0110     // mirrors
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(), // CherenkovRadiator radiator
0120                             m_mirrorSphericalSurface,            // surface
0121                             false                                // bool refractive
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     // complete the radiator volume description; this is the rear side of the container gas volume
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     // sensor modules: search the detector tree for sensors for this sector
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         // get sensor info
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         // - sensor surface position
0160         auto posSensor =
0161             GetVectorFromVariantParameters<dd4hep::Position>(detSensorPars, "pos") / dd4hep::mm;
0162         // - sensor orientation
0163         auto normXdir = GetVectorFromVariantParameters<dd4hep::Direction>(detSensorPars, "normX");
0164         auto normYdir = GetVectorFromVariantParameters<dd4hep::Direction>(detSensorPars, "normY");
0165         auto normZdir = normXdir.Cross(normYdir); // sensor surface normal
0166         // - surface offset, used to convert sensor volume centroid to sensor surface centroid
0167         auto surfaceOffset = normZdir.Unit() * (0.5 * sensorThickness);
0168
0169         // add sensor info to `m_sensor_info` map
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         // create the optical surface
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,                // sector
0180                                                     m_irtPhotonDetector, // CherenkovPhotonDetector
0181                                                     sensorID,            // copy number
0182                                                     m_sensorFlatSurface  // surface
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     } // search for sensors
0191
0192   } // sector loop
0193
0194   // set reference refractive indices // NOTE: numbers may be overridden externally
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   // set refractive index table
0207   SetRefractiveIndexTable();
0208
0209   // define the `cell ID -> pixel position` converter
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 // destructor
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 }