algorithms/reco/InclusiveKinematicsESigma.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 <edm4eic/ReconstructedParticleCollection.h>
0009 #include <edm4hep/MCParticleCollection.h>
0010 #include <edm4hep/Vector3f.h>
0011 #include <fmt/core.h>
0012 #include <cmath>
0013 #include <gsl/pointers>
0014
0015 #include "Beam.h"
0016 #include "Boost.h"
0017 #include "InclusiveKinematicsESigma.h"
0018
0019 using ROOT::Math::PxPyPzEVector;
0020
0021 namespace eicrecon {
0022
0023 void InclusiveKinematicsESigma::init() {}
0024
0025 void InclusiveKinematicsESigma::process(const InclusiveKinematicsESigma::Input& input,
0026                                         const InclusiveKinematicsESigma::Output& output) const {
0027
0028   const auto [mcparts, escat, hfs] = input;
0029   auto [kinematics]                = output;
0030
0031   // Get incoming electron beam
0032   const auto ei_coll = find_first_beam_electron(mcparts);
0033   if (ei_coll.empty()) {
0034     debug("No beam electron found");
0035     return;
0036   }
0037   const PxPyPzEVector ei(round_beam_four_momentum(ei_coll[0].getMomentum(),
0038                                                   m_particleSvc.particle(ei_coll[0].getPDG()).mass,
0039                                                   {-5.0, -10.0, -18.0}, 0.0));
0040
0041   // Get incoming hadron beam
0042   const auto pi_coll = find_first_beam_hadron(mcparts);
0043   if (pi_coll.empty()) {
0044     debug("No beam hadron found");
0045     return;
0046   }
0047   const PxPyPzEVector pi(round_beam_four_momentum(pi_coll[0].getMomentum(),
0048                                                   m_particleSvc.particle(pi_coll[0].getPDG()).mass,
0049                                                   {41.0, 100.0, 275.0}, m_crossingAngle));
0050
0051   // Get boost to colinear frame
0052   auto boost = determine_boost(ei, pi);
0053
0054   // Get electron variables
0055   if (escat->empty()) {
0056     debug("No scattered electron found");
0057     return;
0058   }
0059   auto kf = escat->at(0);
0060   PxPyPzEVector e_lab(kf.getMomentum().x, kf.getMomentum().y, kf.getMomentum().z, kf.getEnergy());
0061   PxPyPzEVector e_boosted = apply_boost(boost, e_lab);
0062   auto pt_e               = e_boosted.Pt();
0063   auto sigma_e            = e_boosted.E() - e_boosted.Pz();
0064
0065   // Get hadronic final state variables
0066   if (hfs->empty()) {
0067     debug("No hadronic final state found");
0068     return;
0069   }
0070   auto sigma_h = hfs->at(0).getSigma();
0071
0072   if (sigma_h <= 0) {
0073     debug("No scattered electron found or sigma zero or negative");
0074     return;
0075   }
0076
0077   auto sigma_tot = sigma_e + sigma_h;
0078
0079   // Calculate kinematic variables
0080   const auto y_e  = 1. - sigma_e / (2. * ei.energy());
0081   const auto Q2_e = (pt_e * pt_e) / (1. - y_e);
0082
0083   const auto y_sig  = sigma_h / sigma_tot;
0084   const auto Q2_sig = (pt_e * pt_e) / (1. - y_sig);
0085   const auto x_sig  = Q2_sig / (4. * ei.energy() * pi.energy() * y_sig);
0086
0087   static const auto m_proton = m_particleSvc.particle(2212).mass;
0088   const auto Q2_esig         = Q2_e;
0089   const auto x_esig          = x_sig;
0090   const auto y_esig          = Q2_esig / (4. * ei.energy() * pi.energy() *
0091                                  x_esig); //equivalent to (2*ei.energy() / sigma_tot)*y_sig
0092   const auto nu_esig         = Q2_esig / (2. * m_proton * x_esig);
0093   const auto W_esig          = sqrt(m_proton * m_proton + 2 * m_proton * nu_esig - Q2_esig);
0094   auto kin                   = kinematics->create(x_esig, Q2_esig, W_esig, y_esig, nu_esig);
0095   kin.setScat(kf);
0096
0097   // Debugging output
0098   debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(),
0099         kin.getNu());
0100 }
0101
0102 } // namespace eicrecon