src/components/InclusiveKinematicsDA.cpp

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