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