geometry/richgeo/IrtGeoPFRICH.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 "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 <RtypesCore.h>
0017 #include <TGeoNode.h>
0018 #include <TRef.h>
0019 #include <TVector3.h>
0020 #include <fmt/core.h>
0021 #include <cmath>
0022 #include <cstdint>
0023 #include <map>
0024 #include <string>
0025 #include <unordered_map>
0026 #include <utility>
0027
0028 #include "services/geometry/richgeo/RichGeo.h"
0029
0030 void richgeo::IrtGeoPFRICH::DD4hep_to_IRT() {
0031
0032   // begin envelope
0033   /* FIXME: have no connection to GEANT G4LogicalVolume pointers; however all is needed
0034    * is to make them unique so that std::map work internally; resort to using integers,
0035    * who cares; material pointer can seemingly be '0', and effective refractive index
0036    * for all radiators will be assigned at the end by hand; FIXME: should assign it on
0037    * per-photon basis, at birth, like standalone GEANT code does;
0038    */
0039   auto vesselZmin     = m_det->constant<double>("PFRICH_zmin") / dd4hep::mm;
0040   auto gasvolMaterial = m_det->constant<std::string>("PFRICH_gasvol_material");
0041   TVector3 normX(1, 0, 0); // normal vectors
0042   TVector3 normY(0, 1, 0);
0043   m_surfEntrance = new FlatSurface(TVector3(0, 0, vesselZmin), normX, normY);
0044   auto* cv       = m_irtDetectorCollection->SetContainerVolume(
0045       m_irtDetector,              // Cherenkov detector
0046       RadiatorName(kGas).c_str(), // name
0047       0,                          // 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   // photon detector
0055   // - FIXME: args (G4Solid,G4Material) inaccessible?
0056   auto cellMask       = uint64_t(std::stoull(m_det->constant<std::string>("PFRICH_cell_mask")));
0057   m_irtPhotonDetector = new CherenkovPhotonDetector(nullptr, nullptr);
0058   m_irtDetector->SetReadoutCellMask(cellMask);
0059   m_irtDetectorCollection->AddPhotonDetector(m_irtDetector,      // Cherenkov detector
0060                                              nullptr,            // G4LogicalVolume (inaccessible?)
0061                                              m_irtPhotonDetector // photon detector
0062   );
0063   m_log->debug("cellMask = {:#X}", cellMask);
0064
0065   // aerogel + filter
0066   /* AddFlatRadiator will create a pair of flat refractive surfaces internally;
0067    * FIXME: should make a small gas gap at the upstream end of the gas volume;
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,                  // Cherenkov detector
0079       RadiatorName(kAerogel).c_str(), // name
0080       0,                              // path
0081       (G4LogicalVolume*)(0x1),        // G4LogicalVolume (inaccessible? use an integer instead)
0082       nullptr,                        // G4RadiatorMaterial
0083       m_aerogelFlatSurface,           // surface
0084       aerogelThickness                // surface thickness
0085   );
0086   auto* filterFlatRadiator = m_irtDetectorCollection->AddFlatRadiator(
0087       m_irtDetector,           // Cherenkov detector
0088       "Filter",                // name
0089       0,                       // path
0090       (G4LogicalVolume*)(0x2), // G4LogicalVolume (inaccessible? use an integer instead)
0091       nullptr,                 // G4RadiatorMaterial
0092       m_filterFlatSurface,     // surface
0093       filterThickness          // surface thickness
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   // sensor modules: search the detector tree for sensors
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       // get sensor info
0110       auto imod = detSensor.id();
0111       // - get sensor centroid position
0112       auto pvSensor  = detSensor.placement();
0113       auto posSensor = (1 / dd4hep::mm) * (m_posRich + pvSensor.position());
0114       // - get sensor surface position
0115       dd4hep::Direction sensorNorm(
0116           0, 0,
0117           1); // FIXME: generalize; this assumes planar layout, with norm along +z axis (toward IP)
0118       auto surfaceOffset    = sensorNorm.Unit() * (0.5 * sensorThickness);
0119       auto posSensorSurface = posSensor + surfaceOffset;
0120       // - add to `m_sensor_info` map
0121       richgeo::Sensor sensor_info;
0122       sensor_info.size             = sensorSize;
0123       sensor_info.surface_centroid = posSensorSurface;
0124       sensor_info.surface_offset   = surfaceOffset;
0125       m_sensor_info.insert({imod, sensor_info});
0126       // - get surface normal and in-plane vectors
0127       double sensorLocalNormX[3] = {1.0, 0.0, 0.0};
0128       double sensorLocalNormY[3] = {0.0, 1.0, 0.0};
0129       double sensorGlobalNormX[3];
0130       double sensorGlobalNormY[3];
0131       pvSensor.ptr()->LocalToMasterVect(
0132           static_cast<const Double_t*>(sensorLocalNormX),
0133           static_cast<Double_t*>(
0134               sensorGlobalNormX)); // ignore vessel transformation, since it is a pure translation
0135       pvSensor.ptr()->LocalToMasterVect(static_cast<const Double_t*>(sensorLocalNormY),
0136                                         static_cast<Double_t*>(sensorGlobalNormY));
0137
0138       // validate sensor position and normal
0139       // - test normal vectors
0140       dd4hep::Direction normXdir;
0141       dd4hep::Direction normYdir;
0142       normXdir.SetCoordinates(static_cast<const Double_t*>(sensorGlobalNormX));
0143       normYdir.SetCoordinates(static_cast<const Double_t*>(sensorGlobalNormY));
0144       auto normZdir =
0145           normXdir.Cross(normYdir); // sensor surface normal, given derived GlobalNormX,Y
0146       auto testOrtho  = normXdir.Dot(normYdir); // should be zero, if normX and normY are orthogonal
0147       auto testRadial = sensorNorm.Cross(normZdir)
0148                             .Mag2(); // should be zero, if sensor surface normal is as expected
0149       if (std::abs(testOrtho) > 1e-6 || std::abs(testRadial) > 1e-6) {
0150         m_log->error(
0151             "sensor normal is wrong: normX.normY = {:f}   |sensorNorm x normZdir|^2 = {:f}",
0152             testOrtho, testRadial);
0153         return;
0154       }
0155
0156       // create the optical surface
0157       m_sensorFlatSurface = new FlatSurface(
0158           TVector3(posSensorSurface.x(), posSensorSurface.y(), posSensorSurface.z()),
0159           TVector3(sensorGlobalNormX), TVector3(sensorGlobalNormY));
0160       m_irtDetector->CreatePhotonDetectorInstance(0,                   // sector
0161                                                   m_irtPhotonDetector, // CherenkovPhotonDetector
0162                                                   imod,                // copy number
0163                                                   m_sensorFlatSurface  // surface
0164       );
0165       m_log->trace("sensor: id={:#08X} pos=({:5.2f}, {:5.2f}, {:5.2f}) normX=({:5.2f}, {:5.2f}, "
0166                    "{:5.2f}) normY=({:5.2f}, {:5.2f}, {:5.2f})",
0167                    imod, posSensorSurface.x(), posSensorSurface.y(), posSensorSurface.z(),
0168                    normXdir.x(), normXdir.y(), normXdir.z(), normYdir.x(), normYdir.y(),
0169                    normYdir.z());
0170
0171       // complete the radiator volume description; this is the rear side of the container gas volume
0172       // Yes, since there are no mirrors in this detector, just close the gas radiator volume by hand (once),
0173       // assuming that all the sensors will be sitting at roughly the same location along the beam line anyway;
0174       if (firstSensor) {
0175         m_irtDetector->GetRadiator(RadiatorName(kGas).c_str())->m_Borders[0].second =
0176             dynamic_cast<ParametricSurface*>(m_sensorFlatSurface);
0177         firstSensor = false;
0178       }
0179
0180     } // if sensor found
0181   }   // search for sensors
0182
0183   // set reference refractive indices // NOTE: numbers may be overridden externally
0184   std::map<const char*, double> rIndices;
0185   rIndices.insert({RadiatorName(kGas).c_str(), 1.0013});
0186   rIndices.insert({RadiatorName(kAerogel).c_str(), 1.0190});
0187   rIndices.insert({"Filter", 1.5017});
0188   for (auto const& [rName, rIndex] : rIndices) {
0189     auto* rad = m_irtDetector->GetRadiator(rName);
0190     if (rad != nullptr) {
0191       rad->SetReferenceRefractiveIndex(rIndex);
0192     }
0193   }
0194
0195   // set refractive index table
0196   SetRefractiveIndexTable();
0197
0198   // define the `cell ID -> pixel position` converter
0199   SetReadoutIDToPositionLambda();
0200 }
0201
0202 // destructor
0203 richgeo::IrtGeoPFRICH::~IrtGeoPFRICH() {
0204   delete m_surfEntrance;
0205   delete m_irtPhotonDetector;
0206   delete m_aerogelFlatSurface;
0207   delete m_filterFlatSurface;
0208   delete m_sensorFlatSurface;
0209 }