src/todo/InclusiveKinematicsElectron.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Wouter Deconinck
0003
0004 #include "Gaudi/Algorithm.h"
0005 #include "GaudiAlg/GaudiTool.h"
0006 #include "GaudiAlg/Producer.h"
0007 #include "GaudiAlg/Transformer.h"
0008 #include "GaudiKernel/RndmGenerators.h"
0009 #include "GaudiKernel/PhysicalConstants.h"
0010 #include <algorithm>
0011 #include <cmath>
0012
0013 #include "JugBase/IParticleSvc.h"
0014 #include "JugBase/DataHandle.h"
0015
0016 #include "JugBase/Utilities/Beam.h"
0017 #include "JugBase/Utilities/Boost.h"
0018
0019 #include "Math/Vector4D.h"
0020 using ROOT::Math::PxPyPzEVector;
0021
0022 // Event Model related classes
0023 #include "edm4hep/MCParticleCollection.h"
0024 #include "edm4eic/MCRecoParticleAssociationCollection.h"
0025 #include "edm4eic/ReconstructedParticleCollection.h"
0026 #include "edm4eic/InclusiveKinematicsCollection.h"
0027
0028 namespace Jug::Reco {
0029
0030 class InclusiveKinematicsElectron : public GaudiAlgorithm {
0031 private:
0032   DataHandle<edm4hep::MCParticleCollection> m_inputMCParticleCollection{
0033     "inputMCParticles",
0034     Gaudi::DataHandle::Reader,
0035     this};
0036   DataHandle<edm4eic::ReconstructedParticleCollection> m_inputParticleCollection{
0037     "inputReconstructedParticles",
0038     Gaudi::DataHandle::Reader,
0039     this};
0040   DataHandle<edm4eic::MCRecoParticleAssociationCollection> m_inputParticleAssociation{
0041     "inputParticleAssociations",
0042     Gaudi::DataHandle::Reader,
0043     this};
0044   DataHandle<edm4eic::InclusiveKinematicsCollection> m_outputInclusiveKinematicsCollection{
0045     "outputInclusiveKinematics",
0046     Gaudi::DataHandle::Writer,
0047     this};
0048
0049   Gaudi::Property<double> m_crossingAngle{this, "crossingAngle", -0.025 * Gaudi::Units::radian};
0050
0051   SmartIF<IParticleSvc> m_pidSvc;
0052   double m_proton{0}, m_neutron{0}, m_electron{0};
0053
0054 public:
0055   InclusiveKinematicsElectron(const std::string& name, ISvcLocator* svcLoc)
0056       : GaudiAlgorithm(name, svcLoc) {
0057     declareProperty("inputMCParticles", m_inputMCParticleCollection, "MCParticles");
0058     declareProperty("inputReconstructedParticles", m_inputParticleCollection, "ReconstructedParticles");
0059     declareProperty("inputParticleAssociations", m_inputParticleAssociation, "MCRecoParticleAssociation");
0060     declareProperty("outputInclusiveKinematics", m_outputInclusiveKinematicsCollection, "InclusiveKinematicsElectron");
0061   }
0062
0063   StatusCode initialize() override {
0064     if (GaudiAlgorithm::initialize().isFailure())
0065       return StatusCode::FAILURE;
0066
0067     m_pidSvc = service("ParticleSvc");
0068     if (!m_pidSvc) {
0069       error() << "Unable to locate Particle Service. "
0070               << "Make sure you have ParticleSvc in the configuration."
0071               << endmsg;
0072       return StatusCode::FAILURE;
0073     }
0074     m_proton = m_pidSvc->particle(2212).mass;
0075     m_neutron = m_pidSvc->particle(2112).mass;
0076     m_electron = m_pidSvc->particle(11).mass;
0077
0078     return StatusCode::SUCCESS;
0079   }
0080
0081   StatusCode execute() override {
0082     // input collections
0083     const auto& mcparts = *(m_inputMCParticleCollection.get());
0084     const auto& rcparts = *(m_inputParticleCollection.get());
0085     const auto& rcassoc = *(m_inputParticleAssociation.get());
0086     // output collection
0087     auto& out_kinematics = *(m_outputInclusiveKinematicsCollection.createAndPut());
0088
0089     // 1. find_if
0090     //const auto mc_first_electron = std::find_if(
0091     //  mcparts.begin(),
0092     //  mcparts.end(),
0093     //  [](const auto& p){ return p.getPDG() == 11; });
0094
0095     // 2a. simple loop over iterator (post-increment)
0096     //auto mc_first_electron = mcparts.end();
0097     //for (auto p = mcparts.begin(); p != mcparts.end(); p++) {
0098     //  if (p->getPDG() == 11) {
0099     //    mc_first_electron = p;
0100     //    break;
0101     //  }
0102     //}
0103     // 2b. simple loop over iterator (pre-increment)
0104     //auto mc_first_electron = mcparts.end();
0105     //for (auto p = mcparts.begin(); p != mcparts.end(); ++p) {
0106     //  if (p->getPDG() == 11) {
0107     //    mc_first_electron = p;
0108     //    break;
0109     //  }
0110     //}
0111
0112     // 3. pre-initialized simple loop
0113     //auto mc_first_electron = mcparts.begin();
0114     //for (; mc_first_electron != mcparts.end(); ++mc_first_electron) {
0115     //  if (mc_first_electron->getPDG() == 11) {
0116     //    break;
0117     //  }
0118     //}
0119
0120     // 4a. iterator equality
0121     //if (mc_first_electron == mcparts.end()) {
0122     //  debug() << "No electron found" << endmsg;
0123     //  return StatusCode::FAILURE;
0124     //}
0125     // 4b. iterator inequality
0126     //if (!(mc_first_electron != mcparts.end())) {
0127     //  debug() << "No electron found" << endmsg;
0128     //  return StatusCode::FAILURE;
0129     //}
0130
0131     // 5. ranges and views
0132     //auto is_electron = [](const auto& p){ return p.getPDG() == 11; };
0133     //for (const auto& e: mcparts | std::views::filter(is_electron)) {
0134     //  break;
0135     //}
0136
0137     // Get incoming electron beam
0138     const auto ei_coll = Jug::Base::Beam::find_first_beam_electron(mcparts);
0139     if (ei_coll.size() == 0) {
0140       if (msgLevel(MSG::DEBUG)) {
0141         debug() << "No beam electron found" << endmsg;
0142       }
0143       return StatusCode::SUCCESS;
0144     }
0145     const PxPyPzEVector ei(
0146       Jug::Base::Beam::round_beam_four_momentum(
0147         ei_coll[0].getMomentum(),
0148         m_electron,
0149         {-5.0, -10.0, -18.0},
0150         0.0)
0151       );
0152
0153     // Get incoming hadron beam
0154     const auto pi_coll = Jug::Base::Beam::find_first_beam_hadron(mcparts);
0155     if (pi_coll.size() == 0) {
0156       if (msgLevel(MSG::DEBUG)) {
0157         debug() << "No beam hadron found" << endmsg;
0158       }
0159       return StatusCode::SUCCESS;
0160     }
0161     const PxPyPzEVector pi(
0162       Jug::Base::Beam::round_beam_four_momentum(
0163         pi_coll[0].getMomentum(),
0164         pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron,
0165         {41.0, 100.0, 275.0},
0166         m_crossingAngle)
0167       );
0168
0169     // Get first scattered electron
0170     const auto ef_coll = Jug::Base::Beam::find_first_scattered_electron(mcparts);
0171     if (ef_coll.size() == 0) {
0172       if (msgLevel(MSG::DEBUG)) {
0173         debug() << "No truth scattered electron found" << endmsg;
0174       }
0175       return StatusCode::SUCCESS;
0176     }
0177     // Associate first scattered electron with reconstructed electrons
0178     //const auto ef_assoc = std::find_if(
0179     //  rcassoc.begin(),
0180     //  rcassoc.end(),
0181     //  [&ef_coll](const auto& a){ return a.getSimID() == ef_coll[0].getObjectID().index; });
0182     auto ef_assoc = rcassoc.begin();
0183     for (; ef_assoc != rcassoc.end(); ++ef_assoc) {
0184       if (ef_assoc->getSimID() == (unsigned) ef_coll[0].getObjectID().index) {
0185         break;
0186       }
0187     }
0188     if (!(ef_assoc != rcassoc.end())) {
0189       if (msgLevel(MSG::DEBUG)) {
0190         debug() << "Truth scattered electron not in reconstructed particles" << endmsg;
0191       }
0192       return StatusCode::SUCCESS;
0193     }
0194     const auto ef_rc{ef_assoc->getRec()};
0195     const auto ef_rc_id{ef_rc.getObjectID().index};
0196
0197     // Loop over reconstructed particles to get outgoing scattered electron
0198     // Use the true scattered electron from the MC information
0199     std::vector<PxPyPzEVector> electrons;
0200     for (const auto& p: rcparts) {
0201       if (p.getObjectID().index == ef_rc_id) {
0202         electrons.emplace_back(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy());
0203         break;
0204       }
0205     }
0206
0207     // If no scattered electron was found
0208     if (electrons.size() == 0) {
0209       if (msgLevel(MSG::DEBUG)) {
0210         debug() << "No scattered electron found" << endmsg;
0211       }
0212       return StatusCode::SUCCESS;
0213     }
0214
0215     // DIS kinematics calculations
0216     const auto ef = electrons.front();
0217     const auto q = ei - ef;
0218     const auto q_dot_pi = q.Dot(pi);
0219     const auto Q2 = -q.Dot(q);
0220     const auto y = q_dot_pi / ei.Dot(pi);
0221     const auto nu = q_dot_pi / m_proton;
0222     const auto x = Q2 / (2. * q_dot_pi);
0223     const auto W = sqrt( + 2.*q_dot_pi - Q2);
0224     auto kin = out_kinematics.create(x, Q2, W, y, nu);
0225     kin.setScat(ef_rc);
0226
0227     // Debugging output
0228     if (msgLevel(MSG::DEBUG)) {
0229       debug() << "pi = " << pi << endmsg;
0230       debug() << "ei = " << ei << endmsg;
0231       debug() << "ef = " << ef << endmsg;
0232       debug() << "q = " << q << endmsg;
0233       debug() << "x,Q2,W,y,nu = "
0234               << kin.getX() << ","
0235               << kin.getQ2() << ","
0236               << kin.getW() << ","
0237               << kin.getY() << ","
0238               << kin.getNu()
0239               << endmsg;
0240     }
0241
0242     return StatusCode::SUCCESS;
0243   }
0244 };
0245
0246 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0247 DECLARE_COMPONENT(InclusiveKinematicsElectron)
0248
0249 } // namespace Jug::Reco