Commit 3199753a authored by Wouter Deconinck's avatar Wouter Deconinck
Browse files

Reenable InclusiveKinematics algorithms

parent bd093139
......@@ -24,6 +24,11 @@ gaudi_add_module(JugRecoPlugins
src/components/ImagingPixelDataSorter.cpp
src/components/ImagingTopoCluster.cpp
src/components/ImagingClusterReco.cpp
src/components/InclusiveKinematicsElectron.cpp
src/components/InclusiveKinematicsDA.cpp
src/components/InclusiveKinematicsJB.cpp
src/components/InclusiveKinematicsSigma.cpp
src/components/InclusiveKinematicseSigma.cpp
src/components/PhotoMultiplierReco.cpp
LINK
Gaudi::GaudiAlgLib Gaudi::GaudiKernel
......
#pragma once
#include "Math/Vector4D.h"
using ROOT::Math::PxPyPzEVector;
#include "edm4hep/MCParticleCollection.h"
#include "eicd/ReconstructedParticleCollection.h"
namespace Jug::Reco::Beam {
template<class collection>
static auto find_first_with_pdg(
const collection& parts,
const std::set<int32_t>& pdg) {
collection c;
for (const auto& p: parts) {
if (pdg.count(p.getPDG()) > 0) {
c.push_back(p);
break;
}
}
return c;
}
template<class collection>
static auto find_first_with_status_pdg(
const collection& parts,
const std::set<int32_t>& status,
const std::set<int32_t>& pdg) {
collection c;
for (const auto& p: parts) {
if (status.count(p.getGeneratorStatus()) > 0 && pdg.count(p.getPDG()) > 0) {
c.push_back(p);
break;
}
}
return c;
}
static auto find_first_beam_electron(const edm4hep::MCParticleCollection& mcparts) {
return find_first_with_status_pdg(mcparts, {4}, {11});
}
static auto find_first_beam_hadron(const edm4hep::MCParticleCollection& mcparts) {
return find_first_with_status_pdg(mcparts, {4}, {2212, 2112});
}
static auto find_first_scattered_electron(const edm4hep::MCParticleCollection& mcparts) {
return find_first_with_status_pdg(mcparts, {1}, {11});
}
static auto find_first_scattered_electron(const eicd::ReconstructedParticleCollection& rcparts) {
return find_first_with_pdg(rcparts, {11});
}
static
PxPyPzEVector
round_beam_four_momentum(
const edm4hep::Vector3f& p_in,
const float mass,
const std::vector<float>& pz_set,
const float crossing_angle = 0.0) {
PxPyPzEVector p_out;
for (const auto& pz : pz_set) {
if (fabs(p_in.z / pz - 1) < 0.1) {
p_out.SetPz(pz);
break;
}
}
p_out.SetPx(p_out.Pz() * sin(crossing_angle));
p_out.SetPz(p_out.Pz() * cos(crossing_angle));
p_out.SetE(std::hypot(p_out.Px(), p_out.Pz(), mass));
return p_out;
}
} // namespace Jug::Reco::Beam
#pragma once
#include "Math/Vector4D.h"
using ROOT::Math::PxPyPzEVector;
#include "Math/LorentzRotation.h"
#include "Math/LorentzVector.h"
#include "Math/RotationX.h"
#include "Math/RotationY.h"
#include "Math/Boost.h"
namespace Jug::Reco::Boost {
using ROOT::Math::LorentzRotation;
static LorentzRotation determine_boost(PxPyPzEVector ei, PxPyPzEVector pi) {
using ROOT::Math::RotationX;
using ROOT::Math::RotationY;
using ROOT::Math::Boost;
// Step 1: Find the needed boosts and rotations from the incoming lepton and hadron beams
// (note, this will give you a perfect boost, in principle you will not know the beam momenta exactly and should use an average)
// Define the Boost to make beams back-to-back
const auto cmBoost = (ei + pi).BoostToCM();
const Boost boost_to_cm(-cmBoost);
// This will boost beams from a center of momentum frame back to (nearly) their original energies
const Boost boost_to_headon(cmBoost); // FIXME
// Boost and rotate the incoming beams to find the proper rotations TLorentzVector
// Boost to COM frame
boost_to_cm(pi);
boost_to_cm(ei);
// Rotate to head-on
RotationY rotAboutY(-1.0*atan2(pi.Px(), pi.Pz())); // Rotate to remove x component of beams
RotationX rotAboutX(+1.0*atan2(pi.Py(), pi.Pz())); // Rotate to remove y component of beams
LorentzRotation tf;
tf *= boost_to_cm;
tf *= rotAboutY;
tf *= rotAboutX;
tf *= boost_to_headon;
return tf;
}
static PxPyPzEVector apply_boost(const LorentzRotation& tf, PxPyPzEVector part) {
// 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)
// Boost and rotate particle 4-momenta into the headon frame
tf(part);
return part;
}
} // namespace Jug::Reco::Boost
......@@ -10,13 +10,15 @@
#include "JugBase/IParticleSvc.h"
#include "JugBase/DataHandle.h"
#include "eicd/VectorXYZT.h"
#include "JugReco/Utilities/Beam.h"
#include "JugReco/Utilities/Boost.h"
#include "TLorentzVector.h"
#include "TVector3.h"
#include "Math/Vector4D.h"
using ROOT::Math::PxPyPzEVector;
// Event Model related classes
#include "edm4hep/MCParticleCollection.h"
#include "eicd/MCRecoParticleAssociationCollection.h"
#include "eicd/ReconstructedParticleCollection.h"
#include "eicd/InclusiveKinematicsCollection.h"
......@@ -32,6 +34,10 @@ private:
"ReconstructedParticles",
Gaudi::DataHandle::Reader,
this};
DataHandle<eicd::MCRecoParticleAssociationCollection> m_inputParticleAssociation{
"MCRecoParticleAssociation",
Gaudi::DataHandle::Reader,
this};
DataHandle<eicd::InclusiveKinematicsCollection> m_outputInclusiveKinematicsCollection{
"InclusiveKinematicsDA",
Gaudi::DataHandle::Writer,
......@@ -40,13 +46,14 @@ private:
Gaudi::Property<double> m_crossingAngle{this, "crossingAngle", -0.025 * Gaudi::Units::radian};
SmartIF<IParticleSvc> m_pidSvc;
double m_proton, m_neutron, m_electron;
double m_proton{0}, m_neutron{0}, m_electron{0};
public:
InclusiveKinematicsDA(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc) {
declareProperty("inputMCParticles", m_inputMCParticleCollection, "MCParticles");
declareProperty("inputParticles", m_inputParticleCollection, "ReconstructedParticles");
declareProperty("inputAssociations", m_inputParticleAssociation, "MCRecoParticleAssociation");
declareProperty("outputData", m_outputInclusiveKinematicsCollection, "InclusiveKinematicsDA");
}
......@@ -68,138 +75,73 @@ public:
return StatusCode::SUCCESS;
}
TLorentzVector apply_boost(eicd::VectorXYZT ei_eicv, eicd::VectorXYZT pi_eicv, TLorentzVector part){
//Step 1: Find the needed boosts and rotations from the incoming lepton and hadron beams
//(note, this will give you a perfect boost, in principle you will not know the beam momenta exactly and should use an average)
// Define the Boost to make beams back-to-back
TLorentzVector ei(ei_eicv.x,ei_eicv.y,ei_eicv.z,ei_eicv.t);
TLorentzVector pi(pi_eicv.x,pi_eicv.y,pi_eicv.z,pi_eicv.t);
TLorentzVector cmBoost = (1./ei.E())*ei + (1./pi.E())*pi;
TLorentzVector boost(-cmBoost.Px(),-cmBoost.Py(),-cmBoost.Pz(),cmBoost.E());
TVector3 b;
b = boost.BoostVector();
TLorentzVector boostBack(0.0,0.0,cmBoost.Pz(),cmBoost.E());
TVector3 bb;
bb = boostBack.BoostVector(); // This will boost beams from a center of momentum frame back to (nearly) their original energies
// Boost and rotate the incoming beams to find the proper rotations TLorentzVector
pi.Boost(b); // Boost to COM frame
ei.Boost(b);
double rotAboutY = -1.0*TMath::ATan2(pi.Px(),pi.Pz()); // Rotate to remove x component of beams
double rotAboutX = 1.0*TMath::ATan2(pi.Py(),pi.Pz()); // Rotate to remove y component of beams
//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)
//Boost and rotate particle 4-momenta into the headon frame
part.Boost(b);
part.RotateY(rotAboutY);
part.RotateX(rotAboutX);
part.Boost(bb);
return part;
}
StatusCode execute() override {
// input collections
const edm4hep::MCParticleCollection& mcparts = *(m_inputMCParticleCollection.get());
const eicd::ReconstructedParticleCollection& parts = *(m_inputParticleCollection.get());
const auto& mcparts = *(m_inputMCParticleCollection.get());
const auto& rcparts = *(m_inputParticleCollection.get());
const auto& rcassoc = *(m_inputParticleAssociation.get());
// output collection
auto& out_kinematics = *(m_outputInclusiveKinematicsCollection.createAndPut());
// Loop over generated particles to get incoming electron beam
eicd::VectorXYZT ei, pi;
bool found_electron = false;
bool found_proton = false;
//Also get the true scattered electron, which will not be included in the sum
//over final-state particles for the JB reconstruction
int32_t mcscatID = -1;
for (const auto& p : mcparts) {
if (p.getGeneratorStatus() == 4 && p.getPDG() == 11) {
// Incoming electron
ei.x = p.getMomentum().x;
ei.y = p.getMomentum().y;
ei.z = p.getMomentum().z;
ei.t = p.getEnergy();
//Should not include true event-by-event smearing of beam particles for reconstruction
//Find a better way to do this...
ei.x = 0.;
ei.y = 0.;
if( fabs(ei.z - 5.0) < 1.0 )
ei.z = -5.0;
else if( fabs(ei.z - 10.0) < 1.0 )
ei.z = -10.0;
else if( fabs(ei.z - 18.0) < 1.0 )
ei.z = -18.0;
ei.t = sqrt( ei.z*ei.z + m_electron*m_electron);
found_electron = true;
}
if (p.getGeneratorStatus() == 4 && (p.getPDG() == 2212 || p.getPDG() == 2112)) {
// Incoming proton
pi.x = p.getMomentum().x;
pi.y = p.getMomentum().y;
pi.z = p.getMomentum().z;
pi.t = p.getEnergy();
//Should not include true event-by-event smearing of beam particles for reconstruction
//Find a better way to do this...
pi.y = 0;
if( fabs(pi.z - 41.0) < 5.0 ){
pi.x = 41.0*sin(m_crossingAngle);
pi.z = 41.0*cos(m_crossingAngle);
}
else if( fabs(pi.z - 100.0) < 5.0 ){
pi.x = 100.0*sin(m_crossingAngle);
pi.z = 100.0*cos(m_crossingAngle);
}
else if( fabs(pi.z - 275.0) < 5.0 ){
pi.x = 275.0*sin(m_crossingAngle);
pi.z = 275.0*cos(m_crossingAngle);
}
pi.t = std::hypot(pi.x, pi.z, (p.getPDG() == 2212) ? m_proton : m_neutron);
found_proton = true;
}
// Index of true Scattered electron. Currently taken as first status==1 electron in HEPMC record.
if (p.getGeneratorStatus() == 1 && p.getPDG() == 11 && mcscatID == -1) {
mcscatID = p.id();
}
if (found_electron && found_proton && mcscatID != -1) {
break;
// Get incoming electron beam
const auto ei_coll = Jug::Reco::Beam::find_first_beam_electron(mcparts);
if (ei_coll.size() == 0) {
if (msgLevel(MSG::DEBUG)) {
debug() << "No beam electron found" << endmsg;
}
return StatusCode::SUCCESS;
}
if (found_electron == false) {
const PxPyPzEVector ei(
Jug::Reco::Beam::round_beam_four_momentum(
ei_coll[0].getMomentum(),
m_electron,
{-5.0, -10.0, -18.0},
0.0)
);
// Get incoming hadron beam
const auto pi_coll = Jug::Reco::Beam::find_first_beam_hadron(mcparts);
if (pi_coll.size() == 0) {
if (msgLevel(MSG::DEBUG)) {
debug() << "No initial electron found" << endmsg;
debug() << "No beam hadron found" << endmsg;
}
return StatusCode::SUCCESS;
}
if (found_proton == false) {
const PxPyPzEVector pi(
Jug::Reco::Beam::round_beam_four_momentum(
pi_coll[0].getMomentum(),
pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron,
{41.0, 100.0, 275.0},
m_crossingAngle)
);
// Get first scattered electron
const auto ef_coll = Jug::Reco::Beam::find_first_scattered_electron(mcparts);
if (ef_coll.size() == 0) {
if (msgLevel(MSG::DEBUG)) {
debug() << "No initial proton found" << endmsg;
debug() << "No truth scattered electron found" << endmsg;
}
return StatusCode::SUCCESS;
}
if (mcscatID == -1) {
// Associate first scattered electron with reconstructed electrons
//const auto ef_assoc = std::find_if(
// rcassoc.begin(),
// rcassoc.end(),
// [&ef_coll](const auto& a){ return a.getSimID() == ef_coll[0].id(); });
auto ef_assoc = rcassoc.begin();
for (; ef_assoc != rcassoc.end(); ++ef_assoc) {
if (ef_assoc->getSimID() == ef_coll[0].id()) {
break;
}
}
if (!(ef_assoc != rcassoc.end())) {
if (msgLevel(MSG::DEBUG)) {
debug() << "No truth scattered electron found" << endmsg;
debug() << "Truth scattered electron not in reconstructed particles" << endmsg;
}
return StatusCode::SUCCESS;
return StatusCode::SUCCESS;
}
const auto ef_rc{ef_assoc->getRec()};
const auto ef_rc_id{ef_rc.id()};
// Loop over reconstructed particles to get all outgoing particles
// -----------------------------------------------------------------
......@@ -214,9 +156,6 @@ public:
// magnitude of the momentum.
// -----------------------------------------------------------------
// Reconstructed Index of scattered electron
eicd::Index scatID;
//Sums in colinear frame
double pxsum = 0;
double pysum = 0;
......@@ -224,29 +163,31 @@ public:
double Esum = 0;
double theta_e = 0;
for(const auto& p : parts){
//Get the scattered electron index and angle
if(p.mcID().value == mcscatID){
TLorentzVector e_lab(p.p().x,p.p().y,p.p().z,p.energy());
TLorentzVector e_boosted = apply_boost(ei,pi,e_lab);
scatID = p.ID();
theta_e = e_boosted.Theta();
}
//Sum over all particles other than scattered electron
else{
//Lorentz vector in lab frame
TLorentzVector hf_lab(p.p().x,p.p().y,p.p().z,p.energy());
//Boost to colinear frame
TLorentzVector hf_boosted = apply_boost(ei,pi,hf_lab);
pxsum += hf_boosted.Px();
pysum += hf_boosted.Py();
pzsum += hf_boosted.Pz();
Esum += hf_boosted.E();
}
// Get boost to colinear frame
auto boost = Jug::Reco::Boost::determine_boost(ei, pi);
for(const auto& p : rcparts) {
// Get the scattered electron index and angle
if (p.id() == ef_rc_id) {
// Lorentz vector in lab frame
PxPyPzEVector e_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy());
// Boost to colinear frame
PxPyPzEVector e_boosted = Jug::Reco::Boost::apply_boost(boost, e_lab);
theta_e = e_boosted.Theta();
// Sum over all particles other than scattered electron
} else {
// Lorentz vector in lab frame
PxPyPzEVector hf_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy());
// Boost to colinear frame
PxPyPzEVector hf_boosted = Jug::Reco::Boost::apply_boost(boost, hf_lab);
pxsum += hf_boosted.Px();
pysum += hf_boosted.Py();
pzsum += hf_boosted.Pz();
Esum += hf_boosted.E();
}
}
// DIS kinematics calculations
......@@ -254,35 +195,34 @@ public:
auto ptsum = sqrt(pxsum*pxsum + pysum*pysum);
auto theta_h = 2.*atan(sigma_h/ptsum);
if (scatID && sigma_h>0) {
auto y_da = tan(theta_h/2.) / ( tan(theta_e/2.) + tan(theta_h/2.) );
auto Q2_da = 4.*ei.energy()*ei.energy() * ( 1. / tan(theta_e/2.) ) * ( 1. / (tan(theta_e/2.) + tan(theta_h/2.)) );
auto x_da = Q2_da / (4.*ei.energy()*pi.energy()*y_da);
auto nu_da = Q2_da / (2.*m_proton*x_da);
auto W_da = sqrt ( m_proton*m_proton + 2*m_proton*nu_da - Q2_da );
auto kin = out_kinematics.create();
kin.Q2(Q2_da);
kin.y(y_da);
kin.nu(nu_da);
kin.x(x_da);
kin.W(W_da);
kin.scatID(scatID); //MC index of scattered electron
// Debugging output
// If no scattered electron was found
if (sigma_h <= 0) {
if (msgLevel(MSG::DEBUG)) {
debug() << "pi = " << pi << endmsg;
debug() << "ei = " << ei << endmsg;
debug() << "x,y,Q2,W,nu = "
<< kin.x() << ","
<< kin.y() << ","
<< kin.Q2() << ","
<< kin.W() << ","
<< kin.nu()
<< endmsg;
debug() << "No scattered electron found" << endmsg;
}
return StatusCode::SUCCESS;
}
// Calculate kinematic variables
const auto y_da = tan(theta_h/2.) / ( tan(theta_e/2.) + tan(theta_h/2.) );
const auto Q2_da = 4.*ei.energy()*ei.energy() * ( 1. / tan(theta_e/2.) ) * ( 1. / (tan(theta_e/2.) + tan(theta_h/2.)) );
const auto x_da = Q2_da / (4.*ei.energy()*pi.energy()*y_da);
const auto nu_da = Q2_da / (2.*m_proton*x_da);
const auto W_da = sqrt(m_proton*m_proton + 2*m_proton*nu_da - Q2_da);
auto kin = out_kinematics.create(x_da, Q2_da, W_da, y_da, nu_da);
kin.setScat(ef_rc);
// Debugging output
if (msgLevel(MSG::DEBUG)) {
debug() << "pi = " << pi << endmsg;
debug() << "ei = " << ei << endmsg;
debug() << "x,Q2,W,y,nu = "
<< kin.getX() << ","
<< kin.getQ2() << ","
<< kin.getW() << ","
<< kin.getY() << ","
<< kin.getNu()
<< endmsg;
}
return StatusCode::SUCCESS;
......
......@@ -5,15 +5,21 @@
#include "GaudiKernel/RndmGenerators.h"
#include "GaudiKernel/PhysicalConstants.h"
#include <algorithm>
#include <ranges>
#include <cmath>
#include "JugBase/IParticleSvc.h"
#include "JugBase/DataHandle.h"
#include "eicd/VectorXYZT.h"
#include "JugReco/Utilities/Beam.h"
#include "JugReco/Utilities/Boost.h"
#include "Math/Vector4D.h"
using ROOT::Math::PxPyPzEVector;
// Event Model related classes
#include "edm4hep/MCParticleCollection.h"
#include "eicd/MCRecoParticleAssociationCollection.h"
#include "eicd/ReconstructedParticleCollection.h"
#include "eicd/InclusiveKinematicsCollection.h"
......@@ -29,6 +35,10 @@ private:
"ReconstructedParticles",
Gaudi::DataHandle::Reader,
this};
DataHandle<eicd::MCRecoParticleAssociationCollection> m_inputParticleAssociation{
"MCRecoParticleAssociation",
Gaudi::DataHandle::Reader,
this};
DataHandle<eicd::InclusiveKinematicsCollection> m_outputInclusiveKinematicsCollection{
"InclusiveKinematicsElectron",
Gaudi::DataHandle::Writer,
......@@ -37,13 +47,14 @@ private:
Gaudi::Property<double> m_crossingAngle{this, "crossingAngle", -0.025 * Gaudi::Units::radian};
SmartIF<IParticleSvc> m_pidSvc;
double m_proton, m_neutron, m_electron;
double m_proton{0}, m_neutron{0}, m_electron{0};
public:
InclusiveKinematicsElectron(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc) {
declareProperty("inputMCParticles", m_inputMCParticleCollection, "MCParticles");
declareProperty("inputParticles", m_inputParticleCollection, "ReconstructedParticles");
declareProperty("inputAssociations", m_inputParticleAssociation, "MCRecoParticleAssociation");
declareProperty("outputData", m_outputInclusiveKinematicsCollection, "InclusiveKinematicsElectron");
}
......@@ -67,138 +78,164 @@ public:
StatusCode execute() override {
// input collections
const edm4hep::MCParticleCollection& mcparts = *(m_inputMCParticleCollection.get());
const eicd::ReconstructedParticleCollection& parts = *(m_inputParticleCollection.get());
const auto& mcparts = *(m_inputMCParticleCollection.get());
const auto& rcparts = *(m_inputParticleCollection.get());
const auto& rcassoc = *(m_inputParticleAssociation.get());
// output collection
auto& out_kinematics = *(m_outputInclusiveKinematicsCollection.createAndPut());
// Loop over generated particles to get incoming electron beam
eicd::VectorXYZT ei, pi;
bool found_electron = false;
bool found_proton = false;
//Also get the true scattered electron, which will not be included in the sum
//over final-state particles for the JB reconstruction
int32_t mcscatID = -1;
for (const auto& p : mcparts) {
if (p.getGeneratorStatus() == 4 && p.getPDG() == 11) {
// Incoming electron
ei.x = p.getMomentum().x;
ei.y = p.getMomentum().y;
ei.z = p.getMomentum().z;
ei.t = p.getEnergy();
//Should not include true event-by-event smearing of beam particles for reconstruction
//Find a better way to do this...
ei.x = 0.;
ei.y = 0.;
if( fabs(ei.z - 5.0) < 1.0 )
ei.z = -5.0;