Skip to content

GitLab

Commit 35fba059 authored by Whitney Armstrong's avatar Whitney Armstrong
Browse files

modified: ../../CMakeLists.txt 

	modified:   MC2DummyParticle.cpp
	modified:   ../../../JugReco/src/components/ParticlesFromTrackFit.cpp
parent 208cefdf
Hide whitespace changes
Inline Side-by-side
Showing with 78 additions and 77 deletions
JugBase/CMakeLists.txt
View file @ 35fba059
......@@ -27,7 +27,7 @@ gaudi_add_library(JugBase
src/Plugins/BFieldUtils.cpp
INCLUDE_DIRS EICD PODIO ROOT $ENV{HOME}/stow/podio/include
LINK_LIBRARIES GaudiAlgLib GaudiKernel ROOT DD4hep::DDG4IO
PUBLIC_HEADERS JugBase)
PUBLIC_HEADERS JugBase )
target_link_libraries(JugBase
podio::podioRootIO
ActsCore ActsPluginDD4hep
......@@ -37,7 +37,7 @@ target_compile_options(JugBase PRIVATE -Wno-suggest-override)
gaudi_add_module(JugBasePlugins
src/components/*.cpp
src/components/ReadTestConsumer.cxx
LINK_LIBRARIES GaudiAlgLib GaudiKernel JugBase ROOT NPDet::DD4podIO ActsCore ActsPluginDD4hep)
LINK_LIBRARIES GaudiAlgLib GaudiKernel JugBase ROOT NPDet::DD4podIO ActsCore ActsPluginDD4hep EICD::eicd)
target_compile_options(JugBasePlugins PRIVATE -Wno-suggest-override)
#gaudi_add_test(ProduceForReadTest
......
JugBase/src/components/MC2DummyParticle.cpp
View file @ 35fba059
......@@ -18,7 +18,9 @@ namespace Jug {
class MC2DummyParticle : public GaudiAlgorithm {
public:
// ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
DataHandle<dd4pod::Geant4ParticleCollection> m_inputHitCollection{"mcparticles", Gaudi::DataHandle::Reader, this};
DataHandle<eic::ReconstructedParticleCollection> m_outputHitCollection{"DummyReconstructedParticles", Gaudi::DataHandle::Writer, this};
MC2DummyParticle(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc)
{
declareProperty("inputCollection", m_inputHitCollection, "mcparticles");
......@@ -37,16 +39,13 @@ namespace Jug {
// output collection
auto out_parts = m_outputHitCollection.createAndPut();
for (const auto& p : *parts) {
double momentum = std::hypot(p.psx(), p.psy(), p.psz());
double energy = std::hypot(momentum, p.mass());
out_parts->push_back({p.pdgID(), energy, {p.psx(),p.psy(),p.psz()}, p.charge(), p.mass()});
eic::ReconstructedParticle rec_part(p.pdgID(), energy, {p.psx(),p.psy(),p.psz()}, (double)p.charge(), p.mass());
out_parts->push_back(rec_part);
}
return StatusCode::SUCCESS;
}
DataHandle<dd4pod::Geant4ParticleCollection> m_inputHitCollection{"mcparticles", Gaudi::DataHandle::Reader, this};
DataHandle<eic::ReconstructedParticleCollection> m_outputHitCollection{"DummyReconstructedParticles", Gaudi::DataHandle::Writer, this};
};
DECLARE_COMPONENT(MC2DummyParticle)
......
JugReco/src/components/ParticlesFromTrackFit.cpp
View file @ 35fba059
......@@ -63,6 +63,8 @@ namespace Jug {
auto rec_parts = m_outputParticles.createAndPut();
auto track_pars = m_outputTrackParameters.createAndPut();
debug() << std::size(*trajectories) << " trajectories " << endmsg;
for(const auto& traj : *trajectories) {
//traj.trajectory().first
const auto& [trackTips, mj] = traj.trajectory();
......@@ -71,76 +73,76 @@ namespace Jug {
continue;
}
// Get the entry index for the single trajectory
auto& trackTip = trackTips.front();
// Collect the trajectory summary info
auto trajState = Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip);
int m_nMeasurements = trajState.nMeasurements;
int m_nStates = trajState.nStates;
// Get the fitted track parameter
bool m_hasFittedParams = false;
if (traj.hasTrackParameters(trackTip)) {
m_hasFittedParams = true;
const auto& boundParam = traj.trackParameters(trackTip);
const auto& parameter = boundParam.parameters();
const auto& covariance = *boundParam.covariance();
debug() << "loc 0 = " << parameter[Acts::eBoundLoc0] << endmsg;
debug() << "loc 1 = " << parameter[Acts::eBoundLoc1] << endmsg;
debug() << "phi = " << parameter[Acts::eBoundPhi] << endmsg;
debug() << "theta = " << parameter[Acts::eBoundTheta] << endmsg;
debug() << "q/p = " << parameter[Acts::eBoundQOverP] << endmsg;
debug() << "p = " << 1.0/parameter[Acts::eBoundQOverP] << endmsg;
debug() << "err phi = " << sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)) << endmsg;
debug() << "err th = " << sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)) << endmsg;
debug() << "err q/p = " << sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP))<< endmsg;
debug() << " chi2 = " << trajState.chi2Sum << endmsg;
eic::TrackParameters pars({
parameter[Acts::eBoundLoc0], parameter[Acts::eBoundLoc1], parameter[Acts::eBoundPhi],
parameter[Acts::eBoundTheta], parameter[Acts::eBoundQOverP],parameter[Acts::eBoundTime],
sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)), sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)),
sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)), sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)),
sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)),sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime))});
track_pars->push_back(pars);
//m_ePHI_fit = parameter[Acts::eBoundPhi];
//m_eTHETA_fit = parameter[Acts::eBoundTheta];
//m_eQOP_fit = parameter[Acts::eBoundQOverP];
//m_eT_fit = parameter[Acts::eBoundTime];
//m_err_eLOC0_fit = sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0));
//m_err_eLOC1_fit = sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1));
//m_err_ePHI_fit = sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi));
//m_err_eTHETA_fit = sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta));
//m_err_eQOP_fit = sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP));
//m_err_eT_fit = sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime));
}
auto tsize = traj.trajectory().first.size();
debug() << "# fitted parameters : " << tsize << endmsg;
if(tsize == 0 ) continue;
traj.trajectory().second.visitBackwards(tsize-1, [&](auto&& trackstate) {
//debug() << trackstate.hasPredicted() << endmsg;
//debug() << trackstate.predicted() << endmsg;
auto params = trackstate.predicted() ;//<< endmsg;
double p0 = (1.0 / params[Acts::eBoundQOverP]) / Acts::UnitConstants::GeV;
debug() << "track predicted p = " << p0 << " GeV" << endmsg;
if ( std::abs(p0) > 500) {
debug() << "skipping" << endmsg;
return;
}
// Get the entry index for the single trajectory
auto& trackTip = trackTips.front();
// Collect the trajectory summary info
auto trajState = Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip);
int m_nMeasurements = trajState.nMeasurements;
int m_nStates = trajState.nStates;
// Get the fitted track parameter
bool m_hasFittedParams = false;
if (traj.hasTrackParameters(trackTip)) {
m_hasFittedParams = true;
const auto& boundParam = traj.trackParameters(trackTip);
const auto& parameter = boundParam.parameters();
const auto& covariance = *boundParam.covariance();
debug() << "loc 0 = " << parameter[Acts::eBoundLoc0] << endmsg;
debug() << "loc 1 = " << parameter[Acts::eBoundLoc1] << endmsg;
debug() << "phi = " << parameter[Acts::eBoundPhi] << endmsg;
debug() << "theta = " << parameter[Acts::eBoundTheta] << endmsg;
debug() << "q/p = " << parameter[Acts::eBoundQOverP] << endmsg;
debug() << "p = " << 1.0/parameter[Acts::eBoundQOverP] << endmsg;
debug() << "err phi = " << sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)) << endmsg;
debug() << "err th = " << sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)) << endmsg;
debug() << "err q/p = " << sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP))<< endmsg;
debug() << " chi2 = " << trajState.chi2Sum << endmsg;
eic::TrackParameters pars({
parameter[Acts::eBoundLoc0], parameter[Acts::eBoundLoc1], parameter[Acts::eBoundPhi],
parameter[Acts::eBoundTheta], parameter[Acts::eBoundQOverP],parameter[Acts::eBoundTime],
sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)), sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)),
sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)), sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)),
sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)),sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime))});
track_pars->push_back(pars);
//m_ePHI_fit = parameter[Acts::eBoundPhi];
//m_eTHETA_fit = parameter[Acts::eBoundTheta];
//m_eQOP_fit = parameter[Acts::eBoundQOverP];
//m_eT_fit = parameter[Acts::eBoundTime];
//m_err_eLOC0_fit = sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0));
//m_err_eLOC1_fit = sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1));
//m_err_ePHI_fit = sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi));
//m_err_eTHETA_fit = sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta));
//m_err_eQOP_fit = sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP));
//m_err_eT_fit = sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime));
}
eic::Particle p({params[Acts::eBoundPhi], params[Acts::eBoundTheta], 1.0 / std::abs(params[Acts::eBoundQOverP]), 0.000511},
{0.0, 0.0, 0.0, params[Acts::eBoundTime]},
(long long)11 * params[Acts::eBoundQOverP] / std::abs(params[Acts::eBoundQOverP]), 0);
//debug() << p << endmsg;
rec_parts->push_back(p);
});
auto tsize = traj.trajectory().first.size();
debug() << "# fitted parameters : " << tsize << endmsg;
if(tsize == 0 ) continue;
traj.trajectory().second.visitBackwards(tsize-1, [&](auto&& trackstate) {
//debug() << trackstate.hasPredicted() << endmsg;
//debug() << trackstate.predicted() << endmsg;
auto params = trackstate.predicted() ;//<< endmsg;
double p0 = (1.0 / params[Acts::eBoundQOverP]) / Acts::UnitConstants::GeV;
debug() << "track predicted p = " << p0 << " GeV" << endmsg;
if ( std::abs(p0) > 500) {
debug() << "skipping" << endmsg;
return;
}
eic::Particle p({params[Acts::eBoundPhi], params[Acts::eBoundTheta], 1.0 / std::abs(params[Acts::eBoundQOverP]), 0.000511},
{0.0, 0.0, 0.0, params[Acts::eBoundTime]},
(long long)11 * params[Acts::eBoundQOverP] / std::abs(params[Acts::eBoundQOverP]), 0);
//debug() << p << endmsg;
rec_parts->push_back(p);
});
}
return StatusCode::SUCCESS;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment