Commit 258894e7 authored by Whitney Armstrong's avatar Whitney Armstrong
Browse files

Adding tracker param init from cluster.

- All charges (-1,0,1) are added to the candidates
- The initial vector points to the EM cluster, set with momentum equal
to the energy where we assume a massless particle
parent 3d3870af
...@@ -33,6 +33,7 @@ gaudi_add_module(JugRecoPlugins ...@@ -33,6 +33,7 @@ gaudi_add_module(JugRecoPlugins
src/components/TrackFindingAlgorithm.cpp src/components/TrackFindingAlgorithm.cpp
src/components/TestACTSLogger.cpp src/components/TestACTSLogger.cpp
src/components/TrackParamTruthInit.cpp src/components/TrackParamTruthInit.cpp
src/components/TrackParamClusterInit.cpp
src/components/SimpleClustering.cpp src/components/SimpleClustering.cpp
src/components/ParticlesFromTrackFit.cpp src/components/ParticlesFromTrackFit.cpp
src/components/CalorimeterIslandCluster.cpp src/components/CalorimeterIslandCluster.cpp
......
...@@ -86,11 +86,9 @@ namespace Jug::Reco { ...@@ -86,11 +86,9 @@ namespace Jug::Reco {
while (continue_clustering) { while (continue_clustering) {
std::vector<eic::CalorimeterHit> cluster_hits; std::vector<eic::CalorimeterHit> cluster_hits;
debug() << " ref_hit" << ref_hit << endmsg; //debug() << " ref_hit" << ref_hit << endmsg;
// auto ref_hit = *(std::max_element(std::begin(hits), std::end(hits), // auto ref_hit = *(std::max_element(std::begin(hits), std::end(hits),
// [](const auto& h1, const auto& h2) { return h1.energy() < h2.energy(); })); // [](const auto& h1, const auto& h2) { return h1.energy() < h2.energy(); }));
// std::partition_copy(std::begin(hits), std::end(hits), std::begin(first_cluster_hits), // std::partition_copy(std::begin(hits), std::end(hits), std::begin(first_cluster_hits),
// std::begin(remaining_hits), // std::begin(remaining_hits),
// [=](const auto& h) { // [=](const auto& h) {
...@@ -122,7 +120,7 @@ namespace Jug::Reco { ...@@ -122,7 +120,7 @@ namespace Jug::Reco {
cluster0.y(cluster0.y() + h.energy() * h.y() / total_energy); cluster0.y(cluster0.y() + h.energy() * h.y() / total_energy);
cluster0.z(cluster0.z() + h.energy() * h.z() / total_energy); cluster0.z(cluster0.z() + h.energy() * h.z() / total_energy);
} }
debug() << " cluster = " << cluster0 << endmsg; //debug() << " cluster = " << cluster0 << endmsg;
clusters.push_back(cluster0); clusters.push_back(cluster0);
......
...@@ -70,7 +70,9 @@ namespace Jug::Reco { ...@@ -70,7 +70,9 @@ namespace Jug::Reco {
} }
m_BField = std::make_shared<Acts::ConstantBField>(Acts::Vector3D{0.0, 0.0, m_geoSvc->centralMagneticField()}); m_BField = std::make_shared<Acts::ConstantBField>(Acts::Vector3D{0.0, 0.0, m_geoSvc->centralMagneticField()});
m_fieldctx = BFieldVariant(m_BField); m_fieldctx = BFieldVariant(m_BField);
m_sourcelinkSelectorCfg = { {Acts::GeometryIdentifier(), {15, 10}},
// chi2 and #sourclinks per surface cutoffs
m_sourcelinkSelectorCfg = { {Acts::GeometryIdentifier(), {3, 3}},
}; };
findTracks = TrackFindingAlgorithm::makeTrackFinderFunction(m_geoSvc->trackingGeometry(), m_BField); findTracks = TrackFindingAlgorithm::makeTrackFinderFunction(m_geoSvc->trackingGeometry(), m_BField);
...@@ -123,7 +125,7 @@ namespace Jug::Reco { ...@@ -123,7 +125,7 @@ namespace Jug::Reco {
debug() << "Track finding failed for truth seed " << iseed << endmsg; debug() << "Track finding failed for truth seed " << iseed << endmsg;
ACTS_WARNING("Track finding failed for truth seed " << iseed << " with error" << result.error()); ACTS_WARNING("Track finding failed for truth seed " << iseed << " with error" << result.error());
// Track finding failed, but still create an empty SimMultiTrajectory // Track finding failed, but still create an empty SimMultiTrajectory
trajectories->push_back(SimMultiTrajectory()); //trajectories->push_back(SimMultiTrajectory());
} }
} }
......
#include <cmath>
// Gaudi
#include "GaudiAlg/GaudiAlgorithm.h"
#include "GaudiKernel/ToolHandle.h"
#include "GaudiAlg/Transformer.h"
#include "GaudiAlg/GaudiTool.h"
#include "GaudiKernel/RndmGenerators.h"
#include "GaudiKernel/Property.h"
#include "JugBase/DataHandle.h"
#include "JugBase/IGeoSvc.h"
#include "JugReco/Track.hpp"
#include "Acts/Utilities/Units.hpp"
#include "eicd/TrackerHitCollection.h"
#include "eicd/ClusterCollection.h"
///// (Reconstructed) track parameters e.g. close to the vertex.
//using TrackParameters = Acts::CurvilinearTrackParameters;
///// Container of reconstructed track states for multiple tracks.
//using TrackParametersContainer = std::vector<TrackParameters>;
///// MultiTrajectory definition
//using Trajectory = Acts::MultiTrajectory<SourceLink>;
///// Container for the truth fitting/finding track(s)
//using TrajectoryContainer = std::vector<SimMultiTrajectory>;
namespace Jug::Reco {
/** Initial Track parameters from MC truth.
*
* TrackParmetersContainer
*/
class TrackParamClusterInit : public GaudiAlgorithm {
public:
using Clusters = eic::ClusterCollection;
DataHandle<Clusters> m_inputClusters{"inputClusters", Gaudi::DataHandle::Reader, this};
DataHandle<TrackParametersContainer> m_outputInitialTrackParameters{"outputInitialTrackParameters",
Gaudi::DataHandle::Writer, this};
public:
TrackParamClusterInit(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc) {
declareProperty("inputClusters", m_inputClusters, "Input clusters");
declareProperty("outputInitialTrackParameters", m_outputInitialTrackParameters, "");
}
StatusCode initialize() override {
if (GaudiAlgorithm::initialize().isFailure())
return StatusCode::FAILURE;
IRndmGenSvc* randSvc = svc<IRndmGenSvc>("RndmGenSvc", true);
if(!randSvc)
return StatusCode::FAILURE;
return StatusCode::SUCCESS;
}
StatusCode execute() override {
// input collection
const eic::ClusterCollection* clusters = m_inputClusters.get();
// Create output collections
auto init_trk_params = m_outputInitialTrackParameters.createAndPut();
for(const auto& c : *clusters) {
using Acts::UnitConstants::GeV;
using Acts::UnitConstants::MeV;
using Acts::UnitConstants::mm;
using Acts::UnitConstants::ns;
double p = c.energy()*GeV;
double len = std::hypot( c.x() , c.y() , c.z() );
// build some track cov matrix
Acts::BoundSymMatrix cov = Acts::BoundSymMatrix::Zero();
cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = 1.0 * mm*1.0 * mm;
cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = 1.0 * mm*1.0 * mm;
cov(Acts::eBoundPhi, Acts::eBoundPhi) = M_PI / 180.0;
cov(Acts::eBoundTheta, Acts::eBoundTheta) = M_PI / 180.0;
cov(Acts::eBoundQOverP, Acts::eBoundQOverP) = 1.0 / (0.05*0.05*p*p);
cov(Acts::eBoundTime, Acts::eBoundTime) = Acts::UnitConstants::ns;
// add all charges to the track candidate...
init_trk_params->emplace_back(Acts::Vector4D(0 * mm, 0 * mm, 0 * mm, 0),
Acts::Vector3D(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, -1,
std::make_optional(cov));
init_trk_params->emplace_back(Acts::Vector4D(0 * mm, 0 * mm, 0 * mm, 0),
Acts::Vector3D(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, 1,
std::make_optional(cov));
init_trk_params->emplace_back(Acts::Vector4D(0 * mm, 0 * mm, 0 * mm, 0),
Acts::Vector3D(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, 0,
std::make_optional(cov));
debug() << "Invoke track finding seeded by truth particle with p = " << p/GeV << " GeV" << endmsg;
}
return StatusCode::SUCCESS;
}
};
DECLARE_COMPONENT(TrackParamClusterInit)
} // namespace Jug::reco
...@@ -83,10 +83,10 @@ namespace Jug { ...@@ -83,10 +83,10 @@ namespace Jug {
//eic::TrackerHit hit; //eic::TrackerHit hit;
eic::TrackerHit hit((long long)ahit.cellID(), eic::TrackerHit hit((long long)ahit.cellID(),
(float)ahit.time()/1000, // ps (float)ahit.time()/1000, // ps
(float)ahit.charge()/ 1.0e6, // MeV (float)ahit.charge()/ 1.0e6, // GeV
(float)0.0, (float)0.0,
{pos.x()/dd4hep::mm, pos.y()/dd4hep::mm,pos.z()/dd4hep::mm}, {pos.x()/dd4hep::mm, pos.y()/dd4hep::mm,pos.z()/dd4hep::mm},
{dim[0],dim[1],0.0}); {dim[0]/dd4hep::mm,dim[1]/dd4hep::mm,0.0});
rec_hits->push_back(hit); rec_hits->push_back(hit);
} }
return StatusCode::SUCCESS; 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