algorithms/fardetectors/MatrixTransferStatic.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 - 2025 Alex Jentsch, Wouter Deconinck, Sylvester Joosten, David Lawrence, Simon Gardner
0003 //
0004 // This converted from: https://eicweb.phy.anl.gov/EIC/juggler/-/blob/master/JugReco/src/components/FarForwardParticles.cpp
0005
0006 #include "MatrixTransferStatic.h"
0007
0008 #include <DD4hep/Alignments.h>
0009 #include <DD4hep/DetElement.h>
0010 #include <DD4hep/Objects.h>
0011 #include <DD4hep/VolumeManager.h>
0012 #include <Evaluator/DD4hepUnits.h>
0013 #include <Math/GenVector/Cartesian3D.h>
0014 #include <Math/GenVector/DisplacementVector3D.h>
0015 #include <edm4hep/Vector3d.h>
0016 #include <edm4hep/Vector3f.h>
0017 #include <edm4hep/utils/vector_utils.h>
0018 #include <fmt/core.h>
0019 #include <cmath>
0020
0021 #include <gsl/pointers>
0022 #include <stdexcept>
0023 #include <vector>
0024
0025 #include "algorithms/fardetectors/MatrixTransferStaticConfig.h"
0026
0027 void eicrecon::MatrixTransferStatic::init() {}
0028
0029 void eicrecon::MatrixTransferStatic::process(const MatrixTransferStatic::Input& input,
0030                                              const MatrixTransferStatic::Output& output) const {
0031
0032   const auto [mcparts, rechits] = input;
0033   auto [outputParticles]        = output;
0034
0035   std::vector<std::vector<double>> aX;
0036   std::vector<std::vector<double>> aY;
0037
0038   //----- Define constants here ------
0039   double aXinv[2][2] = {{0.0, 0.0}, {0.0, 0.0}};
0040   double aYinv[2][2] = {{0.0, 0.0}, {0.0, 0.0}};
0041
0042   double nomMomentum = NAN;
0043   double local_x_offset{0.0};
0044   double local_y_offset{0.0};
0045   double local_x_slope_offset{0.0};
0046   double local_y_slope_offset{0.0};
0047
0048   double numBeamProtons  = 0;
0049   double runningMomentum = 0.0;
0050
0051   for (const auto& p : *mcparts) {
0052     if (mcparts->size() == 1 && p.getPDG() == 2212) { //proton particle gun for testing
0053       runningMomentum = p.getMomentum().z;
0054       numBeamProtons++;
0055     }
0056     if (p.getGeneratorStatus() == 4 && p.getPDG() == 2212) { //look for "beam" proton
0057       runningMomentum += p.getMomentum().z;
0058       numBeamProtons++;
0059     }
0060     if (p.getGeneratorStatus() == 4 &&
0061         p.getPDG() == 2112) { //look for "beam" neutron (for deuterons)
0062       runningMomentum += p.getMomentum().z;
0063       numBeamProtons++;
0064     }
0065   }
0066
0067   if (numBeamProtons == 0) {
0068     if (m_cfg.requireBeamProton) {
0069       critical("No beam protons found");
0070       throw std::runtime_error("No beam protons found");
0071     }
0072     return;
0073   }
0074
0075   nomMomentum = runningMomentum / numBeamProtons;
0076
0077   double nomMomentumError = 0.05;
0078
0079   //This is a temporary solution to get the beam energy information
0080   //needed to select the correct matrix
0081
0082   bool matrix_found = false;
0083   for (const MatrixConfig& matrix_config : m_cfg.matrix_configs) {
0084     if (std::abs(matrix_config.nomMomentum - nomMomentum) / matrix_config.nomMomentum <
0085         nomMomentumError) {
0086       if (matrix_found) {
0087         error("Conflicting matrix values matching momentum {}", nomMomentum);
0088       }
0089       matrix_found = true;
0090
0091       aX = matrix_config.aX;
0092       aY = matrix_config.aY;
0093
0094       local_x_offset       = matrix_config.local_x_offset;
0095       local_y_offset       = matrix_config.local_y_offset;
0096       local_x_slope_offset = matrix_config.local_x_slope_offset;
0097       local_y_slope_offset = matrix_config.local_y_slope_offset;
0098     }
0099   }
0100   if (not matrix_found) {
0101     if (m_cfg.requireMatchingMatrix) {
0102       critical("No matrix found with matching beam momentum");
0103       throw std::runtime_error("No matrix found with matching beam momentum");
0104     }
0105     return;
0106   }
0107
0108   double determinant = aX[0][0] * aX[1][1] - aX[0][1] * aX[1][0];
0109
0110   if (determinant == 0) {
0111     error("Reco matrix determinant = 0! Matrix cannot be inverted! Double-check matrix!");
0112     return;
0113   }
0114
0115   aXinv[0][0] = aX[1][1] / determinant;
0116   aXinv[0][1] = -aX[0][1] / determinant;
0117   aXinv[1][0] = -aX[1][0] / determinant;
0118   aXinv[1][1] = aX[0][0] / determinant;
0119
0120   determinant = aY[0][0] * aY[1][1] - aY[0][1] * aY[1][0];
0121
0122   if (determinant == 0) {
0123     error("Reco matrix determinant = 0! Matrix cannot be inverted! Double-check matrix!");
0124     return;
0125   }
0126
0127   aYinv[0][0] = aY[1][1] / determinant;
0128   aYinv[0][1] = -aY[0][1] / determinant;
0129   aYinv[1][0] = -aY[1][0] / determinant;
0130   aYinv[1][1] = aY[0][0] / determinant;
0131
0132   //---- begin Reconstruction code ----
0133
0134   edm4hep::Vector3f goodHit[2] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
0135
0136   bool goodHit1 = false;
0137   bool goodHit2 = false;
0138
0139   int numGoodHits1 = 0;
0140   int numGoodHits2 = 0;
0141
0142   trace("size of RP hit array = {}", rechits->size());
0143
0144   for (const auto& h : *rechits) {
0145
0146     auto cellID = h.getCellID();
0147     // The actual hit position in Global Coordinates
0148     auto gpos = m_converter->position(cellID);
0149     // local positions
0150     auto volman = m_detector->volumeManager();
0151     auto local  = volman.lookupDetElement(cellID);
0152
0153     auto pos0 = local.nominal().worldToLocal(
0154         dd4hep::Position(gpos.x(), gpos.y(), gpos.z())); // hit position in local coordinates
0155
0156     // convert into mm
0157     gpos = gpos / dd4hep::mm;
0158     pos0 = pos0 / dd4hep::mm;
0159
0160     trace("gpos.z() = {}, pos0.z() = {}, E_dep = {}", gpos.z(), pos0.z(), h.getEdep());
0161
0162     if (gpos.z() > m_cfg.hit2minZ && gpos.z() < m_cfg.hit2maxZ) {
0163
0164       trace("[gpos.x(), gpos.y(), gpos.z()] = {}, {}, {};  E_dep = {} MeV", gpos.x(), gpos.y(),
0165             gpos.z(), h.getEdep() * 1000);
0166       numGoodHits2++;
0167       goodHit[1].x = gpos.x(); //pos0.x() - pos0 is local coordinates, gpos is global
0168       goodHit[1].y =
0169           gpos.y(); //pos0.y() - temporarily changing to global to solve the local coordinate issue
0170       goodHit[1].z = gpos.z(); //         - which is unique to the Roman pots situation
0171       goodHit2     = numGoodHits2 == 1;
0172     }
0173     if (gpos.z() > m_cfg.hit1minZ && gpos.z() < m_cfg.hit1maxZ) {
0174
0175       trace("[gpos.x(), gpos.y(), gpos.z()] = {}, {}, {};  E_dep = {} MeV", gpos.x(), gpos.y(),
0176             gpos.z(), h.getEdep() * 1000);
0177       numGoodHits1++;
0178       goodHit[0].x = gpos.x(); //pos0.x()
0179       goodHit[0].y = gpos.y(); //pos0.y()
0180       goodHit[0].z = gpos.z();
0181       goodHit1     = numGoodHits1 == 1;
0182     }
0183   }
0184
0185   // NB:
0186   // This is a "dumb" algorithm - I am just checking the basic thing works with a simple single-proton test.
0187   // Will need to update and modify for generic # of hits for more complicated final-states.
0188
0189   if (goodHit1 && goodHit2) {
0190
0191     // extract hit, subtract orbit offset – this is to get the hits in the coordinate system of the orbit
0192     // trajectory
0193     double XL[2] = {goodHit[0].x - local_x_offset, goodHit[1].x - local_x_offset};
0194     double YL[2] = {goodHit[0].y - local_y_offset, goodHit[1].y - local_y_offset};
0195
0196     double base = goodHit[1].z - goodHit[0].z;
0197
0198     if (base == 0) {
0199       info("Detector separation = 0! Cannot calculate slope!");
0200     } else {
0201
0202       double Xip[2] = {0.0, 0.0};
0203       double Xrp[2] = {XL[1], ((XL[1] - XL[0]) / (base)) / dd4hep::mrad - local_x_slope_offset};
0204       double Yip[2] = {0.0, 0.0};
0205       double Yrp[2] = {YL[1], ((YL[1] - YL[0]) / (base)) / dd4hep::mrad - local_y_slope_offset};
0206
0207       // use the hit information and calculated slope at the RP + the transfer matrix inverse to calculate the
0208       // Polar Angle and deltaP at the IP
0209
0210       for (unsigned i0 = 0; i0 < 2; i0++) {
0211         for (unsigned i1 = 0; i1 < 2; i1++) {
0212           Xip[i0] += aXinv[i0][i1] * Xrp[i1];
0213           Yip[i0] += aYinv[i0][i1] * Yrp[i1];
0214         }
0215       }
0216
0217       Yip[1] = Yrp[0] / aY[0][1];
0218
0219       // convert polar angles to radians
0220       double rsx = Xip[1] * dd4hep::mrad;
0221       double rsy = Yip[1] * dd4hep::mrad;
0222
0223       // calculate momentum magnitude from measured deltaP – using thin lens optics.
0224       double p    = nomMomentum * (1 + 0.01 * Xip[0]);
0225       double norm = std::sqrt(1.0 + rsx * rsx + rsy * rsy);
0226
0227       edm4hep::Vector3f prec = {static_cast<float>(p * rsx / norm),
0228                                 static_cast<float>(p * rsy / norm), static_cast<float>(p / norm)};
0229       auto refPoint          = goodHit[0];
0230
0231       trace("RP Reco Momentum ---> px = {},  py = {}, pz = {}", prec.x, prec.y, prec.z);
0232
0233       //----- end reconstruction code ------
0234
0235       edm4eic::MutableReconstructedParticle reconTrack;
0236       reconTrack.setType(0);
0237       reconTrack.setMomentum(prec);
0238       reconTrack.setEnergy(
0239           std::hypot(edm4hep::utils::magnitude(reconTrack.getMomentum()), m_cfg.partMass));
0240       reconTrack.setReferencePoint(refPoint);
0241       reconTrack.setCharge(m_cfg.partCharge);
0242       reconTrack.setMass(m_cfg.partMass);
0243       reconTrack.setGoodnessOfPID(1.);
0244       reconTrack.setPDG(m_cfg.partPDG);
0245       //reconTrack.covMatrix(); // @TODO: Errors
0246       outputParticles->push_back(reconTrack);
0247     }
0248   } // end enough hits for matrix reco
0249 }