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