src/components/FarForwardParticlesOMD.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Alex Jentsch, Sylvester Joosten, Wouter Deconinck
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 <k4FWCore/DataHandle.h>
0015
0016 // Event Model related classes
0017 #include "edm4eic/ReconstructedParticleCollection.h"
0018 #include "edm4eic/TrackerHitCollection.h"
0019 #include <edm4hep/utils/vector_utils.h>
0020
0021 namespace Jug::Reco {
0022
0023 class FarForwardParticlesOMD : public GaudiAlgorithm {
0024 private:
0025   DataHandle<edm4eic::TrackerHitCollection> m_inputHitCollection{"FarForwardTrackerHits", Gaudi::DataHandle::Reader, this};
0026   DataHandle<edm4eic::ReconstructedParticleCollection> m_outputParticles{"outputParticles", Gaudi::DataHandle::Writer,
0027                                                                      this};
0028
0029   //----- Define constants here ------
0030
0031   Gaudi::Property<double> local_x_offset_station_1{this, "localXOffsetSta1", -762.5006104};
0032   Gaudi::Property<double> local_x_offset_station_2{this, "localXOffsetSta2", -881.9621277};
0033   Gaudi::Property<double> local_x_slope_offset{this, "localXSlopeOffset", -59.73075865};
0034   Gaudi::Property<double> local_y_slope_offset{this, "localYSlopeOffset", 0.0012755};
0035   Gaudi::Property<double> crossingAngle{this, "crossingAngle", -0.025};
0036   Gaudi::Property<double> nomMomentum{this, "beamMomentum", 137.5}; // This number is set to 50% maximum beam momentum
0037
0038   const double aXOMD[2][2] = {{1.6229248, 12.9519653}, {-2.86056525, 0.1830292}};
0039   const double aYOMD[2][2] = {{0.0000185, -28.599739}, {0.00000925, -2.8795791}};
0040
0041   double aXOMDinv[2][2] = {{0.0, 0.0}, {0.0, 0.0}};
0042   double aYOMDinv[2][2] = {{0.0, 0.0}, {0.0, 0.0}};
0043
0044 public:
0045   FarForwardParticlesOMD(const std::string& name, ISvcLocator* svcLoc)
0046       : GaudiAlgorithm(name, svcLoc) {
0047     declareProperty("inputCollection", m_inputHitCollection, "FarForwardTrackerHits");
0048     declareProperty("outputCollection", m_outputParticles, "ReconstructedParticles");
0049   }
0050
0051   StatusCode initialize() override {
0052     if (GaudiAlgorithm::initialize().isFailure()) {
0053       return StatusCode::FAILURE;
0054     }
0055     double det = aXOMD[0][0] * aXOMD[1][1] - aXOMD[0][1] * aXOMD[1][0];
0056
0057     if (det == 0) {
0058       error() << "Reco matrix determinant = 0!"
0059               << "Matrix cannot be inverted! Double-check matrix!" << endmsg;
0060       return StatusCode::FAILURE;
0061     }
0062
0063     aXOMDinv[0][0] = aXOMD[1][1] / det;
0064     aXOMDinv[0][1] = -aXOMD[0][1] / det;
0065     aXOMDinv[1][0] = -aXOMD[1][0] / det;
0066     aXOMDinv[1][1] = aXOMD[0][0] / det;
0067
0068     det            = aYOMD[0][0] * aYOMD[1][1] - aYOMD[0][1] * aYOMD[1][0];
0069     aYOMDinv[0][0] = aYOMD[1][1] / det;
0070     aYOMDinv[0][1] = -aYOMD[0][1] / det;
0071     aYOMDinv[1][0] = -aYOMD[1][0] / det;
0072     aYOMDinv[1][1] = aYOMD[0][0] / det;
0073
0074     return StatusCode::SUCCESS;
0075   }
0076
0077   StatusCode execute() override {
0078     const edm4eic::TrackerHitCollection* rawhits = m_inputHitCollection.get();
0079     auto& rc                                 = *(m_outputParticles.createAndPut());
0080
0081     // for (const auto& part : mc) {
0082     //    if (part.genStatus() > 1) {
0083     //        if (msgLevel(MSG::DEBUG)) {
0084     //            debug() << "ignoring particle with genStatus = " << part.genStatus() << endmsg;
0085     //        }
0086     //        continue;
0087     //    }
0088
0089     //---- begin Roman Pot Reconstruction code ----
0090
0091     int eventReset = 0; // counter for IDing at least one hit per layer
0092     std::vector<double> hitx;
0093     std::vector<double> hity;
0094     std::vector<double> hitz;
0095
0096     for (const auto& h : *rawhits) {
0097
0098       // The actual hit position:
0099       auto pos0 = h.getPosition();
0100       // auto mom0 = h.momentum;
0101       // auto pidCode = h.g4ID;
0102       auto eDep = h.getEdep();
0103
0104       if (eDep < 0.00001) {
0105         continue;
0106       }
0107
0108       if (eventReset < 2) {
0109         hitx.push_back(pos0.x - local_x_offset_station_2);
0110       } // use station 2 for both offsets since it is used for the reference orbit
0111       else {
0112         hitx.push_back(pos0.x - local_x_offset_station_2);
0113       }
0114
0115       hity.push_back(pos0.y);
0116       hitz.push_back(pos0.z);
0117
0118       eventReset++;
0119     }
0120
0121     // NB:
0122     // This is a "dumb" algorithm - I am just checking the basic thing works with a simple single-proton test.
0123     // Will need to update and modify for generic # of hits for more complicated final-states.
0124
0125     if (eventReset == 4) {
0126
0127       // extract hit, subtract orbit offset – this is to get the hits in the coordinate system of the orbit trajectory
0128       double XL[2] = {hitx[0], hitx[2]};
0129       double YL[2] = {hity[0], hity[2]};
0130
0131       double base = hitz[2] - hitz[0];
0132
0133       if (base == 0) {
0134         warning() << "Detector separation = 0!"
0135                 << "Cannot calculate slope!" << endmsg;
0136         return StatusCode::SUCCESS;
0137       }
0138
0139       double Xomd[2] = {XL[1], (1000 * (XL[1] - XL[0]) / (base)) - local_x_slope_offset};
0140       double Xip[2]  = {0.0, 0.0};
0141       double Yomd[2] = {YL[1], (1000 * (YL[1] - YL[0]) / (base)) - local_y_slope_offset};
0142       double Yip[2]  = {0.0, 0.0};
0143
0144       // use the hit information and calculated slope at the RP + the transfer matrix inverse to calculate the Polar
0145       // Angle and deltaP at the IP
0146
0147       for (unsigned i0 = 0; i0 < 2; i0++) {
0148         for (unsigned i1 = 0; i1 < 2; i1++) {
0149           Xip[i0] += aXOMDinv[i0][i1] * Xomd[i1];
0150           Yip[i0] += aYOMDinv[i0][i1] * Yomd[i1];
0151         }
0152       }
0153
0154       // convert polar angles to radians
0155       double rsx = Xip[1] / 1000.;
0156       double rsy = Yip[1] / 1000.;
0157
0158       // calculate momentum magnitude from measured deltaP – using thin lens optics.
0159       double p    = nomMomentum * (1 + 0.01 * Xip[0]);
0160       double norm = std::sqrt(1.0 + rsx * rsx + rsy * rsy);
0161
0162       const float prec[3] = {static_cast<float>(p * rsx / norm), static_cast<float>(p * rsy / norm),
0163                              static_cast<float>(p / norm)};
0164
0165       //----- end RP reconstruction code ------
0166
0167       edm4eic::MutableReconstructedParticle rpTrack;
0168       rpTrack.setType(0);
0169       rpTrack.setMomentum({prec});
0170       rpTrack.setEnergy(std::hypot(edm4hep::utils::magnitude(rpTrack.getMomentum()), .938272));
0171       rpTrack.setReferencePoint({0, 0, 0});
0172       rpTrack.setCharge(1);
0173       rpTrack.setMass(.938272);
0174       rpTrack.setGoodnessOfPID(1.);
0175       rpTrack.setPDG(2122);
0176       //rpTrack.covMatrix(); // @TODO: Errors
0177       rc->push_back(rpTrack);
0178
0179     } // end enough hits for matrix reco
0180
0181     return StatusCode::SUCCESS;
0182   }
0183 };
0184
0185 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0186 DECLARE_COMPONENT(FarForwardParticlesOMD)
0187
0188 } // namespace Jug::Reco