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