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 <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   // 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 = 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,              // Cherenkov detector
0045         RadiatorName(kGas).c_str(), // name
0046         isec,                       // path
0047         (G4LogicalVolume*)(0x0),    // G4LogicalVolume (inaccessible? use an integer instead)
0048         nullptr,                    // G4RadiatorMaterial (inaccessible?)
0049         m_surfEntrance              // surface
0050         );
0051     cv->SetAlternativeMaterialName(gasvolMaterial.c_str());
0052   }
0053
0054   // photon detector
0055   // - FIXME: args (G4Solid,G4Material) inaccessible?
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,      // 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       );
0117     m_mirrorSphericalSurface = new SphericalSurface(TVector3(mirrorCenter.x(), mirrorCenter.y(), mirrorCenter.z()), mirrorRadius);
0118     m_mirrorOpticalBoundary  = new OpticalBoundary(
0119         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) rad->m_Borders[isec].second = m_mirrorSphericalSurface;
0134     else throw std::runtime_error("Gas radiator not built in IrtGeo");
0135
0136     // sensor modules: search the detector tree for sensors for this sector
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         // get sensor info
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         // - sensor surface position
0152         auto posSensor = GetVectorFromVariantParameters<dd4hep::Position>(detSensorPars, "pos") / dd4hep::mm;
0153         // - sensor orientation
0154         auto normXdir = GetVectorFromVariantParameters<dd4hep::Direction>(detSensorPars, "normX");
0155         auto normYdir = GetVectorFromVariantParameters<dd4hep::Direction>(detSensorPars, "normY");
0156         auto normZdir = normXdir.Cross(normYdir); // sensor surface normal
0157         // - surface offset, used to convert sensor volume centroid to sensor surface centroid
0158         auto surfaceOffset = normZdir.Unit() * (0.5*sensorThickness);
0159
0160         // add sensor info to `m_sensor_info` map
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         // create the optical surface
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,                // sector
0174             m_irtPhotonDetector, // CherenkovPhotonDetector
0175             sensorID,            // copy number
0176             m_sensorFlatSurface  // surface
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     } // search for sensors
0187
0188   } // sector loop
0189
0190   // set reference refractive indices // NOTE: numbers may be overridden externally
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   // set refractive index table
0202   SetRefractiveIndexTable();
0203
0204   // define the `cell ID -> pixel position` converter
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 // destructor
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 }