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