algorithms/reco/HadronicFinalState.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2024 Tyler Kutz
0003
0004 #include <edm4eic/EDM4eicVersion.h>
0005 #if EDM4EIC_VERSION_MAJOR >= 6
0006
0007 #include <Math/GenVector/LorentzVector.h>
0008 #include <Math/GenVector/PxPyPzE4D.h>
0009 #include <Math/Vector4Dfwd.h>
0010 #include <edm4eic/HadronicFinalStateCollection.h>
0011 #include <edm4eic/MCRecoParticleAssociationCollection.h>
0012 #include <edm4eic/ReconstructedParticleCollection.h>
0013 #include <edm4hep/MCParticleCollection.h>
0014 #include <edm4hep/Vector3f.h>
0015 #include <fmt/core.h>
0016 #include <podio/ObjectID.h>
0017 #include <cmath>
0018 #include <gsl/pointers>
0019 #include <vector>
0020
0021 #include "Beam.h"
0022 #include "Boost.h"
0023 #include "HadronicFinalState.h"
0024
0025 using ROOT::Math::PxPyPzEVector;
0026
0027 namespace eicrecon {
0028
0029   void HadronicFinalState::init(std::shared_ptr<spdlog::logger>& logger) {
0030     m_log = logger;
0031     // m_pidSvc = service("ParticleSvc");
0032     // if (!m_pidSvc) {
0033     //   m_log->debug("Unable to locate Particle Service. "
0034     //     "Make sure you have ParticleSvc in the configuration.");
0035     // }
0036   }
0037
0038   void HadronicFinalState::process(
0039       const HadronicFinalState::Input& input,
0040       const HadronicFinalState::Output& output) const {
0041
0042     const auto [mcparts, rcparts, rcassoc] = input;
0043     auto [hadronicfinalstate] = output;
0044
0045     // Get incoming electron beam
0046     const auto ei_coll = find_first_beam_electron(mcparts);
0047     if (ei_coll.size() == 0) {
0048       m_log->debug("No beam electron found");
0049       return;
0050     }
0051     const PxPyPzEVector ei(
0052       round_beam_four_momentum(
0053         ei_coll[0].getMomentum(),
0054         m_electron,
0055         {-5.0, -10.0, -18.0},
0056         0.0)
0057       );
0058
0059     // Get incoming hadron beam
0060     const auto pi_coll = find_first_beam_hadron(mcparts);
0061     if (pi_coll.size() == 0) {
0062       m_log->debug("No beam hadron found");
0063       return;
0064     }
0065     const PxPyPzEVector pi(
0066       round_beam_four_momentum(
0067         pi_coll[0].getMomentum(),
0068         pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron,
0069         {41.0, 100.0, 275.0},
0070         m_crossingAngle)
0071       );
0072
0073     // Get first scattered electron
0074     const auto ef_coll = find_first_scattered_electron(mcparts);
0075     if (ef_coll.size() == 0) {
0076       m_log->debug("No truth scattered electron found");
0077       return;
0078     }
0079     // Associate first scattered electron with reconstructed electrons
0080     //const auto ef_assoc = std::find_if(
0081     //  rcassoc->begin(),
0082     //  rcassoc->end(),
0083     //  [&ef_coll](const auto& a){ return a.getSim().getObjectID() == ef_coll[0].getObjectID(); });
0084     auto ef_assoc = rcassoc->begin();
0085     for (; ef_assoc != rcassoc->end(); ++ef_assoc) {
0086       if (ef_assoc->getSim().getObjectID() == ef_coll[0].getObjectID()) {
0087         break;
0088       }
0089     }
0090     if (!(ef_assoc != rcassoc->end())) {
0091       m_log->debug("Truth scattered electron not in reconstructed particles");
0092       return;
0093     }
0094     const auto ef_rc{ef_assoc->getRec()};
0095     const auto ef_rc_id{ef_rc.getObjectID().index};
0096
0097     // Sums in colinear frame
0098     double pxsum = 0;
0099     double pysum = 0;
0100     double pzsum = 0;
0101     double Esum = 0;
0102
0103     // Get boost to colinear frame
0104     auto boost = determine_boost(ei, pi);
0105
0106     auto hfs = hadronicfinalstate->create(0., 0., 0.);
0107
0108     for (const auto& p: *rcparts) {
0109
0110       bool isHadron = true;
0111       // Check if it's the scattered electron
0112       if (p.getObjectID().index == ef_rc_id) isHadron = false;
0113       // Check for non-hadron PDG codes
0114       if (p.getPDG() == 11) isHadron = false;
0115       if (p.getPDG() == 22) isHadron = false;
0116       if (p.getPDG() == 13) isHadron = false;
0117       // If it's the scattered electron or not a hadron, skip
0118       if(!isHadron) continue;
0119
0120       // Lorentz vector in lab frame
0121       PxPyPzEVector hf_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy());
0122       // Boost to colinear frame
0123       PxPyPzEVector hf_boosted = apply_boost(boost, hf_lab);
0124
0125       pxsum += hf_boosted.Px();
0126       pysum += hf_boosted.Py();
0127       pzsum += hf_boosted.Pz();
0128       Esum += hf_boosted.E();
0129
0130       hfs.addToHadrons(p);
0131
0132     }
0133
0134     // Hadronic final state calculations
0135     auto sigma = Esum - pzsum;
0136     auto pT = sqrt(pxsum*pxsum + pysum*pysum);
0137     auto gamma = (pT*pT - sigma*sigma)/(pT*pT + sigma*sigma);
0138
0139     hfs.setSigma(sigma);
0140     hfs.setPT(pT);
0141     hfs.setGamma(gamma);
0142
0143     // Sigma zero or negative
0144     if (sigma <= 0) {
0145       m_log->debug("Sigma zero or negative");
0146       return;
0147     }
0148
0149     m_log->debug("sigma_h, pT_h, gamma_h = {},{},{}", hfs.getSigma(), hfs.getPT(), hfs.getGamma());
0150
0151   }
0152
0153 } // namespace Jug::Reco
0154 #endif