algorithms/reco/InclusiveKinematicsSigma.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Wouter Deconinck, Barak Schmookler
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 <cmath>
0011 #include <tuple>
0012
0013 #include "Beam.h"
0014 #include "Boost.h"
0015 #include "InclusiveKinematicsSigma.h"
0016
0017 using ROOT::Math::PxPyPzEVector;
0018
0019 namespace eicrecon {
0020
0021 void InclusiveKinematicsSigma::init() {}
0022
0023 void InclusiveKinematicsSigma::process(const InclusiveKinematicsSigma::Input& input,
0024                                        const InclusiveKinematicsSigma::Output& output) const {
0025
0026   const auto [mc_beam_electrons, mc_beam_protons, escat, hfs] = input;
0027   auto [out_kinematics]                                       = output;
0028
0029   // Get first (should be only) beam electron
0030   if (mc_beam_electrons->empty()) {
0031     debug("No beam electron found");
0032     return;
0033   }
0034   const auto& ei_particle = (*mc_beam_electrons)[0];
0035   const PxPyPzEVector ei(round_beam_four_momentum(ei_particle.getMomentum(),
0036                                                   m_particleSvc.particle(ei_particle.getPDG()).mass,
0037                                                   electron_beam_pz_set, 0.0));
0038
0039   // Get first (should be only) beam proton
0040   if (mc_beam_protons->empty()) {
0041     debug("No beam hadron found");
0042     return;
0043   }
0044   const auto& pi_particle = (*mc_beam_protons)[0];
0045   const PxPyPzEVector pi(round_beam_four_momentum(pi_particle.getMomentum(),
0046                                                   m_particleSvc.particle(pi_particle.getPDG()).mass,
0047                                                   hadron_beam_pz_set, m_crossingAngle));
0048
0049   // Get boost to colinear frame
0050   auto boost = determine_boost(ei, pi);
0051
0052   // Get electron variables
0053   if (escat->empty()) {
0054     debug("No scattered electron found");
0055     return;
0056   }
0057   auto kf = escat->at(0);
0058   PxPyPzEVector e_lab(kf.getMomentum().x, kf.getMomentum().y, kf.getMomentum().z, kf.getEnergy());
0059   PxPyPzEVector e_boosted = apply_boost(boost, e_lab);
0060   auto pt_e               = e_boosted.Pt();
0061   auto sigma_e            = e_boosted.E() - e_boosted.Pz();
0062
0063   // Get hadronic final state variables
0064   if (hfs->empty()) {
0065     debug("No hadronic final state found");
0066     return;
0067   }
0068   auto sigma_h = hfs->at(0).getSigma();
0069
0070   if (sigma_h <= 0) {
0071     debug("No scattered electron found or sigma zero or negative");
0072     return;
0073   }
0074
0075   auto sigma_tot = sigma_e + sigma_h;
0076
0077   // Calculate kinematic variables
0078   static const auto m_proton = m_particleSvc.particle(2212).mass;
0079   const auto y_sig           = sigma_h / sigma_tot;
0080   const auto Q2_sig          = (pt_e * pt_e) / (1. - y_sig);
0081   const auto x_sig           = Q2_sig / (4. * ei.energy() * pi.energy() * y_sig);
0082   const auto nu_sig          = Q2_sig / (2. * m_proton * x_sig);
0083   const auto W_sig           = sqrt(m_proton * m_proton + 2 * m_proton * nu_sig - Q2_sig);
0084   auto kin                   = out_kinematics->create(x_sig, Q2_sig, W_sig, y_sig, nu_sig);
0085   kin.setScat(kf);
0086
0087   debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(),
0088         kin.getNu());
0089 }
0090
0091 } // namespace eicrecon