truth/src/MC2SmearedParticle.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Sylvester Joosten, Whitney Armstrong, Wouter Deconinck
0003
0004 #include <algorithms/truth/MC2SmearedParticle.h>
0005
0006 #include <cmath>
0007 #include <edm4hep/utils/vector_utils.h>
0008
0009 namespace algorithms::truth {
0010
0011 void MC2SmearedParticle::init() {
0012   ; // do nothing
0013 }
0014
0015 void MC2SmearedParticle::process(const MC2SmearedParticle::Input& input,
0016                                  const MC2SmearedParticle::Output& output) const {
0017   const auto [parts] = input;
0018   auto [out_parts]   = output;
0019
0020   for (const auto& p : *parts) {
0021     if (p.getGeneratorStatus() > 1) {
0022       if (aboveDebugThreshold()) {
0023         debug() << "ignoring particle with generatorStatus = " << p.getGeneratorStatus() << endmsg;
0024       }
0025       continue;
0026     }
0027
0028     // for now just use total momentum smearing as this is the largest effect,
0029     // ideally we should also smear the angles but this should be good enough
0030     // for now.
0031     const auto pvec     = p.getMomentum();
0032     const auto pgen     = std::hypot(pvec.x, pvec.y, pvec.z);
0033     const auto momentum = pgen * m_rng.gaussian<double>(0., m_smearing);
0034     // make sure we keep energy consistent
0035     using MomType = decltype(edm4eic::ReconstructedParticle().getMomentum().x);
0036     const MomType energy =
0037         std::sqrt(p.getEnergy() * p.getEnergy() - pgen * pgen + momentum * momentum);
0038     const MomType px = p.getMomentum().x * momentum / pgen;
0039     const MomType py = p.getMomentum().y * momentum / pgen;
0040     const MomType pz = p.getMomentum().z * momentum / pgen;
0041
0042     const MomType dpx = m_smearing * px;
0043     const MomType dpy = m_smearing * py;
0044     const MomType dpz = m_smearing * pz;
0045     const MomType dE  = m_smearing * energy;
0046     // ignore covariance for now
0047     // @TODO: vertex smearing
0048     const MomType vx = p.getVertex().x;
0049     const MomType vy = p.getVertex().y;
0050     const MomType vz = p.getVertex().z;
0051
0052     auto rec_part = out_parts->create();
0053     rec_part.setType(-1); // @TODO: determine type codes
0054     rec_part.setEnergy(energy);
0055     rec_part.setMomentum({px, py, pz});
0056     rec_part.setReferencePoint({vx, vy, vz}); // @FIXME: probably not what we want?
0057     rec_part.setCharge(p.getCharge());
0058     rec_part.setMass(p.getMass());
0059     rec_part.setGoodnessOfPID(1); // Perfect PID
0060     rec_part.setCovMatrix({dpx, dpy, dpz, dE});
0061     rec_part.setPDG(p.getPDG());
0062   }
0063 }
0064
0065 } // namespace algorithms::truth
0066