algorithms/reco/InclusiveKinematicsJB.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/HadronicFinalStateCollection.h>
0008 #include <edm4eic/InclusiveKinematicsCollection.h>
0009 #include <cmath>
0010 #include <tuple>
0011
0012 #include "Beam.h"
0013 #include "InclusiveKinematicsJB.h"
0014
0015 using ROOT::Math::PxPyPzEVector;
0016
0017 namespace eicrecon {
0018
0019 void InclusiveKinematicsJB::init() {}
0020
0021 void InclusiveKinematicsJB::process(const InclusiveKinematicsJB::Input& input,
0022                                     const InclusiveKinematicsJB::Output& output) const {
0023
0024   const auto [mc_beam_electrons, mc_beam_protons, escat, hfs] = input;
0025   auto [out_kinematics]                                       = output;
0026
0027   // Get first (should be only) beam electron
0028   if (mc_beam_electrons->empty()) {
0029     debug("No beam electron found");
0030     return;
0031   }
0032   const auto& ei_particle = (*mc_beam_electrons)[0];
0033   const PxPyPzEVector ei(round_beam_four_momentum(ei_particle.getMomentum(),
0034                                                   m_particleSvc.particle(ei_particle.getPDG()).mass,
0035                                                   electron_beam_pz_set, 0.0));
0036
0037   // Get first (should be only) beam proton
0038   if (mc_beam_protons->empty()) {
0039     debug("No beam hadron found");
0040     return;
0041   }
0042   const auto& pi_particle = (*mc_beam_protons)[0];
0043   const PxPyPzEVector pi(round_beam_four_momentum(pi_particle.getMomentum(),
0044                                                   m_particleSvc.particle(pi_particle.getPDG()).mass,
0045                                                   hadron_beam_pz_set, m_crossingAngle));
0046
0047   // Get hadronic final state variables
0048   if (hfs->empty()) {
0049     debug("No hadronic final state found");
0050     return;
0051   }
0052   auto sigma_h = hfs->at(0).getSigma();
0053   auto ptsum   = hfs->at(0).getPT();
0054
0055   // Sigma zero or negative
0056   if (sigma_h <= 0) {
0057     debug("Sigma zero or negative");
0058     return;
0059   }
0060
0061   // Calculate kinematic variables
0062   static const auto m_proton = m_particleSvc.particle(2212).mass;
0063   const auto y_jb            = sigma_h / (2. * ei.energy());
0064   if (y_jb >= 1) {
0065     // y > 0 is mathematically guaranteed by sigma_h > 0, but y < 1 is not
0066     debug("InclusiveKinematicsJB: event with y >= 1 skipped");
0067     return;
0068   }
0069   const auto Q2_jb = ptsum * ptsum / (1. - y_jb);
0070   const auto x_jb  = Q2_jb / (4. * ei.energy() * pi.energy() * y_jb);
0071   if (x_jb >= 1) {
0072     // x > 0 is mathematically guaranteed by 0 < y < 1, but x < 1 is not
0073     debug("InclusiveKinematicsJB: event with x >= 1 skipped");
0074     return;
0075   }
0076   const auto nu_jb = Q2_jb / (2. * m_proton * x_jb);
0077   const auto W_jb  = sqrt(m_proton * m_proton + 2 * m_proton * nu_jb - Q2_jb);
0078   auto kin         = out_kinematics->create(x_jb, Q2_jb, W_jb, y_jb, nu_jb);
0079   if (escat->empty()) {
0080     debug("No scattered electron found");
0081   } else {
0082     kin.setScat(escat->at(0));
0083   }
0084
0085   debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(),
0086         kin.getNu());
0087 }
0088
0089 } // namespace eicrecon