src/components/InclusiveKinematicsDA.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/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 InclusiveKinematicsDA : 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   InclusiveKinematicsDA(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, "InclusiveKinematicsDA");
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 pxsum = 0;
0161     double pysum = 0;
0162     double pzsum = 0;
0163     double Esum = 0;
0164     double theta_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         theta_e = e_boosted.Theta();
0178
0179       // Sum over all particles other than scattered electron
0180       } else {
0181         // Lorentz vector in lab frame
0182         PxPyPzEVector hf_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy());
0183         // Boost to colinear frame
0184         PxPyPzEVector hf_boosted = Jug::Base::Boost::apply_boost(boost, hf_lab);
0185
0186         pxsum += hf_boosted.Px();
0187         pysum += hf_boosted.Py();
0188         pzsum += hf_boosted.Pz();
0189         Esum += hf_boosted.E();
0190       }
0191     }
0192
0193     // DIS kinematics calculations
0194     auto sigma_h = Esum - pzsum;
0195     auto ptsum = sqrt(pxsum*pxsum + pysum*pysum);
0196     auto theta_h = 2.*atan(sigma_h/ptsum);
0197
0198     // If no scattered electron was found
0199     if (sigma_h <= 0) {
0200       if (msgLevel(MSG::DEBUG)) {
0201         debug() << "No scattered electron found" << endmsg;
0202       }
0203       return StatusCode::SUCCESS;
0204     }
0205
0206     // Calculate kinematic variables
0207     const auto y_da = tan(theta_h/2.) / ( tan(theta_e/2.) + tan(theta_h/2.) );
0208     const auto Q2_da = 4.*ei.energy()*ei.energy() * ( 1. / tan(theta_e/2.) ) * ( 1. / (tan(theta_e/2.) + tan(theta_h/2.)) );
0209     const auto x_da = Q2_da / (4.*ei.energy()*pi.energy()*y_da);
0210     const auto nu_da = Q2_da / (2.*m_proton*x_da);
0211     const auto W_da = sqrt(m_proton*m_proton + 2*m_proton*nu_da - Q2_da);
0212     auto kin = out_kinematics.create(x_da, Q2_da, W_da, y_da, nu_da);
0213     kin.setScat(ef_rc);
0214
0215     // Debugging output
0216     if (msgLevel(MSG::DEBUG)) {
0217       debug() << "pi = " << pi << endmsg;
0218       debug() << "ei = " << ei << endmsg;
0219       debug() << "x,Q2,W,y,nu = "
0220               << kin.getX() << ","
0221               << kin.getQ2() << ","
0222               << kin.getW() << ","
0223               << kin.getY() << ","
0224               << kin.getNu()
0225               << endmsg;
0226     }
0227
0228     return StatusCode::SUCCESS;
0229   }
0230 };
0231
0232 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0233 DECLARE_COMPONENT(InclusiveKinematicsDA)
0234
0235 } // namespace Jug::Reco