src/components/InclusiveKinematicsElectron.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Wouter Deconinck
0003
0004 #include "Gaudi/Algorithm.h"
0005 #include "GaudiKernel/RndmGenerators.h"
0006 #include "GaudiKernel/PhysicalConstants.h"
0007 #include <algorithm>
0008 #include <cmath>
0009
0010 #include "JugBase/IParticleSvc.h"
0011 #include <k4FWCore/DataHandle.h>
0012
0013 #include "JugBase/Utilities/Beam.h"
0014 #include "JugBase/Utilities/Boost.h"
0015
0016 #include "Math/Vector4D.h"
0017 using ROOT::Math::PxPyPzEVector;
0018
0019 // Event Model related classes
0020 #include "edm4hep/MCParticleCollection.h"
0021 #include "edm4eic/MCRecoParticleAssociationCollection.h"
0022 #include "edm4eic/ReconstructedParticleCollection.h"
0023 #include "edm4eic/InclusiveKinematicsCollection.h"
0024
0025 namespace Jug::Reco {
0026
0027 class InclusiveKinematicsElectron : public Gaudi::Algorithm {
0028 private:
0029   mutable DataHandle<edm4hep::MCParticleCollection> m_inputMCParticleCollection{
0030     "inputMCParticles",
0031     Gaudi::DataHandle::Reader,
0032     this};
0033   mutable DataHandle<edm4eic::ReconstructedParticleCollection> m_inputParticleCollection{
0034     "inputReconstructedParticles",
0035     Gaudi::DataHandle::Reader,
0036     this};
0037   mutable DataHandle<edm4eic::MCRecoParticleAssociationCollection> m_inputParticleAssociation{
0038     "inputParticleAssociations",
0039     Gaudi::DataHandle::Reader,
0040     this};
0041   mutable DataHandle<edm4eic::InclusiveKinematicsCollection> m_outputInclusiveKinematicsCollection{
0042     "outputInclusiveKinematics",
0043     Gaudi::DataHandle::Writer,
0044     this};
0045
0046   Gaudi::Property<double> m_crossingAngle{this, "crossingAngle", -0.025 * Gaudi::Units::radian};
0047
0048   SmartIF<IParticleSvc> m_pidSvc;
0049   double m_proton{0}, m_neutron{0}, m_electron{0};
0050
0051 public:
0052   InclusiveKinematicsElectron(const std::string& name, ISvcLocator* svcLoc)
0053       : Gaudi::Algorithm(name, svcLoc) {
0054     declareProperty("inputMCParticles", m_inputMCParticleCollection, "MCParticles");
0055     declareProperty("inputReconstructedParticles", m_inputParticleCollection, "ReconstructedParticles");
0056     declareProperty("inputParticleAssociations", m_inputParticleAssociation, "MCRecoParticleAssociation");
0057     declareProperty("outputInclusiveKinematics", m_outputInclusiveKinematicsCollection, "InclusiveKinematicsElectron");
0058   }
0059
0060   StatusCode initialize() override {
0061     if (Gaudi::Algorithm::initialize().isFailure())
0062       return StatusCode::FAILURE;
0063
0064     m_pidSvc = service("ParticleSvc");
0065     if (!m_pidSvc) {
0066       error() << "Unable to locate Particle Service. "
0067               << "Make sure you have ParticleSvc in the configuration."
0068               << endmsg;
0069       return StatusCode::FAILURE;
0070     }
0071     m_proton = m_pidSvc->particle(2212).mass;
0072     m_neutron = m_pidSvc->particle(2112).mass;
0073     m_electron = m_pidSvc->particle(11).mass;
0074
0075     return StatusCode::SUCCESS;
0076   }
0077
0078   StatusCode execute(const EventContext&) const override {
0079     // input collections
0080     const auto& mcparts = *(m_inputMCParticleCollection.get());
0081     const auto& rcparts = *(m_inputParticleCollection.get());
0082     const auto& rcassoc = *(m_inputParticleAssociation.get());
0083     // output collection
0084     auto& out_kinematics = *(m_outputInclusiveKinematicsCollection.createAndPut());
0085
0086     // 1. find_if
0087     //const auto mc_first_electron = std::find_if(
0088     //  mcparts.begin(),
0089     //  mcparts.end(),
0090     //  [](const auto& p){ return p.getPDG() == 11; });
0091
0092     // 2a. simple loop over iterator (post-increment)
0093     //auto mc_first_electron = mcparts.end();
0094     //for (auto p = mcparts.begin(); p != mcparts.end(); p++) {
0095     //  if (p->getPDG() == 11) {
0096     //    mc_first_electron = p;
0097     //    break;
0098     //  }
0099     //}
0100     // 2b. simple loop over iterator (pre-increment)
0101     //auto mc_first_electron = mcparts.end();
0102     //for (auto p = mcparts.begin(); p != mcparts.end(); ++p) {
0103     //  if (p->getPDG() == 11) {
0104     //    mc_first_electron = p;
0105     //    break;
0106     //  }
0107     //}
0108
0109     // 3. pre-initialized simple loop
0110     //auto mc_first_electron = mcparts.begin();
0111     //for (; mc_first_electron != mcparts.end(); ++mc_first_electron) {
0112     //  if (mc_first_electron->getPDG() == 11) {
0113     //    break;
0114     //  }
0115     //}
0116
0117     // 4a. iterator equality
0118     //if (mc_first_electron == mcparts.end()) {
0119     //  debug() << "No electron found" << endmsg;
0120     //  return StatusCode::FAILURE;
0121     //}
0122     // 4b. iterator inequality
0123     //if (!(mc_first_electron != mcparts.end())) {
0124     //  debug() << "No electron found" << endmsg;
0125     //  return StatusCode::FAILURE;
0126     //}
0127
0128     // 5. ranges and views
0129     //auto is_electron = [](const auto& p){ return p.getPDG() == 11; };
0130     //for (const auto& e: mcparts | std::views::filter(is_electron)) {
0131     //  break;
0132     //}
0133
0134     // Get incoming electron beam
0135     const auto ei_coll = Jug::Base::Beam::find_first_beam_electron(mcparts);
0136     if (ei_coll.size() == 0) {
0137       if (msgLevel(MSG::DEBUG)) {
0138         debug() << "No beam electron found" << endmsg;
0139       }
0140       return StatusCode::SUCCESS;
0141     }
0142     const PxPyPzEVector ei(
0143       Jug::Base::Beam::round_beam_four_momentum(
0144         ei_coll[0].getMomentum(),
0145         m_electron,
0146         {-5.0, -10.0, -18.0},
0147         0.0)
0148       );
0149
0150     // Get incoming hadron beam
0151     const auto pi_coll = Jug::Base::Beam::find_first_beam_hadron(mcparts);
0152     if (pi_coll.size() == 0) {
0153       if (msgLevel(MSG::DEBUG)) {
0154         debug() << "No beam hadron found" << endmsg;
0155       }
0156       return StatusCode::SUCCESS;
0157     }
0158     const PxPyPzEVector pi(
0159       Jug::Base::Beam::round_beam_four_momentum(
0160         pi_coll[0].getMomentum(),
0161         pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron,
0162         {41.0, 100.0, 275.0},
0163         m_crossingAngle)
0164       );
0165
0166     // Get first scattered electron
0167     const auto ef_coll = Jug::Base::Beam::find_first_scattered_electron(mcparts);
0168     if (ef_coll.size() == 0) {
0169       if (msgLevel(MSG::DEBUG)) {
0170         debug() << "No truth scattered electron found" << endmsg;
0171       }
0172       return StatusCode::SUCCESS;
0173     }
0174     // Associate first scattered electron with reconstructed electrons
0175     //const auto ef_assoc = std::find_if(
0176     //  rcassoc.begin(),
0177     //  rcassoc.end(),
0178     //  [&ef_coll](const auto& a){ return a.getSimID() == ef_coll[0].getObjectID().index; });
0179     auto ef_assoc = rcassoc.begin();
0180     for (; ef_assoc != rcassoc.end(); ++ef_assoc) {
0181       if (ef_assoc->getSimID() == (unsigned) ef_coll[0].getObjectID().index) {
0182         break;
0183       }
0184     }
0185     if (!(ef_assoc != rcassoc.end())) {
0186       if (msgLevel(MSG::DEBUG)) {
0187         debug() << "Truth scattered electron not in reconstructed particles" << endmsg;
0188       }
0189       return StatusCode::SUCCESS;
0190     }
0191     const auto ef_rc{ef_assoc->getRec()};
0192     const auto ef_rc_id{ef_rc.getObjectID().index};
0193
0194     // Loop over reconstructed particles to get outgoing scattered electron
0195     // Use the true scattered electron from the MC information
0196     std::vector<PxPyPzEVector> electrons;
0197     for (const auto& p: rcparts) {
0198       if (p.getObjectID().index == ef_rc_id) {
0199         electrons.emplace_back(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy());
0200         break;
0201       }
0202     }
0203
0204     // If no scattered electron was found
0205     if (electrons.size() == 0) {
0206       if (msgLevel(MSG::DEBUG)) {
0207         debug() << "No scattered electron found" << endmsg;
0208       }
0209       return StatusCode::SUCCESS;
0210     }
0211
0212     // DIS kinematics calculations
0213     const auto ef = electrons.front();
0214     const auto q = ei - ef;
0215     const auto q_dot_pi = q.Dot(pi);
0216     const auto Q2 = -q.Dot(q);
0217     const auto y = q_dot_pi / ei.Dot(pi);
0218     const auto nu = q_dot_pi / m_proton;
0219     const auto x = Q2 / (2. * q_dot_pi);
0220     const auto W = sqrt( + 2.*q_dot_pi - Q2);
0221     auto kin = out_kinematics.create(x, Q2, W, y, nu);
0222     kin.setScat(ef_rc);
0223
0224     // Debugging output
0225     if (msgLevel(MSG::DEBUG)) {
0226       debug() << "pi = " << pi << endmsg;
0227       debug() << "ei = " << ei << endmsg;
0228       debug() << "ef = " << ef << endmsg;
0229       debug() << "q = " << q << endmsg;
0230       debug() << "x,Q2,W,y,nu = "
0231               << kin.getX() << ","
0232               << kin.getQ2() << ","
0233               << kin.getW() << ","
0234               << kin.getY() << ","
0235               << kin.getNu()
0236               << endmsg;
0237     }
0238
0239     return StatusCode::SUCCESS;
0240   }
0241 };
0242
0243 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0244 DECLARE_COMPONENT(InclusiveKinematicsElectron)
0245
0246 } // namespace Jug::Reco