algorithms/reco/InclusiveKinematicsJB.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Wouter Deconinck
0003
0004 #include <edm4eic/EDM4eicVersion.h>
0005 #if EDM4EIC_VERSION_MAJOR >= 6
0006
0007 #include <Math/GenVector/LorentzVector.h>
0008 #include <Math/GenVector/PxPyPzE4D.h>
0009 #include <Math/Vector4Dfwd.h>
0010 #include <edm4eic/HadronicFinalStateCollection.h>
0011 #include <edm4eic/InclusiveKinematicsCollection.h>
0012 #include <fmt/core.h>
0013 #include <cmath>
0014 #include <gsl/pointers>
0015
0016 #include "Beam.h"
0017 #include "InclusiveKinematicsJB.h"
0018
0019 using ROOT::Math::PxPyPzEVector;
0020
0021 namespace eicrecon {
0022
0023   void InclusiveKinematicsJB::init() {
0024     // m_pidSvc = service("ParticleSvc");
0025     // if (!m_pidSvc) {
0026     //   debug("Unable to locate Particle Service. "
0027     //     "Make sure you have ParticleSvc in the configuration.");
0028     // }
0029   }
0030
0031   void InclusiveKinematicsJB::process(
0032       const InclusiveKinematicsJB::Input& input,
0033       const InclusiveKinematicsJB::Output& output) const {
0034
0035     const auto [mcparts, escat, hfs] = input;
0036     auto [kinematics] = output;
0037
0038     // Get incoming electron beam
0039     const auto ei_coll = find_first_beam_electron(mcparts);
0040     if (ei_coll.size() == 0) {
0041       debug("No beam electron found");
0042       return;
0043     }
0044     const PxPyPzEVector ei(
0045       round_beam_four_momentum(
0046         ei_coll[0].getMomentum(),
0047         m_electron,
0048         {-5.0, -10.0, -18.0},
0049         0.0)
0050       );
0051
0052     // Get incoming hadron beam
0053     const auto pi_coll = find_first_beam_hadron(mcparts);
0054     if (pi_coll.size() == 0) {
0055       debug("No beam hadron found");
0056       return;
0057     }
0058     const PxPyPzEVector pi(
0059       round_beam_four_momentum(
0060         pi_coll[0].getMomentum(),
0061         pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron,
0062         {41.0, 100.0, 275.0},
0063         m_crossingAngle)
0064       );
0065
0066     // Get hadronic final state variables
0067     auto sigma_h = hfs->at(0).getSigma();
0068     auto ptsum = hfs->at(0).getPT();
0069     auto gamma_h = hfs->at(0).getGamma();
0070
0071     // Sigma zero or negative
0072     if (sigma_h <= 0) {
0073       debug("Sigma zero or negative");
0074       return;
0075     }
0076
0077     // Calculate kinematic variables
0078     const auto y_jb = sigma_h / (2.*ei.energy());
0079     const auto Q2_jb = ptsum*ptsum / (1. - y_jb);
0080     const auto x_jb = Q2_jb / (4.*ei.energy()*pi.energy()*y_jb);
0081     const auto nu_jb = Q2_jb / (2.*m_proton*x_jb);
0082     const auto W_jb = sqrt(m_proton*m_proton + 2*m_proton*nu_jb - Q2_jb);
0083     auto kin = kinematics->create(x_jb, Q2_jb, W_jb, y_jb, nu_jb);
0084     kin.setScat(escat->at(0));
0085
0086     debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(),
0087             kin.getQ2(), kin.getW(), kin.getY(), kin.getNu());
0088   }
0089
0090 }
0091 #endif