Skip to content
Snippets Groups Projects
TrackerSourceLinker.cpp 9.3 KiB
Newer Older
  • Learn to ignore specific revisions
  • // This file is part of the Acts project.
    //
    // Copyright (C) 2019 CERN for the benefit of the Acts project
    //
    // This Source Code Form is subject to the terms of the Mozilla Public
    // License, v. 2.0. If a copy of the MPL was not distributed with this
    // file, You can obtain one at http://mozilla.org/MPL/2.0/.
    
    //#include "HitSmearing.hpp"
    
    #include "JugReco/GeometryContainers.hpp"
    
    // Gaudi
    #include "GaudiAlg/GaudiAlgorithm.h"
    #include "GaudiKernel/ToolHandle.h"
    #include "GaudiAlg/Transformer.h"
    #include "GaudiAlg/GaudiTool.h"
    #include "GaudiKernel/RndmGenerators.h"
    #include "GaudiKernel/Property.h"
    
    #include "JugBase/DataHandle.h"
    #include "JugBase/IGeoSvc.h"
    
    #include "DDRec/CellIDPositionConverter.h"
    #include "DDRec/SurfaceManager.h"
    #include "DDRec/Surface.h"
    
    #include "DD4hep/Volumes.h"
    #include "DD4hep/DD4hepUnits.h"
    
    
    #include "Acts/Utilities/Units.hpp"
    
    #include "Acts/Utilities/Definitions.hpp"
    #include "Acts/Geometry/TrackingGeometry.hpp"
    #include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp"
    
    #include "JugReco/SourceLinks.h"
    
    #include "eicd/TrackerHitCollection.h"
    
    namespace Jug::Reco {
    
      class TrackerSourceLinker : public GaudiAlgorithm {
      public:
        DataHandle<eic::TrackerHitCollection>    m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
        DataHandle<SourceLinkContainer> m_outputSourceLinks{"outputSourceLinks", Gaudi::DataHandle::Writer, this};
        /// Pointer to the geometry service
        SmartIF<IGeoSvc> m_geoSvc;
    
        /// Lookup container for hit surfaces that generate smeared hits
        std::unordered_map<uint64_t, const Acts::Surface*> m_surfaces;
    
      public:
        TrackerSourceLinker(const std::string& name, ISvcLocator* svcLoc)
            : GaudiAlgorithm(name, svcLoc) {
          declareProperty("inputHitCollection", m_inputHitCollection, "");
          declareProperty("outputSourceLinks", m_outputSourceLinks, "");
        }
    
        StatusCode initialize() override {
          if (GaudiAlgorithm::initialize().isFailure())
            return StatusCode::FAILURE;
          m_geoSvc = service("GeoSvc");
          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;
          }
          debug() << "visiting all the surfaces  " << endmsg;
          m_geoSvc->trackingGeometry()->visitSurfaces([this](const Acts::Surface* surface) {
            // for now we just require a valid surface
            if (not surface) {
              return;
            }
            auto det_element = dynamic_cast<const Acts::DD4hepDetectorElement*>(surface->associatedDetectorElement());
            if(!det_element) {
              debug() << "invalid det_element!!! " << endmsg;
              return;
            }
            this->m_surfaces.insert_or_assign(det_element->identifier(), surface);
          });
    
          return StatusCode::SUCCESS;
        }
    
        StatusCode execute() override {
          // input collection
          const eic::TrackerHitCollection* hits = m_inputHitCollection.get();
          // Create output collections
          auto source_links = m_outputSourceLinks.createAndPut();
          // setup local covariance
          // TODO add support for per volume/layer/module settings
    
    
          debug() << (*hits).size() << " hits " << endmsg;
    
          for(const auto& ahit : *hits) {
    
            Acts::BoundMatrix cov           = Acts::BoundMatrix::Zero();
    
            cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = ahit.covMatrix(0)*Acts::UnitConstants::mm*ahit.covMatrix(0)*Acts::UnitConstants::mm;
            cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = ahit.covMatrix(1)*Acts::UnitConstants::mm*ahit.covMatrix(1)*Acts::UnitConstants::mm;
    
            auto vol_ctx = m_geoSvc->cellIDPositionConverter()->findContext(ahit.cellID());
    
            auto vol_id = vol_ctx->identifier;
    
            debug() << " hit          : \n" <<  ahit << endmsg;
            //debug() << "cell ID : " << ahit.cellID() << endmsg;
            //debug() << " position : (" <<  ahit.position(0) << ", " <<  ahit.position(1) << ", "<<  ahit.position(2) << ") " << endmsg;
            debug() << " vol_id       : " <<  vol_id << endmsg;
            debug() << " placment pos : " << vol_ctx->volumePlacement().position() << endmsg;
    
    
            const auto is = m_surfaces.find(vol_id);
            if (is == m_surfaces.end()) {
    
              debug() << " vol_id (" <<  vol_id << ")  not found in m_surfaces." <<endmsg;
    
              continue;
            }
            const Acts::Surface* surface = is->second;
    
    
            debug() << " surface center : " << surface->center(Acts::GeometryContext()) << endmsg;
    
            // transform global position into local coordinates
            Acts::Vector2D pos(0, 0);
            // geometry context contains nothing here
    
            pos = surface->globalToLocal(Acts::GeometryContext(), {ahit.position(0), ahit.position(1), ahit.position(2)}, {0, 0, 0}).value();//, pos);
    
    
            //// smear truth to create local measurement
            Acts::BoundVector loc = Acts::BoundVector::Zero();
    
            loc[Acts::eBoundLoc0]     = pos[0] ;//+ m_cfg.sigmaLoc0 * stdNormal(rng);
            loc[Acts::eBoundLoc0]     = pos[1] ;//+ m_cfg.sigmaLoc1 * stdNormal(rng);
    
            debug() << "loc : (" << loc[0] << ", " << loc[1] << ")" << endmsg;
    
    
            // create source link at the end of the container
            //auto it = source_links->emplace_hint(source_links->end(), *surface, hit, 2, loc, cov);
            auto it = source_links->emplace_hint(source_links->end(), *surface,  2, loc, cov);
            // ensure hits and links share the same order to prevent ugly surprises
            if (std::next(it) != source_links->end()) {
              error() << "The hit ordering broke. Run for your life." << endmsg;
              return StatusCode::FAILURE;
            }
    
            //std::array<double,3> posarr; pos.GetCoordinates(posarr);
            //std::array<double,3> dimarr; dim.GetCoordinates(posarr);
            //eic::TrackerHit hit;
            //eic::TrackerHit hit((long long)ahit.cellID0(), (long long)ahit.cellID(), (long long)ahit.time(),
            //                    (float)ahit.charge() / 10000.0, (float)0.0, {{pos.x(), pos.y(),pos.z()}},{{dim[0],dim[1],0.0}});
            //rec_hits->push_back(hit);
          }
          return StatusCode::SUCCESS;
        }
    
       };
      DECLARE_COMPONENT(TrackerSourceLinker)
    
    } // namespace Jug::reco
    
    //HitSmearing::HitSmearing(const Config& cfg, Acts::Logging::Level lvl)
    //    : BareAlgorithm("HitSmearing", lvl), m_cfg(cfg) {
    //  if (m_cfg.inputSimulatedHits.empty()) {
    //    throw std::invalid_argument("Missing input simulated hits collection");
    //  }
    //  if (m_cfg.outputSourceLinks.empty()) {
    //    throw std::invalid_argument("Missing output source links collection");
    //  }
    //  if ((m_cfg.sigmaLoc0 < 0) or (m_cfg.sigmaLoc1 < 0)) {
    //    throw std::invalid_argument("Invalid resolution setting");
    //  }
    //  if (not m_cfg.trackingGeometry) {
    //    throw std::invalid_argument("Missing tracking geometry");
    //  }
    //  if (!m_cfg.randomNumbers) {
    //    throw std::invalid_argument("Missing random numbers tool");
    //  }
    //  // fill the surface map to allow lookup by geometry id only
    //  m_cfg.trackingGeometry->visitSurfaces([this](const Acts::Surface* surface) {
    //    // for now we just require a valid surface
    //    if (not surface) {
    //      return;
    //    }
    //    this->m_surfaces.insert_or_assign(surface->associatedDetectorElement()->identifier(), surface);
    //  });
    //     }
    //}
    
    //ProcessCode HitSmearing::execute(const AlgorithmContext& ctx) const {
    //  // setup input and output containers
    //  const auto& hits =
    //      ctx.eventStore.get<SimHitContainer>(m_cfg.inputSimulatedHits);
    //  SimSourceLinkContainer sourceLinks;
    //  sourceLinks.reserve(hits.size());
    //
    //  // setup random number generator
    //  auto rng = m_cfg.randomNumbers->spawnGenerator(ctx);
    //  std::normal_distribution<double> stdNormal(0.0, 1.0);
    //
    //  // setup local covariance
    //  // TODO add support for per volume/layer/module settings
    //  Acts::BoundMatrix cov = Acts::BoundMatrix::Zero();
    //  cov(Acts::eLOC_0, Acts::eLOC_0) = m_cfg.sigmaLoc0 * m_cfg.sigmaLoc0;
    //  cov(Acts::eLOC_1, Acts::eLOC_1) = m_cfg.sigmaLoc1 * m_cfg.sigmaLoc1;
    //
    //  for (auto&& [moduleGeoId, moduleHits] : groupByModule(hits)) {
    //    // check if we should create hits for this surface
    //    const auto is = m_surfaces.find(moduleGeoId);
    //    if (is == m_surfaces.end()) {
    //      continue;
    //    }
    //
    //    // smear all truth hits for this module
    //    const Acts::Surface* surface = is->second;
    //    for (const auto& hit : moduleHits) {
    //      // transform global position into local coordinates
    //      Acts::Vector2D pos(0, 0);
    //      surface->globalToLocal(ctx.geoContext, hit.position(),
    //                             hit.unitDirection(), pos);
    //
    //      // smear truth to create local measurement
    //      Acts::BoundVector loc = Acts::BoundVector::Zero();
    //      loc[Acts::eLOC_0] = pos[0] + m_cfg.sigmaLoc0 * stdNormal(rng);
    //      loc[Acts::eLOC_1] = pos[1] + m_cfg.sigmaLoc1 * stdNormal(rng);
    //
    //      // create source link at the end of the container
    //      auto it = sourceLinks.emplace_hint(sourceLinks.end(), *surface, hit, 2,
    //                                         loc, cov);
    //      // ensure hits and links share the same order to prevent ugly surprises
    //      if (std::next(it) != sourceLinks.end()) {
    //        ACTS_FATAL("The hit ordering broke. Run for your life.");
    //        return ProcessCode::ABORT;
    //      }
    //    }
    //  }
    //
    //  ctx.eventStore.add(m_cfg.outputSourceLinks, std::move(sourceLinks));
    //  return ProcessCode::SUCCESS;
    //}