diff --git a/CMakeLists.txt b/CMakeLists.txt index 0b8a4bfd306c0f231fe6061885b58986b773f116..f501726de156894144f59ed0d61294db6d0643f5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -28,6 +28,8 @@ find_package(ROOT COMPONENTS MathCore GenVector Geom REQUIRED) find_package(DD4hep COMPONENTS DDG4 DDG4IO DDRec REQUIRED) find_package(Geant4) +find_library(genfit2 genfit2 /usr/local/lib REQUIRED) + find_package(GaudiProject) gaudi_project(Juggler v1r0 diff --git a/JugBase/JugBase/IGeoSvc.h b/JugBase/JugBase/IGeoSvc.h index a6b4bd5a79317c795bcb047452db5ffcf0c6db87..01b7509613c6fbe60d2352223b67cc6b07aa8bb8 100644 --- a/JugBase/JugBase/IGeoSvc.h +++ b/JugBase/JugBase/IGeoSvc.h @@ -10,22 +10,26 @@ #define IGEOSVC_H #include "GaudiKernel/IService.h" +#include <unordered_map> namespace dd4hep { -class Detector; -class DetElement; -namespace rec { -class CellIDPositionConverter; -} -} + class Detector; + class DetElement; + namespace rec { + class CellIDPositionConverter; + } +} // namespace dd4hep namespace Acts { -class TrackingGeometry; + class TrackingGeometry; + class Surface; } class G4VUserDetectorConstruction; class GAUDI_API IGeoSvc : virtual public IService { +public: + using VolumeSurfaceMap = std::unordered_map<uint64_t, const Acts::Surface*>; public: /// InterfaceID @@ -39,6 +43,7 @@ public: virtual double centralMagneticField() const = 0; // receive Geant4 Geometry //virtual G4VUserDetectorConstruction* getGeant4Geo() = 0; + virtual const VolumeSurfaceMap& surfaceMap() const = 0; virtual ~IGeoSvc() {} }; diff --git a/JugBase/src/components/GeoSvc.cpp b/JugBase/src/components/GeoSvc.cpp index 50cb8134302862ff8208ef4f7e7460efcb6721bb..3d60994b7fd2fa0949d011c589165b5efd830845 100644 --- a/JugBase/src/components/GeoSvc.cpp +++ b/JugBase/src/components/GeoSvc.cpp @@ -29,13 +29,14 @@ void draw_surfaces(std::shared_ptr<const Acts::TrackingGeometry> trk_geo, const trk_geo->visitSurfaces([&](const Acts::Surface* surface) { // for now we just require a valid surface if (not surface) { + std::cout << " Not a surface \n"; return; } surfaces.push_back(surface); }); std::ofstream os; os.open(fname); - os << std::fixed << std::setprecision(4); + os << std::fixed << std::setprecision(6); size_t nVtx = 0; for (const auto& srfx : surfaces) { const PlaneSurface* srf = dynamic_cast<const PlaneSurface*>(srfx); @@ -126,8 +127,24 @@ StatusCode GeoSvc::initialize() { m_trackingGeo = std::move(Acts::convertDD4hepDetector(m_dd4hepgeo->world(), Acts::Logging::VERBOSE, Acts::equidistant, Acts::equidistant, Acts::equidistant)); if (m_trackingGeo) { - draw_surfaces(m_trackingGeo, "Derp.obj"); + draw_surfaces(m_trackingGeo, "tracking_geometry.obj"); } + + debug() << "visiting all the surfaces " << endmsg; + m_trackingGeo->visitSurfaces([this](const Acts::Surface* surface) { + // for now we just require a valid surface + if (not surface) { + return; + } + auto det_element = + dynamic_cast<const Acts::DD4hepDetectorElement*>(surface->associatedDetectorElement()); + if (!det_element) { + debug() << "invalid det_element!!! " << endmsg; + return; + } + this->m_surfaces.insert_or_assign(det_element->identifier(), surface); + }); + return StatusCode::SUCCESS; } diff --git a/JugBase/src/components/GeoSvc.h b/JugBase/src/components/GeoSvc.h index 2efeeec3c284fb31087fab16d1b2f0c42aa808c7..afc6fd830736a38e2a5410e5cddd58df54e46e52 100644 --- a/JugBase/src/components/GeoSvc.h +++ b/JugBase/src/components/GeoSvc.h @@ -11,10 +11,13 @@ // Interface #include "JugBase/IGeoSvc.h" -#include "Acts/Definitions/Units.hpp" -#include "DD4hep/DD4hepUnits.h" -//using namespace Acts::UnitLiterals; +// ACTS +#include "Acts/Definitions/Units.hpp" +#include "Acts/Surfaces/Surface.hpp" +#include "Acts/Definitions/Common.hpp" +#include "Acts/Geometry/TrackingGeometry.hpp" +#include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp" // Gaudi #include "GaudiKernel/MsgStream.h" @@ -26,6 +29,7 @@ #include "DDRec/CellIDPositionConverter.h" #include "DDRec/SurfaceManager.h" #include "DDRec/Surface.h" +#include "DD4hep/DD4hepUnits.h" // Create a test context //#define CHECK_ROTATION_ANGLE(t, a, tolerance) \ @@ -37,8 +41,29 @@ //using SrfVec = std::vector<std::shared_ptr<const Surface>>; void draw_surfaces(std::shared_ptr<const Acts::TrackingGeometry> trk_geo, const std::string& fname); - class GeoSvc : public extends<Service, IGeoSvc> { +public: + using VolumeSurfaceMap = std::unordered_map<uint64_t, const Acts::Surface*>; + +private: + /// ACTS Tracking Geometry + std::shared_ptr<const Acts::TrackingGeometry> m_trackingGeo; + + /// Lookup container for hit surfaces that generate smeared hits + VolumeSurfaceMap m_surfaces; + + /// Pointer to the interface to the DD4hep geometry + dd4hep::Detector* m_dd4hepgeo; + + std::shared_ptr<const dd4hep::rec::CellIDPositionConverter> m_cellid_converter = + nullptr; //(*(m_geoSvc->detector())); + + /// XML-files with the detector description + Gaudi::Property<std::vector<std::string>> m_xmlFileNames{ + this, "detectors", {}, "Detector descriptions XML-files"}; + + /// output + MsgStream m_log; public: GeoSvc(const std::string& name, ISvcLocator* svc); @@ -59,41 +84,33 @@ public: */ virtual dd4hep::DetElement getDD4HepGeo() override; - virtual std::shared_ptr<const dd4hep::rec::CellIDPositionConverter> cellIDPositionConverter() const { + virtual std::shared_ptr<const dd4hep::rec::CellIDPositionConverter> + cellIDPositionConverter() const + { return m_cellid_converter; } /** Get the main dd4hep Detector. * Returns the pointer to the main dd4hep detector class. */ - virtual dd4hep::Detector* detector() override ; + virtual dd4hep::Detector* detector() override; /** Gets the ACTS tracking geometry. */ virtual std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry() const; - virtual double centralMagneticField() const { - return m_dd4hepgeo->field().magneticField({0,0,0}).z()*(Acts::UnitConstants::T/dd4hep::tesla); + virtual double centralMagneticField() const + { + return m_dd4hepgeo->field().magneticField({0, 0, 0}).z() * + (Acts::UnitConstants::T / dd4hep::tesla); } - -private: - - /// ACTS Tracking Geometry - std::shared_ptr<const Acts::TrackingGeometry> m_trackingGeo; - - /// Pointer to the interface to the DD4hep geometry - dd4hep::Detector* m_dd4hepgeo; - - std::shared_ptr<const dd4hep::rec::CellIDPositionConverter> m_cellid_converter = nullptr;//(*(m_geoSvc->detector())); - - /// XML-files with the detector description - Gaudi::Property<std::vector<std::string>> m_xmlFileNames{this, "detectors", {}, "Detector descriptions XML-files"}; - - /// output - MsgStream m_log; + virtual const VolumeSurfaceMap& surfaceMap() const { return m_surfaces; } }; -inline std::shared_ptr<const Acts::TrackingGeometry> GeoSvc::trackingGeometry() const { return m_trackingGeo; } +inline std::shared_ptr<const Acts::TrackingGeometry> GeoSvc::trackingGeometry() const +{ + return m_trackingGeo; +} -#endif // GEOSVC_H +#endif // GEOSVC_H diff --git a/JugTrack/CMakeLists.txt b/JugTrack/CMakeLists.txt index 53425a5901895fbd8fa6bdb27152bf1336ea8605..d8da234fb51a2de9f2ff0180f980bbd104108cf2 100644 --- a/JugTrack/CMakeLists.txt +++ b/JugTrack/CMakeLists.txt @@ -13,6 +13,8 @@ find_package(Acts REQUIRED COMPONENTS Core PluginIdentification PluginTGeo Plugi find_package(ROOT COMPONENTS RIO Tree Core REQUIRED) find_package(DD4hep COMPONENTS DDG4 DDG4IO DDRec REQUIRED) +find_library(genfit2 genfit2 /usr/local/lib REQUIRED) + # this declaration will not be needed in the future gaudi_depends_on_subdirs(GaudiAlg GaudiKernel) #gaudi_install_scripts() @@ -28,6 +30,8 @@ gaudi_depends_on_subdirs(GaudiAlg GaudiKernel) # ) gaudi_add_module(JugTrackPlugins + src/components/GenFitTrackFitter.cpp + src/components/SingleTrackSourceLinker.cpp src/components/TrackerSourceLinker.cpp src/components/Tracker2SourceLinker.cpp src/components/TrackerSourcesLinker.cpp @@ -40,7 +44,7 @@ gaudi_add_module(JugTrackPlugins src/components/TrackParamClusterInit.cpp src/components/TrackParamVertexClusterInit.cpp src/components/ParticlesFromTrackFit.cpp - LINK_LIBRARIES GaudiAlgLib GaudiKernel JugBase ROOT NPDet::DD4podIO EICD::eicd DDRec Acts) + LINK_LIBRARIES GaudiAlgLib GaudiKernel JugBase ROOT NPDet::DD4podIO EICD::eicd DDRec Acts genfit2) target_compile_options(JugTrackPlugins PRIVATE -Wno-suggest-override) diff --git a/JugTrack/JugTrack/ProtoTrack.hpp b/JugTrack/JugTrack/ProtoTrack.hpp index 12fdb56f9d5c92902ed040adcdb1a63f820823c1..8b2b7e280c53e0bca792003f9d4cd18641b9d016 100644 --- a/JugTrack/JugTrack/ProtoTrack.hpp +++ b/JugTrack/JugTrack/ProtoTrack.hpp @@ -11,7 +11,7 @@ #include <cstddef> #include <vector> -namespace FW { +namespace Jug { /// A proto track is a collection of hits identified by their indices. using ProtoTrack = std::vector<size_t>; diff --git a/JugTrack/JugTrack/SimHit.hpp b/JugTrack/JugTrack/SimHit.hpp deleted file mode 100644 index c17eebebbc5f6c06b9e20897250408dfb08d610c..0000000000000000000000000000000000000000 --- a/JugTrack/JugTrack/SimHit.hpp +++ /dev/null @@ -1,20 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2017 CERN for the benefit of the Acts project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#pragma once - -//#include "ACTFW/EventData/GeometryContainers.hpp" -//#include "ActsFatras/EventData/Hit.hpp" -// -//namespace FW { -// -//using SimHit = ::ActsFatras::Hit; -///// Store hits ordered by geometry identifier. -//using SimHitContainer = GeometryIdMultiset<::ActsFatras::Hit>; -// -//} // end of namespace FW diff --git a/JugTrack/JugTrack/SimIdentifier.hpp b/JugTrack/JugTrack/SimIdentifier.hpp deleted file mode 100644 index 1024479e460e881ec2fa12ee24105feb4236cd9c..0000000000000000000000000000000000000000 --- a/JugTrack/JugTrack/SimIdentifier.hpp +++ /dev/null @@ -1,81 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2016-2018 CERN for the benefit of the Acts project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#pragma once - -//#include "Acts/EventData/MeasurementHelpers.hpp" -// -//#include <cstddef> -//#include <vector> -// -//namespace FW { -// -///// A hit identifier with additional truth information. -///// -///// In an addition to a unique identifier, an list of additional indices is -///// stored. These can refer e.g. to particle indices or truth hit indices. -///// Using indices instead of pointers allows more flexibility, i.e. we are not -///// fixed to a specific object type when using e.g. pointers, and is more -///// robust since no requirements on stable memory locations of the pointed-to -///// objects are necessary (as would be the case for pointers). -//class SimIdentifier : public Acts::MinimalSourceLink { -// public: -// using Value = uint64_t; -// using Difference = int64_t; -// -// /// Constructor from encoded identifier value. -// /// -// /// @param value is the identifier value -// explicit SimIdentifier(Value value) : m_value(value) {} -// /// Constructor from encoded identifier value and truth information. -// /// -// /// @param value is the identifier value -// /// @param indices -// SimIdentifier(Value value, std::vector<std::size_t> indices) -// : m_value(value), m_indices(std::move(indices)) {} -// -// // Explicitely defaulted constructors and assignment operators -// SimIdentifier() = default; -// SimIdentifier(const SimIdentifier&) = default; -// SimIdentifier(SimIdentifier&&) = default; -// SimIdentifier& operator=(const SimIdentifier&) = default; -// SimIdentifier& operator=(SimIdentifier&&) = default; -// -// /// Assign from an identifier value. -// SimIdentifier& operator=(Value value); -// /// Cast to an identifier value. -// operator Value() const { return m_value; } -// -// /// Explicit access the underlying identifier value. -// Value value() const { return m_value; } -// /// Access all associated truth indices. -// const std::vector<std::size_t>& indices() const { return m_indices; } -// -// /// Attach a truth index to the identifier. -// void addIndex(std::size_t index) { m_indices.push_back(index); } -// -// private: -// /// The stored identifier value. -// Value m_value = 0u; -// /// Associated truth indices. -// std::vector<std::size_t> m_indices; -// -// friend constexpr bool operator<(const SimIdentifier& lhs, -// const SimIdentifier& rhs) { -// return lhs.m_value < rhs.m_value; -// } -// friend bool operator==(const SimIdentifier& lhs, const SimIdentifier& rhs) { -// return lhs.m_value == rhs.m_value; -// } -//}; -// -//} // end of namespace FW -// -//using identifier_type = ::FW::SimIdentifier::Value; -//using identifier_diff = ::FW::SimIdentifier::Difference; -//using Identifier = ::FW::SimIdentifier; diff --git a/JugTrack/JugTrack/SimParticle.hpp b/JugTrack/JugTrack/SimParticle.hpp deleted file mode 100644 index 64f9b9a8aaff5f1455bac202533c60c2bc1179dc..0000000000000000000000000000000000000000 --- a/JugTrack/JugTrack/SimParticle.hpp +++ /dev/null @@ -1,40 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2018 CERN for the benefit of the Acts project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#pragma once - -//#include "ActsFatras/EventData/Particle.hpp" -// -//#include <boost/container/flat_set.hpp> -// -//namespace FW { -//namespace detail { -//struct CompareParticleId { -// using is_transparent = void; -// constexpr bool operator()(const ActsFatras::Particle& lhs, -// const ActsFatras::Particle& rhs) const { -// return lhs.particleId() < rhs.particleId(); -// } -// constexpr bool operator()(ActsFatras::Barcode lhs, -// const ActsFatras::Particle& rhs) const { -// return lhs < rhs.particleId(); -// } -// constexpr bool operator()(const ActsFatras::Particle& lhs, -// ActsFatras::Barcode rhs) const { -// return lhs.particleId() < rhs; -// } -//}; -//} // namespace detail -// -//using SimParticle = ::ActsFatras::Particle; -///// Store particles ordered by particle identifier. -//using SimParticleContainer = -// ::boost::container::flat_set<::ActsFatras::Particle, -// detail::CompareParticleId>; -// -//} // end of namespace FW diff --git a/JugTrack/JugTrack/SimSourceLink.hpp b/JugTrack/JugTrack/SimSourceLink.hpp deleted file mode 100644 index 46e1c9a8c3d34b01e8898a15b5ff4fdbbfbed67e..0000000000000000000000000000000000000000 --- a/JugTrack/JugTrack/SimSourceLink.hpp +++ /dev/null @@ -1,85 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2016-2019 CERN for the benefit of the Acts project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#pragma once - -////#include "ACTFW/EventData/GeometryContainers.hpp" -//#include "Acts/EventData/Measurement.hpp" -//#include "ActsFatras/EventData/Hit.hpp" -// -//#include <stdexcept> -//#include <string> -// -//namespace FW { -// -///// Source link class for simulation in the acts-framework. -///// -///// The source link stores the measuremts, surface, and the associated simulated -///// truth hit. -///// -///// @todo Allow multiple truth hits e.g. for merged hits. -//class SimSourceLink { -// public: -// SimSourceLink(const Acts::Surface& surface, const ActsFatras::Hit& truthHit, -// size_t dim, Acts::BoundVector values, Acts::BoundMatrix cov) -// : m_values(values), -// m_cov(cov), -// m_dim(dim), -// m_geometryId(truthHit.geometryId()), -// m_surface(&surface), -// m_truthHit(&truthHit) {} -// /// Must be default_constructible to satisfy SourceLinkConcept. -// SimSourceLink() = default; -// SimSourceLink(SimSourceLink&&) = default; -// SimSourceLink(const SimSourceLink&) = default; -// SimSourceLink& operator=(SimSourceLink&&) = default; -// SimSourceLink& operator=(const SimSourceLink&) = default; -// -// constexpr Acts::GeometryIdentifier geometryId() const { return m_geometryId; } -// constexpr const Acts::Surface& referenceSurface() const { return *m_surface; } -// constexpr const ActsFatras::Hit& truthHit() const { return *m_truthHit; } -// -// Acts::FittableMeasurement<SimSourceLink> operator*() const { -// if (m_dim == 0) { -// throw std::runtime_error("Cannot create dim 0 measurement"); -// } else if (m_dim == 1) { -// return Acts::Measurement<SimSourceLink, Acts::BoundParametersIndices, -// Acts::eBoundLoc0>{ -// m_surface->getSharedPtr(), *this, m_cov.topLeftCorner<1, 1>(), -// m_values[0]}; -// } else if (m_dim == 2) { -// return Acts::Measurement<SimSourceLink, Acts::BoundParametersIndices, -// Acts::eBoundLoc0, Acts::eBoundLoc1>{ -// m_surface->getSharedPtr(), *this, m_cov.topLeftCorner<2, 2>(), -// m_values[0], m_values[1]}; -// } else { -// throw std::runtime_error("Dim " + std::to_string(m_dim) + -// " currently not supported."); -// } -// } -// -// private: -// Acts::BoundVector m_values; -// Acts::BoundMatrix m_cov; -// size_t m_dim = 0u; -// // store geo id copy to avoid indirection via truth hit -// Acts::GeometryIdentifier m_geometryId; -// // need to store pointers to make the object copyable -// const Acts::Surface* m_surface; -// const ActsFatras::Hit* m_truthHit; -// -// friend constexpr bool operator==(const SimSourceLink& lhs, -// const SimSourceLink& rhs) { -// return lhs.m_truthHit == rhs.m_truthHit; -// } -//}; -// -///// Store source links ordered by geometry identifier. -//using SimSourceLinkContainer = GeometryIdMultiset<SimSourceLink>; -// -//} // end of namespace FW diff --git a/JugTrack/JugTrack/SimVertex.hpp b/JugTrack/JugTrack/SimVertex.hpp deleted file mode 100644 index 8095d3d3dd1a432a27d50551c205388e52753153..0000000000000000000000000000000000000000 --- a/JugTrack/JugTrack/SimVertex.hpp +++ /dev/null @@ -1,56 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2018 CERN for the benefit of the Acts project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#pragma once - -//#include "Acts/Utilities/Definitions.hpp" -//#include "ActsFatras/EventData/Particle.hpp" -//#include "ActsFatras/EventData/ProcessType.hpp" -// -//#include <vector> -// -//namespace FW { -// -///// A simultated vertex e.g. from a physics process. -//struct SimVertex { -// using Scalar = double; -// using Vector4 = Acts::ActsVector<Scalar, 4>; -// -// /// The vertex four-position. -// Vector4 position4 = Vector4::Zero(); -// /// The vertex process type. -// ActsFatras::ProcessType process = ActsFatras::ProcessType::eUndefined; -// /// The incoming particles into the vertex -// std::vector<ActsFatras::Particle> incoming = {}; -// /// The outgoing particles from the vertex -// std::vector<ActsFatras::Particle> outgoing = {}; -// -// /// Construct the vertex from a position and optional process type. -// /// -// /// @param position4_ the vertex four-position -// /// @param process_ the process type that generated this vertex -// /// -// /// Associated particles are left empty by default and must be filled by the -// /// user after construction. -// SimVertex(const Vector4& position4_, ActsFatras::ProcessType process_ = -// ActsFatras::ProcessType::eUndefined) -// : position4(position4_), process(process_) {} -// // explicitely default rule-of-five. -// SimVertex() = default; -// SimVertex(const SimVertex&) = default; -// SimVertex(SimVertex&&) = default; -// SimVertex& operator=(const SimVertex&) = default; -// SimVertex& operator=(SimVertex&&) = default; -// -// /// The vertex three-position. -// auto position() const { return position4.head<3>(); } -// /// The vertex time. -// Scalar time() const { return position4[3]; } -//}; -// -//} // namespace FW diff --git a/JugTrack/src/components/GenFitTrackFitter.cpp b/JugTrack/src/components/GenFitTrackFitter.cpp new file mode 100644 index 0000000000000000000000000000000000000000..d526094c0312d70022fa0549d5b6f29c1722160e --- /dev/null +++ b/JugTrack/src/components/GenFitTrackFitter.cpp @@ -0,0 +1,164 @@ +#include "GenFitTrackFitter.h" +// Gaudi +#include "GaudiAlg/GaudiAlgorithm.h" +#include "GaudiKernel/ToolHandle.h" +#include "GaudiAlg/Transformer.h" +#include "GaudiAlg/GaudiTool.h" +#include "GaudiKernel/RndmGenerators.h" +#include "Gaudi/Property.h" + +#include "DDRec/CellIDPositionConverter.h" +#include "DDRec/SurfaceManager.h" +#include "DDRec/Surface.h" + +#include "JugBase/DataHandle.h" +#include "JugBase/IGeoSvc.h" +#include "JugTrack/GeometryContainers.hpp" +#include "JugTrack/IndexSourceLink.hpp" +#include "JugTrack/Track.hpp" +#include "JugTrack/BField.h" +#include "JugTrack/Measurement.hpp" + +#include "eicd/TrackerHitCollection.h" + +#include <functional> +#include <stdexcept> +#include <vector> +#include <random> +#include <stdexcept> + +//# genfit +#include "ConstField.h" +#include "DAF.h" +#include "Exception.h" +#include "FieldManager.h" +#include "KalmanFitterRefTrack.h" +#include "StateOnPlane.h" +#include "Track.h" +#include "TrackPoint.h" +#include "MaterialEffects.h" +#include "RKTrackRep.h" +#include "TGeoMaterialInterface.h" +//#include <EventDisplay.h> +#include "PlanarMeasurement.h" +#include "HelixTrackModel.h" +//#include "MeasurementCreator.h" +#include "WireMeasurement.h" + +#include "TDatabasePDG.h" +#include "TEveManager.h" +#include "TGeoManager.h" +#include "TRandom.h" +#include "TVector3.h" +#include <vector> + +namespace Jug::Reco { + + GenFitTrackFitter::GenFitTrackFitter(const std::string& name, ISvcLocator* svcLoc) + : GaudiAlgorithm(name, svcLoc) + { + declareProperty("inputHitCollection", m_inputHitCollection, ""); + declareProperty("initialTrackParameters", m_initialTrackParameters, ""); + declareProperty("inputProtoTracks", m_inputProtoTracks, ""); + declareProperty("foundTracks", m_foundTracks, ""); + declareProperty("outputTrajectories", m_outputTrajectories, ""); + } + + StatusCode GenFitTrackFitter::initialize() + { + if (GaudiAlgorithm::initialize().isFailure()) + return StatusCode::FAILURE; + m_geoSvc = service("GeoSvc"); + if (!m_geoSvc) { + error() << "Unable to locate Geometry Service. " + << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." + << endmsg; + return StatusCode::FAILURE; + } + + genfit::FieldManager::getInstance()->init(new genfit::ConstField( + 0., 0., m_geoSvc->centralMagneticField() * 10.0)); // gentfit uses kilo-Gauss + genfit::MaterialEffects::getInstance()->init(new genfit::TGeoMaterialInterface()); + + return StatusCode::SUCCESS; + } + + StatusCode GenFitTrackFitter::execute() + { + // Read input data + const eic::TrackerHitCollection* hits = m_inputHitCollection.get(); + const TrackParametersContainer* initialParameters = m_initialTrackParameters.get(); + const ProtoTrackContainer* protoTracks = m_inputProtoTracks.get(); + + const auto& single_track_param = (*initialParameters)[0]; + // TrajectoryContainer trajectories; + //auto trajectories = m_outputTrajectories.createAndPut(); + + debug() << "Single track mom : " << single_track_param.absoluteMomentum() << endmsg; + if( hits->size() < 2 ) { + return StatusCode::SUCCESS; + } + // init fitter + genfit::AbsKalmanFitter* fitter = new genfit::KalmanFitterRefTrack(); + + // particle pdg code; pion hypothesis + const int pdg = 11; + + // start values for the fit, e.g. from pattern recognition + TVector3 pos(0, 0, 0); + TVector3 mom(single_track_param.momentum()[0], single_track_param.momentum()[1], single_track_param.momentum()[2]); + + // trackrep + genfit::AbsTrackRep* rep = new genfit::RKTrackRep(pdg); + + // create track + genfit::Track fitTrack(rep, pos, mom); + + debug() << (*hits).size() << " hits " << endmsg; + int ihit = 0; + for (const auto& ahit : *hits) { + + //auto volman = m_geoSvc->detector()->volumeManager(); + //auto alignment = volman.lookupDetElement(vol_id).nominal(); + //auto local_position = alignment.worldToLocal(global_position); + TMatrixDSym hitCov(2); + hitCov.UnitMatrix(); + hitCov(0,0) = ahit.covsym_xx()*ahit.covsym_xx(); + hitCov(1,1) = ahit.covsym_yy()*ahit.covsym_yy(); + + TVector3 point = {ahit.x(), ahit.y(), ahit.z()}; + TVector3 u_dir = {1, 0, 0}; + u_dir.SetPhi(point.Phi()); + TVector3 v_dir = {0, 0, 1}; + + // add some planar hits to track with coordinates I just made up + TVectorD hitCoords(2); + hitCoords[0] = 0; + hitCoords[1] = 0; + auto measurement = new genfit::PlanarMeasurement(hitCoords, hitCov, ahit.cellID(), ihit, nullptr); + + measurement->setPlane(genfit::SharedPlanePtr(new genfit::DetPlane(point, u_dir, v_dir)), + ihit); + fitTrack.insertPoint(new genfit::TrackPoint(measurement, &fitTrack)); + + ihit++; + } + //check + fitTrack.checkConsistency(); + + // do the fit + fitter->processTrack(&fitTrack); + + // print fit result + fitTrack.getFittedState().Print(); + + //check + fitTrack.checkConsistency(); + + delete fitter; + + return StatusCode::SUCCESS; + } + + DECLARE_COMPONENT(GenFitTrackFitter) +} // namespace Jug::Reco diff --git a/JugTrack/src/components/GenFitTrackFitter.h b/JugTrack/src/components/GenFitTrackFitter.h new file mode 100644 index 0000000000000000000000000000000000000000..82cc181b0bb9e9fd6bc31741685dd3ea9a3ab866 --- /dev/null +++ b/JugTrack/src/components/GenFitTrackFitter.h @@ -0,0 +1,57 @@ +#ifndef JUGGLER_JUGRECO_GenFitTrackFitter_HH +#define JUGGLER_JUGRECO_GenFitTrackFitter_HH + +#include "JugTrack/GeometryContainers.hpp" + +// Gaudi +#include "GaudiAlg/GaudiAlgorithm.h" +#include "GaudiKernel/ToolHandle.h" +#include "Gaudi/Property.h" + +#include "JugBase/DataHandle.h" +#include "JugBase/IGeoSvc.h" + +#include <functional> +#include <stdexcept> +#include <vector> + +#include "JugTrack/IndexSourceLink.hpp" +#include "JugTrack/Track.hpp" +#include "JugTrack/BField.h" +#include "JugTrack/Measurement.hpp" +#include "JugTrack/Trajectories.hpp" +#include "JugTrack/ProtoTrack.hpp" + +#include "eicd/TrackerHitCollection.h" + +#include <random> +#include <stdexcept> + +namespace Jug::Reco { + + class GenFitTrackFitter : public GaudiAlgorithm { + public: + + public: + DataHandle<eic::TrackerHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; + DataHandle<TrackParametersContainer> m_initialTrackParameters{"initialTrackParameters", Gaudi::DataHandle::Reader, this}; + DataHandle<ProtoTrackContainer> m_inputProtoTracks{"inputProtoTracks", Gaudi::DataHandle::Reader, this}; + DataHandle<TrajectoriesContainer> m_foundTracks{"foundTracks", Gaudi::DataHandle::Reader, this}; + DataHandle<TrajectoriesContainer> m_outputTrajectories{"outputTrajectories", Gaudi::DataHandle::Writer, this}; + + SmartIF<IGeoSvc> m_geoSvc; + //Acts::GeometryContext m_geoctx; + //Acts::CalibrationContext m_calibctx; + //Acts::MagneticFieldContext m_fieldctx; + + + GenFitTrackFitter(const std::string& name, ISvcLocator* svcLoc); + + StatusCode initialize() override; + StatusCode execute() override; + }; + + +} // namespace Jug::Reco + +#endif diff --git a/JugTrack/src/components/ParticlesFromTrackFit.cpp b/JugTrack/src/components/ParticlesFromTrackFit.cpp index a2811bbfd43c151c22e16b91f437e3ebf3f1b75d..ab3eb76b606ba296a3be221a60212be1b9dfc338 100644 --- a/JugTrack/src/components/ParticlesFromTrackFit.cpp +++ b/JugTrack/src/components/ParticlesFromTrackFit.cpp @@ -70,7 +70,6 @@ namespace Jug { for (size_t itraj = 0; itraj < trajectories->size(); ++itraj) { const auto& traj = (*trajectories)[itraj]; - // Get the entry index for the single trajectory // The trajectory entry indices and the multiTrajectory const auto& mj = traj.multiTrajectory(); diff --git a/JugTrack/src/components/SingleTrackSourceLinker.cpp b/JugTrack/src/components/SingleTrackSourceLinker.cpp new file mode 100644 index 0000000000000000000000000000000000000000..356034d364d4f41f50075bbc2a0a691e05bf1a72 --- /dev/null +++ b/JugTrack/src/components/SingleTrackSourceLinker.cpp @@ -0,0 +1,182 @@ +#include "JugTrack/GeometryContainers.hpp" + +// Gaudi +#include "GaudiAlg/GaudiAlgorithm.h" +#include "GaudiKernel/ToolHandle.h" +#include "GaudiAlg/Transformer.h" +#include "GaudiAlg/GaudiTool.h" +#include "GaudiKernel/RndmGenerators.h" +#include "Gaudi/Property.h" + +#include "JugBase/DataHandle.h" +#include "JugBase/IGeoSvc.h" + +#include "DDRec/CellIDPositionConverter.h" +#include "DDRec/SurfaceManager.h" +#include "DDRec/Surface.h" +#include "DD4hep/Volumes.h" +#include "DD4hep/DD4hepUnits.h" + + +#include "Acts/Surfaces/Surface.hpp" +#include "Acts/Definitions/Units.hpp" +#include "Acts/Definitions/Common.hpp" +#include "Acts/Geometry/TrackingGeometry.hpp" +#include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp" + +#include "JugTrack/Index.hpp" +#include "JugTrack/IndexSourceLink.hpp" +#include "JugTrack/Measurement.hpp" +#include "JugTrack/ProtoTrack.hpp" + +#include "eicd/TrackerHitCollection.h" + +namespace Jug::Reco { + + /** Single Track source Linker and proto tracks. + * This algorithm assumes only single track events. + * + */ + class SingleTrackSourceLinker : public GaudiAlgorithm { + public: + DataHandle<eic::TrackerHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; + DataHandle<IndexSourceLinkContainer> m_outputSourceLinks{"outputSourceLinks", Gaudi::DataHandle::Writer, this}; + DataHandle<MeasurementContainer> m_outputMeasurements{"outputMeasurements", Gaudi::DataHandle::Writer, this}; + DataHandle<ProtoTrackContainer> m_outputProtoTracks{"outputProtoTracks", Gaudi::DataHandle::Writer, this}; + /// Pointer to the geometry service + SmartIF<IGeoSvc> m_geoSvc; + + /// Lookup container for hit surfaces that generate smeared hits + std::unordered_map<uint64_t, const Acts::Surface*> m_surfaces; + + public: + SingleTrackSourceLinker(const std::string& name, ISvcLocator* svcLoc) + : GaudiAlgorithm(name, svcLoc) { + declareProperty("inputHitCollection", m_inputHitCollection, ""); + declareProperty("outputSourceLinks", m_outputSourceLinks, ""); + declareProperty("outputMeasurements", m_outputMeasurements, ""); + declareProperty("outputProtoTracks", m_outputProtoTracks, ""); + } + + StatusCode initialize() override { + if (GaudiAlgorithm::initialize().isFailure()) + return StatusCode::FAILURE; + m_geoSvc = service("GeoSvc"); + if (!m_geoSvc) { + error() << "Unable to locate Geometry Service. " + << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." + << endmsg; + return StatusCode::FAILURE; + } + return StatusCode::SUCCESS; + } + + StatusCode execute() override { + // input collection + const eic::TrackerHitCollection* hits = m_inputHitCollection.get(); + // Create output collections + auto sourceLinks = m_outputSourceLinks.createAndPut(); + auto measurements = m_outputMeasurements.createAndPut(); + auto protoTracks = m_outputProtoTracks.createAndPut(); + // IndexMultimap<ActsFatras::Barcode> hitParticlesMap; + // IndexMultimap<Index> hitSimHitsMap; + sourceLinks->reserve(hits->size()); + measurements->reserve(hits->size()); + + + // assume single track --> one ProtoTrack + ProtoTrack track; + track.reserve((*hits).size()); + + debug() << (*hits).size() << " hits " << endmsg; + int ihit = 0; + for(const auto& ahit : *hits) { + + track.emplace_back(ihit); + + auto vol_ctx = m_geoSvc->cellIDPositionConverter()->findContext(ahit.cellID()); + auto vol_id = vol_ctx->identifier; + const auto is = m_geoSvc->surfaceMap().find(vol_id); + if (is == m_surfaces.end()) { + debug() << " vol_id (" << vol_id << ") not found in m_surfaces!!!!" <<endmsg; + continue; + } + const Acts::Surface* surface = is->second; + + // NOTE + // Here is where the important hit and tracking geometry is connected. + // A "Measurement" is constructed to for each hit which makes the connection to + // the tracking surface and covariance matrix + + dd4hep::Position global_position(ahit.x(), ahit.y(), ahit.z()); + + auto volman = m_geoSvc->detector()->volumeManager(); + auto alignment = volman.lookupDetElement(vol_id).nominal(); + auto local_position = alignment.worldToLocal(global_position); + + debug() << "===== Debugging hit =====" << endmsg; + debug() << "global pos (" << global_position.x() << "," << global_position.y() << "," + << global_position.z() << ")" << endmsg; + //debug() << "local pos (" << local_position.x() << "," << local_position.y() << "," + // << local_position.z() << ")" << endmsg; + auto acts_pos = surface->globalToLocal(Acts::GeometryContext(), {ahit.x(), ahit.y(), ahit.z()}, {0, 0, 0}).value();//, pos); + debug() << " ACTS local position : (" << acts_pos[0] << "," << acts_pos[1] << "," << acts_pos[2] << ")"<< endmsg; + // construct the vector of measured parameters (2d position in this case) + Acts::Vector2 pos(acts_pos.x(), acts_pos.y()); + + // construct the covariance matrix + Acts::SymMatrix2 cov = Acts::SymMatrix2::Zero(); + cov(0,0) = ahit.covsym_xx()*Acts::UnitConstants::mm*ahit.covsym_xx()*Acts::UnitConstants::mm; + cov(1,1) = ahit.covsym_yy()*Acts::UnitConstants::mm*ahit.covsym_yy()*Acts::UnitConstants::mm; + + // Above we only consider the two position coordinates the comment below shows how to add time + // which we will probably want to try later. + // + //Acts::SymMatrix3 cov; + //cov << 0.05, 0., 0., 0., 0.05, 0., 0., 0., 900. * Acts::UnitConstants::ps * Acts::UnitConstants::ps; + //Acts::Vector3 par(localX, localY, simHit.time()); + + Index hitIdx = ihit; + IndexSourceLink sourceLink(vol_id, ihit); + auto meas = Acts::makeMeasurement(sourceLink, pos, cov, + Acts::eBoundLoc0, Acts::eBoundLoc1);//, Acts::eBoundTime); + + // TODO: check that local to global is the same in dd4hep and acts. + // + //debug() << "ACTS surface center : " << surface->center(Acts::GeometryContext()) << endmsg; + // transform global position into local coordinates + // geometry context contains nothing here + //Acts::Vector2 loc = Acts::Vector2::Zero(); + //loc[Acts::eBoundLoc0] = pos[0] ;//+ m_cfg.sigmaLoc0 * stdNormal(rng); + //loc[Acts::eBoundLoc1] = pos[1] ;//+ m_cfg.sigmaLoc1 * stdNormal(rng); + ////debug() << "loc : (" << loc[0] << ", " << loc[1] << ")" << endmsg; + + ////local position + ////auto loc = {ahit.x(), ahit.y(), ahit.z()} - vol_ctx->volumePlacement().position() + ////debug() << " hit : \n" << ahit << endmsg; + ////debug() << " cell ID : " << ahit.cellID() << endmsg; + ////debug() << " position : (" << ahit.position(0) << ", " << ahit.position(1) << ", "<< ahit.position(2) << ") " << endmsg; + ////debug() << " vol_id : " << vol_id << endmsg; + ////debug() << " placment pos : " << vol_ctx->volumePlacement().position() << endmsg; + + //// the measurement container is unordered and the index under which the + //// measurement will be stored is known before adding it. + //auto meas = Acts::makeMeasurement(sourceLink, loc, cov, Acts::eBoundLoc0, Acts::eBoundLoc1); + + // add to output containers. since the input is already geometry-order, + // new elements in geometry containers can just be appended at the end. + sourceLinks->emplace_hint(sourceLinks->end(), std::move(sourceLink)); + measurements->emplace_back(std::move(meas)); + + ihit++; + } + // add proto track to the output collection + protoTracks->emplace_back(std::move(track)); + return StatusCode::SUCCESS; + } + + }; + DECLARE_COMPONENT(SingleTrackSourceLinker) + +} // namespace Jug::reco + diff --git a/JugTrack/src/components/TrackFittingAlgorithm.cpp b/JugTrack/src/components/TrackFittingAlgorithm.cpp index 5e12f2f08b08d395cee51ab0ca50277c69a9d2ab..c16bbf9f2780229d68bf58ec2f409c0f0312d58e 100644 --- a/JugTrack/src/components/TrackFittingAlgorithm.cpp +++ b/JugTrack/src/components/TrackFittingAlgorithm.cpp @@ -53,6 +53,7 @@ namespace Jug::Reco { declareProperty("inputSourceLinks", m_inputSourceLinks, ""); declareProperty("initialTrackParameters", m_initialTrackParameters, ""); declareProperty("inputMeasurements", m_inputMeasurements, ""); + declareProperty("inputProtoTracks", m_inputProtoTracks, ""); declareProperty("foundTracks", m_foundTracks, ""); declareProperty("outputTrajectories", m_outputTrajectories, ""); } @@ -80,46 +81,79 @@ namespace Jug::Reco { StatusCode TrackFittingAlgorithm::execute() { // Read input data - const IndexSourceLinkContainer* src_links = m_inputSourceLinks.get(); - const TrackParametersContainer* init_trk_params = m_initialTrackParameters.get(); - const MeasurementContainer* measurements = m_inputMeasurements.get(); + const IndexSourceLinkContainer* sourceLinks = m_inputSourceLinks.get(); + const TrackParametersContainer* initialParameters = m_initialTrackParameters.get(); + const MeasurementContainer* measurements = m_inputMeasurements.get(); + const ProtoTrackContainer* protoTracks = m_inputProtoTracks.get(); + ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("TrackFittingAlgorithm Logger", Acts::Logging::INFO)); + + // Consistency cross checks + if (protoTracks->size() != initialParameters->size()) { + ACTS_FATAL("Inconsistent number of proto tracks and initial parameters"); + return StatusCode::FAILURE; + } // TrajectoryContainer trajectories; auto trajectories = m_outputTrajectories.createAndPut(); - trajectories->reserve(init_trk_params->size()); + trajectories->reserve(initialParameters->size()); - //// Construct a perigee surface as the target surface + // Construct a perigee surface as the target surface auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3{0., 0., 0.}); - ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("TrackFittingAlgorithm Logger", Acts::Logging::INFO)); - std::vector<IndexSourceLink> trackSourceLinks; + //Acts::KalmanFitterOptions<MeasurementCalibrator, Acts::VoidOutlierFinder> kfOptions( + // ctx.geoContext, ctx.magFieldContext, ctx.calibContext, MeasurementCalibrator(measurements), + // Acts::VoidOutlierFinder(), Acts::LoggerWrapper{logger()}, Acts::PropagatorPlainOptions(), + // &(*pSurface)); + + // kfOptions.multipleScattering = m_cfg.multipleScattering; + // kfOptions.energyLoss = m_cfg.energyLoss; + // Acts::KalmanFitterOptions<Acts::VoidOutlierFinder> kfOptions( + // m_geoctx, m_fieldctx, m_calibctx, Acts::VoidOutlierFinder(), + // Acts::LoggerWrapper{logger()}, Acts::PropagatorPlainOptions(), &(*pSurface)); + // // Set the KalmanFitter options + + Acts::PropagatorPlainOptions pOptions; + pOptions.maxSteps = 10000; + + Acts::KalmanFitterOptions<MeasurementCalibrator, Acts::VoidOutlierFinder> kfOptions( + m_geoctx, m_fieldctx, m_calibctx, MeasurementCalibrator(*measurements), + Acts::VoidOutlierFinder(), Acts::LoggerWrapper{logger()}, pOptions, &(*pSurface)); + + // used for processing the data + std::vector<IndexSourceLink> trackSourceLinks; + std::vector<const Acts::Surface*> surfSequence; + + debug() << "initialParams size: " << initialParameters->size() << endmsg; + debug() << "measurements size: " << measurements->size() << endmsg; // Perform the track finding for each starting parameter // @TODO: use seeds from track seeding algorithm as starting parameter - // initial track params and proto tracks might likely have the same size. - for (std::size_t iseed = 0; iseed < init_trk_params->size(); ++iseed) { + // initial track params and proto tracks might likely have the same size. + //for (std::size_t iseed = 0; iseed < init_trk_params->size(); ++iseed) { + for (std::size_t itrack = 0; itrack < (*protoTracks).size(); ++itrack) { - // this will eventually be per-track. for now we assume all src links are from the same single track! - for (auto& lnk : (*src_links)) { - // auto sourceLink = sourceLinks.nth(hitIndex); - trackSourceLinks.push_back(lnk); - } + const auto& protoTrack = (*protoTracks)[itrack]; + const auto& initialParams = (*initialParameters)[itrack]; + debug() << "protoTrack size: " << protoTrack.size() << endmsg; - const auto& initialParams = (*init_trk_params)[iseed]; - Acts::PropagatorPlainOptions pOptions; - pOptions.maxSteps = 10000; + trackSourceLinks.clear(); + trackSourceLinks.reserve(protoTrack.size()); - //Acts::KalmanFitterOptions<Acts::VoidOutlierFinder> kfOptions( - // m_geoctx, m_fieldctx, m_calibctx, Acts::VoidOutlierFinder(), - // Acts::LoggerWrapper{logger()}, Acts::PropagatorPlainOptions(), &(*pSurface)); - // // Set the KalmanFitter options - Acts::KalmanFitterOptions<MeasurementCalibrator, Acts::VoidOutlierFinder> kfOptions( - m_geoctx, m_fieldctx, m_calibctx, - MeasurementCalibrator(*measurements), Acts::VoidOutlierFinder(), - Acts::LoggerWrapper{logger()}, Acts::PropagatorPlainOptions(), &(*pSurface)); + for (auto hitIndex : protoTrack) { + debug() << " hit index = " << hitIndex << endmsg; + auto sourceLink = sourceLinks->nth(hitIndex); + auto geoId = sourceLink->geometryId(); + if (sourceLink == sourceLinks->end()) { + ACTS_FATAL("Proto track " << itrack << " contains invalid hit index" + << hitIndex); + return StatusCode::FAILURE; + } + trackSourceLinks.push_back(*sourceLink); + //surfSequence.push_back(m_cfg.trackingGeometry->findSurface(geoId)); + } - debug() << "Invoke track fitting ... " << iseed << endmsg; + debug() << "Invoke track fitting ... " << itrack << endmsg; auto result = fitTrack(trackSourceLinks, initialParams, kfOptions); debug() << "fitting done." << endmsg; // if (result.ok()) { @@ -159,7 +193,7 @@ namespace Jug::Reco { trajectories->emplace_back(std::move(fitOutput.fittedStates), std::move(trackTips), std::move(indexedParams)); } else { - ACTS_WARNING("Fit failed for track " << iseed << " with error" << result.error()); + ACTS_WARNING("Fit failed for track " << itrack << " with error" << result.error()); // Fit failed. Add an empty result so the output container has // the same number of entries as the input. trajectories->push_back(Trajectories()); diff --git a/JugTrack/src/components/TrackFittingAlgorithm.h b/JugTrack/src/components/TrackFittingAlgorithm.h index 567dc4f408bd30631edb877f7a9a4fc886f940bc..a9593a3dbdf3212eeec4f05416f18e51667ba7ec 100644 --- a/JugTrack/src/components/TrackFittingAlgorithm.h +++ b/JugTrack/src/components/TrackFittingAlgorithm.h @@ -29,6 +29,7 @@ #include "JugTrack/BField.h" #include "JugTrack/Measurement.hpp" #include "JugTrack/Trajectories.hpp" +#include "JugTrack/ProtoTrack.hpp" #include "eicd/TrackerHitCollection.h" @@ -73,11 +74,12 @@ namespace Jug::Reco { const std::vector<IndexSourceLink>&, const TrackParameters&, const TrackFitterOptions&)>; public: - DataHandle<IndexSourceLinkContainer> m_inputSourceLinks{"inputSourceLinks", Gaudi::DataHandle::Reader, this}; + DataHandle<IndexSourceLinkContainer> m_inputSourceLinks{"inputSourceLinks", Gaudi::DataHandle::Reader, this}; DataHandle<TrackParametersContainer> m_initialTrackParameters{"initialTrackParameters", Gaudi::DataHandle::Reader, this}; DataHandle<MeasurementContainer> m_inputMeasurements{"inputMeasurements", Gaudi::DataHandle::Reader, this}; - DataHandle<TrajectoriesContainer> m_foundTracks{"foundTracks", Gaudi::DataHandle::Reader, this}; - DataHandle<TrajectoriesContainer> m_outputTrajectories{"outputTrajectories", Gaudi::DataHandle::Writer, this}; + DataHandle<ProtoTrackContainer> m_inputProtoTracks{"inputProtoTracks", Gaudi::DataHandle::Reader, this}; + DataHandle<TrajectoriesContainer> m_foundTracks{"foundTracks", Gaudi::DataHandle::Reader, this}; + DataHandle<TrajectoriesContainer> m_outputTrajectories{"outputTrajectories", Gaudi::DataHandle::Writer, this}; FitterFunction m_trackFittingFunc; SmartIF<IGeoSvc> m_geoSvc; @@ -112,10 +114,10 @@ namespace Jug::Reco { //, const std::vector<const Acts::Surface*>& surfSequence) const; }; - inline TrackFittingAlgorithm::FitterResult TrackFittingAlgorithm::fitTrack( - const std::vector<IndexSourceLink>& sourceLinks, const TrackParameters& initialParameters, - const Acts::KalmanFitterOptions<MeasurementCalibrator, Acts::VoidOutlierFinder>& options) - const + inline TrackFittingAlgorithm::FitterResult + TrackFittingAlgorithm::fitTrack(const std::vector<IndexSourceLink>& sourceLinks, + const TrackParameters& initialParameters, + const TrackFitterOptions& options) const { // const std::vector<const Acts::Surface*>& surfSequence) const // if (m_cfg.directNavigation) { diff --git a/JugTrack/src/components/TrackFittingFunction.cpp b/JugTrack/src/components/TrackFittingFunction.cpp index 53abe7ac138c90c2a2a2e0303bf411fb7e749eae..584e9c392dbeea78c7b0964a145b6baa013b5ce0 100644 --- a/JugTrack/src/components/TrackFittingFunction.cpp +++ b/JugTrack/src/components/TrackFittingFunction.cpp @@ -18,6 +18,15 @@ namespace { + using Updater = Acts::GainMatrixUpdater; + using Smoother = Acts::GainMatrixSmoother; + using Stepper = Acts::EigenStepper<>; + using Propagator = Acts::Propagator<Stepper, Acts::Navigator>; + using Fitter = Acts::KalmanFitter<Propagator, Updater, Smoother>; + using DirectPropagator = Acts::Propagator<Stepper, Acts::DirectNavigator>; + using DirectFitter = Acts::KalmanFitter<DirectPropagator, Updater, Smoother>; + + template <typename track_fitter_t> struct TrackFitterFunctionImpl { track_fitter_t trackFitter; diff --git a/JugTrack/src/components/TrackParamTruthInit.cpp b/JugTrack/src/components/TrackParamTruthInit.cpp index 55dca3869b7e7ef5cc8304ca974150cdccd795d1..292801ced6a761de72ba6131cc4158f64053f625 100644 --- a/JugTrack/src/components/TrackParamTruthInit.cpp +++ b/JugTrack/src/components/TrackParamTruthInit.cpp @@ -83,19 +83,19 @@ namespace Jug::Reco { // build some track cov matrix Acts::BoundSymMatrix cov = Acts::BoundSymMatrix::Zero(); - cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = 25*um*25*um; - cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = 100*um*100*um; - cov(Acts::eBoundPhi, Acts::eBoundPhi) = 0.005*0.005; - cov(Acts::eBoundTheta, Acts::eBoundTheta) = 0.001*0.001; + cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = 1000*um*1000*um; + cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = 1000*um*1000*um; + cov(Acts::eBoundPhi, Acts::eBoundPhi) = 0.05*0.05; + cov(Acts::eBoundTheta, Acts::eBoundTheta) = 0.01*0.01; cov(Acts::eBoundQOverP, Acts::eBoundQOverP) = (0.1*0.1) / (GeV*GeV); - cov(Acts::eBoundTime, Acts::eBoundTime) = 1.0e9*ns*1.0e9*ns; + cov(Acts::eBoundTime, Acts::eBoundTime) = 10.0e9*ns*10.0e9*ns; Acts::BoundVector params; params(Acts::eBoundLoc0) = 0.0 * mm ; // cylinder radius params(Acts::eBoundLoc1) = 0.0 * mm ; // cylinder length params(Acts::eBoundPhi) = momentum.Phi(); params(Acts::eBoundTheta) = momentum.Theta(); - params(Acts::eBoundQOverP) = 1/p; + params(Acts::eBoundQOverP) = -1/p; params(Acts::eBoundTime) = part.time() * ns; /// \todo create or find better particle data interface. // get the particle charge @@ -110,7 +110,7 @@ namespace Jug::Reco { Acts::Vector3{part.vsx() * mm, part.vsy() * mm, part.vsz() * mm}); //params(Acts::eBoundQOverP) = charge/p; - init_trk_params->push_back({pSurface, params, charge}); + init_trk_params->push_back({pSurface, params, charge,cov}); // std::make_optional(std::move(cov)) debug() << "Invoke track finding seeded by truth particle with p = " << p/GeV << " GeV" << endmsg;