src/components/InclusiveKinematicseSigma.cpp

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