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