src/components/InclusiveKinematicsSigma.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Wouter Deconinck, Barak Schmookler
0003
0004 #include "GaudiAlg/GaudiAlgorithm.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 <k4FWCore/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 InclusiveKinematicsSigma : 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   InclusiveKinematicsSigma(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, "InclusiveKinematicsSigma");
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     // Get incoming electron beam
0090     const auto ei_coll = Jug::Base::Beam::find_first_beam_electron(mcparts);
0091     if (ei_coll.size() == 0) {
0092       if (msgLevel(MSG::DEBUG)) {
0093         debug() << "No beam electron found" << endmsg;
0094       }
0095       return StatusCode::SUCCESS;
0096     }
0097     const PxPyPzEVector ei(
0098       Jug::Base::Beam::round_beam_four_momentum(
0099         ei_coll[0].getMomentum(),
0100         m_electron,
0101         {-5.0, -10.0, -18.0},
0102         0.0)
0103       );
0104
0105     // Get incoming hadron beam
0106     const auto pi_coll = Jug::Base::Beam::find_first_beam_hadron(mcparts);
0107     if (pi_coll.size() == 0) {
0108       if (msgLevel(MSG::DEBUG)) {
0109         debug() << "No beam hadron found" << endmsg;
0110       }
0111       return StatusCode::SUCCESS;
0112     }
0113     const PxPyPzEVector pi(
0114       Jug::Base::Beam::round_beam_four_momentum(
0115         pi_coll[0].getMomentum(),
0116         pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron,
0117         {41.0, 100.0, 275.0},
0118         m_crossingAngle)
0119       );
0120
0121     // Get first scattered electron
0122     const auto ef_coll = Jug::Base::Beam::find_first_scattered_electron(mcparts);
0123     if (ef_coll.size() == 0) {
0124       if (msgLevel(MSG::DEBUG)) {
0125         debug() << "No truth scattered electron found" << endmsg;
0126       }
0127       return StatusCode::SUCCESS;
0128     }
0129     // Associate first scattered electron with reconstructed electrons
0130     //const auto ef_assoc = std::find_if(
0131     //  rcassoc.begin(),
0132     //  rcassoc.end(),
0133     //  [&ef_coll](const auto& a){ return a.getSimID() == ef_coll[0].getObjectID().index; });
0134     auto ef_assoc = rcassoc.begin();
0135     for (; ef_assoc != rcassoc.end(); ++ef_assoc) {
0136       if (ef_assoc->getSimID() == (unsigned) ef_coll[0].getObjectID().index) {
0137         break;
0138       }
0139     }
0140     if (!(ef_assoc != rcassoc.end())) {
0141       if (msgLevel(MSG::DEBUG)) {
0142         debug() << "Truth scattered electron not in reconstructed particles" << endmsg;
0143       }
0144       return StatusCode::SUCCESS;
0145     }
0146     const auto ef_rc{ef_assoc->getRec()};
0147     const auto ef_rc_id{ef_rc.getObjectID().index};
0148
0149     // Loop over reconstructed particles to get all outgoing particles
0150     // -----------------------------------------------------------------
0151     // Right now, everything is taken from Reconstructed particles branches.
0152     //
0153     // This means the tracking detector is used for charged particles to caculate the momentum,
0154     // and the magnitude of this momentum plus the true PID to calculate the energy.
0155     // No requirement is made that these particles produce a hit in any other detector
0156     //
0157     // Using the Reconstructed particles branches also means that the reconstruction for neutrals is done using the
0158     // calorimeter(s) information for the energy and angles, and then using this energy and the true PID to get the
0159     // magnitude of the momentum.
0160     // -----------------------------------------------------------------
0161
0162     // Sums in colinear frame
0163     double pzsum = 0;
0164     double Esum = 0;
0165
0166     double pt_e = 0;
0167     double sigma_e = 0;
0168
0169     // Get boost to colinear frame
0170     auto boost = Jug::Base::Boost::determine_boost(ei, pi);
0171
0172     for(const auto& p: rcparts) {
0173       // Get the scattered electron index and angle
0174       if (p.getObjectID().index == ef_rc_id) {
0175         // Lorentz vector in lab frame
0176         PxPyPzEVector e_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy());
0177         // Boost to colinear frame
0178         PxPyPzEVector e_boosted = Jug::Base::Boost::apply_boost(boost, e_lab);
0179
0180         pt_e = e_boosted.Pt();
0181         sigma_e = e_boosted.E() - e_boosted.Pz();
0182
0183       // Sum over all particles other than scattered electron
0184       } else{
0185         // Lorentz vector in lab frame
0186         PxPyPzEVector hf_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy());
0187         // Boost to colinear frame
0188         PxPyPzEVector hf_boosted = Jug::Base::Boost::apply_boost(boost, hf_lab);
0189         pzsum += hf_boosted.Pz();
0190         Esum += hf_boosted.E();
0191       }
0192     }
0193
0194     // DIS kinematics calculations
0195     auto sigma_h = Esum - pzsum;
0196     auto sigma_tot = sigma_e + sigma_h;
0197
0198     if (sigma_h <= 0) {
0199       if (msgLevel(MSG::DEBUG)) {
0200         debug() << "No scattered electron found or sigma zero or negative" << endmsg;
0201       }
0202       return StatusCode::SUCCESS;
0203     }
0204
0205     // Calculate kinematic variables
0206     const auto y_sig = sigma_h / sigma_tot;
0207     const auto Q2_sig = (pt_e*pt_e) / (1. - y_sig);
0208     const auto x_sig = Q2_sig / (4.*ei.energy()*pi.energy()*y_sig);
0209     const auto nu_sig = Q2_sig / (2.*m_proton*x_sig);
0210     const auto W_sig = sqrt(m_proton*m_proton + 2*m_proton*nu_sig - Q2_sig);
0211     auto kin = out_kinematics.create(x_sig, Q2_sig, W_sig, y_sig, nu_sig);
0212     kin.setScat(ef_rc);
0213
0214     // Debugging output
0215     if (msgLevel(MSG::DEBUG)) {
0216       debug() << "pi = " << pi << endmsg;
0217       debug() << "ei = " << ei << endmsg;
0218       debug() << "x,Q2,W,y,nu = "
0219               << kin.getX() << ","
0220               << kin.getQ2() << ","
0221               << kin.getW() << ","
0222               << kin.getY() << ","
0223               << kin.getNu()
0224               << endmsg;
0225     }
0226
0227     return StatusCode::SUCCESS;
0228   }
0229 };
0230
0231 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0232 DECLARE_COMPONENT(InclusiveKinematicsSigma)
0233
0234 } // namespace Jug::Reco