algorithms/reco/InclusiveKinematicsElectron.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022, 2023 Wouter Deconinck, Tooba Ali
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 <cmath>
0012 #include <tuple>
0013 #include <vector>
0014
0015 #include "Beam.h"
0016 #include "InclusiveKinematicsElectron.h"
0017
0018 using ROOT::Math::PxPyPzEVector;
0019
0020 namespace eicrecon {
0021
0022 void InclusiveKinematicsElectron::init() {}
0023
0024 void InclusiveKinematicsElectron::process(const InclusiveKinematicsElectron::Input& input,
0025                                           const InclusiveKinematicsElectron::Output& output) const {
0026
0027   const auto [mc_beam_electrons, mc_beam_protons, escat, hfs] = input;
0028   auto [kinematics]                                           = output;
0029
0030   // Get first (should be only) beam electron
0031   if (mc_beam_electrons->empty()) {
0032     debug("No beam electron found");
0033     return;
0034   }
0035   const auto& ei_particle = (*mc_beam_electrons)[0];
0036   const PxPyPzEVector ei(round_beam_four_momentum(ei_particle.getMomentum(),
0037                                                   m_particleSvc.particle(ei_particle.getPDG()).mass,
0038                                                   electron_beam_pz_set, 0.0));
0039
0040   // Get first (should be only) beam proton
0041   if (mc_beam_protons->empty()) {
0042     debug("No beam hadron found");
0043     return;
0044   }
0045   const auto& pi_particle = (*mc_beam_protons)[0];
0046   const PxPyPzEVector pi(round_beam_four_momentum(pi_particle.getMomentum(),
0047                                                   m_particleSvc.particle(pi_particle.getPDG()).mass,
0048                                                   hadron_beam_pz_set, m_crossingAngle));
0049
0050   // Get scattered electron
0051   std::vector<PxPyPzEVector> electrons;
0052   for (const auto& p : *escat) {
0053     electrons.emplace_back(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy());
0054     break;
0055   }
0056
0057   // If no scattered electron was found
0058   if (electrons.empty()) {
0059     debug("No scattered electron found");
0060     return;
0061   }
0062
0063   // DIS kinematics calculations
0064   static const auto m_proton = m_particleSvc.particle(2212).mass;
0065   const auto ef              = electrons.front();
0066   const auto q               = ei - ef;
0067   const auto q_dot_pi        = q.Dot(pi);
0068   const auto Q2              = -q.Dot(q);
0069   const auto y               = q_dot_pi / ei.Dot(pi);
0070   const auto nu              = q_dot_pi / m_proton;
0071   const auto x               = Q2 / (2. * q_dot_pi);
0072   const auto W               = sqrt(m_proton * m_proton + 2. * q_dot_pi - Q2);
0073   auto kin                   = kinematics->create(x, Q2, W, y, nu);
0074   kin.setScat(escat->at(0));
0075
0076   debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(),
0077         kin.getNu());
0078 }
0079
0080 } // namespace eicrecon