File indexing completed on 2024-09-28 07:03:48
0001
0002
0003
0004 #include "GaudiAlg/GaudiAlgorithm.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 <k4FWCore/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
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
0084 const auto& mcparts = *(m_inputMCParticleCollection.get());
0085 const auto& rcparts = *(m_inputParticleCollection.get());
0086 const auto& rcassoc = *(m_inputParticleAssociation.get());
0087
0088 auto& out_kinematics = *(m_outputInclusiveKinematicsCollection.createAndPut());
0089
0090
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
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
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
0131
0132
0133
0134
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
0151
0152
0153
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163
0164 double pxsum = 0;
0165 double pysum = 0;
0166 double pzsum = 0;
0167 double Esum = 0;
0168
0169
0170 auto boost = Jug::Base::Boost::determine_boost(ei, pi);
0171
0172 for (const auto& p: rcparts) {
0173
0174 if (p.getObjectID().index == ef_rc_id) {
0175
0176
0177 } else {
0178
0179 PxPyPzEVector hf_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy());
0180
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
0191 auto sigma_h = Esum - pzsum;
0192 auto ptsum = sqrt(pxsum*pxsum + pysum*pysum);
0193
0194
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
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
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
0229 DECLARE_COMPONENT(InclusiveKinematicsJB)
0230
0231 }