src/components/TrackParamImagingClusterInit.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Whitney Armstrong, Wouter Deconinck, Sylvester Joosten
0003
0004 #include <cmath>
0005 // Gaudi
0006 #include "Gaudi/Algorithm.h"
0007 #include "GaudiKernel/ToolHandle.h"
0008 #include "GaudiKernel/RndmGenerators.h"
0009 #include "Gaudi/Property.h"
0010
0011 #include <k4FWCore/DataHandle.h>
0012 #include <k4Interface/IGeoSvc.h>
0013 #include "ActsExamples/EventData/Track.hpp"
0014 #include "Acts/Definitions/Units.hpp"
0015 #include "Acts/Definitions/Common.hpp"
0016
0017 #include "edm4eic/TrackerHitCollection.h"
0018 #include "edm4eic/ClusterCollection.h"
0019 #include "edm4hep/utils/vector_utils.h"
0020
0021 #include "Acts/Surfaces/PerigeeSurface.hpp"
0022
0023   ///// (Reconstructed) track parameters e.g. close to the vertex.
0024   //using TrackParameters = Acts::CurvilinearTrackParameters;
0025
0026   ///// Container of reconstructed track states for multiple tracks.
0027   //using TrackParametersContainer = std::vector<TrackParameters>;
0028
0029   ///// MultiTrajectory definition
0030   //using Trajectory = Acts::MultiTrajectory<SourceLink>;
0031
0032   ///// Container for the truth fitting/finding track(s)
0033   //using TrajectoryContainer = std::vector<SimMultiTrajectory>;
0034
0035 namespace Jug::Reco {
0036
0037   /** Initial Track parameters from MC truth.
0038    *
0039    *  TrackParmetersContainer
0040    *
0041    *  \ingroup tracking
0042    */
0043   class TrackParamImagingClusterInit : public Gaudi::Algorithm {
0044   private:
0045     mutable DataHandle<edm4eic::ClusterCollection>          m_inputClusters{"inputClusters", Gaudi::DataHandle::Reader, this};
0046     mutable DataHandle<ActsExamples::TrackParametersContainer> m_outputInitialTrackParameters{"outputInitialTrackParameters",
0047                                                                         Gaudi::DataHandle::Writer, this};
0048
0049   public:
0050     TrackParamImagingClusterInit(const std::string& name, ISvcLocator* svcLoc)
0051         : Gaudi::Algorithm(name, svcLoc) {
0052       declareProperty("inputClusters", m_inputClusters, "Input clusters");
0053       declareProperty("outputInitialTrackParameters", m_outputInitialTrackParameters, "");
0054     }
0055
0056     StatusCode initialize() override {
0057       if (Gaudi::Algorithm::initialize().isFailure()) {
0058         return StatusCode::FAILURE;
0059       }
0060       IRndmGenSvc* randSvc = Gaudi::svcLocator()->service<IRndmGenSvc>("RndmGenSvc", true);
0061       if (randSvc == nullptr) {
0062         return StatusCode::FAILURE;
0063       }
0064       return StatusCode::SUCCESS;
0065     }
0066
0067     StatusCode execute(const EventContext&) const override {
0068       // input collection
0069       const auto* const clusters = m_inputClusters.get();
0070       // Create output collections
0071       auto* init_trk_params = m_outputInitialTrackParameters.createAndPut();
0072
0073       for(const auto& c : *clusters) {
0074
0075         using Acts::UnitConstants::GeV;
0076         using Acts::UnitConstants::MeV;
0077         using Acts::UnitConstants::mm;
0078         using Acts::UnitConstants::ns;
0079
0080         const double p = c.getEnergy()*GeV;
0081         // FIXME hardcoded value
0082         if( p < 0.1*GeV) {
0083           continue;
0084         }
0085         const double theta = edm4hep::utils::anglePolar(c.getPosition());
0086         const double phi = edm4hep::utils::angleAzimuthal(c.getPosition());
0087
0088         Acts::BoundVector  params;
0089         params(Acts::eBoundLoc0)   = 0.0 * mm ;
0090         params(Acts::eBoundLoc1)   = 0.0 * mm ;
0091         params(Acts::eBoundPhi)    = phi;
0092         params(Acts::eBoundTheta)  = theta;
0093         params(Acts::eBoundQOverP) = 1/p;
0094         params(Acts::eBoundTime)   = 0 * ns;
0095
0096         auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3{0,0,0});
0097
0098         debug() << "Invoke track finding seeded by truth particle with p = " << p/GeV  << " GeV" << endmsg;
0099
0100         // add both charges to the track candidate...
0101         init_trk_params->push_back({pSurface, params, {}, Acts::ParticleHypothesis::pion()});
0102
0103         Acts::BoundVector  params2;
0104         params2(Acts::eBoundLoc0)   = 0.0 * mm ;
0105         params2(Acts::eBoundLoc1)   = 0.0 * mm ;
0106         params2(Acts::eBoundPhi)    = phi;
0107         params2(Acts::eBoundTheta)  = theta;
0108         params2(Acts::eBoundQOverP) = -1/p;
0109         params2(Acts::eBoundTime)   = 0 * ns;
0110         init_trk_params->push_back({pSurface, params2, {}, Acts::ParticleHypothesis::pion()});
0111
0112       }
0113       return StatusCode::SUCCESS;
0114     }
0115   };
0116   // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0117   DECLARE_COMPONENT(TrackParamImagingClusterInit)
0118
0119 } // namespace Jug::reco