src/todo/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 "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 "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 InclusiveKinematicseSigma : 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   InclusiveKinematicseSigma(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, "InclusiveKinematicseSigma");
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
0168     double pt_e = 0;
0169     double sigma_e = 0;
0170
0171     // Get boost to colinear frame
0172     auto boost = Jug::Base::Boost::determine_boost(ei, pi);
0173
0174     for (const auto& p: rcparts) {
0175       // Get the scattered electron index and angle
0176       if (p.getObjectID().index == ef_rc_id) {
0177         // Lorentz vector in lab frame
0178         PxPyPzEVector e_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy());
0179         // Boost to colinear frame
0180         PxPyPzEVector e_boosted = Jug::Base::Boost::apply_boost(boost, e_lab);
0181
0182         pt_e = e_boosted.Pt();
0183         sigma_e = e_boosted.E() - e_boosted.Pz();
0184
0185       // Sum over all particles other than scattered electron
0186       } else{
0187         // Lorentz vector in lab frame
0188         PxPyPzEVector hf_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy());
0189         // Boost to colinear frame
0190         PxPyPzEVector hf_boosted = Jug::Base::Boost::apply_boost(boost, hf_lab);
0191
0192         pxsum += hf_boosted.Px();
0193         pysum += hf_boosted.Py();
0194         pzsum += hf_boosted.Pz();
0195         Esum += hf_boosted.E();
0196       }
0197     }
0198
0199     // DIS kinematics calculations
0200     auto sigma_h = Esum - pzsum;
0201     auto sigma_tot = sigma_e + sigma_h;
0202
0203     // If no scattered electron was found
0204     if (sigma_h <= 0) {
0205       if (msgLevel(MSG::DEBUG)) {
0206         debug() << "No scattered electron found" << endmsg;
0207       }
0208       return StatusCode::SUCCESS;
0209     }
0210
0211     // Calculate kinematic variables
0212     const auto y_e = 1. - sigma_e / (2.*ei.energy());
0213     const auto Q2_e = (pt_e*pt_e) / (1. - y_e);
0214
0215     const auto y_sig = sigma_h / sigma_tot;
0216     const auto Q2_sig = (pt_e*pt_e) / (1. - y_sig);
0217     const auto x_sig = Q2_sig / (4.*ei.energy()*pi.energy()*y_sig);
0218
0219     const auto Q2_esig = Q2_e;
0220     const auto x_esig = x_sig;
0221     const auto y_esig = Q2_esig / (4.*ei.energy()*pi.energy()*x_esig); //equivalent to (2*ei.energy() / sigma_tot)*y_sig
0222     const auto nu_esig = Q2_esig / (2.*m_proton*x_esig);
0223     const auto W_esig = sqrt(m_proton*m_proton + 2*m_proton*nu_esig - Q2_esig);
0224     auto kin = out_kinematics.create(x_esig, Q2_esig, W_esig, y_esig, nu_esig);
0225     kin.setScat(ef_rc);
0226
0227     // Debugging output
0228     if (msgLevel(MSG::DEBUG)) {
0229       debug() << "pi = " << pi << endmsg;
0230       debug() << "ei = " << ei << endmsg;
0231       debug() << "x,Q2,W,y,nu = "
0232               << kin.getX() << ","
0233               << kin.getQ2() << ","
0234               << kin.getW() << ","
0235               << kin.getY() << ","
0236               << kin.getNu()
0237               << endmsg;
0238     }
0239
0240     return StatusCode::SUCCESS;
0241   }
0242 };
0243
0244 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0245 DECLARE_COMPONENT(InclusiveKinematicseSigma)
0246
0247 } // namespace Jug::Reco