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 <Parsers/Primitives.h>
0018 #include <TRef.h>
0019 #include <fmt/core.h>
0020 #include <fmt/format.h>
0021 #include <cstdint>
0022 #include <map>
0023 #include <stdexcept>
0024 #include <string>
0025 #include <unordered_map>
0026 #include <utility>
0027
0028 void richgeo::IrtGeoDRICH::DD4hep_to_IRT() {
0029
0030   // begin envelope
0031   /* FIXME: have no connection to GEANT G4LogicalVolume pointers; however all is needed
0032    * is to make them unique so that std::map work internally; resort to using integers,
0033    * who cares; material pointer can seemingly be '0', and effective refractive index
0034    * for all radiators will be assigned at the end by hand; FIXME: should assign it on
0035    * per-photon basis, at birth, like standalone GEANT code does;
0036    */
0037   auto nSectors              = m_det->constant<int>("DRICH_num_sectors");
0038   auto vesselZmin            = m_det->constant<double>("DRICH_zmin") / dd4hep::mm;
0039   auto vesselWindowThickness = m_det->constant<double>("DRICH_window_thickness") / dd4hep::mm;
0040   auto gasvolMaterial        = m_det->constant<std::string>("DRICH_gasvol_material");
0041   TVector3 normX(1, 0, 0); // normal vectors
0042   TVector3 normY(0, -1, 0);
0043   m_surfEntrance =
0044       new FlatSurface(TVector3(0, 0, vesselZmin + vesselWindowThickness), normX, normY);
0045   for (int isec = 0; isec < nSectors; isec++) {
0046     auto* cv = m_irtDetectorCollection->SetContainerVolume(
0047         m_irtDetector,              // Cherenkov detector
0048         RadiatorName(kGas).c_str(), // name
0049         isec,                       // path
0050         (G4LogicalVolume*)nullptr,  // G4LogicalVolume (inaccessible? use an integer instead)
0051         nullptr,                    // G4RadiatorMaterial (inaccessible?)
0052         m_surfEntrance              // surface
0053     );
0054     cv->SetAlternativeMaterialName(gasvolMaterial.c_str());
0055   }
0056
0057   // photon detector
0058   // - FIXME: args (G4Solid,G4Material) inaccessible?
0059   auto cellMask       = uint64_t(std::stoull(m_det->constant<std::string>("DRICH_cell_mask")));
0060   m_irtPhotonDetector = new CherenkovPhotonDetector(nullptr, nullptr);
0061   m_irtDetector->SetReadoutCellMask(cellMask);
0062   m_irtDetectorCollection->AddPhotonDetector(m_irtDetector,      // Cherenkov detector
0063                                              nullptr,            // G4LogicalVolume (inaccessible?)
0064                                              m_irtPhotonDetector // photon detector
0065   );
0066   m_log->debug("cellMask = {:#X}", cellMask);
0067
0068   // aerogel + filter
0069   /* AddFlatRadiator will create a pair of flat refractive surfaces internally;
0070    * FIXME: should make a small gas gap at the upstream end of the gas volume;
0071    * FIXME: do we need a sector loop?
0072    * FIXME: airgap radiator?
0073    */
0074   auto aerogelZpos      = m_det->constant<double>("DRICH_aerogel_zpos") / dd4hep::mm;
0075   auto aerogelThickness = m_det->constant<double>("DRICH_aerogel_thickness") / dd4hep::mm;
0076   auto aerogelMaterial  = m_det->constant<std::string>("DRICH_aerogel_material");
0077   auto filterZpos       = m_det->constant<double>("DRICH_filter_zpos") / dd4hep::mm;
0078   auto filterThickness  = m_det->constant<double>("DRICH_filter_thickness") / dd4hep::mm;
0079   auto filterMaterial   = m_det->constant<std::string>("DRICH_filter_material");
0080   m_aerogelFlatSurface  = new FlatSurface(TVector3(0, 0, aerogelZpos), normX, normY);
0081   m_filterFlatSurface   = new FlatSurface(TVector3(0, 0, filterZpos), normX, normY);
0082   for (int isec = 0; isec < nSectors; isec++) {
0083     auto* aerogelFlatRadiator = m_irtDetectorCollection->AddFlatRadiator(
0084         m_irtDetector,                  // Cherenkov detector
0085         RadiatorName(kAerogel).c_str(), // name
0086         isec,                           // path
0087         (G4LogicalVolume*)(0x1),        // G4LogicalVolume (inaccessible? use an integer instead)
0088         nullptr,                        // G4RadiatorMaterial
0089         m_aerogelFlatSurface,           // surface
0090         aerogelThickness                // surface thickness
0091     );
0092     auto* filterFlatRadiator = m_irtDetectorCollection->AddFlatRadiator(
0093         m_irtDetector,           // Cherenkov detector
0094         "Filter",                // name
0095         isec,                    // path
0096         (G4LogicalVolume*)(0x2), // G4LogicalVolume (inaccessible? use an integer instead)
0097         nullptr,                 // G4RadiatorMaterial
0098         m_filterFlatSurface,     // surface
0099         filterThickness          // surface thickness
0100     );
0101     aerogelFlatRadiator->SetAlternativeMaterialName(aerogelMaterial.c_str());
0102     filterFlatRadiator->SetAlternativeMaterialName(filterMaterial.c_str());
0103   }
0104   m_log->debug("aerogelZpos = {:f} mm", aerogelZpos);
0105   m_log->debug("filterZpos  = {:f} mm", filterZpos);
0106   m_log->debug("aerogel thickness = {:f} mm", aerogelThickness);
0107   m_log->debug("filter thickness  = {:f} mm", filterThickness);
0108
0109   // sector loop
0110   for (int isec = 0; isec < nSectors; isec++) {
0111     std::string secName = "sec" + std::to_string(isec);
0112     // mirrors
0113     auto mirrorRadius = m_det->constant<double>("DRICH_mirror_radius") / dd4hep::mm;
0114     dd4hep::Position mirrorCenter(
0115         m_det->constant<double>("DRICH_mirror_center_x_" + secName) / dd4hep::mm,
0116         m_det->constant<double>("DRICH_mirror_center_y_" + secName) / dd4hep::mm,
0117         m_det->constant<double>("DRICH_mirror_center_z_" + secName) / dd4hep::mm);
0118     m_mirrorSphericalSurface = new SphericalSurface(
0119         TVector3(mirrorCenter.x(), mirrorCenter.y(), mirrorCenter.z()), mirrorRadius);
0120     m_mirrorOpticalBoundary =
0121         new OpticalBoundary(m_irtDetector->GetContainerVolume(), // CherenkovRadiator radiator
0122                             m_mirrorSphericalSurface,            // surface
0123                             false                                // bool refractive
0124         );
0125     m_irtDetector->AddOpticalBoundary(isec, m_mirrorOpticalBoundary);
0126     m_log->debug("");
0127     m_log->debug("  SECTOR {:d} MIRROR:", isec);
0128     m_log->debug("    mirror x = {:f} mm", mirrorCenter.x());
0129     m_log->debug("    mirror y = {:f} mm", mirrorCenter.y());
0130     m_log->debug("    mirror z = {:f} mm", mirrorCenter.z());
0131     m_log->debug("    mirror R = {:f} mm", mirrorRadius);
0132
0133     // complete the radiator volume description; this is the rear side of the container gas volume
0134     auto* rad = m_irtDetector->GetRadiator(RadiatorName(kGas).c_str());
0135     if (rad != nullptr) {
0136       rad->m_Borders[isec].second = m_mirrorSphericalSurface;
0137     } else {
0138       throw std::runtime_error("Gas radiator not built in IrtGeo");
0139     }
0140
0141     // sensor modules: search the detector tree for sensors for this sector
0142     m_log->trace("  SENSORS:");
0143     m_log->trace(
0144         "--------------------------------------------------------------------------------------");
0145     m_log->trace(
0146         "name ID sector   pos_x pos_y pos_z   normX_x normX_y normX_z   normY_x normY_y normY_z");
0147     m_log->trace(
0148         "--------------------------------------------------------------------------------------");
0149     auto sensorThickness = m_det->constant<double>("DRICH_sensor_thickness") / dd4hep::mm;
0150     auto sensorSize      = m_det->constant<double>("DRICH_sensor_size") / dd4hep::mm;
0151     for (auto const& [de_name, detSensor] : m_detRich.children()) {
0152       if (de_name.find("sensor_de_" + secName) != std::string::npos) {
0153
0154         // get sensor info
0155         const auto sensorID       = detSensor.id();
0156         auto* const detSensorPars = detSensor.extension<dd4hep::rec::VariantParameters>(true);
0157         if (detSensorPars == nullptr) {
0158           throw std::runtime_error(
0159               fmt::format("sensor '{}' does not have VariantParameters", de_name));
0160         }
0161         // - sensor surface position
0162         auto posSensor =
0163             GetVectorFromVariantParameters<dd4hep::Position>(detSensorPars, "pos") / dd4hep::mm;
0164         // - sensor orientation
0165         auto normXdir = GetVectorFromVariantParameters<dd4hep::Direction>(detSensorPars, "normX");
0166         auto normYdir = GetVectorFromVariantParameters<dd4hep::Direction>(detSensorPars, "normY");
0167         auto normZdir = normXdir.Cross(normYdir); // sensor surface normal
0168         // - surface offset, used to convert sensor volume centroid to sensor surface centroid
0169         auto surfaceOffset = normZdir.Unit() * (0.5 * sensorThickness);
0170
0171         // add sensor info to `m_sensor_info` map
0172         richgeo::Sensor sensor_info;
0173         sensor_info.size             = sensorSize;
0174         sensor_info.surface_centroid = posSensor;
0175         sensor_info.surface_offset   = surfaceOffset;
0176         m_sensor_info.insert({sensorID, sensor_info});
0177         // create the optical surface
0178         m_sensorFlatSurface = new FlatSurface(TVector3(posSensor.x(), posSensor.y(), posSensor.z()),
0179                                               TVector3(normXdir.x(), normXdir.y(), normXdir.z()),
0180                                               TVector3(normYdir.x(), normYdir.y(), normYdir.z()));
0181         m_irtDetector->CreatePhotonDetectorInstance(isec,                // sector
0182                                                     m_irtPhotonDetector, // CherenkovPhotonDetector
0183                                                     sensorID,            // copy number
0184                                                     m_sensorFlatSurface  // surface
0185         );
0186         m_log->trace("{} {:#X} {}   {:5.2f} {:5.2f} {:5.2f}   {:5.2f} {:5.2f} {:5.2f}   {:5.2f} "
0187                      "{:5.2f} {:5.2f}",
0188                      de_name, sensorID, isec, posSensor.x(), posSensor.y(), posSensor.z(),
0189                      normXdir.x(), normXdir.y(), normXdir.z(), normYdir.x(), normYdir.y(),
0190                      normYdir.z());
0191       }
0192     } // search for sensors
0193
0194   } // sector loop
0195
0196   // set reference refractive indices // NOTE: numbers may be overridden externally
0197   std::map<const std::string, double> rIndices;
0198   rIndices.insert({RadiatorName(kGas), 1.00076});
0199   rIndices.insert({RadiatorName(kAerogel), 1.0190});
0200   rIndices.insert({"Filter", 1.5017});
0201   for (auto const& [rName, rIndex] : rIndices) {
0202     auto* rad = m_irtDetector->GetRadiator(rName.c_str());
0203     if (rad != nullptr) {
0204       rad->SetReferenceRefractiveIndex(rIndex);
0205     }
0206   }
0207
0208   // set refractive index table
0209   SetRefractiveIndexTable();
0210
0211   // define the `cell ID -> pixel position` converter
0212   SetReadoutIDToPositionLambda();
0213 }
0214 TVector3 richgeo::IrtGeoDRICH::GetSensorSurfaceNorm(CellIDType id) {
0215   TVector3 sensorNorm;
0216   auto cellMask       = uint64_t(std::stoull(m_det->constant<std::string>("DRICH_cell_mask")));
0217   auto sensor_info    = this->m_sensor_info;
0218   auto sID            = id & cellMask;
0219   auto sensor_info_it = sensor_info.find(sID);
0220   if (sensor_info_it != sensor_info.end()) {
0221     auto sensor_obj = sensor_info_it->second;
0222     auto normZdir   = sensor_obj.surface_offset.Unit();
0223     sensorNorm.SetX(static_cast<double>(normZdir.x()));
0224     sensorNorm.SetY(static_cast<double>(normZdir.y()));
0225     sensorNorm.SetZ(static_cast<double>(normZdir.z()));
0226   } else {
0227     m_log->error("Cannot find sensor {} in IrtGeoDRICH::GetSensorSurface", id);
0228     throw std::runtime_error("sensor not found in IrtGeoDRIC::GetSensorSurfaceNormal");
0229   }
0230   return sensorNorm;
0231 }
0232 // destructor
0233 richgeo::IrtGeoDRICH::~IrtGeoDRICH() {
0234   delete m_surfEntrance;
0235   delete m_irtPhotonDetector;
0236   delete m_aerogelFlatSurface;
0237   delete m_filterFlatSurface;
0238   delete m_mirrorSphericalSurface;
0239   delete m_mirrorOpticalBoundary;
0240   delete m_sensorFlatSurface;
0241 }