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 <fmt/core.h>
0010 #include <cmath>
0011 #include <gsl/pointers>
0012 #include <vector>
0013
0014 #include "Beam.h"
0015 #include "InclusiveKinematicsJB.h"
0016
0017 using ROOT::Math::PxPyPzEVector;
0018
0019 namespace eicrecon {
0020
0021 void InclusiveKinematicsJB::init() {}
0022
0023 void InclusiveKinematicsJB::process(const InclusiveKinematicsJB::Input& input,
0024                                     const InclusiveKinematicsJB::Output& output) const {
0025
0026   const auto [mcparts, escat, hfs] = input;
0027   auto [kinematics]                = output;
0028
0029   // Get incoming electron beam
0030   const auto ei_coll = find_first_beam_electron(mcparts);
0031   if (ei_coll.empty()) {
0032     debug("No beam electron found");
0033     return;
0034   }
0035   const PxPyPzEVector ei(round_beam_four_momentum(ei_coll[0].getMomentum(),
0036                                                   m_particleSvc.particle(ei_coll[0].getPDG()).mass,
0037                                                   {-5.0, -10.0, -18.0}, 0.0));
0038
0039   // Get incoming hadron beam
0040   const auto pi_coll = find_first_beam_hadron(mcparts);
0041   if (pi_coll.empty()) {
0042     debug("No beam hadron found");
0043     return;
0044   }
0045   const PxPyPzEVector pi(round_beam_four_momentum(pi_coll[0].getMomentum(),
0046                                                   m_particleSvc.particle(pi_coll[0].getPDG()).mass,
0047                                                   {41.0, 100.0, 275.0}, m_crossingAngle));
0048
0049   // Get hadronic final state variables
0050   if (hfs->empty()) {
0051     debug("No hadronic final state found");
0052     return;
0053   }
0054   auto sigma_h = hfs->at(0).getSigma();
0055   auto ptsum   = hfs->at(0).getPT();
0056
0057   // Sigma zero or negative
0058   if (sigma_h <= 0) {
0059     debug("Sigma zero or negative");
0060     return;
0061   }
0062
0063   // Calculate kinematic variables
0064   static const auto m_proton = m_particleSvc.particle(2212).mass;
0065   const auto y_jb            = sigma_h / (2. * ei.energy());
0066   const auto Q2_jb           = ptsum * ptsum / (1. - y_jb);
0067   const auto x_jb            = Q2_jb / (4. * ei.energy() * pi.energy() * y_jb);
0068   const auto nu_jb           = Q2_jb / (2. * m_proton * x_jb);
0069   const auto W_jb            = sqrt(m_proton * m_proton + 2 * m_proton * nu_jb - Q2_jb);
0070   auto kin                   = kinematics->create(x_jb, Q2_jb, W_jb, y_jb, nu_jb);
0071   if (escat->empty()) {
0072     debug("No scattered electron found");
0073   } else {
0074     kin.setScat(escat->at(0));
0075   }
0076
0077   debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(),
0078         kin.getNu());
0079 }
0080
0081 } // namespace eicrecon