algorithms/reco/InclusiveKinematicsDA.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Wouter Deconinck
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
0014 #include "Beam.h"
0015 #include "Boost.h"
0016 #include "InclusiveKinematicsDA.h"
0017
0018 using ROOT::Math::PxPyPzEVector;
0019
0020 namespace eicrecon {
0021
0022 void InclusiveKinematicsDA::init() {}
0023
0024 void InclusiveKinematicsDA::process(const InclusiveKinematicsDA::Input& input,
0025                                     const InclusiveKinematicsDA::Output& output) const {
0026
0027   const auto [mc_beam_electrons, mc_beam_protons, escat, hfs] = input;
0028   auto [out_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 boost to colinear frame
0051   auto boost = determine_boost(ei, pi);
0052
0053   // Get electron angle
0054   if (escat->empty()) {
0055     debug("No scattered electron found");
0056     return;
0057   }
0058   auto kf = escat->at(0);
0059   PxPyPzEVector e_lab(kf.getMomentum().x, kf.getMomentum().y, kf.getMomentum().z, kf.getEnergy());
0060   PxPyPzEVector e_boosted = apply_boost(boost, e_lab);
0061   auto theta_e            = e_boosted.Theta();
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   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   static const auto m_proton = m_particleSvc.particle(2212).mass;
0079   const auto y_da            = tan(gamma_h / 2.) / (tan(theta_e / 2.) + tan(gamma_h / 2.));
0080   const auto Q2_da           = 4. * ei.energy() * ei.energy() * (1. / tan(theta_e / 2.)) *
0081                                (1. / (tan(theta_e / 2.) + tan(gamma_h / 2.)));
0082   const auto x_da            = Q2_da / (4. * ei.energy() * pi.energy() * y_da);
0083   const auto nu_da           = Q2_da / (2. * m_proton * x_da);
0084   const auto W_da            = sqrt(m_proton * m_proton + 2 * m_proton * nu_da - Q2_da);
0085   auto kin                   = out_kinematics->create(x_da, Q2_da, W_da, y_da, nu_da);
0086   kin.setScat(kf);
0087
0088   debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(),
0089         kin.getNu());
0090 }
0091
0092 } // namespace eicrecon