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