src/components/PhotoRingClusters.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Chao Peng, Whitney Armstrong, Sylvester Joosten
0003
0004 /*  Clustering Algorithm for Ring Imaging Cherenkov (RICH) events
0005  *
0006  *  Author: Chao Peng (ANL)
0007  *  Date: 10/04/2020
0008  *
0009  */
0010
0011 #include <algorithm>
0012
0013 #include "Gaudi/Property.h"
0014 #include "Gaudi/Algorithm.h"
0015 #include "GaudiKernel/PhysicalConstants.h"
0016 #include "GaudiKernel/RndmGenerators.h"
0017 #include "GaudiKernel/ToolHandle.h"
0018
0019 #include "DDRec/CellIDPositionConverter.h"
0020 #include "DDRec/Surface.h"
0021 #include "DDRec/SurfaceManager.h"
0022
0023 #include <k4FWCore/DataHandle.h>
0024 #include <k4Interface/IGeoSvc.h>
0025
0026 // Event Model related classes
0027 #include "FuzzyKClusters.h"
0028 #include "edm4eic/PMTHitCollection.h"
0029 #include "edm4eic/RingImageCollection.h"
0030
0031 using namespace Gaudi::Units;
0032 using namespace Eigen;
0033
0034 namespace Jug::Reco {
0035
0036 /**  Clustering Algorithm for Ring Imaging Cherenkov (RICH) events.
0037  *
0038  * \ingroup reco
0039  */
0040 class PhotoRingClusters : public Gaudi::Algorithm {
0041 private:
0042   mutable DataHandle<edm4eic::PMTHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
0043   mutable DataHandle<edm4eic::RingImageCollection> m_outputClusterCollection{"outputClusterCollection", Gaudi::DataHandle::Writer,
0044                                                                   this};
0045   // @TODO
0046   // A more realistic way is to have tracker info as the input to determine how much clusters should be found
0047   Gaudi::Property<int> m_nRings{this, "nRings", 1};
0048   Gaudi::Property<int> m_nIters{this, "nIters", 1000};
0049   Gaudi::Property<double> m_q{this, "q", 2.0};
0050   Gaudi::Property<double> m_eps{this, "epsilon", 1e-4};
0051   Gaudi::Property<double> m_minNpe{this, "minNpe", 0.5};
0052   // Pointer to the geometry service
0053   SmartIF<IGeoSvc> m_geoSvc;
0054
0055 public:
0056   // ill-formed: using Gaudi::Algorithm::GaudiAlgorithm;
0057   PhotoRingClusters(const std::string& name, ISvcLocator* svcLoc) : Gaudi::Algorithm(name, svcLoc) {
0058     declareProperty("inputHitCollection", m_inputHitCollection, "");
0059     declareProperty("outputClusterCollection", m_outputClusterCollection, "");
0060   }
0061
0062   StatusCode initialize() override {
0063     if (Gaudi::Algorithm::initialize().isFailure()) {
0064       return StatusCode::FAILURE;
0065     }
0066     m_geoSvc = service("GeoSvc");
0067     if (!m_geoSvc) {
0068       error() << "Unable to locate Geometry Service. "
0069               << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
0070       return StatusCode::FAILURE;
0071     }
0072     return StatusCode::SUCCESS;
0073   }
0074
0075   StatusCode execute(const EventContext&) const override {
0076     // input collections
0077     const auto& rawhits = *m_inputHitCollection.get();
0078     // Create output collections
0079     auto& clusters = *m_outputClusterCollection.createAndPut();
0080
0081     // algorithm
0082     auto alg = fkc::KRings();
0083
0084     // fill data
0085     MatrixXd data(rawhits.size(), 2);
0086     for (int i = 0; i < data.rows(); ++i) {
0087       if (rawhits[i].getNpe() > m_minNpe) {
0088         data.row(i) << rawhits[i].getLocal().x, rawhits[i].getLocal().y;
0089       }
0090     }
0091
0092     // clustering
0093     auto res = alg.Fit(data, m_nRings, m_q, m_eps, m_nIters);
0094
0095     // local position
0096     // @TODO: Many fields in RingImage not filled, need to assess
0097     //        if those are in fact needed
0098     for (int i = 0; i < res.rows(); ++i) {
0099       auto cl = clusters.create();
0100       cl.setPosition({static_cast<float>(res(i, 0)), static_cast<float>(res(i, 1)), 0});
0101       // @TODO: positionError() not set
0102       // @TODO: theta() not set
0103       // @TODO: thetaError() not set
0104       cl.setRadius(res(i, 2));
0105       // @TODO: radiusError not set
0106     }
0107
0108     return StatusCode::SUCCESS;
0109   }
0110 };
0111
0112 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0113 DECLARE_COMPONENT(PhotoRingClusters)
0114
0115 } // namespace Jug::Reco