geometry/richgeo/ReadoutGeo.cc

0001 // Copyright (C) 2023, Christopher Dilks, Luigi Dello Stritto
0002 // Subject to the terms in the LICENSE file found in the top-level directory.
0003 //
0004 //
0005
0006 #include "ReadoutGeo.h"
0007
0008 #include <DD4hep/Alignments.h>
0009 #include <DD4hep/Fields.h>
0010 #include <DD4hep/IDDescriptor.h>
0011 #include <DD4hep/Readout.h>
0012 #include <DD4hep/VolumeManager.h>
0013 #include <DD4hep/Volumes.h>
0014 #include <Evaluator/DD4hepUnits.h>
0015 #include <Math/GenVector/Cartesian3D.h>
0016 #include <Math/GenVector/DisplacementVector3D.h>
0017 #include <TGeoMatrix.h>
0018 #include <ctype.h>
0019 #include <fmt/core.h>
0020 #include <algorithm>
0021 #include <cmath>
0022 #include <map>
0023
0024 #include "services/geometry/richgeo/RichGeo.h"
0025
0026 // constructor
0027 richgeo::ReadoutGeo::ReadoutGeo(std::string detName_, gsl::not_null<const dd4hep::Detector*> det_, gsl::not_null<const dd4hep::rec::CellIDPositionConverter*> conv_, std::shared_ptr<spdlog::logger> log_)
0028   : m_detName(detName_), m_det(det_), m_conv(conv_), m_log(log_)
0029 {
0030   // capitalize m_detName
0031   std::transform(m_detName.begin(), m_detName.end(), m_detName.begin(), ::toupper);
0032
0033   // random number generators
0034   m_random.SetSeed(1); // default seed
0035
0036   // default (empty) cellID looper
0037   m_loopCellIDs = [] (std::function<void(CellIDType)> lambda) { return; };
0038
0039   // default (empty) cellID rng generator
0040   m_rngCellIDs = [] (std::function<void(CellIDType)> lambda, float p) { return; };
0041
0042   // common objects
0043   m_readoutCoder = m_det->readout(m_detName+"Hits").idSpec().decoder();
0044   m_detRich      = m_det->detector(m_detName);
0045   m_systemID     = m_detRich.id();
0046
0047   // dRICH readout --------------------------------------------------------------------
0048   if(m_detName=="DRICH") {
0049
0050     // get constants from geometry
0051     m_num_sec           = m_det->constant<int>("DRICH_num_sectors");
0052     m_num_pdus          = m_det->constant<int>("DRICH_num_pdus");
0053     m_num_sipms_per_pdu = std::pow(m_det->constant<int>("DRICH_pdu_num_sensors"), 2);
0054     m_num_px            = m_det->constant<int>("DRICH_num_px");
0055     m_pixel_size        = m_det->constant<double>("DRICH_pixel_size") / dd4hep::mm;
0056
0057     // define cellID looper
0058     m_loopCellIDs = [this] (std::function<void(CellIDType)> lambda) {
0059       m_log->trace("call VisitAllReadoutPixels for systemID = {} = {}", m_systemID, m_detName);
0060
0061       // loop over sensors (for all sectors)
0062       for(auto const& [deName, detSensor] : m_detRich.children()) {
0063         if(deName.find("sensor_de_sec")!=std::string::npos) {
0064
0065           // decode `sensorID` to module number and sector number
0066           auto sensorID = detSensor.id();
0067           auto ipdu     = m_readoutCoder->get(sensorID, "pdu");
0068           auto isipm    = m_readoutCoder->get(sensorID, "sipm");
0069           auto isec     = m_readoutCoder->get(sensorID, "sector");
0070           // m_log->trace("  module: sensorID={:#018X} => ipdu={:<6} isipm={:<6} isec={:<2} name={}", sensorID, ipdu, isipm, isec, deName);
0071
0072           // loop over xy-segmentation
0073           for (int x = 0; x < m_num_px; x++) {
0074             for (int y = 0; y < m_num_px; y++) {
0075
0076               auto cellID = cellIDEncoding(isec, ipdu, isipm, x, y);
0077
0078               // then execute the user's lambda function
0079               lambda(cellID);
0080             }
0081           } // end xy-segmentation loop
0082         }
0083       } // end sensor loop (for all sectors)
0084     }; // end definition of m_loopCellIDs
0085
0086     // define k random cell IDs generator
0087     m_rngCellIDs = [this] (std::function<void(CellIDType)> lambda, float p) {
0088       m_log->trace("call RngReadoutPixels for systemID = {} = {}", m_systemID, m_detName);
0089
0090       int k = p * m_num_sec * m_num_pdus * m_num_sipms_per_pdu * m_num_px * m_num_px;
0091
0092       for (int i = 0; i < k; i++) {
0093         int isec = m_random.Uniform(0., m_num_sec);
0094         int ipdu = m_random.Uniform(0., m_num_pdus);
0095         int isipm = m_random.Uniform(0., m_num_sipms_per_pdu);
0096         int x = m_random.Uniform(0., m_num_px);
0097         int y = m_random.Uniform(0., m_num_px);
0098
0099         auto cellID = cellIDEncoding(isec, ipdu, isipm, x, y);
0100
0101         lambda(cellID);
0102       }
0103     };
0104
0105   }
0106
0107   // pfRICH readout --------------------------------------------------------------------
0108   else if(m_detName=="PFRICH") {
0109     m_log->error("TODO: pfRICH readout bindings have not yet been implemented");
0110   }
0111
0112   // ------------------------------------------------------------------------------------------------
0113   else m_log->error("ReadoutGeo is not defined for detector '{}'",m_detName);
0114
0115 }
0116
0117
0118 // pixel gap mask
0119 // FIXME: generalize; this assumes the segmentation is `CartesianGridXY`
0120 bool richgeo::ReadoutGeo::PixelGapMask(CellIDType cellID, dd4hep::Position pos_hit_global) {
0121   auto pos_pixel_global = m_conv->position(cellID);
0122   auto pos_pixel_local  = GetSensorLocalPosition(cellID, pos_pixel_global);
0123   auto pos_hit_local    = GetSensorLocalPosition(cellID, pos_hit_global);
0124   return ! (
0125       std::abs( pos_hit_local.x()/dd4hep::mm - pos_pixel_local.x()/dd4hep::mm ) > m_pixel_size/2 ||
0126       std::abs( pos_hit_local.y()/dd4hep::mm - pos_pixel_local.y()/dd4hep::mm ) > m_pixel_size/2
0127       );
0128 }
0129
0130
0131 // transform global position `pos` to sensor `cellID` frame position
0132 // IMPORTANT NOTE: this has only been tested for the dRICH; if you use it, test it carefully...
0133 dd4hep::Position richgeo::ReadoutGeo::GetSensorLocalPosition(CellIDType cellID, dd4hep::Position pos) {
0134
0135   // get the VolumeManagerContext for this sensitive detector
0136   auto context = m_conv->findContext(cellID);
0137
0138   // transformation vector buffers
0139   double xyz_l[3], xyz_e[3], xyz_g[3];
0140   double pv_g[3], pv_l[3];
0141
0142   // get sensor position w.r.t. its parent
0143   auto sensor_elem = context->element;
0144   sensor_elem.placement().position().GetCoordinates(xyz_l);
0145
0146   // convert sensor position to global position (cf. `CellIDPositionConverter::positionNominal()`)
0147   const auto& volToElement = context->toElement();
0148   volToElement.LocalToMaster(xyz_l, xyz_e);
0149   const auto& elementToGlobal = sensor_elem.nominal().worldTransformation();
0150   elementToGlobal.LocalToMaster(xyz_e, xyz_g);
0151   dd4hep::Position pos_sensor;
0152   pos_sensor.SetCoordinates(xyz_g);
0153
0154   // get the position vector of `pos` w.r.t. the sensor position `pos_sensor`
0155   dd4hep::Direction pos_pv = pos - pos_sensor;
0156
0157   // then transform it to the sensor's local frame
0158   pos_pv.GetCoordinates(pv_g);
0159   volToElement.MasterToLocalVect(pv_g, pv_l);
0160   dd4hep::Position pos_transformed;
0161   pos_transformed.SetCoordinates(pv_l);
0162
0163   // trace log
0164   /*
0165   if(m_log->level() <= spdlog::level::trace) {
0166     m_log->trace("pixel hit on cellID={:#018x}",cellID);
0167     auto print_pos = [&] (std::string name, dd4hep::Position p) {
0168       m_log->trace("  {:>30} x={:.2f} y={:.2f} z={:.2f} [mm]: ", name, p.x()/dd4hep::mm,  p.y()/dd4hep::mm,  p.z()/dd4hep::mm);
0169     };
0170     print_pos("input position",  pos);
0171     print_pos("sensor position", pos_sensor);
0172     print_pos("output position", pos_transformed);
0173     // auto dim = m_conv->cellDimensions(cellID);
0174     // for (std::size_t j = 0; j < std::size(dim); ++j)
0175     //   m_log->trace("   - dimension {:<5} size: {:.2}",  j, dim[j]);
0176   }
0177   */
0178
0179   return pos_transformed;
0180 }