Commit 0a8cdc20 authored by Whitney Armstrong's avatar Whitney Armstrong
Browse files

Merging in new genfit testing and acts changes.

parent 93d2422c
......@@ -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
......
......@@ -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() {}
};
......
......@@ -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;
}
......
......@@ -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
......@@ -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)
......
......@@ -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>;
......
// 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
// 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;
// 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
// 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
// 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
#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();