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