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