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
0025   void InclusiveKinematicsJB::process(
0026       const InclusiveKinematicsJB::Input& input,
0027       const InclusiveKinematicsJB::Output& output) const {
0028
0029     const auto [mcparts, escat, hfs] = input;
0030     auto [kinematics] = output;
0031
0032     // Get incoming electron beam
0033     const auto ei_coll = find_first_beam_electron(mcparts);
0034     if (ei_coll.size() == 0) {
0035       debug("No beam electron found");
0036       return;
0037     }
0038     const PxPyPzEVector ei(
0039       round_beam_four_momentum(
0040         ei_coll[0].getMomentum(),
0041         m_particleSvc.particle(ei_coll[0].getPDG()).mass,
0042         {-5.0, -10.0, -18.0},
0043         0.0)
0044       );
0045
0046     // Get incoming hadron beam
0047     const auto pi_coll = find_first_beam_hadron(mcparts);
0048     if (pi_coll.size() == 0) {
0049       debug("No beam hadron found");
0050       return;
0051     }
0052     const PxPyPzEVector pi(
0053       round_beam_four_momentum(
0054         pi_coll[0].getMomentum(),
0055         m_particleSvc.particle(pi_coll[0].getPDG()).mass,
0056         {41.0, 100.0, 275.0},
0057         m_crossingAngle)
0058       );
0059
0060     // Get hadronic final state variables
0061     auto sigma_h = hfs->at(0).getSigma();
0062     auto ptsum = hfs->at(0).getPT();
0063     auto gamma_h = hfs->at(0).getGamma();
0064
0065     // Sigma zero or negative
0066     if (sigma_h <= 0) {
0067       debug("Sigma zero or negative");
0068       return;
0069     }
0070
0071     // Calculate kinematic variables
0072     static const auto m_proton = m_particleSvc.particle(2212).mass;
0073     const auto y_jb = sigma_h / (2.*ei.energy());
0074     const auto Q2_jb = ptsum*ptsum / (1. - y_jb);
0075     const auto x_jb = Q2_jb / (4.*ei.energy()*pi.energy()*y_jb);
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 = kinematics->create(x_jb, Q2_jb, W_jb, y_jb, nu_jb);
0079     kin.setScat(escat->at(0));
0080
0081     debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(),
0082             kin.getQ2(), kin.getW(), kin.getY(), kin.getNu());
0083   }
0084
0085 }
0086 #endif