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 <fmt/core.h>
0012 #include <cmath>
0013 #include <gsl/pointers>
0014 #include <vector>
0015
0016 #include "Beam.h"
0017 #include "Boost.h"
0018 #include "InclusiveKinematicsDA.h"
0019
0020 using ROOT::Math::PxPyPzEVector;
0021
0022 namespace eicrecon {
0023
0024 void InclusiveKinematicsDA::init() {}
0025
0026 void InclusiveKinematicsDA::process(const InclusiveKinematicsDA::Input& input,
0027                                     const InclusiveKinematicsDA::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.empty()) {
0035     debug("No beam electron found");
0036     return;
0037   }
0038   const PxPyPzEVector ei(round_beam_four_momentum(ei_coll[0].getMomentum(),
0039                                                   m_particleSvc.particle(ei_coll[0].getPDG()).mass,
0040                                                   {-5.0, -10.0, -18.0}, 0.0));
0041
0042   // Get incoming hadron beam
0043   const auto pi_coll = find_first_beam_hadron(mcparts);
0044   if (pi_coll.empty()) {
0045     debug("No beam hadron found");
0046     return;
0047   }
0048   const PxPyPzEVector pi(round_beam_four_momentum(pi_coll[0].getMomentum(),
0049                                                   m_particleSvc.particle(pi_coll[0].getPDG()).mass,
0050                                                   {41.0, 100.0, 275.0}, m_crossingAngle));
0051
0052   // Get boost to colinear frame
0053   auto boost = determine_boost(ei, pi);
0054
0055   // Get electron angle
0056   if (escat->empty()) {
0057     debug("No scattered electron found");
0058     return;
0059   }
0060   auto kf = escat->at(0);
0061   PxPyPzEVector e_lab(kf.getMomentum().x, kf.getMomentum().y, kf.getMomentum().z, kf.getEnergy());
0062   PxPyPzEVector e_boosted = apply_boost(boost, e_lab);
0063   auto theta_e            = e_boosted.Theta();
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   auto gamma_h = hfs->at(0).getGamma();
0072
0073   // Sigma zero or negative
0074   if (sigma_h <= 0) {
0075     debug("Sigma zero or negative");
0076     return;
0077   }
0078
0079   // Calculate kinematic variables
0080   static const auto m_proton = m_particleSvc.particle(2212).mass;
0081   const auto y_da            = tan(gamma_h / 2.) / (tan(theta_e / 2.) + tan(gamma_h / 2.));
0082   const auto Q2_da           = 4. * ei.energy() * ei.energy() * (1. / tan(theta_e / 2.)) *
0083                      (1. / (tan(theta_e / 2.) + tan(gamma_h / 2.)));
0084   const auto x_da  = Q2_da / (4. * ei.energy() * pi.energy() * y_da);
0085   const auto nu_da = Q2_da / (2. * m_proton * x_da);
0086   const auto W_da  = sqrt(m_proton * m_proton + 2 * m_proton * nu_da - Q2_da);
0087   auto kin         = kinematics->create(x_da, Q2_da, W_da, y_da, nu_da);
0088   kin.setScat(kf);
0089
0090   debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(),
0091         kin.getNu());
0092 }
0093
0094 } // namespace eicrecon