src/todo/ImagingPixelReco.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Chao Peng, Sylvester Joosten, Chao, Whitney Armstrong, Wouter Deconinck, Jihee Kim
0003
0004 // Reconstruct digitized outputs of ImagingCalorimeter
0005 // It converts digitized ADC/TDC values to energy/time, and looks for geometrical information of the
0006 // readout pixels Author: Chao Peng Date: 06/02/2021
0007
0008 #include <algorithm>
0009 #include <bitset>
0010
0011 #include "Gaudi/Property.h"
0012 #include "GaudiAlg/GaudiAlgorithm.h"
0013 #include "GaudiAlg/GaudiTool.h"
0014 #include "GaudiAlg/Transformer.h"
0015 #include "GaudiKernel/PhysicalConstants.h"
0016 #include "GaudiKernel/RndmGenerators.h"
0017 #include "GaudiKernel/ToolHandle.h"
0018
0019 #include "DDRec/CellIDPositionConverter.h"
0020 #include "DDRec/Surface.h"
0021 #include "DDRec/SurfaceManager.h"
0022
0023 #include "JugBase/DataHandle.h"
0024 #include "JugBase/IGeoSvc.h"
0025
0026 // Event Model related classes
0027 #include "edm4eic/CalorimeterHitCollection.h"
0028 #include "edm4eic/RawCalorimeterHitCollection.h"
0029
0030 using namespace Gaudi::Units;
0031
0032 namespace Jug::Reco {
0033
0034 /** Imaging calorimeter pixel hit reconstruction.
0035  *
0036  * Reconstruct digitized outputs of ImagingCalorimeter
0037  * It converts digitized ADC/TDC values to energy/time, and looks for geometrical information of the
0038  *
0039  * \ingroup reco
0040  */
0041 class ImagingPixelReco : public GaudiAlgorithm {
0042 private:
0043   // geometry service
0044   Gaudi::Property<std::string> m_geoSvcName{this, "geoServiceName", "GeoSvc"};
0045   Gaudi::Property<std::string> m_readout{this, "readoutClass", ""};
0046   Gaudi::Property<std::string> m_layerField{this, "layerField", "layer"};
0047   Gaudi::Property<std::string> m_sectorField{this, "sectorField", "sector"};
0048   // length unit (from dd4hep geometry service)
0049   Gaudi::Property<double> m_lUnit{this, "lengthUnit", dd4hep::mm};
0050   // digitization parameters
0051   Gaudi::Property<unsigned int> m_capADC{this, "capacityADC", 8096};
0052   Gaudi::Property<unsigned int> m_pedMeanADC{this, "pedestalMean", 400};
0053   Gaudi::Property<double> m_dyRangeADC{this, "dynamicRangeADC", 100 * MeV};
0054   Gaudi::Property<double> m_pedSigmaADC{this, "pedestalSigma", 3.2};
0055   Gaudi::Property<double> m_thresholdADC{this, "thresholdFactor", 3.0};
0056   // Calibration!
0057   Gaudi::Property<double> m_sampFrac{this, "samplingFraction", 1.0};
0058
0059   // unitless counterparts for the input parameters
0060   double dyRangeADC{0};
0061
0062   // hits containers
0063   DataHandle<edm4eic::RawCalorimeterHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader,
0064                                                                     this};
0065   DataHandle<edm4eic::CalorimeterHitCollection> m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer,
0066                                                                   this};
0067
0068   // Pointer to the geometry service
0069   SmartIF<IGeoSvc> m_geoSvc;
0070   // visit readout fields
0071   dd4hep::BitFieldCoder* id_dec;
0072   size_t sector_idx{0}, layer_idx{0};
0073
0074 public:
0075   ImagingPixelReco(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {
0076     declareProperty("inputHitCollection", m_inputHitCollection, "");
0077     declareProperty("outputHitCollection", m_outputHitCollection, "");
0078   }
0079
0080   StatusCode initialize() override {
0081     if (GaudiAlgorithm::initialize().isFailure()) {
0082       return StatusCode::FAILURE;
0083     }
0084     m_geoSvc = service(m_geoSvcName);
0085     if (!m_geoSvc) {
0086       error() << "Unable to locate Geometry Service. "
0087               << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
0088       return StatusCode::FAILURE;
0089     }
0090
0091     if (m_readout.value().empty()) {
0092       error() << "readoutClass is not provided, it is needed to know the fields in readout ids" << endmsg;
0093       return StatusCode::FAILURE;
0094     }
0095
0096     try {
0097       id_dec     = m_geoSvc->detector()->readout(m_readout).idSpec().decoder();
0098       sector_idx = id_dec->index(m_sectorField);
0099       layer_idx  = id_dec->index(m_layerField);
0100     } catch (...) {
0101       error() << "Failed to load ID decoder for " << m_readout << endmsg;
0102       return StatusCode::FAILURE;
0103     }
0104
0105     // unitless conversion
0106     dyRangeADC = m_dyRangeADC.value() / GeV;
0107
0108     return StatusCode::SUCCESS;
0109   }
0110
0111   StatusCode execute() override {
0112     // input collections
0113     const auto& rawhits = *m_inputHitCollection.get();
0114     // Create output collections
0115     auto& hits = *m_outputHitCollection.createAndPut();
0116
0117     // energy time reconstruction
0118     for (const auto& rh : rawhits) {
0119
0120       #pragma GCC diagnostic push
0121       #pragma GCC diagnostic error "-Wsign-conversion"
0122
0123       // did not pass the threshold
0124       if (rh.getAmplitude() < m_pedMeanADC + m_thresholdADC * m_pedSigmaADC) {
0125         continue;
0126       }
0127       const double energy =
0128         (((signed)rh.getAmplitude() - (signed)m_pedMeanADC)) / (double)m_capADC * dyRangeADC / m_sampFrac; // convert ADC -> energy
0129       const double time = rh.getTimeStamp() * 1.e-6;                                       // ns
0130
0131       #pragma GCC diagnostic pop
0132
0133       const auto id     = rh.getCellID();
0134       // @TODO remove
0135       const int lid = (int)id_dec->get(id, layer_idx);
0136       const int sid = (int)id_dec->get(id, sector_idx);
0137
0138       // global positions
0139       const auto gpos = m_geoSvc->cellIDPositionConverter()->position(id);
0140       // local positions
0141       const auto volman = m_geoSvc->detector()->volumeManager();
0142       // TODO remove
0143       const auto alignment = volman.lookupDetElement(id).nominal();
0144       const auto pos       = alignment.worldToLocal(dd4hep::Position(gpos.x(), gpos.y(), gpos.z()));
0145
0146
0147       // create const vectors for passing to hit initializer list
0148       const decltype(edm4eic::CalorimeterHitData::position) position(
0149         gpos.x() / m_lUnit, gpos.y() / m_lUnit, gpos.z() / m_lUnit
0150       );
0151       const decltype(edm4eic::CalorimeterHitData::local) local(
0152         pos.x() / m_lUnit, pos.y() / m_lUnit, pos.z() / m_lUnit
0153       );
0154
0155       hits.push_back(edm4eic::CalorimeterHit{id,                         // cellID
0156                                           static_cast<float>(energy), // energy
0157                                           0,                          // energyError
0158                                           static_cast<float>(time),   // time
0159                                           0,                          // timeError TODO
0160                                           position,                   // global pos
0161                                           {0, 0, 0}, // @TODO: add dimension
0162                                           sid,lid,
0163                                           local});                    // local pos
0164     }
0165     return StatusCode::SUCCESS;
0166   }
0167 };
0168
0169 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0170 DECLARE_COMPONENT(ImagingPixelReco)
0171
0172 } // namespace Jug::Reco