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