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