src/components/ConformalXYPeakProtoTracks.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Whitney Armstrong, Sylvester Joosten
0003
0004 #include <cmath>
0005 // Gaudi
0006 #include "Gaudi/Property.h"
0007 #include "Gaudi/Algorithm.h"
0008 #include "GaudiKernel/ToolHandle.h"
0009
0010 #include <k4FWCore/DataHandle.h>
0011 #include <k4Interface/IGeoSvc.h>
0012 #include "ActsExamples/EventData/ProtoTrack.hpp"
0013 #include "ActsExamples/EventData/Track.hpp"
0014
0015 #include "TH1F.h"
0016 #include "Math/Vector3D.h"
0017 #include "Math/Vector2D.h"
0018
0019 #include "edm4eic/TrackerHitCollection.h"
0020
0021 namespace Jug::Reco {
0022
0023 /** Conformal XY hits.
0024  *
0025  *  Conformal mapping which turns circles to lines.
0026  *
0027  *  \ingroup tracking
0028  */
0029 class ConformalXYPeakProtoTracks : public Gaudi::Algorithm {
0030 private:
0031   mutable DataHandle<edm4eic::TrackerHitCollection> m_inputTrackerHits{"inputTrackerHits", Gaudi::DataHandle::Reader, this};
0032   mutable DataHandle<ActsExamples::ProtoTrackContainer> m_outputProtoTracks{"outputProtoTracks", Gaudi::DataHandle::Writer, this};
0033   mutable DataHandle<int> m_nProtoTracks{"nProtoTracks", Gaudi::DataHandle::Writer, this};
0034
0035   Gaudi::Property<int> m_nPhiBins{this, "nPhiBins", 100};
0036
0037   using ConformalHit = ROOT::Math::XYVector;
0038
0039 public:
0040   ConformalXYPeakProtoTracks(const std::string& name, ISvcLocator* svcLoc) : Gaudi::Algorithm(name, svcLoc) {
0041     declareProperty("inputTrackerHits", m_inputTrackerHits, "tracker hits whose indices are used in proto-tracks");
0042     declareProperty("outputProtoTracks", m_outputProtoTracks, "grouped hit indicies");
0043     declareProperty("nProtoTracks", m_nProtoTracks, "number of proto tracks");
0044   }
0045
0046   StatusCode initialize() override {
0047     if (Gaudi::Algorithm::initialize().isFailure()) {
0048       return StatusCode::FAILURE;
0049     }
0050     return StatusCode::SUCCESS;
0051   }
0052
0053   StatusCode execute(const EventContext&) const override {
0054     // input collection
0055     const edm4eic::TrackerHitCollection* hits = m_inputTrackerHits.get();
0056     // Create output collections
0057     auto* proto_tracks = m_outputProtoTracks.createAndPut();
0058     int n_proto_tracks = 0;
0059
0060     std::vector<ConformalHit> conformal_hits;
0061
0062     ConformalHit ref_hit(0.0,0.0); // future versions will improve on this.
0063
0064     TH1F h_phi("h_phi",";phi",m_nPhiBins.value(),-M_PI,M_PI);
0065
0066     // 1. conformal XY transform hits
0067     // 2. fill histogram with phi
0068     for(const auto& ahit : *hits) {
0069       double xc = ahit.getPosition().x - ref_hit.x();
0070       double yc = ahit.getPosition().y - ref_hit.y();
0071       double r = std::hypot(xc, yc);
0072       conformal_hits.emplace_back(2.0*xc/r,2.0*yc/r);
0073       h_phi.Fill(conformal_hits.back().phi());
0074     }
0075     // 3. Get location of maxima
0076     std::vector<int> max_bins;
0077     while(max_bins.size() < 100) {
0078       int    max_bin = h_phi.GetMaximumBin();
0079       double max_val = h_phi.GetMaximum();
0080       if(max_val < 3)  {
0081         break;
0082       }
0083       max_bins.push_back(max_bin);
0084       h_phi.SetBinContent(max_bin, 0.0); // zero bin and continue
0085     }
0086     n_proto_tracks = max_bins.size();
0087     if (msgLevel(MSG::DEBUG)) {
0088       debug() << " Found " << n_proto_tracks << " proto tracks." << endmsg;
0089     }
0090     // 4. Group hits peaked in phi
0091     for(auto b : max_bins) {
0092       ActsExamples::ProtoTrack proto_track; // this is just a std::vector<int>
0093       for(size_t ihit = 0 ; ihit< hits->size() ; ihit++) {
0094         double phi = conformal_hits[ihit].phi();
0095         double bin_phi = h_phi.GetXaxis()->GetBinCenter(b);
0096         double bin_width = h_phi.GetXaxis()->GetBinWidth(b); /// \todo make bin width an algo parameter
0097         if (std::abs(phi - bin_phi) < bin_width/2.0) {
0098           proto_track.push_back(ihit);
0099         }
0100       }
0101       if(proto_track.size()>3) {
0102         proto_tracks->push_back(proto_track);
0103       }
0104     }
0105     m_nProtoTracks.put(new int{n_proto_tracks});
0106     // 5. profit!
0107
0108     return StatusCode::SUCCESS;
0109   }
0110 };
0111 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0112 DECLARE_COMPONENT(ConformalXYPeakProtoTracks)
0113
0114 } // namespace Jug::Reco