File indexing completed on 2024-09-28 07:03:04
0001
0002
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
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
0031 std::transform(m_detName.begin(), m_detName.end(), m_detName.begin(), ::toupper);
0032
0033
0034 m_random.SetSeed(1);
0035
0036
0037 m_loopCellIDs = [] (std::function<void(CellIDType)> lambda) { return; };
0038
0039
0040 m_rngCellIDs = [] (std::function<void(CellIDType)> lambda, float p) { return; };
0041
0042
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
0048 if(m_detName=="DRICH") {
0049
0050
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
0058 m_loopCellIDs = [this] (std::function<void(CellIDType)> lambda) {
0059 m_log->trace("call VisitAllReadoutPixels for systemID = {} = {}", m_systemID, m_detName);
0060
0061
0062 for(auto const& [deName, detSensor] : m_detRich.children()) {
0063 if(deName.find("sensor_de_sec")!=std::string::npos) {
0064
0065
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
0071
0072
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
0079 lambda(cellID);
0080 }
0081 }
0082 }
0083 }
0084 };
0085
0086
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
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
0119
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
0132
0133 dd4hep::Position richgeo::ReadoutGeo::GetSensorLocalPosition(CellIDType cellID, dd4hep::Position pos) {
0134
0135
0136 auto context = m_conv->findContext(cellID);
0137
0138
0139 double xyz_l[3], xyz_e[3], xyz_g[3];
0140 double pv_g[3], pv_l[3];
0141
0142
0143 auto sensor_elem = context->element;
0144 sensor_elem.placement().position().GetCoordinates(xyz_l);
0145
0146
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
0155 dd4hep::Direction pos_pv = pos - pos_sensor;
0156
0157
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
0164
0165
0166
0167
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179 return pos_transformed;
0180 }