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