src/todo/InclusiveKinematicsTruth.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Wouter Deconinck
0003
0004 #include "Gaudi/Algorithm.h"
0005 #include "GaudiAlg/GaudiTool.h"
0006 #include "GaudiAlg/Producer.h"
0007 #include "GaudiAlg/Transformer.h"
0008 #include "GaudiKernel/RndmGenerators.h"
0009 #include "GaudiKernel/PhysicalConstants.h"
0010 #include <algorithm>
0011 #include <cmath>
0012
0013 #include "JugBase/IParticleSvc.h"
0014 #include "JugBase/DataHandle.h"
0015
0016 #include "JugBase/Utilities/Beam.h"
0017
0018 #include "Math/Vector4D.h"
0019 using ROOT::Math::PxPyPzEVector;
0020
0021 // Event Model related classes
0022 #include "edm4hep/MCParticleCollection.h"
0023 #include "edm4eic/InclusiveKinematicsCollection.h"
0024
0025 #include "edm4hep/utils/vector_utils.h"
0026
0027 namespace Jug::Fast {
0028
0029 class InclusiveKinematicsTruth : public GaudiAlgorithm {
0030 private:
0031   DataHandle<edm4hep::MCParticleCollection> m_inputMCParticleCollection{
0032     "inputMCParticles",
0033     Gaudi::DataHandle::Reader,
0034     this};
0035   DataHandle<edm4eic::InclusiveKinematicsCollection> m_outputInclusiveKinematicsCollection{
0036     "outputInclusiveKinematics",
0037     Gaudi::DataHandle::Writer,
0038     this};
0039
0040   SmartIF<IParticleSvc> m_pidSvc;
0041   double m_proton{0};
0042   double m_neutron{0};
0043   double m_electron{0};
0044
0045 public:
0046   InclusiveKinematicsTruth(const std::string& name, ISvcLocator* svcLoc)
0047       : GaudiAlgorithm(name, svcLoc) {
0048     declareProperty("inputMCParticles", m_inputMCParticleCollection, "MCParticles");
0049     declareProperty("outputInclusiveKinematics", m_outputInclusiveKinematicsCollection, "InclusiveKinematicsTruth");
0050   }
0051
0052   StatusCode initialize() override {
0053     if (GaudiAlgorithm::initialize().isFailure()) {
0054       return StatusCode::FAILURE;
0055     }
0056
0057     m_pidSvc = service("ParticleSvc");
0058     if (!m_pidSvc) {
0059       error() << "Unable to locate Particle Service. "
0060               << "Make sure you have ParticleSvc in the configuration."
0061               << endmsg;
0062       return StatusCode::FAILURE;
0063     }
0064     m_proton = m_pidSvc->particle(2212).mass;
0065     m_neutron = m_pidSvc->particle(2112).mass;
0066     m_electron = m_pidSvc->particle(11).mass;
0067
0068     return StatusCode::SUCCESS;
0069   }
0070
0071   StatusCode execute() override {
0072     // input collections
0073     const auto& mcparts = *(m_inputMCParticleCollection.get());
0074     // output collection
0075     auto& out_kinematics = *(m_outputInclusiveKinematicsCollection.createAndPut());
0076
0077     // Loop over generated particles to get incoming electron and proton beams
0078     // and the scattered electron. In the presence of QED radition on the incoming
0079     // or outgoing electron line, the vertex kinematics will be different than the
0080     // kinematics calculated using the scattered electron as done here.
0081     // Also need to update for CC events.
0082
0083     // Get incoming electron beam
0084     const auto ei_coll = Jug::Base::Beam::find_first_beam_electron(mcparts);
0085     if (ei_coll.size() == 0) {
0086       if (msgLevel(MSG::DEBUG)) {
0087         debug() << "No beam electron found" << endmsg;
0088       }
0089       return StatusCode::SUCCESS;
0090     }
0091     const auto ei_p = ei_coll[0].getMomentum();
0092     const auto ei_p_mag = edm4hep::utils::magnitude(ei_p);
0093     const auto ei_mass = m_electron;
0094     const PxPyPzEVector ei(ei_p.x, ei_p.y, ei_p.z, std::hypot(ei_p_mag, ei_mass));
0095
0096     // Get incoming hadron beam
0097     const auto pi_coll = Jug::Base::Beam::find_first_beam_hadron(mcparts);
0098     if (pi_coll.size() == 0) {
0099       if (msgLevel(MSG::DEBUG)) {
0100         debug() << "No beam hadron found" << endmsg;
0101       }
0102       return StatusCode::SUCCESS;
0103     }
0104     const auto pi_p = pi_coll[0].getMomentum();
0105     const auto pi_p_mag = edm4hep::utils::magnitude(pi_p);
0106     const auto pi_mass = pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron;
0107     const PxPyPzEVector pi(pi_p.x, pi_p.y, pi_p.z, std::hypot(pi_p_mag, pi_mass));
0108
0109     // Get first scattered electron
0110     // Scattered electron. Currently taken as first status==1 electron in HEPMC record,
0111     // which seems to be correct based on a cursory glance at the Pythia8 output. In the future,
0112     // it may be better to trace back each final-state electron and see which one originates from
0113     // the beam.
0114     const auto ef_coll = Jug::Base::Beam::find_first_scattered_electron(mcparts);
0115     if (ef_coll.size() == 0) {
0116       if (msgLevel(MSG::DEBUG)) {
0117         debug() << "No truth scattered electron found" << endmsg;
0118       }
0119       return StatusCode::SUCCESS;
0120     }
0121     const auto ef_p = ef_coll[0].getMomentum();
0122     const auto ef_p_mag = edm4hep::utils::magnitude(ef_p);
0123     const auto ef_mass = m_electron;
0124     const PxPyPzEVector ef(ef_p.x, ef_p.y, ef_p.z, std::hypot(ef_p_mag, ef_mass));
0125
0126     // DIS kinematics calculations
0127     const auto q = ei - ef;
0128     const auto q_dot_pi = q.Dot(pi);
0129     const auto Q2 = -q.Dot(q);
0130     const auto y = q_dot_pi / ei.Dot(pi);
0131     const auto nu = q_dot_pi / pi_mass;
0132     const auto x = Q2 / (2.*q_dot_pi);
0133     const auto W = sqrt(pi_mass*pi_mass + 2.*q_dot_pi - Q2);
0134     const auto kin = out_kinematics.create(x, Q2, W, y, nu);
0135
0136     // Debugging output
0137     if (msgLevel(MSG::DEBUG)) {
0138       debug() << "pi = " << pi << endmsg;
0139       debug() << "ei = " << ei << endmsg;
0140       debug() << "ef = " << ef << endmsg;
0141       debug() << "q = " << q << endmsg;
0142       debug() << "x,Q2,W,y,nu = "
0143               << kin.getX() << ","
0144               << kin.getQ2() << ","
0145               << kin.getW() << ","
0146               << kin.getY() << ","
0147               << kin.getNu()
0148               << endmsg;
0149     }
0150
0151     return StatusCode::SUCCESS;
0152   }
0153 };
0154
0155 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0156 DECLARE_COMPONENT(InclusiveKinematicsTruth)
0157
0158 } // namespace Jug::Fast