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 <fmt/core.h>
0012 #include <cmath>
0013 #include <gsl/pointers>
0014 #include <vector>
0015
0016 #include "Beam.h"
0017 #include "InclusiveKinematicsElectron.h"
0018
0019 using ROOT::Math::PxPyPzEVector;
0020
0021 namespace eicrecon {
0022
0023 void InclusiveKinematicsElectron::init() {}
0024
0025 void InclusiveKinematicsElectron::process(const InclusiveKinematicsElectron::Input& input,
0026                                           const InclusiveKinematicsElectron::Output& output) const {
0027
0028   const auto [mcparts, escat, hfs] = input;
0029   auto [kinematics]                = output;
0030
0031   // 1. find_if
0032   //const auto mc_first_electron = std::find_if(
0033   //  mcparts.begin(),
0034   //  mcparts.end(),
0035   //  [](const auto& p){ return p.getPDG() == 11; });
0036
0037   // 2a. simple loop over iterator (post-increment)
0038   //auto mc_first_electron = mcparts.end();
0039   //for (auto p = mcparts.begin(); p != mcparts.end(); p++) {
0040   //  if (p.getPDG() == 11) {
0041   //    mc_first_electron = p;
0042   //    break;
0043   //  }
0044   //}
0045   // 2b. simple loop over iterator (pre-increment)
0046   //auto mc_first_electron = mcparts.end();
0047   //for (auto p = mcparts.begin(); p != mcparts.end(); ++p) {
0048   //  if (p.getPDG() == 11) {
0049   //    mc_first_electron = p;
0050   //    break;
0051   //  }
0052   //}
0053
0054   // 3. pre-initialized simple loop
0055   //auto mc_first_electron = mcparts.begin();
0056   //for (; mc_first_electron != mcparts.end(); ++mc_first_electron) {
0057   //  if (mc_first_electron.getPDG() == 11) {
0058   //    break;
0059   //  }
0060   //}
0061
0062   // 4a. iterator equality
0063   //if (mc_first_electron == mcparts.end()) {
0064   //  debug() << "No electron found" << endmsg;
0065   //  return StatusCode::FAILURE;
0066   //}
0067   // 4b. iterator inequality
0068   //if (!(mc_first_electron != mcparts.end())) {
0069   //  debug() << "No electron found" << endmsg;
0070   //  return StatusCode::FAILURE;
0071   //}
0072
0073   // 5. ranges and views
0074   //auto is_electron = [](const auto& p){ return p.getPDG() == 11; };
0075   //for (const auto& e: mcparts | std::views::filter(is_electron)) {
0076   //  break;
0077   //}
0078
0079   // Get incoming electron beam
0080   const auto ei_coll = find_first_beam_electron(mcparts);
0081   if (ei_coll.empty()) {
0082     debug("No beam electron found");
0083     return;
0084   }
0085   const PxPyPzEVector ei(round_beam_four_momentum(ei_coll[0].getMomentum(),
0086                                                   m_particleSvc.particle(ei_coll[0].getPDG()).mass,
0087                                                   {-5.0, -10.0, -18.0}, 0.0));
0088
0089   // Get incoming hadron beam
0090   const auto pi_coll = find_first_beam_hadron(mcparts);
0091   if (pi_coll.empty()) {
0092     debug("No beam hadron found");
0093     return;
0094   }
0095   const PxPyPzEVector pi(round_beam_four_momentum(pi_coll[0].getMomentum(),
0096                                                   m_particleSvc.particle(pi_coll[0].getPDG()).mass,
0097                                                   {41.0, 100.0, 275.0}, m_crossingAngle));
0098
0099   // Get scattered electron
0100   std::vector<PxPyPzEVector> electrons;
0101   for (const auto& p : *escat) {
0102     electrons.emplace_back(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy());
0103     break;
0104   }
0105
0106   // If no scattered electron was found
0107   if (electrons.empty()) {
0108     debug("No scattered electron found");
0109     return;
0110   }
0111
0112   // DIS kinematics calculations
0113   static const auto m_proton = m_particleSvc.particle(2212).mass;
0114   const auto ef              = electrons.front();
0115   const auto q               = ei - ef;
0116   const auto q_dot_pi        = q.Dot(pi);
0117   const auto Q2              = -q.Dot(q);
0118   const auto y               = q_dot_pi / ei.Dot(pi);
0119   const auto nu              = q_dot_pi / m_proton;
0120   const auto x               = Q2 / (2. * q_dot_pi);
0121   const auto W               = sqrt(m_proton * m_proton + 2. * q_dot_pi - Q2);
0122   auto kin                   = kinematics->create(x, Q2, W, y, nu);
0123   kin.setScat(escat->at(0));
0124
0125   debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(),
0126         kin.getNu());
0127 }
0128
0129 } // namespace eicrecon