src/components/InclusiveKinematicsJB.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Wouter Deconinck
0003
0004 #include "Gaudi/Algorithm.h"
0005 #include "GaudiKernel/PhysicalConstants.h"
0006 #include "GaudiKernel/RndmGenerators.h"
0007 #include <algorithm>
0008 #include <cmath>
0009
0010 #include "JugBase/IParticleSvc.h"
0011 #include <k4FWCore/DataHandle.h>
0012
0013 #include "JugBase/Utilities/Beam.h"
0014 #include "JugBase/Utilities/Boost.h"
0015
0016 #include "Math/Vector4D.h"
0017 using ROOT::Math::PxPyPzEVector;
0018
0019 // Event Model related classes
0020 #include "edm4hep/MCParticleCollection.h"
0021 #include "edm4eic/MCRecoParticleAssociationCollection.h"
0022 #include "edm4eic/ReconstructedParticleCollection.h"
0023 #include "edm4eic/InclusiveKinematicsCollection.h"
0024
0025 namespace Jug::Reco {
0026
0027 class InclusiveKinematicsJB : public Gaudi::Algorithm {
0028 private:
0029   mutable DataHandle<edm4hep::MCParticleCollection> m_inputMCParticleCollection{
0030     "inputMCParticles",
0031     Gaudi::DataHandle::Reader,
0032     this};
0033   mutable DataHandle<edm4eic::ReconstructedParticleCollection> m_inputParticleCollection{
0034     "inputReconstructedParticles",
0035     Gaudi::DataHandle::Reader,
0036     this};
0037   mutable DataHandle<edm4eic::MCRecoParticleAssociationCollection> m_inputParticleAssociation{
0038     "inputParticleAssociations",
0039     Gaudi::DataHandle::Reader,
0040     this};
0041   mutable DataHandle<edm4eic::InclusiveKinematicsCollection> m_outputInclusiveKinematicsCollection{
0042     "outputInclusiveKinematics",
0043     Gaudi::DataHandle::Writer,
0044     this};
0045
0046   Gaudi::Property<double> m_crossingAngle{this, "crossingAngle", -0.025 * Gaudi::Units::radian};
0047
0048   SmartIF<IParticleSvc> m_pidSvc;
0049   double m_proton{0}, m_neutron{0}, m_electron{0};
0050
0051 public:
0052   InclusiveKinematicsJB(const std::string& name, ISvcLocator* svcLoc)
0053       : Gaudi::Algorithm(name, svcLoc) {
0054     declareProperty("inputMCParticles", m_inputMCParticleCollection, "MCParticles");
0055     declareProperty("inputReconstructedParticles", m_inputParticleCollection, "ReconstructedParticles");
0056     declareProperty("inputParticleAssociations", m_inputParticleAssociation, "MCRecoParticleAssociation");
0057     declareProperty("outputInclusiveKinematics", m_outputInclusiveKinematicsCollection, "InclusiveKinematicsJB");
0058   }
0059
0060   StatusCode initialize() override {
0061     if (Gaudi::Algorithm::initialize().isFailure())
0062       return StatusCode::FAILURE;
0063
0064     m_pidSvc = service("ParticleSvc");
0065     if (!m_pidSvc) {
0066       error() << "Unable to locate Particle Service. "
0067               << "Make sure you have ParticleSvc in the configuration."
0068               << endmsg;
0069       return StatusCode::FAILURE;
0070     }
0071     m_proton = m_pidSvc->particle(2212).mass;
0072     m_neutron = m_pidSvc->particle(2112).mass;
0073     m_electron = m_pidSvc->particle(11).mass;
0074
0075
0076     return StatusCode::SUCCESS;
0077   }
0078
0079   StatusCode execute(const EventContext&) const override {
0080     // input collections
0081     const auto& mcparts = *(m_inputMCParticleCollection.get());
0082     const auto& rcparts = *(m_inputParticleCollection.get());
0083     const auto& rcassoc = *(m_inputParticleAssociation.get());
0084     // output collection
0085     auto& out_kinematics = *(m_outputInclusiveKinematicsCollection.createAndPut());
0086
0087     // Get incoming electron beam
0088     const auto ei_coll = Jug::Base::Beam::find_first_beam_electron(mcparts);
0089     if (ei_coll.size() == 0) {
0090       if (msgLevel(MSG::DEBUG)) {
0091         debug() << "No beam electron found" << endmsg;
0092       }
0093       return StatusCode::SUCCESS;
0094     }
0095     const PxPyPzEVector ei(
0096       Jug::Base::Beam::round_beam_four_momentum(
0097         ei_coll[0].getMomentum(),
0098         m_electron,
0099         {-5.0, -10.0, -18.0},
0100         0.0)
0101       );
0102
0103     // Get incoming hadron beam
0104     const auto pi_coll = Jug::Base::Beam::find_first_beam_hadron(mcparts);
0105     if (pi_coll.size() == 0) {
0106       if (msgLevel(MSG::DEBUG)) {
0107         debug() << "No beam hadron found" << endmsg;
0108       }
0109       return StatusCode::SUCCESS;
0110     }
0111     const PxPyPzEVector pi(
0112       Jug::Base::Beam::round_beam_four_momentum(
0113         pi_coll[0].getMomentum(),
0114         pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron,
0115         {41.0, 100.0, 275.0},
0116         m_crossingAngle)
0117       );
0118
0119     // Get first scattered electron
0120     const auto ef_coll = Jug::Base::Beam::find_first_scattered_electron(mcparts);
0121     if (ef_coll.size() == 0) {
0122       if (msgLevel(MSG::DEBUG)) {
0123         debug() << "No truth scattered electron found" << endmsg;
0124       }
0125       return StatusCode::SUCCESS;
0126     }
0127     // Associate first scattered electron with reconstructed electrons
0128     //const auto ef_assoc = std::find_if(
0129     //  rcassoc.begin(),
0130     //  rcassoc.end(),
0131     //  [&ef_coll](const auto& a){ return a.getSimID() == ef_coll[0].getObjectID().index; });
0132     auto ef_assoc = rcassoc.begin();
0133     for (; ef_assoc != rcassoc.end(); ++ef_assoc) {
0134       if (ef_assoc->getSimID() == (unsigned) ef_coll[0].getObjectID().index) {
0135         break;
0136       }
0137     }
0138     if (!(ef_assoc != rcassoc.end())) {
0139       if (msgLevel(MSG::DEBUG)) {
0140         debug() << "Truth scattered electron not in reconstructed particles" << endmsg;
0141       }
0142       return StatusCode::SUCCESS;
0143     }
0144     const auto ef_rc{ef_assoc->getRec()};
0145     const auto ef_rc_id{ef_rc.getObjectID().index};
0146
0147     // Loop over reconstructed particles to get all outgoing particles other than the scattered electron
0148     // -----------------------------------------------------------------
0149     // Right now, everything is taken from Reconstructed particles branches.
0150     //
0151     // This means the tracking detector is used for charged particles to caculate the momentum,
0152     // and the magnitude of this momentum plus the true PID to calculate the energy.
0153     // No requirement is made that these particles produce a hit in any other detector
0154     //
0155     // Using the Reconstructed particles branches also means that the reconstruction for neutrals is done using the
0156     // calorimeter(s) information for the energy and angles, and then using this energy and the true PID to get the
0157     // magnitude of the momentum.
0158     // -----------------------------------------------------------------
0159
0160     // Sums in colinear frame
0161     double pxsum = 0;
0162     double pysum = 0;
0163     double pzsum = 0;
0164     double Esum = 0;
0165
0166     // Get boost to colinear frame
0167     auto boost = Jug::Base::Boost::determine_boost(ei, pi);
0168
0169     for (const auto& p: rcparts) {
0170       // Get the scattered electron index and angle
0171       if (p.getObjectID().index == ef_rc_id) {
0172
0173       // Sum over all particles other than scattered electron
0174       } else {
0175         // Lorentz vector in lab frame
0176         PxPyPzEVector hf_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy());
0177         // Boost to colinear frame
0178         PxPyPzEVector hf_boosted = Jug::Base::Boost::apply_boost(boost, hf_lab);
0179
0180         pxsum += hf_boosted.Px();
0181         pysum += hf_boosted.Py();
0182         pzsum += hf_boosted.Pz();
0183         Esum += hf_boosted.E();
0184       }
0185     }
0186
0187     // DIS kinematics calculations
0188     auto sigma_h = Esum - pzsum;
0189     auto ptsum = sqrt(pxsum*pxsum + pysum*pysum);
0190
0191     // Sigma zero or negative
0192     if (sigma_h <= 0) {
0193       if (msgLevel(MSG::DEBUG)) {
0194         debug() << "Sigma zero or negative" << endmsg;
0195       }
0196       return StatusCode::SUCCESS;
0197     }
0198
0199     // Calculate kinematic variables
0200     const auto y_jb = sigma_h / (2.*ei.energy());
0201     const auto Q2_jb = ptsum*ptsum / (1. - y_jb);
0202     const auto x_jb = Q2_jb / (4.*ei.energy()*pi.energy()*y_jb);
0203     const auto nu_jb = Q2_jb / (2.*m_proton*x_jb);
0204     const auto W_jb = sqrt(m_proton*m_proton + 2*m_proton*nu_jb - Q2_jb);
0205     auto kin = out_kinematics.create(x_jb, Q2_jb, W_jb, y_jb, nu_jb);
0206     kin.setScat(ef_rc);
0207
0208     // Debugging output
0209     if (msgLevel(MSG::DEBUG)) {
0210       debug() << "pi = " << pi << endmsg;
0211       debug() << "ei = " << ei << endmsg;
0212       debug() << "x,Q2,W,y,nu = "
0213               << kin.getX() << ","
0214               << kin.getQ2() << ","
0215               << kin.getW() << ","
0216               << kin.getY() << ","
0217               << kin.getNu()
0218               << endmsg;
0219     }
0220
0221     return StatusCode::SUCCESS;
0222   }
0223 };
0224
0225 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0226 DECLARE_COMPONENT(InclusiveKinematicsJB)
0227
0228 } // namespace Jug::Reco