Commit 1723867b authored by Wouter Deconinck's avatar Wouter Deconinck
Browse files

Implement InclusiveKinematicsTruth with updated functionality in Beam.h

parent 3199753a
...@@ -6,10 +6,10 @@ using ROOT::Math::PxPyPzEVector; ...@@ -6,10 +6,10 @@ using ROOT::Math::PxPyPzEVector;
#include "edm4hep/MCParticleCollection.h" #include "edm4hep/MCParticleCollection.h"
#include "eicd/ReconstructedParticleCollection.h" #include "eicd/ReconstructedParticleCollection.h"
namespace Jug::Reco::Beam { namespace Jug::Base::Beam {
template<class collection> template<class collection>
static auto find_first_with_pdg( auto find_first_with_pdg(
const collection& parts, const collection& parts,
const std::set<int32_t>& pdg) { const std::set<int32_t>& pdg) {
collection c; collection c;
...@@ -23,7 +23,7 @@ namespace Jug::Reco::Beam { ...@@ -23,7 +23,7 @@ namespace Jug::Reco::Beam {
} }
template<class collection> template<class collection>
static auto find_first_with_status_pdg( auto find_first_with_status_pdg(
const collection& parts, const collection& parts,
const std::set<int32_t>& status, const std::set<int32_t>& status,
const std::set<int32_t>& pdg) { const std::set<int32_t>& pdg) {
...@@ -37,23 +37,23 @@ namespace Jug::Reco::Beam { ...@@ -37,23 +37,23 @@ namespace Jug::Reco::Beam {
return c; return c;
} }
static auto find_first_beam_electron(const edm4hep::MCParticleCollection& mcparts) { inline auto find_first_beam_electron(const edm4hep::MCParticleCollection& mcparts) {
return find_first_with_status_pdg(mcparts, {4}, {11}); return find_first_with_status_pdg(mcparts, {4}, {11});
} }
static auto find_first_beam_hadron(const edm4hep::MCParticleCollection& mcparts) { inline auto find_first_beam_hadron(const edm4hep::MCParticleCollection& mcparts) {
return find_first_with_status_pdg(mcparts, {4}, {2212, 2112}); return find_first_with_status_pdg(mcparts, {4}, {2212, 2112});
} }
static auto find_first_scattered_electron(const edm4hep::MCParticleCollection& mcparts) { inline auto find_first_scattered_electron(const edm4hep::MCParticleCollection& mcparts) {
return find_first_with_status_pdg(mcparts, {1}, {11}); return find_first_with_status_pdg(mcparts, {1}, {11});
} }
static auto find_first_scattered_electron(const eicd::ReconstructedParticleCollection& rcparts) { inline auto find_first_scattered_electron(const eicd::ReconstructedParticleCollection& rcparts) {
return find_first_with_pdg(rcparts, {11}); return find_first_with_pdg(rcparts, {11});
} }
static inline
PxPyPzEVector PxPyPzEVector
round_beam_four_momentum( round_beam_four_momentum(
const edm4hep::Vector3f& p_in, const edm4hep::Vector3f& p_in,
...@@ -73,4 +73,4 @@ namespace Jug::Reco::Beam { ...@@ -73,4 +73,4 @@ namespace Jug::Reco::Beam {
return p_out; return p_out;
} }
} // namespace Jug::Reco::Beam } // namespace Jug::Base::Beam
...@@ -9,11 +9,11 @@ using ROOT::Math::PxPyPzEVector; ...@@ -9,11 +9,11 @@ using ROOT::Math::PxPyPzEVector;
#include "Math/RotationY.h" #include "Math/RotationY.h"
#include "Math/Boost.h" #include "Math/Boost.h"
namespace Jug::Reco::Boost { namespace Jug::Base::Boost {
using ROOT::Math::LorentzRotation; using ROOT::Math::LorentzRotation;
static LorentzRotation determine_boost(PxPyPzEVector ei, PxPyPzEVector pi) { inline LorentzRotation determine_boost(PxPyPzEVector ei, PxPyPzEVector pi) {
using ROOT::Math::RotationX; using ROOT::Math::RotationX;
using ROOT::Math::RotationY; using ROOT::Math::RotationY;
...@@ -47,7 +47,7 @@ namespace Jug::Reco::Boost { ...@@ -47,7 +47,7 @@ namespace Jug::Reco::Boost {
return tf; return tf;
} }
static PxPyPzEVector apply_boost(const LorentzRotation& tf, PxPyPzEVector part) { inline PxPyPzEVector apply_boost(const LorentzRotation& tf, PxPyPzEVector part) {
// Step 2: Apply boosts and rotations to any particle 4-vector // Step 2: Apply boosts and rotations to any particle 4-vector
// (here too, choices will have to be made as to what the 4-vector is for reconstructed particles) // (here too, choices will have to be made as to what the 4-vector is for reconstructed particles)
...@@ -57,4 +57,4 @@ namespace Jug::Reco::Boost { ...@@ -57,4 +57,4 @@ namespace Jug::Reco::Boost {
return part; return part;
} }
} // namespace Jug::Reco::Boost } // namespace Jug::Base::Boost
...@@ -3,35 +3,47 @@ ...@@ -3,35 +3,47 @@
#include "GaudiAlg/Producer.h" #include "GaudiAlg/Producer.h"
#include "GaudiAlg/Transformer.h" #include "GaudiAlg/Transformer.h"
#include "GaudiKernel/RndmGenerators.h" #include "GaudiKernel/RndmGenerators.h"
#include "GaudiKernel/PhysicalConstants.h"
#include <algorithm> #include <algorithm>
#include <ranges>
#include <cmath> #include <cmath>
#include "JugBase/IParticleSvc.h" #include "JugBase/IParticleSvc.h"
#include "JugBase/DataHandle.h" #include "JugBase/DataHandle.h"
#include <CLHEP/Vector/LorentzVector.h> #include "JugBase/Utilities/Beam.h"
#include "Math/Vector4D.h"
using ROOT::Math::PxPyPzEVector;
// Event Model related classes // Event Model related classes
#include "edm4hep/MCParticleCollection.h" #include "edm4hep/MCParticleCollection.h"
#include "eicd/InclusiveKinematicsCollection.h" #include "eicd/InclusiveKinematicsCollection.h"
#include "eicd/vector_utils.h"
namespace Jug::Fast { namespace Jug::Fast {
class InclusiveKinematicsTruth : public GaudiAlgorithm { class InclusiveKinematicsTruth : public GaudiAlgorithm {
private: private:
DataHandle<edm4hep::MCParticleCollection> m_inputParticleCollection{"MCParticles", Gaudi::DataHandle::Reader, DataHandle<edm4hep::MCParticleCollection> m_inputMCParticleCollection{
this}; "inputMCParticles",
DataHandle<eicd::InclusiveKinematicsCollection> m_outputInclusiveKinematicsCollection{"InclusiveKinematicsTruth", Gaudi::DataHandle::Reader,
Gaudi::DataHandle::Writer, this}; this};
DataHandle<eicd::InclusiveKinematicsCollection> m_outputInclusiveKinematicsCollection{
"InclusiveKinematicsTruth",
Gaudi::DataHandle::Writer,
this};
SmartIF<IParticleSvc> m_pidSvc; SmartIF<IParticleSvc> m_pidSvc;
double m_proton{0}; double m_proton{0};
double m_neutron{0}; double m_neutron{0};
double m_electron{0};
public: public:
InclusiveKinematicsTruth(const std::string& name, ISvcLocator* svcLoc) InclusiveKinematicsTruth(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc) { : GaudiAlgorithm(name, svcLoc) {
declareProperty("inputMCParticles", m_inputParticleCollection, "MCParticles"); declareProperty("inputMCParticles", m_inputMCParticleCollection, "MCParticles");
declareProperty("outputData", m_outputInclusiveKinematicsCollection, "InclusiveKinematicsTruth"); declareProperty("outputData", m_outputInclusiveKinematicsCollection, "InclusiveKinematicsTruth");
} }
...@@ -49,13 +61,14 @@ public: ...@@ -49,13 +61,14 @@ public:
} }
m_proton = m_pidSvc->particle(2212).mass; m_proton = m_pidSvc->particle(2212).mass;
m_neutron = m_pidSvc->particle(2112).mass; m_neutron = m_pidSvc->particle(2112).mass;
m_electron = m_pidSvc->particle(11).mass;
return StatusCode::SUCCESS; return StatusCode::SUCCESS;
} }
StatusCode execute() override { StatusCode execute() override {
// input collection // input collections
const edm4hep::MCParticleCollection& parts = *(m_inputParticleCollection.get()); const auto& mcparts = *(m_inputMCParticleCollection.get());
// output collection // output collection
auto& out_kinematics = *(m_outputInclusiveKinematicsCollection.createAndPut()); auto& out_kinematics = *(m_outputInclusiveKinematicsCollection.createAndPut());
...@@ -64,85 +77,59 @@ public: ...@@ -64,85 +77,59 @@ public:
// or outgoing electron line, the vertex kinematics will be different than the // or outgoing electron line, the vertex kinematics will be different than the
// kinematics calculated using the scattered electron as done here. // kinematics calculated using the scattered electron as done here.
// Also need to update for CC events. // Also need to update for CC events.
auto ei = CLHEP::HepLorentzVector(0., 0., 0., 0.);
auto pi = CLHEP::HepLorentzVector(0., 0., 0., 0.);
auto ef = CLHEP::HepLorentzVector(0., 0., 0., 0.);
bool ebeam_found = false;
bool pbeam_found = false;
int32_t mcscatID = -1;
for (const auto& p : parts) {
if (p.getGeneratorStatus() == 4 && p.getPDG() == 11) { // Incoming electron
ei.setPx(p.getMomentum().x);
ei.setPy(p.getMomentum().y);
ei.setPz(p.getMomentum().z);
ei.setE(p.getEnergy());
ebeam_found = true;
}
else if (p.getGeneratorStatus() == 4 && p.getPDG() == 2212) { // Incoming proton
pi.setPx(p.getMomentum().x);
pi.setPy(p.getMomentum().y);
pi.setPz(p.getMomentum().z);
pi.setE(p.getEnergy());
pbeam_found = true;
}
else if (p.getGeneratorStatus() == 4 && p.getPDG() == 2112) { // Incoming neutron
pi.setPx(p.getMomentum().x);
pi.setPy(p.getMomentum().y);
pi.setPz(p.getMomentum().z);
pi.setE(p.getEnergy());
pbeam_found = true;
}
// Scattered electron. Currently taken as first status==1 electron in HEPMC record,
// which seems to be correct based on a cursory glance at the Pythia8 output. In the future,
// it may be better to trace back each final-state electron and see which one originates from
// the beam.
else if (p.getGeneratorStatus() == 1 && p.getPDG() == 11 && mcscatID == -1) {
ef.setPx(p.getMomentum().x);
ef.setPy(p.getMomentum().y);
ef.setPz(p.getMomentum().z);
ef.setE(p.getEnergy());
mcscatID = p.id();
}
if (ebeam_found && pbeam_found && mcscatID != -1) {
// all done!
break;
}
}
// Not all particles found // Get incoming electron beam
if (ebeam_found == false) { const auto ei_coll = Jug::Base::Beam::find_first_beam_electron(mcparts);
if (ei_coll.size() == 0) {
if (msgLevel(MSG::DEBUG)) { if (msgLevel(MSG::DEBUG)) {
debug() << "No initial electron found" << endmsg; debug() << "No beam electron found" << endmsg;
} }
return StatusCode::SUCCESS; return StatusCode::SUCCESS;
} }
if (pbeam_found == false) { const auto ei_p = ei_coll[0].getMomentum();
const auto ei_p_mag = eicd::magnitude(ei_p);
const auto ei_mass = m_electron;
const PxPyPzEVector ei(ei_p.x, ei_p.y, ei_p.z, std::hypot(ei_p_mag, ei_mass));
// Get incoming hadron beam
const auto pi_coll = Jug::Base::Beam::find_first_beam_hadron(mcparts);
if (pi_coll.size() == 0) {
if (msgLevel(MSG::DEBUG)) { if (msgLevel(MSG::DEBUG)) {
debug() << "No initial proton found" << endmsg; debug() << "No beam hadron found" << endmsg;
} }
return StatusCode::SUCCESS; return StatusCode::SUCCESS;
} }
if (mcscatID == -1) { const auto pi_p = pi_coll[0].getMomentum();
const auto pi_p_mag = eicd::magnitude(pi_p);
const auto pi_mass = pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron;
const PxPyPzEVector pi(pi_p.x, pi_p.y, pi_p.z, std::hypot(pi_p_mag, pi_mass));
// Get first scattered electron
// Scattered electron. Currently taken as first status==1 electron in HEPMC record,
// which seems to be correct based on a cursory glance at the Pythia8 output. In the future,
// it may be better to trace back each final-state electron and see which one originates from
// the beam.
const auto ef_coll = Jug::Base::Beam::find_first_scattered_electron(mcparts);
if (ef_coll.size() == 0) {
if (msgLevel(MSG::DEBUG)) { if (msgLevel(MSG::DEBUG)) {
debug() << "No scattered electron found" << endmsg; debug() << "No truth scattered electron found" << endmsg;
} }
return StatusCode::SUCCESS; return StatusCode::SUCCESS;
} }
const auto ef_p = ef_coll[0].getMomentum();
const auto ef_p_mag = eicd::magnitude(ef_p);
const auto ef_mass = m_electron;
const PxPyPzEVector ef(ef_p.x, ef_p.y, ef_p.z, std::hypot(ef_p_mag, ef_mass));
// DIS kinematics calculations // DIS kinematics calculations
auto kin = out_kinematics.create();
const auto q = ei - ef; const auto q = ei - ef;
kin.setQ2(-1. * q.m2()); const auto q_dot_pi = q.Dot(pi);
kin.setY((q * pi) / (ei * pi)); const auto Q2 = -q.Dot(q);
kin.setNu(q * pi / m_proton); const auto y = q_dot_pi / ei.Dot(pi);
kin.setX(kin.getQ2() / (2. * q * pi)); const auto nu = q_dot_pi / pi_mass;
kin.setW(sqrt((pi + q).m2())); const auto x = Q2 / (2.*q_dot_pi);
// kin.scat(????); @TODO: this is now set as a OneToOneRelation to ReconstructedParticle, const auto W = sqrt(pi_mass*pi_mass + 2.*q_dot_pi - Q2);
// which breaks for this algorithm that takes raw MCParticles const auto kin = out_kinematics.create(x, Q2, W, y, nu);
// Debugging output // Debugging output
if (msgLevel(MSG::DEBUG)) { if (msgLevel(MSG::DEBUG)) {
...@@ -150,11 +137,11 @@ public: ...@@ -150,11 +137,11 @@ public:
debug() << "ei = " << ei << endmsg; debug() << "ei = " << ei << endmsg;
debug() << "ef = " << ef << endmsg; debug() << "ef = " << ef << endmsg;
debug() << "q = " << q << endmsg; debug() << "q = " << q << endmsg;
debug() << "x,y,Q2,W,nu = " debug() << "x,Q2,W,y,nu = "
<< kin.getX() << "," << kin.getX() << ","
<< kin.getY() << ","
<< kin.getQ2() << "," << kin.getQ2() << ","
<< kin.getW() << "," << kin.getW() << ","
<< kin.getY() << ","
<< kin.getNu() << kin.getNu()
<< endmsg; << endmsg;
} }
......
...@@ -10,8 +10,8 @@ ...@@ -10,8 +10,8 @@
#include "JugBase/IParticleSvc.h" #include "JugBase/IParticleSvc.h"
#include "JugBase/DataHandle.h" #include "JugBase/DataHandle.h"
#include "JugReco/Utilities/Beam.h" #include "JugBase/Utilities/Beam.h"
#include "JugReco/Utilities/Boost.h" #include "JugBase/Utilities/Boost.h"
#include "Math/Vector4D.h" #include "Math/Vector4D.h"
using ROOT::Math::PxPyPzEVector; using ROOT::Math::PxPyPzEVector;
...@@ -84,7 +84,7 @@ public: ...@@ -84,7 +84,7 @@ public:
auto& out_kinematics = *(m_outputInclusiveKinematicsCollection.createAndPut()); auto& out_kinematics = *(m_outputInclusiveKinematicsCollection.createAndPut());
// Get incoming electron beam // Get incoming electron beam
const auto ei_coll = Jug::Reco::Beam::find_first_beam_electron(mcparts); const auto ei_coll = Jug::Base::Beam::find_first_beam_electron(mcparts);
if (ei_coll.size() == 0) { if (ei_coll.size() == 0) {
if (msgLevel(MSG::DEBUG)) { if (msgLevel(MSG::DEBUG)) {
debug() << "No beam electron found" << endmsg; debug() << "No beam electron found" << endmsg;
...@@ -92,7 +92,7 @@ public: ...@@ -92,7 +92,7 @@ public:
return StatusCode::SUCCESS; return StatusCode::SUCCESS;
} }
const PxPyPzEVector ei( const PxPyPzEVector ei(
Jug::Reco::Beam::round_beam_four_momentum( Jug::Base::Beam::round_beam_four_momentum(
ei_coll[0].getMomentum(), ei_coll[0].getMomentum(),
m_electron, m_electron,
{-5.0, -10.0, -18.0}, {-5.0, -10.0, -18.0},
...@@ -100,7 +100,7 @@ public: ...@@ -100,7 +100,7 @@ public:
); );
// Get incoming hadron beam // Get incoming hadron beam
const auto pi_coll = Jug::Reco::Beam::find_first_beam_hadron(mcparts); const auto pi_coll = Jug::Base::Beam::find_first_beam_hadron(mcparts);
if (pi_coll.size() == 0) { if (pi_coll.size() == 0) {
if (msgLevel(MSG::DEBUG)) { if (msgLevel(MSG::DEBUG)) {
debug() << "No beam hadron found" << endmsg; debug() << "No beam hadron found" << endmsg;
...@@ -108,7 +108,7 @@ public: ...@@ -108,7 +108,7 @@ public:
return StatusCode::SUCCESS; return StatusCode::SUCCESS;
} }
const PxPyPzEVector pi( const PxPyPzEVector pi(
Jug::Reco::Beam::round_beam_four_momentum( Jug::Base::Beam::round_beam_four_momentum(
pi_coll[0].getMomentum(), pi_coll[0].getMomentum(),
pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron, pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron,
{41.0, 100.0, 275.0}, {41.0, 100.0, 275.0},
...@@ -116,7 +116,7 @@ public: ...@@ -116,7 +116,7 @@ public:
); );
// Get first scattered electron // Get first scattered electron
const auto ef_coll = Jug::Reco::Beam::find_first_scattered_electron(mcparts); const auto ef_coll = Jug::Base::Beam::find_first_scattered_electron(mcparts);
if (ef_coll.size() == 0) { if (ef_coll.size() == 0) {
if (msgLevel(MSG::DEBUG)) { if (msgLevel(MSG::DEBUG)) {
debug() << "No truth scattered electron found" << endmsg; debug() << "No truth scattered electron found" << endmsg;
...@@ -164,7 +164,7 @@ public: ...@@ -164,7 +164,7 @@ public:
double theta_e = 0; double theta_e = 0;
// Get boost to colinear frame // Get boost to colinear frame
auto boost = Jug::Reco::Boost::determine_boost(ei, pi); auto boost = Jug::Base::Boost::determine_boost(ei, pi);
for(const auto& p : rcparts) { for(const auto& p : rcparts) {
// Get the scattered electron index and angle // Get the scattered electron index and angle
...@@ -172,7 +172,7 @@ public: ...@@ -172,7 +172,7 @@ public:
// Lorentz vector in lab frame // Lorentz vector in lab frame
PxPyPzEVector e_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); PxPyPzEVector e_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy());
// Boost to colinear frame // Boost to colinear frame
PxPyPzEVector e_boosted = Jug::Reco::Boost::apply_boost(boost, e_lab); PxPyPzEVector e_boosted = Jug::Base::Boost::apply_boost(boost, e_lab);
theta_e = e_boosted.Theta(); theta_e = e_boosted.Theta();
...@@ -181,7 +181,7 @@ public: ...@@ -181,7 +181,7 @@ public:
// Lorentz vector in lab frame // Lorentz vector in lab frame
PxPyPzEVector hf_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); PxPyPzEVector hf_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy());
// Boost to colinear frame // Boost to colinear frame
PxPyPzEVector hf_boosted = Jug::Reco::Boost::apply_boost(boost, hf_lab); PxPyPzEVector hf_boosted = Jug::Base::Boost::apply_boost(boost, hf_lab);
pxsum += hf_boosted.Px(); pxsum += hf_boosted.Px();
pysum += hf_boosted.Py(); pysum += hf_boosted.Py();
......
...@@ -11,8 +11,8 @@ ...@@ -11,8 +11,8 @@
#include "JugBase/IParticleSvc.h" #include "JugBase/IParticleSvc.h"
#include "JugBase/DataHandle.h" #include "JugBase/DataHandle.h"
#include "JugReco/Utilities/Beam.h" #include "JugBase/Utilities/Beam.h"
#include "JugReco/Utilities/Boost.h" #include "JugBase/Utilities/Boost.h"
#include "Math/Vector4D.h" #include "Math/Vector4D.h"
using ROOT::Math::PxPyPzEVector; using ROOT::Math::PxPyPzEVector;
...@@ -133,7 +133,7 @@ public: ...@@ -133,7 +133,7 @@ public:
//} //}
// Get incoming electron beam // Get incoming electron beam
const auto ei_coll = Jug::Reco::Beam::find_first_beam_electron(mcparts); const auto ei_coll = Jug::Base::Beam::find_first_beam_electron(mcparts);
if (ei_coll.size() == 0) { if (ei_coll.size() == 0) {
if (msgLevel(MSG::DEBUG)) { if (msgLevel(MSG::DEBUG)) {
debug() << "No beam electron found" << endmsg; debug() << "No beam electron found" << endmsg;
...@@ -141,7 +141,7 @@ public: ...@@ -141,7 +141,7 @@ public:
return StatusCode::SUCCESS; return StatusCode::SUCCESS;
} }
const PxPyPzEVector ei( const PxPyPzEVector ei(
Jug::Reco::Beam::round_beam_four_momentum( Jug::Base::Beam::round_beam_four_momentum(
ei_coll[0].getMomentum(), ei_coll[0].getMomentum(),
m_electron, m_electron,
{-5.0, -10.0, -18.0}, {-5.0, -10.0, -18.0},
...@@ -149,7 +149,7 @@ public: ...@@ -149,7 +149,7 @@ public:
); );
// Get incoming hadron beam // Get incoming hadron beam
const auto pi_coll = Jug::Reco::Beam::find_first_beam_hadron(mcparts); const auto pi_coll = Jug::Base::Beam::find_first_beam_hadron(mcparts);
if (pi_coll.size() == 0) { if (pi_coll.size() == 0) {
if (msgLevel(MSG::DEBUG)) { if (msgLevel(MSG::DEBUG)) {
debug() << "No beam hadron found" << endmsg; debug() << "No beam hadron found" << endmsg;
...@@ -157,7 +157,7 @@ public: ...@@ -157,7 +157,7 @@ public:
return StatusCode::SUCCESS; return StatusCode::SUCCESS;
} }
const PxPyPzEVector pi( const PxPyPzEVector pi(
Jug::Reco::Beam::round_beam_four_momentum( Jug::Base::Beam::round_beam_four_momentum(
pi_coll[0].getMomentum(), pi_coll[0].getMomentum(),
pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron, pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron,
{41.0, 100.0, 275.0}, {41.0, 100.0, 275.0},
...@@ -165,7 +165,7 @@ public: ...@@ -165,7 +165,7 @@ public:
); );
// Get first scattered electron // Get first scattered electron
const auto ef_coll = Jug::Reco::Beam::find_first_scattered_electron(mcparts); const auto ef_coll = Jug::Base::Beam::find_first_scattered_electron(mcparts);
if (ef_coll.size() == 0) { if (ef_coll.size() == 0) {
if (msgLevel(MSG