algorithms/reco/InclusiveKinematicsTruth.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Wouter Deconinck
0003
0004 #include <Math/GenVector/LorentzVector.h>
0005 #include <Math/GenVector/PxPyPzE4D.h>
0006 #include <Math/Vector4Dfwd.h>
0007 #include <edm4eic/InclusiveKinematicsCollection.h>
0008 #include <edm4hep/MCParticleCollection.h>
0009 #include <edm4hep/Vector3f.h>
0010 #include <edm4hep/utils/vector_utils.h>
0011 #include <fmt/core.h>
0012 #include <cmath>
0013 #include <gsl/pointers>
0014
0015 #include "Beam.h"
0016 #include "InclusiveKinematicsTruth.h"
0017
0018 using ROOT::Math::PxPyPzEVector;
0019
0020 namespace eicrecon {
0021
0022   void InclusiveKinematicsTruth::init() { }
0023
0024   void InclusiveKinematicsTruth::process(
0025       const InclusiveKinematicsTruth::Input& input,
0026       const InclusiveKinematicsTruth::Output& output) const {
0027
0028     const auto [mcparts] = input;
0029     auto [kinematics] = output;
0030
0031     // Loop over generated particles to get incoming electron and proton beams
0032     // and the scattered electron. In the presence of QED radition on the incoming
0033     // or outgoing electron line, the vertex kinematics will be different than the
0034     // kinematics calculated using the scattered electron as done here.
0035     // Also need to update for CC events.
0036
0037     // Get incoming electron beam
0038     const auto ei_coll = find_first_beam_electron(mcparts);
0039     if (ei_coll.size() == 0) {
0040       debug("No beam electron found");
0041       return;
0042     }
0043     const auto ei_p = ei_coll[0].getMomentum();
0044     const auto ei_p_mag = edm4hep::utils::magnitude(ei_p);
0045     static const auto ei_mass = m_particleSvc.particle(11).mass;
0046     const PxPyPzEVector ei(ei_p.x, ei_p.y, ei_p.z, std::hypot(ei_p_mag, ei_mass));
0047
0048     // Get incoming hadron beam
0049     const auto pi_coll = find_first_beam_hadron(mcparts);
0050     if (pi_coll.size() == 0) {
0051       debug("No beam hadron found");
0052       return;
0053     }
0054     const auto pi_p = pi_coll[0].getMomentum();
0055     const auto pi_p_mag = edm4hep::utils::magnitude(pi_p);
0056     const auto pi_mass = m_particleSvc.particle(pi_coll[0].getPDG()).mass;
0057     const PxPyPzEVector pi(pi_p.x, pi_p.y, pi_p.z, std::hypot(pi_p_mag, pi_mass));
0058
0059     // Get first scattered electron
0060     // Scattered electron. Currently taken as first status==1 electron in HEPMC record,
0061     // which seems to be correct based on a cursory glance at the Pythia8 output. In the future,
0062     // it may be better to trace back each final-state electron and see which one originates from
0063     // the beam.
0064     const auto ef_coll = find_first_scattered_electron(mcparts);
0065     if (ef_coll.size() == 0) {
0066       debug("No truth scattered electron found");
0067       return;
0068     }
0069     const auto ef_p = ef_coll[0].getMomentum();
0070     const auto ef_p_mag = edm4hep::utils::magnitude(ef_p);
0071     static const auto ef_mass = m_particleSvc.particle(11).mass;
0072     const PxPyPzEVector ef(ef_p.x, ef_p.y, ef_p.z, std::hypot(ef_p_mag, ef_mass));
0073
0074     // DIS kinematics calculations
0075     const auto q = ei - ef;
0076     const auto q_dot_pi = q.Dot(pi);
0077     const auto Q2 = -q.Dot(q);
0078     const auto y = q_dot_pi / ei.Dot(pi);
0079     const auto nu = q_dot_pi / pi_mass;
0080     const auto x = Q2 / (2.*q_dot_pi);
0081     const auto W = sqrt(pi_mass*pi_mass + 2.*q_dot_pi - Q2);
0082     auto kin = kinematics->create(x, Q2, W, y, nu);
0083
0084     debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(),
0085             kin.getQ2(), kin.getW(), kin.getY(), kin.getNu());
0086   }
0087
0088 } // namespace Jug::Reco