algorithms/fardetectors/PolynomialMatrixReconstruction.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2025 Alex Jentsch
0003 //
0004
0005 #include "algorithms/fardetectors/PolynomialMatrixReconstruction.h"
0006
0007 #include <DD4hep/Alignments.h>
0008 #include <DD4hep/DetElement.h>
0009 #include <DD4hep/Objects.h>
0010 #include <DD4hep/VolumeManager.h>
0011 #include <Evaluator/DD4hepUnits.h>
0012 #include <Math/GenVector/Cartesian3D.h>
0013 #include <Math/GenVector/DisplacementVector3D.h>
0014 #include <TGraph2D.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 #include <filesystem>
0021 #include <gsl/pointers>
0022 #include <memory>
0023 #include <stdexcept>
0024 #include <vector>
0025
0026 #include "algorithms/fardetectors/PolynomialMatrixReconstructionConfig.h"
0027
0028 namespace eicrecon {
0029
0030 void eicrecon::PolynomialMatrixReconstruction::init() {}
0031
0032 void eicrecon::PolynomialMatrixReconstruction::process(
0033     const PolynomialMatrixReconstruction::Input& input,
0034     const PolynomialMatrixReconstruction::Output& output) const {
0035
0036   const auto [mcparts, rechits] = input;
0037   auto [outputParticles]        = output;
0038
0039   std::vector<std::vector<double>> aX;
0040   std::vector<std::vector<double>> aY;
0041
0042   //----- Define constants here ------
0043   double aXRP[2][2] = {{0.0, 0.0}, {0.0, 0.0}};
0044   double aYRP[2][2] = {{0.0, 0.0}, {0.0, 0.0}};
0045
0046   double aXInv[2][2] = {{0.0, 0.0}, {0.0, 0.0}};
0047   double aYInv[2][2] = {{0.0, 0.0}, {0.0, 0.0}};
0048
0049   double nomMomentum          = NAN;
0050   double local_x_offset       = 0.0;
0051   double local_y_offset       = 0.0;
0052   double local_x_slope_offset = 0.0;
0053   double local_y_slope_offset = 0.0;
0054
0055   double numBeamProtons  = 0;
0056   double runningMomentum = 0.0;
0057
0058   for (const auto& p : *mcparts) {
0059     if (mcparts->size() == 1 && p.getPDG() == 2212) { //proton particle gun for testing
0060       runningMomentum = p.getMomentum().z;
0061       numBeamProtons++;
0062     }
0063     if (p.getGeneratorStatus() == 4 && p.getPDG() == 2212) { //look for "beam" proton
0064       runningMomentum += p.getMomentum().z;
0065       numBeamProtons++;
0066     }
0067     if (p.getGeneratorStatus() == 4 &&
0068         p.getPDG() == 2112) { //look for "beam" neutron (for deuterons)
0069       runningMomentum += p.getMomentum().z;
0070       numBeamProtons++;
0071     }
0072     if (p.getGeneratorStatus() == 1 && p.getPDG() == 2212) {
0073       double xMom_noXAngle = sin(m_cfg.crossingAngle) * p.getMomentum().z +
0074                              cos(m_cfg.crossingAngle) * p.getMomentum().x;
0075
0076       trace("Final-state proton MCParticle momentum --> px = {}, py = {}, pz = {}, theta_x_IP = "
0077             "{}, theta_y_IP = {}",
0078             xMom_noXAngle, p.getMomentum().y, p.getMomentum().z, xMom_noXAngle / p.getMomentum().z,
0079             p.getMomentum().y / p.getMomentum().z);
0080     }
0081   }
0082
0083   if (numBeamProtons == 0) {
0084     if (m_cfg.requireBeamProton) {
0085       error("No beam protons found");
0086       throw std::runtime_error("No beam protons found");
0087     }
0088     return;
0089   }
0090
0091   nomMomentum = runningMomentum / numBeamProtons;
0092
0093   trace("nomMomentum extracted from beam protons --> nomMomentum = {}", nomMomentum);
0094
0095   double nomMomentumError = 0.05;
0096
0097   //method for polynomial method of reconstruction
0098   //
0099   // 1) use the hit information to find the correct X_L value needed to evaluate the matrix
0100   // 2) evaluate the correct matrix for x and y reconstruction
0101   // 3) invert the matrices and extract the needed deltaP/p and theta_x,y values
0102   // 4) use 3) to calculate the momentum vector
0103   // 5) save output
0104   //
0105   // matrices are as follows
0106   //
0107   // | a0  a1 |   | theta_x_ip |   | x_rp       |
0108   // |        | x |            | = |            |
0109   // | b0  b1 |   | deltap/p   |   | theta_x_rp |
0110   //
0111   // | c0  c1 |   | y_ip       |   | y_rp       |
0112   // |        | x |            | = |            |
0113   // | d0  d1 |   | theta_y_ip |   | theta_y_rp |
0114   //
0115   // the "y_ip" part is a bit extraneous is more for some level of testing - we generally "assume" the vectors come from
0116   // (x, y, z) = (0, 0, 0), because we don't really have a choice.
0117   //
0118   // We will test how vertex smearing contributes to overall worsening of the resolution with the generally "correct" matrix being used
0119   // allowing us to isolate the various smearing effects.
0120   //
0121
0122   bool valid_energy_found = false;
0123   for (const PolynomialMatrixConfig& poly_matrix_configs : m_cfg.poly_matrix_configs) {
0124     if (std::abs(poly_matrix_configs.nomMomentum - nomMomentum) / poly_matrix_configs.nomMomentum <
0125         nomMomentumError) {
0126       if (valid_energy_found) {
0127         error("Conflicting matrix values matching momentum {}", nomMomentum);
0128       }
0129       valid_energy_found = true;
0130       nomMomentum        = poly_matrix_configs.nomMomentum;
0131     }
0132   }
0133   if (not valid_energy_found) {
0134     if (m_cfg.requireValidBeamEnergy) {
0135       error("No tune beam energy found - cannot acquire lookup table");
0136       throw std::runtime_error("No valid beam energy found, cannot reconstruct momentum");
0137     }
0138     return;
0139   }
0140
0141   //xL table filled here from LUT -- Graph2D used for nice interpolation functionality and simple loading of LUT file
0142
0143   thread_local std::string filename(
0144       fmt::format("calibrations/RP_60_xL_100_beamEnergy_{:.0f}.xL.lut", nomMomentum));
0145   thread_local std::unique_ptr<TGraph2D> xLGraph{nullptr};
0146   if (xLGraph == nullptr) {
0147     if (std::filesystem::exists(filename)) {
0148       xLGraph = std::make_unique<TGraph2D>(filename.c_str(), "%lf %lf %lf");
0149     } else {
0150       error("Cannot find lookup xL table for {}", nomMomentum);
0151       throw std::runtime_error("Cannot find xL lookup table from calibrations -- cannot proceed");
0152     }
0153   }
0154
0155   trace("filename for lookup --> {}", filename);
0156
0157   //important to ensure interoplation works correctly -- do not remove -- not available until ROOT v6.36, will need to add back in later
0158   //xLGraph->RemoveDuplicates();
0159
0160   //---- begin Reconstruction code ----
0161
0162   edm4hep::Vector3f goodHit[2] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
0163
0164   bool goodHit1 = false;
0165   bool goodHit2 = false;
0166
0167   int numGoodHits1 = 0;
0168   int numGoodHits2 = 0;
0169
0170   trace("size of RP hit array = {}", rechits->size());
0171
0172   //
0173   //we need the hit information FIRST in order to start the reconstruction procedure
0174   //
0175
0176   for (const auto& h : *rechits) {
0177
0178     auto cellID = h.getCellID();
0179     // The actual hit position in Global Coordinates
0180     auto gpos = m_converter->position(cellID);
0181     // local positions
0182     auto volman = m_detector->volumeManager();
0183     auto local  = volman.lookupDetElement(cellID);
0184
0185     auto pos0 = local.nominal().worldToLocal(
0186         dd4hep::Position(gpos.x(), gpos.y(), gpos.z())); // hit position in local coordinates
0187
0188     // convert into mm
0189     gpos = gpos / dd4hep::mm;
0190     pos0 = pos0 / dd4hep::mm;
0191
0192     trace("gpos.z() = {}, gpos.x() = {}, gpos.y() = {}, E_dep = {}", gpos.z(), gpos.x(), gpos.y(),
0193           h.getEdep());
0194
0195     if (gpos.z() > m_cfg.hit2minZ && gpos.z() < m_cfg.hit2maxZ) {
0196
0197       trace("[gpos.x(), gpos.y(), gpos.z()] = {}, {}, {};  E_dep = {} MeV", gpos.x(), gpos.y(),
0198             gpos.z(), h.getEdep() * 1000);
0199       numGoodHits2++;
0200       goodHit[1].x = gpos.x(); //pos0.x() - pos0 is local coordinates, gpos is global
0201       goodHit[1].y =
0202           gpos.y(); //pos0.y() - temporarily changing to global to solve the local coordinate issue
0203       goodHit[1].z = gpos.z(); //         - which is unique to the Roman pots situation
0204       if (numGoodHits2 == 1) {
0205         goodHit2 = true;
0206       } else
0207         goodHit2 = false;
0208     }
0209     if (gpos.z() > m_cfg.hit1minZ && gpos.z() < m_cfg.hit1maxZ) {
0210
0211       trace("[gpos.x(), gpos.y(), gpos.z()] = {}, {}, {};  E_dep = {} MeV", gpos.x(), gpos.y(),
0212             gpos.z(), h.getEdep() * 1000);
0213       numGoodHits1++;
0214       goodHit[0].x = gpos.x(); //pos0.x()
0215       goodHit[0].y = gpos.y(); //pos0.y()
0216       goodHit[0].z = gpos.z();
0217       if (numGoodHits1 == 1) {
0218         goodHit1 = true;
0219       } else
0220         goodHit1 = false;
0221     }
0222   }
0223
0224   // This algorithm uses single hits from the first and last layer of the RP
0225   // Better to do a linear fit with all 4 layers to improve resolution and to
0226   // better handle both multi-particle final-states and cases where large numbers
0227   // of secondaries are produced from spallation on detector material
0228
0229   if (goodHit1 && goodHit2) {
0230
0231     double base = goodHit[1].z - goodHit[0].z;
0232
0233     if (base == 0) {
0234       info("Detector separation = 0! Cannot calculate slope!");
0235       //delete avoids overwriting of TGraph2D and memory leak warning
0236       // cannot load the table in init stage because the beam energy is
0237       // extracted event by event
0238       // can be fixed with beam energy meta data being added in the npsim
0239       // output tree
0240       throw std::runtime_error("Detector separation = 0! Cannot calculate slope!");
0241     }
0242
0243     double XL[2] = {goodHit[0].x, goodHit[1].x};
0244     double YL[2] = {goodHit[0].y, goodHit[1].y};
0245
0246     double Xrp[2] = {XL[1], ((XL[1] - XL[0]) / (base)) / dd4hep::mrad};
0247     double Yrp[2] = {YL[1], ((YL[1] - YL[0]) / (base)) / dd4hep::mrad};
0248
0249     //First, we use the hit information to extract the x_L value needed for the matrix.
0250     //This x_L evaluation only uses the x_rp and theta_x_rp values (2D lookup).
0251
0252     //100 converts xL from decimal to percentage -- it's how the fits were done.
0253     double extracted_xL_value = 100 * xLGraph->Interpolate(Xrp[0], Xrp[1]);
0254
0255     trace("RP extracted x_L ---> x_rp = {}, theta_x_rp = {}, x_L = {}", Xrp[0], Xrp[1],
0256           extracted_xL_value);
0257
0258     if (extracted_xL_value == 0) {
0259       error("Extracted X_L value is 0 --> cannot calculate matrix");
0260       throw std::runtime_error("Extracted X_L value is 0 --> cannot calculate matrix");
0261     }
0262
0263     local_x_offset       = calculateOffsetFromXL(0, extracted_xL_value, nomMomentum);
0264     local_x_slope_offset = calculateOffsetFromXL(1, extracted_xL_value, nomMomentum);
0265     local_y_offset       = calculateOffsetFromXL(2, extracted_xL_value, nomMomentum);
0266     local_y_slope_offset = calculateOffsetFromXL(3, extracted_xL_value, nomMomentum);
0267
0268     trace("RP offsets ---> local_x_offset = {},  local_x_slope_offset = {}, local_y_offset = {}, "
0269           "local_y_slope_offset = {}",
0270           local_x_offset, local_x_slope_offset, local_y_offset, local_y_slope_offset);
0271
0272     // extract hit, subtract orbit offset – this is to get the hits in the coordinate system of the orbit
0273     // trajectory --> then we can get the matrix parameters
0274
0275     Xrp[0] = Xrp[0] - local_x_offset;
0276     Xrp[1] = Xrp[1] - local_x_slope_offset;
0277     Yrp[0] = Yrp[0] - local_y_offset;
0278     Yrp[1] = Yrp[1] - local_y_slope_offset;
0279
0280     aXRP[0][0] = calculateMatrixValueFromXL(0, extracted_xL_value, nomMomentum);
0281     aXRP[0][1] = calculateMatrixValueFromXL(1, extracted_xL_value, nomMomentum);
0282     aXRP[1][0] = calculateMatrixValueFromXL(2, extracted_xL_value, nomMomentum);
0283     aXRP[1][1] = calculateMatrixValueFromXL(3, extracted_xL_value, nomMomentum);
0284
0285     aYRP[0][0] = calculateMatrixValueFromXL(4, extracted_xL_value, nomMomentum);
0286     aYRP[0][1] = calculateMatrixValueFromXL(5, extracted_xL_value, nomMomentum);
0287     aYRP[1][0] = calculateMatrixValueFromXL(6, extracted_xL_value, nomMomentum);
0288     aYRP[1][1] = calculateMatrixValueFromXL(7, extracted_xL_value, nomMomentum);
0289
0290     trace("matrix values --> x[0][0] = {}, x[0][1] = {}, x[1][0] = {}, x[1][1] = {}", aXRP[0][0],
0291           aXRP[0][1], aXRP[1][0], aXRP[1][1]);
0292     trace("matrix values --> y[0][0] = {}, y[0][1] = {}, y[1][0] = {}, y[1][1] = {}", aYRP[0][0],
0293           aYRP[0][1], aYRP[1][0], aYRP[1][1]);
0294
0295     double determinant_x = aXRP[0][0] * aXRP[1][1] - aXRP[0][1] * aXRP[1][0];
0296     double determinant_y = aYRP[0][0] * aYRP[1][1] - aYRP[0][1] * aYRP[1][0];
0297
0298     if (determinant_x == 0 || determinant_y == 0) {
0299       error("Reco matrix determinant = 0! Matrix cannot be inverted! Double-check matrix!");
0300       throw std::runtime_error(
0301           "Reco matrix determinant = 0! Matrix cannot be inverted! Double-check matrix!");
0302     }
0303
0304     aXInv[0][0] = (aXRP[1][1] / determinant_x);
0305     aXInv[0][1] = -1 * (aXRP[0][1] / determinant_x);
0306     aXInv[1][0] = -1 * (aXRP[1][0] / determinant_x);
0307     aXInv[1][1] = (aXRP[0][0] / determinant_x);
0308
0309     aYInv[0][0] = (aYRP[1][1] / determinant_y);
0310     aYInv[0][1] = -1 * (aYRP[0][1] / determinant_y);
0311     aYInv[1][0] = -1 * (aYRP[1][0] / determinant_y);
0312     aYInv[1][1] = (aYRP[0][0] / determinant_y);
0313
0314     double thetax_ip        = aXInv[0][0] * Xrp[0] + aXInv[0][1] * Xrp[1];
0315     double delta_p_over_p   = aXInv[1][0] * Xrp[0] + aXInv[1][1] * Xrp[1];
0316     double thetay_ip        = aYInv[1][0] * Yrp[0] + aYInv[1][1] * Yrp[1];
0317     double thetay_ip_STATIC = Yrp[0] / aYRP[0][1];
0318     double thetay_AVG =
0319         0.5 * (thetay_ip +
0320                thetay_ip_STATIC); //approximation to solve issue with small numbers in y-matrix
0321
0322     trace(
0323         "theta_x_IP_proper = {}, thetay_IP_proper = {}, thetay_IP_STATIC = {}, thetay_IP_AVG = {}",
0324         thetax_ip / 1000, thetay_ip / 1000, thetay_ip_STATIC / 1000, thetay_AVG / 1000);
0325
0326     double rsx = thetax_ip * dd4hep::mrad;
0327     double rsy = thetay_AVG * dd4hep::mrad; //using approximation here for the theta_y
0328
0329     // calculate momentum magnitude from measured deltaP – using thin lens optics.
0330     double momOrbit = (extracted_xL_value / 100.0) * nomMomentum;
0331
0332     double p    = momOrbit * (1 + 0.01 * delta_p_over_p);
0333     double norm = std::sqrt(1.0 + rsx * rsx + rsy * rsy);
0334
0335     edm4hep::Vector3f prec = {static_cast<float>(p * rsx / norm),
0336                               static_cast<float>(p * rsy / norm), static_cast<float>(p / norm)};
0337     auto refPoint          = goodHit[0];
0338
0339     trace("RP polynomial Reco Momentum ---> px = {},  py = {}, pz = {}", prec.x, prec.y, prec.z);
0340
0341     //----- end reconstruction code ------
0342
0343     edm4eic::MutableReconstructedParticle reconTrack;
0344     reconTrack.setType(0);
0345     reconTrack.setMomentum(prec);
0346     reconTrack.setEnergy(
0347         std::hypot(edm4hep::utils::magnitude(reconTrack.getMomentum()), m_cfg.partMass));
0348     reconTrack.setReferencePoint(refPoint);
0349     reconTrack.setCharge(m_cfg.partCharge);
0350     reconTrack.setMass(m_cfg.partMass);
0351     reconTrack.setGoodnessOfPID(1.);
0352     reconTrack.setPDG(m_cfg.partPDG);
0353     //reconTrack.covMatrix(); // @TODO: Errors
0354     outputParticles->push_back(reconTrack);
0355
0356   } // end enough hits for matrix reco
0357
0358 } //end ::process
0359
0360 double PolynomialMatrixReconstruction::calculateOffsetFromXL(int whichOffset, double x_L,
0361                                                              double beamEnergy) const {
0362
0363   if (whichOffset >= 4) {
0364     throw std::runtime_error(fmt::format("Bad offset index {}", whichOffset));
0365   }
0366
0367   double offset_value_and_par[4][3];
0368
0369   //because of the optics configuration, these 6 parameters are zero
0370   //but this may not always be true, so we leave these parameters here
0371   //to account for x-y coupling with different optics or for the OMD
0372
0373   offset_value_and_par[2][0] = 0.0;
0374   offset_value_and_par[2][1] = 0.0;
0375   offset_value_and_par[2][2] = 0.0;
0376   offset_value_and_par[3][0] = 0.0;
0377   offset_value_and_par[3][1] = 0.0;
0378   offset_value_and_par[3][2] = 0.0;
0379
0380   if (beamEnergy == 275) { //275 GeV
0381     offset_value_and_par[0][0] = -1916.507316;
0382     offset_value_and_par[0][1] = 11.100930;
0383     offset_value_and_par[0][2] = -0.040338;
0384     offset_value_and_par[1][0] = -84.913431;
0385     offset_value_and_par[1][1] = 0.614101;
0386     offset_value_and_par[1][2] = -0.002200;
0387   } else if (beamEnergy == 100) { //100 GeV
0388     offset_value_and_par[0][0] = -1839.212116;
0389     offset_value_and_par[0][1] = 9.573747;
0390     offset_value_and_par[0][2] = -0.032770;
0391     offset_value_and_par[1][0] = -80.659006;
0392     offset_value_and_par[1][1] = 0.530685;
0393     offset_value_and_par[1][2] = -0.001790;
0394   } else if (beamEnergy == 130) { //130 GeV
0395     offset_value_and_par[0][0] = -1875.115783;
0396     offset_value_and_par[0][1] = 10.291006;
0397     offset_value_and_par[0][2] = -0.036341;
0398     offset_value_and_par[1][0] = -82.577631;
0399     offset_value_and_par[1][1] = 0.567990;
0400     offset_value_and_par[1][2] = -0.001972;
0401   } else if (beamEnergy == 41) { //41 GeV
0402     offset_value_and_par[0][0] = -1754.718344;
0403     offset_value_and_par[0][1] = 7.767054;
0404     offset_value_and_par[0][2] = -0.023105;
0405     offset_value_and_par[1][0] = -73.956525;
0406     offset_value_and_par[1][1] = 0.391292;
0407     offset_value_and_par[1][2] = -0.001063;
0408   } else
0409     throw std::runtime_error(fmt::format("Unknown beamEnergy {}", beamEnergy));
0410
0411   return (offset_value_and_par[whichOffset][0] + offset_value_and_par[whichOffset][1] * x_L +
0412           offset_value_and_par[whichOffset][2] * x_L * x_L);
0413 }
0414
0415 double PolynomialMatrixReconstruction::calculateMatrixValueFromXL(int whichElement, double x_L,
0416                                                                   double beamEnergy) const {
0417
0418   double matrix_value_and_par[8][3];
0419
0420   if ((whichElement < 0) || (whichElement > 7)) {
0421     throw std::runtime_error(fmt::format("Bad Array index(ces)", whichElement));
0422   }
0423
0424   if (beamEnergy == 275) { //275 GeV
0425     matrix_value_and_par[0][0] = -94.860289;
0426     matrix_value_and_par[0][1] = 2.373011;
0427     matrix_value_and_par[0][2] = -0.011253;
0428     matrix_value_and_par[1][0] = 4.029769;
0429     matrix_value_and_par[1][1] = 0.011753;
0430     matrix_value_and_par[1][2] = -0.000198;
0431     matrix_value_and_par[2][0] = -8.278798;
0432     matrix_value_and_par[2][1] = 0.157343;
0433     matrix_value_and_par[2][2] = -0.000728;
0434     matrix_value_and_par[3][0] = 0.217521;
0435     matrix_value_and_par[3][1] = 0.000788;
0436     matrix_value_and_par[3][2] = -0.000011;
0437     matrix_value_and_par[4][0] = -0.224777;
0438     matrix_value_and_par[4][1] = 0.003369;
0439     matrix_value_and_par[4][2] = -0.000015;
0440     matrix_value_and_par[5][0] = -162.881282;
0441     matrix_value_and_par[5][1] = 2.959217;
0442     matrix_value_and_par[5][2] = -0.013032;
0443     matrix_value_and_par[6][0] = -0.006100;
0444     matrix_value_and_par[6][1] = 0.000045;
0445     matrix_value_and_par[6][2] = 0.000000;
0446     matrix_value_and_par[7][0] = -7.781919;
0447     matrix_value_and_par[7][1] = 0.138049;
0448     matrix_value_and_par[7][2] = -0.000618;
0449   } else if (beamEnergy == 100) { //100 GeV
0450     matrix_value_and_par[0][0] = -145.825102;
0451     matrix_value_and_par[0][1] = 3.099320;
0452     matrix_value_and_par[0][2] = -0.014369;
0453     matrix_value_and_par[1][0] = 1.600321;
0454     matrix_value_and_par[1][1] = 0.055970;
0455     matrix_value_and_par[1][2] = -0.000407;
0456     matrix_value_and_par[2][0] = -10.874678;
0457     matrix_value_and_par[2][1] = 0.193778;
0458     matrix_value_and_par[2][2] = -0.000883;
0459     matrix_value_and_par[3][0] = 0.188975;
0460     matrix_value_and_par[3][1] = 0.000745;
0461     matrix_value_and_par[3][2] = -0.000008;
0462     matrix_value_and_par[4][0] = -0.341312;
0463     matrix_value_and_par[4][1] = 0.005630;
0464     matrix_value_and_par[4][2] = -0.000026;
0465     matrix_value_and_par[5][0] = -193.516983;
0466     matrix_value_and_par[5][1] = 3.532786;
0467     matrix_value_and_par[5][2] = -0.015695;
0468     matrix_value_and_par[6][0] = -0.007973;
0469     matrix_value_and_par[6][1] = 0.000064;
0470     matrix_value_and_par[6][2] = 0.000000;
0471     matrix_value_and_par[7][0] = -7.287608;
0472     matrix_value_and_par[7][1] = 0.118922;
0473     matrix_value_and_par[7][2] = -0.000468;
0474   } else if (beamEnergy == 130) { //130 GeV
0475     matrix_value_and_par[0][0] = -124.028256;
0476     matrix_value_and_par[0][1] = 2.790119;
0477     matrix_value_and_par[0][2] = -0.013056;
0478     matrix_value_and_par[1][0] = 2.985763;
0479     matrix_value_and_par[1][1] = 0.029403;
0480     matrix_value_and_par[1][2] = -0.000275;
0481     matrix_value_and_par[2][0] = -9.921622;
0482     matrix_value_and_par[2][1] = 0.181781;
0483     matrix_value_and_par[2][2] = -0.000838;
0484     matrix_value_and_par[3][0] = 0.322642;
0485     matrix_value_and_par[3][1] = -0.002093;
0486     matrix_value_and_par[3][2] = 0.000007;
0487     matrix_value_and_par[4][0] = -0.291535;
0488     matrix_value_and_par[4][1] = 0.004711;
0489     matrix_value_and_par[4][2] = -0.000022;
0490     matrix_value_and_par[5][0] = -173.506851;
0491     matrix_value_and_par[5][1] = 3.174727;
0492     matrix_value_and_par[5][2] = -0.014032;
0493     matrix_value_and_par[6][0] = -0.006733;
0494     matrix_value_and_par[6][1] = 0.000052;
0495     matrix_value_and_par[6][2] = 0.000000;
0496     matrix_value_and_par[7][0] = -7.030784;
0497     matrix_value_and_par[7][1] = 0.118149;
0498     matrix_value_and_par[7][2] = -0.000485;
0499   } else if (beamEnergy == 41) { //41 GeV
0500     matrix_value_and_par[0][0] = -174.583113;
0501     matrix_value_and_par[0][1] = 3.599373;
0502     matrix_value_and_par[0][2] = -0.016738;
0503     matrix_value_and_par[1][0] = -2.768064;
0504     matrix_value_and_par[1][1] = 0.143176;
0505     matrix_value_and_par[1][2] = -0.000847;
0506     matrix_value_and_par[2][0] = -12.409090;
0507     matrix_value_and_par[2][1] = 0.218275;
0508     matrix_value_and_par[2][2] = -0.000994;
0509     matrix_value_and_par[3][0] = -0.135417;
0510     matrix_value_and_par[3][1] = 0.006663;
0511     matrix_value_and_par[3][2] = -0.000036;
0512     matrix_value_and_par[4][0] = -0.299041;
0513     matrix_value_and_par[4][1] = 0.004490;
0514     matrix_value_and_par[4][2] = -0.000020;
0515     matrix_value_and_par[5][0] = -175.144897;
0516     matrix_value_and_par[5][1] = 3.222149;
0517     matrix_value_and_par[5][2] = -0.014292;
0518     matrix_value_and_par[6][0] = -0.007376;
0519     matrix_value_and_par[6][1] = 0.000052;
0520     matrix_value_and_par[6][2] = 0.000000;
0521     matrix_value_and_par[7][0] = -8.443702;
0522     matrix_value_and_par[7][1] = 0.155002;
0523     matrix_value_and_par[7][2] = -0.000708;
0524   } else
0525     throw std::runtime_error(fmt::format("Unknown beamEnergy {}", beamEnergy));
0526
0527   return (matrix_value_and_par[whichElement][0] + matrix_value_and_par[whichElement][1] * x_L +
0528           matrix_value_and_par[whichElement][2] * x_L * x_L);
0529 }
0530
0531 } //namespace eicrecon