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
0020 #include "Beam.h"
0021 #include "Boost.h"
0022 #include "HadronicFinalState.h"
0023
0024 using ROOT::Math::PxPyPzEVector;
0025
0026 namespace eicrecon {
0027
0028 void HadronicFinalState::init() {}
0029
0030 void HadronicFinalState::process(const HadronicFinalState::Input& input,
0031                                  const HadronicFinalState::Output& output) const {
0032
0033   const auto [mcparts, rcparts, rcassoc] = input;
0034   auto [hadronicfinalstate]              = output;
0035
0036   // Get incoming electron beam
0037   const auto ei_coll = find_first_beam_electron(mcparts);
0038   if (ei_coll.empty()) {
0039     debug("No beam electron found");
0040     return;
0041   }
0042   const PxPyPzEVector ei(round_beam_four_momentum(ei_coll[0].getMomentum(),
0043                                                   m_particleSvc.particle(ei_coll[0].getPDG()).mass,
0044                                                   {-5.0, -10.0, -18.0}, 0.0));
0045
0046   // Get incoming hadron beam
0047   const auto pi_coll = find_first_beam_hadron(mcparts);
0048   if (pi_coll.empty()) {
0049     debug("No beam hadron found");
0050     return;
0051   }
0052   const PxPyPzEVector pi(round_beam_four_momentum(pi_coll[0].getMomentum(),
0053                                                   m_particleSvc.particle(pi_coll[0].getPDG()).mass,
0054                                                   {41.0, 100.0, 275.0}, m_crossingAngle));
0055
0056   // Get first scattered electron
0057   const auto ef_coll = find_first_scattered_electron(mcparts);
0058   if (ef_coll.empty()) {
0059     debug("No truth scattered electron found");
0060     return;
0061   }
0062   // Associate first scattered electron with reconstructed electrons
0063   //const auto ef_assoc = std::find_if(
0064   //  rcassoc->begin(),
0065   //  rcassoc->end(),
0066   //  [&ef_coll](const auto& a){ return a.getSim().getObjectID() == ef_coll[0].getObjectID(); });
0067   auto ef_assoc = rcassoc->begin();
0068   for (; ef_assoc != rcassoc->end(); ++ef_assoc) {
0069     if (ef_assoc->getSim().getObjectID() == ef_coll[0].getObjectID()) {
0070       break;
0071     }
0072   }
0073   if (!(ef_assoc != rcassoc->end())) {
0074     debug("Truth scattered electron not in reconstructed particles");
0075     return;
0076   }
0077   const auto ef_rc{ef_assoc->getRec()};
0078   const auto ef_rc_id{ef_rc.getObjectID().index};
0079
0080   // Sums in colinear frame
0081   double pxsum = 0;
0082   double pysum = 0;
0083   double pzsum = 0;
0084   double Esum  = 0;
0085
0086   // Get boost to colinear frame
0087   auto boost = determine_boost(ei, pi);
0088
0089   auto hfs = hadronicfinalstate->create(0., 0., 0.);
0090
0091   for (const auto& p : *rcparts) {
0092     // Check if it's the scattered electron
0093     if (p.getObjectID().index != ef_rc_id) {
0094       // Lorentz vector in lab frame
0095       PxPyPzEVector hf_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy());
0096       // Boost to colinear frame
0097       PxPyPzEVector hf_boosted = apply_boost(boost, hf_lab);
0098
0099       pxsum += hf_boosted.Px();
0100       pysum += hf_boosted.Py();
0101       pzsum += hf_boosted.Pz();
0102       Esum += hf_boosted.E();
0103
0104       hfs.addToHadrons(p);
0105     }
0106   }
0107
0108   // Hadronic final state calculations
0109   auto sigma = Esum - pzsum;
0110   auto pT    = sqrt(pxsum * pxsum + pysum * pysum);
0111   auto gamma = acos((pT * pT - sigma * sigma) / (pT * pT + sigma * sigma));
0112
0113   hfs.setSigma(sigma);
0114   hfs.setPT(pT);
0115   hfs.setGamma(gamma);
0116
0117   // Sigma zero or negative
0118   if (sigma <= 0) {
0119     debug("Sigma zero or negative");
0120     return;
0121   }
0122
0123   debug("sigma_h, pT_h, gamma_h = {},{},{}", hfs.getSigma(), hfs.getPT(), hfs.getGamma());
0124 }
0125
0126 } // namespace eicrecon
0127 #endif