Skip to content
Snippets Groups Projects
ImagingPixelReco.cpp 5.46 KiB
Newer Older
  • Learn to ignore specific revisions
  • Chao Peng's avatar
    Chao Peng committed
    // Reconstruct digitized outputs of ImagingCalorimeter
    
    Sylvester Joosten's avatar
    Sylvester Joosten committed
    // It converts digitized ADC/TDC values to energy/time, and looks for geometrical information of the
    // readout pixels Author: Chao Peng Date: 06/02/2021
    
    Chao Peng's avatar
    Chao Peng committed
    
    #include <algorithm>
    #include <bitset>
    
    #include "Gaudi/Property.h"
    #include "GaudiAlg/GaudiAlgorithm.h"
    #include "GaudiAlg/GaudiTool.h"
    #include "GaudiAlg/Transformer.h"
    #include "GaudiKernel/PhysicalConstants.h"
    #include "GaudiKernel/RndmGenerators.h"
    #include "GaudiKernel/ToolHandle.h"
    
    #include "DDRec/CellIDPositionConverter.h"
    #include "DDRec/Surface.h"
    #include "DDRec/SurfaceManager.h"
    
    // FCCSW
    #include "JugBase/DataHandle.h"
    #include "JugBase/IGeoSvc.h"
    
    // Event Model related classes
    #include "eicd/ImagingPixelCollection.h"
    #include "eicd/RawCalorimeterHitCollection.h"
    
    using namespace Gaudi::Units;
    
    namespace Jug::Reco {
    
    
    Sylvester Joosten's avatar
    Sylvester Joosten committed
      class ImagingPixelReco : public GaudiAlgorithm {
      public:
    
    Chao Peng's avatar
    Chao Peng committed
        // geometry service
    
    Sylvester Joosten's avatar
    Sylvester Joosten committed
        Gaudi::Property<std::string> m_geoSvcName{this, "geoServiceName", "GeoSvc"};
        Gaudi::Property<std::string> m_readout{this, "readoutClass", ""};
        Gaudi::Property<std::string> m_layerField{this, "layerField", "layer"};
        Gaudi::Property<std::string> m_sectorField{this, "sectorField", "sector"};
    
    Chao Peng's avatar
    Chao Peng committed
        // length unit (from dd4hep geometry service)
    
    Sylvester Joosten's avatar
    Sylvester Joosten committed
        Gaudi::Property<double> m_lUnit{this, "lengthUnit", dd4hep::mm};
    
    Chao Peng's avatar
    Chao Peng committed
        // digitization parameters
    
    Sylvester Joosten's avatar
    Sylvester Joosten committed
        Gaudi::Property<int>    m_capADC{this, "capacityADC", 8096};
        Gaudi::Property<int>    m_pedMeanADC{this, "pedestalMean", 400};
        Gaudi::Property<double> m_dyRangeADC{this, "dynamicRangeADC", 100 * MeV};
        Gaudi::Property<double> m_pedSigmaADC{this, "pedestalSigma", 3.2};
        Gaudi::Property<double> m_thresholdADC{this, "thresholdFactor", 3.0};
    
    Chao Peng's avatar
    Chao Peng committed
        // hits containers
    
    Sylvester Joosten's avatar
    Sylvester Joosten committed
        DataHandle<eic::RawCalorimeterHitCollection> m_inputHitCollection{
            "inputHitCollection", Gaudi::DataHandle::Reader, this};
        DataHandle<eic::ImagingPixelCollection> m_outputHitCollection{"outputHitCollection",
                                                                      Gaudi::DataHandle::Writer, this};
    
    Chao Peng's avatar
    Chao Peng committed
    
        // Pointer to the geometry service
        SmartIF<IGeoSvc> m_geoSvc;
        // visit readout fields
    
    Sylvester Joosten's avatar
    Sylvester Joosten committed
        dd4hep::BitFieldCoder* id_dec;
        size_t                 sector_idx, layer_idx;
    
    Chao Peng's avatar
    Chao Peng committed
    
        ImagingPixelReco(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc)
        {
    
    Sylvester Joosten's avatar
    Sylvester Joosten committed
          declareProperty("inputHitCollection", m_inputHitCollection, "");
          declareProperty("outputHitCollection", m_outputHitCollection, "");
    
    Chao Peng's avatar
    Chao Peng committed
        }
    
        StatusCode initialize() override
        {
    
    Sylvester Joosten's avatar
    Sylvester Joosten committed
          if (GaudiAlgorithm::initialize().isFailure()) {
            return StatusCode::FAILURE;
          }
          m_geoSvc = service(m_geoSvcName);
          if (!m_geoSvc) {
            error() << "Unable to locate Geometry Service. "
                    << "Make sure you have GeoSvc and SimSvc in the right order in the configuration."
                    << endmsg;
            return StatusCode::FAILURE;
          }
    
          if (m_readout.value().empty()) {
            error() << "readoutClass is not provided, it is needed to know the fields in readout ids"
                    << endmsg;
            return StatusCode::FAILURE;
          }
    
          try {
            id_dec     = m_geoSvc->detector()->readout(m_readout).idSpec().decoder();
            sector_idx = id_dec->index(m_sectorField);
            layer_idx  = id_dec->index(m_layerField);
          } catch (...) {
            error() << "Failed to load ID decoder for " << m_readout << endmsg;
            return StatusCode::FAILURE;
          }
    
          return StatusCode::SUCCESS;
    
    Chao Peng's avatar
    Chao Peng committed
        }
    
        StatusCode execute() override
        {
    
    Sylvester Joosten's avatar
    Sylvester Joosten committed
          // input collections
          const auto& rawhits = *m_inputHitCollection.get();
          // Create output collections
          auto& hits = *m_outputHitCollection.createAndPut();
    
          // energy time reconstruction
          for (const auto& rh : rawhits) {
            // did not pass the threshold
            if ((rh.amplitude() - m_pedMeanADC) < m_thresholdADC * m_pedSigmaADC) {
              continue;
    
    Chao Peng's avatar
    Chao Peng committed
            }
    
    Sylvester Joosten's avatar
    Sylvester Joosten committed
            double edep = (rh.amplitude() - m_pedMeanADC) / (double)m_capADC *
                          m_dyRangeADC;   // convert ADC -> energy
            double time = rh.timeStamp(); // ns
            auto   id   = rh.cellID();
            int    lid  = (int)id_dec->get(id, layer_idx);
            int    sid  = (int)id_dec->get(id, sector_idx);
    
            // global positions
            auto gpos = m_geoSvc->cellIDPositionConverter()->position(id);
            // local positions
            auto volman    = m_geoSvc->detector()->volumeManager();
            auto alignment = volman.lookupDetElement(id).nominal();
            auto pos       = alignment.worldToLocal(dd4hep::Position(gpos.x(), gpos.y(), gpos.z()));
            // polar coordinates
            double r   = std::sqrt(gpos.x() * gpos.x() + gpos.y() * gpos.y() + gpos.z() * gpos.z());
            double th  = std::acos(gpos.z() / r);
            double eta = -std::log(std::tan(th / 2.));
            double phi = std::atan2(gpos.y(), gpos.x());
    
            hits.push_back(eic::ImagingPixel{
                -1,
                lid,
                sid,
                -1, // cluster id, layer id, sector id, hit id
                edep,
                time,
                eta, // edep, time, pseudo-rapidity
                {pos.x() / m_lUnit, pos.y() / m_lUnit, pos.z() / m_lUnit},    // local pos
                {gpos.x() / m_lUnit, gpos.y() / m_lUnit, gpos.z() / m_lUnit}, // global pos
                {r, th, phi}                                                  // polar global pos
            });
          }
          return StatusCode::SUCCESS;
    
    Chao Peng's avatar
    Chao Peng committed
        }
    
    Sylvester Joosten's avatar
    Sylvester Joosten committed
      };
    
    Chao Peng's avatar
    Chao Peng committed
    
    
    Sylvester Joosten's avatar
    Sylvester Joosten committed
      DECLARE_COMPONENT(ImagingPixelReco)
    
    Chao Peng's avatar
    Chao Peng committed
    
    } // namespace Jug::Reco