src/components/FarForwardParticles.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Alex Jentsch, Wouter Deconinck, Sylvester Joosten
0003
0004 #include <algorithm>
0005 #include <cmath>
0006 #include <fmt/format.h>
0007
0008 #include "Gaudi/Algorithm.h"
0009 #include "GaudiKernel/RndmGenerators.h"
0010
0011 #include "DDRec/CellIDPositionConverter.h"
0012 #include "DDRec/Surface.h"
0013 #include "DDRec/SurfaceManager.h"
0014
0015 #include <k4FWCore/DataHandle.h>
0016 #include <k4Interface/IGeoSvc.h>
0017
0018 // Event Model related classes
0019 #include "edm4eic/ReconstructedParticleCollection.h"
0020 #include "edm4eic/TrackerHitCollection.h"
0021 #include <edm4hep/utils/vector_utils.h>
0022
0023 namespace Jug::Reco {
0024
0025 class FarForwardParticles : public Gaudi::Algorithm {
0026 private:
0027   mutable DataHandle<edm4eic::TrackerHitCollection> m_inputHitCollection{"FarForwardTrackerHits", Gaudi::DataHandle::Reader, this};
0028   mutable DataHandle<edm4eic::ReconstructedParticleCollection> m_outputParticles{"outputParticles", Gaudi::DataHandle::Writer,
0029                                                                      this};
0030
0031   //----- Define constants here ------
0032
0033   Gaudi::Property<double> local_x_offset_station_1{this, "localXOffsetSta1", -833.3878326};
0034   Gaudi::Property<double> local_x_offset_station_2{this, "localXOffsetSta2", -924.342804};
0035   Gaudi::Property<double> local_x_slope_offset{this, "localXSlopeOffset", -0.00622147};
0036   Gaudi::Property<double> local_y_slope_offset{this, "localYSlopeOffset", -0.0451035};
0037   Gaudi::Property<double> crossingAngle{this, "crossingAngle", -0.025};
0038   Gaudi::Property<double> nomMomentum{this, "beamMomentum", 275.0};
0039
0040   Gaudi::Property<std::string> m_geoSvcName{this, "geoServiceName", "GeoSvc"};
0041   Gaudi::Property<std::string> m_readout{this, "readoutClass", ""};
0042   Gaudi::Property<std::string> m_layerField{this, "layerField", ""};
0043   Gaudi::Property<std::string> m_sectorField{this, "sectorField", ""};
0044   SmartIF<IGeoSvc> m_geoSvc;
0045   std::shared_ptr<const dd4hep::rec::CellIDPositionConverter> m_converter;
0046   dd4hep::BitFieldCoder* id_dec = nullptr;
0047   size_t sector_idx{0}, layer_idx{0};
0048
0049   Gaudi::Property<std::string> m_localDetElement{this, "localDetElement", ""};
0050   Gaudi::Property<std::vector<std::string>> u_localDetFields{this, "localDetFields", {}};
0051   mutable dd4hep::DetElement local;
0052   size_t local_mask = ~0;
0053
0054   const double aXRP[2][2] = {{2.102403743, 29.11067626}, {0.186640381, 0.192604619}};
0055   const double aYRP[2][2] = {{0.0000159900, 3.94082098}, {0.0000079946, -0.1402995}};
0056
0057   double aXRPinv[2][2] = {{0.0, 0.0}, {0.0, 0.0}};
0058   double aYRPinv[2][2] = {{0.0, 0.0}, {0.0, 0.0}};
0059
0060 public:
0061   FarForwardParticles(const std::string& name, ISvcLocator* svcLoc) : Gaudi::Algorithm(name, svcLoc) {
0062     declareProperty("inputCollection", m_inputHitCollection, "FarForwardTrackerHits");
0063     declareProperty("outputCollection", m_outputParticles, "ReconstructedParticles");
0064   }
0065
0066   // See Wouter's example to extract local coordinates CalorimeterHitReco.cpp
0067   // includes DDRec/CellIDPositionConverter.here
0068   // See tutorial
0069   // auto converter = m_GeoSvc ....
0070   // https://eicweb.phy.anl.gov/EIC/juggler/-/blob/master/JugReco/src/components/CalorimeterHitReco.cpp - line 200
0071   // include the Eigen libraries, used in ACTS, for the linear algebra.
0072
0073   StatusCode initialize() override {
0074     if (Gaudi::Algorithm::initialize().isFailure()) {
0075       return StatusCode::FAILURE;
0076     }
0077     m_geoSvc = service(m_geoSvcName);
0078     if (!m_geoSvc) {
0079       error() << "Unable to locate Geometry Service. "
0080               << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
0081       return StatusCode::FAILURE;
0082     }
0083     m_converter = std::make_shared<const dd4hep::rec::CellIDPositionConverter>(*(m_geoSvc->getDetector()));
0084
0085     // do not get the layer/sector ID if no readout class provided
0086     if (m_readout.value().empty()) {
0087       return StatusCode::SUCCESS;
0088     }
0089
0090     auto id_spec = m_geoSvc->getDetector()->readout(m_readout).idSpec();
0091     try {
0092       id_dec = id_spec.decoder();
0093       if (!m_sectorField.value().empty()) {
0094         sector_idx = id_dec->index(m_sectorField);
0095         info() << "Find sector field " << m_sectorField.value() << ", index = " << sector_idx << endmsg;
0096       }
0097       if (!m_layerField.value().empty()) {
0098         layer_idx = id_dec->index(m_layerField);
0099         info() << "Find layer field " << m_layerField.value() << ", index = " << sector_idx << endmsg;
0100       }
0101     } catch (...) {
0102       error() << "Failed to load ID decoder for " << m_readout << endmsg;
0103       return StatusCode::FAILURE;
0104     }
0105
0106     // local detector name has higher priority
0107     if (!m_localDetElement.value().empty()) {
0108       try {
0109         local = m_geoSvc->getDetector()->detector(m_localDetElement.value());
0110         info() << "Local coordinate system from DetElement " << m_localDetElement.value() << endmsg;
0111       } catch (...) {
0112         error() << "Failed to locate local coordinate system from DetElement " << m_localDetElement.value() << endmsg;
0113         return StatusCode::FAILURE;
0114       }
0115       // or get from fields
0116     } else {
0117       std::vector<std::pair<std::string, int>> fields;
0118       for (auto& f : u_localDetFields.value()) {
0119         fields.emplace_back(f, 0);
0120       }
0121       local_mask = id_spec.get_mask(fields);
0122       // use all fields if nothing provided
0123       if (fields.empty()) {
0124         local_mask = ~0;
0125       }
0126       // info() << fmt::format("Local DetElement mask {:#064b} from fields [{}]", local_mask,
0127       //                      fmt::join(fields, ", "))
0128       //        << endmsg;
0129     }
0130
0131     double det = aXRP[0][0] * aXRP[1][1] - aXRP[0][1] * aXRP[1][0];
0132
0133     if (det == 0) {
0134       error() << "Reco matrix determinant = 0!"
0135               << "Matrix cannot be inverted! Double-check matrix!" << endmsg;
0136       return StatusCode::FAILURE;
0137     }
0138
0139     aXRPinv[0][0] = aXRP[1][1] / det;
0140     aXRPinv[0][1] = -aXRP[0][1] / det;
0141     aXRPinv[1][0] = -aXRP[1][0] / det;
0142     aXRPinv[1][1] = aXRP[0][0] / det;
0143
0144     det           = aYRP[0][0] * aYRP[1][1] - aYRP[0][1] * aYRP[1][0];
0145     aYRPinv[0][0] = aYRP[1][1] / det;
0146     aYRPinv[0][1] = -aYRP[0][1] / det;
0147     aYRPinv[1][0] = -aYRP[1][0] / det;
0148     aYRPinv[1][1] = aYRP[0][0] / det;
0149
0150     return StatusCode::SUCCESS;
0151   }
0152
0153   StatusCode execute(const EventContext&) const override {
0154     const edm4eic::TrackerHitCollection* rawhits = m_inputHitCollection.get();
0155     auto& rc                                 = *(m_outputParticles.createAndPut());
0156
0157     auto converter = m_converter;
0158
0159     // for (const auto& part : mc) {
0160     //    if (part.genStatus() > 1) {
0161     //        if (msgLevel(MSG::DEBUG)) {
0162     //            debug() << "ignoring particle with genStatus = " << part.genStatus() << endmsg;
0163     //        }
0164     //        continue;
0165     //    }
0166
0167     //---- begin Roman Pot Reconstruction code ----
0168
0169     int eventReset = 0; // counter for IDing at least one hit per layer
0170     std::vector<double> hitx;
0171     std::vector<double> hity;
0172     std::vector<double> hitz;
0173
0174     for (const auto& h : *rawhits) {
0175
0176       auto cellID = h.getCellID();
0177       // The actual hit position in Global Coordinates
0178       // auto pos0 = h.position();
0179
0180       auto gpos = converter->position(cellID);
0181       // local positions
0182       if (m_localDetElement.value().empty()) {
0183         auto volman = m_geoSvc->getDetector()->volumeManager();
0184         local       = volman.lookupDetElement(cellID);
0185       }
0186       auto pos0 = local.nominal().worldToLocal(
0187           dd4hep::Position(gpos.x(), gpos.y(), gpos.z())); // hit position in local coordinates
0188
0189       // auto mom0 = h.momentum;
0190       // auto pidCode = h.g4ID;
0191       auto eDep = h.getEdep();
0192
0193       if (eDep < 0.00001) {
0194         continue;
0195       }
0196
0197       if (eventReset < 2) {
0198         hitx.push_back(pos0.x()); // - local_x_offset_station_2);
0199       }                           // use station 2 for both offsets since it is used for the reference orbit
0200       else {
0201         hitx.push_back(pos0.x()); // - local_x_offset_station_2);
0202       }
0203
0204       hity.push_back(pos0.y());
0205       hitz.push_back(pos0.z());
0206
0207       eventReset++;
0208     }
0209
0210     // NB:
0211     // This is a "dumb" algorithm - I am just checking the basic thing works with a simple single-proton test.
0212     // Will need to update and modify for generic # of hits for more complicated final-states.
0213
0214     if (eventReset == 4) {
0215
0216       // extract hit, subtract orbit offset – this is to get the hits in the coordinate system of the orbit
0217       // trajectory
0218       double XL[2] = {hitx[0], hitx[2]};
0219       double YL[2] = {hity[0], hity[2]};
0220
0221       double base = hitz[2] - hitz[0];
0222
0223       if (base == 0) {
0224         warning() << "Detector separation = 0!"
0225                   << "Cannot calculate slope!" << endmsg;
0226         return StatusCode::SUCCESS;
0227       }
0228
0229       double Xip[2] = {0.0, 0.0};
0230       double Xrp[2] = {XL[1], (1000 * (XL[1] - XL[0]) / (base)) - local_x_slope_offset}; //- _SX0RP_;
0231       double Yip[2] = {0.0, 0.0};
0232       double Yrp[2] = {YL[1], (1000 * (YL[1] - YL[0]) / (base)) - local_y_slope_offset}; //- _SY0RP_;
0233
0234       // use the hit information and calculated slope at the RP + the transfer matrix inverse to calculate the
0235       // Polar Angle and deltaP at the IP
0236
0237       for (unsigned i0 = 0; i0 < 2; i0++) {
0238         for (unsigned i1 = 0; i1 < 2; i1++) {
0239           Xip[i0] += aXRPinv[i0][i1] * Xrp[i1];
0240           Yip[i0] += aYRPinv[i0][i1] * Yrp[i1];
0241         }
0242       }
0243
0244       // convert polar angles to radians
0245       double rsx = Xip[1] / 1000.;
0246       double rsy = Yip[1] / 1000.;
0247
0248       // calculate momentum magnitude from measured deltaP – using thin lens optics.
0249       double p    = nomMomentum * (1 + 0.01 * Xip[0]);
0250       double norm = std::sqrt(1.0 + rsx * rsx + rsy * rsy);
0251
0252       const float prec[3] = {static_cast<float>(p * rsx / norm), static_cast<float>(p * rsy / norm),
0253                              static_cast<float>(p / norm)};
0254
0255       //----- end RP reconstruction code ------
0256
0257       edm4eic::MutableReconstructedParticle rpTrack;
0258       rpTrack.setType(0);
0259       rpTrack.setMomentum({prec});
0260       rpTrack.setEnergy(std::hypot(edm4hep::utils::magnitude(rpTrack.getMomentum()), .938272));
0261       rpTrack.setReferencePoint({0, 0, 0});
0262       rpTrack.setCharge(1);
0263       rpTrack.setMass(.938272);
0264       rpTrack.setGoodnessOfPID(1.);
0265       rpTrack.setPDG(2122);
0266       //rpTrack.covMatrix(); // @TODO: Errors
0267       rc->push_back(rpTrack);
0268
0269     } // end enough hits for matrix reco
0270
0271     return StatusCode::SUCCESS;
0272   }
0273 };
0274
0275 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0276 DECLARE_COMPONENT(FarForwardParticles)
0277
0278 } // namespace Jug::Reco