diff --git a/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.cxx b/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.cxx index ff163e105d51b02267ce35e3fdfe43737c5a184c..e91386ff5f5f229b0fb2c5e0a52474c4155a013c 100644 --- a/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.cxx +++ b/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.cxx @@ -507,16 +507,19 @@ SegmentFitAlg::selectFits(const FitMap& fits) const auto content = *it; ATH_MSG_VERBOSE("clusters: " << content->clusterMask.count() << " / chi2: " << content->fitChi2) ; } - ClusterSet usedClusters { clusterInfo::nClusters }; while (info.size() > 0) { auto selected = info.front(); ATH_MSG_VERBOSE("Selected nclust = " << selected->clusterMask.count() << " with chi2 = " << selected->fitChi2); selectedFits[selected->clusterMask] = selected; - usedClusters |= selected->clusterMask; - info.remove_if([&](std::shared_ptr<fitInfo> p) {return (p->clusterMask & ~usedClusters) != p->clusterMask;}); + // remove fits for which the fraction of shared hits is larger than m_sharedFraction or the number of hits is + // smaller than m_minClustersPerFit + info.remove_if([&](std::shared_ptr<fitInfo> p) {return + ((p->clusterMask & selected->clusterMask).count() > m_sharedHitFraction * selected->clusterMask.count()) + || (p->clusterMask.count() < m_minClustersPerFit) ;} + ); } m_numberOfSelectedFits += selectedFits.size(); diff --git a/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.h b/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.h index eb41946cd0a9094a24919eaf6824705094afa234..757eb1a22752acbe72e0da7a0c6c7a89878bf76a 100644 --- a/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.h +++ b/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.h @@ -414,6 +414,8 @@ class SegmentFitAlg : public AthReentrantAlgorithm DoubleProperty m_yResidualCut { this, "ResidualCut", 3 * Gaudi::Units::mm, "Cluster y tolerance for compatibility with extrapolated fit." }; DoubleProperty m_tanThetaCut { this, "TanThetaCut", 0.0, "Maximum accepted tangent of track angle from beam axis; 0 means no cut."}; DoubleProperty m_reducedChi2Cut { this, "ReducedChi2Cut", 10.0, "Maximum accepted chi^2 per degree of freedom for final fits; 0 means no cut." }; + DoubleProperty m_sharedHitFraction { this, "SharedHitFraction", -1., "Fraction of hits which are allowed to be shared between two fits." }; + IntegerProperty m_minClustersPerFit { this, "MinClustersPerFit", 4, "Minimum number of clusters a fit has to have." }; mutable std::atomic<int> m_numberOfEvents{0}; mutable std::atomic<int> m_numberExcessOccupancy{0}; diff --git a/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt b/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt index 32633498acaf4d4d2dd866883d557336d02a6abd..d7801a45f716915b6d2c82773551e51d73991307 100755 --- a/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt +++ b/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt @@ -6,8 +6,6 @@ find_package( CLHEP ) find_package( Eigen ) find_package( Boost ) find_package( Acts COMPONENTS Core ) -set( Acts_DIR /home/tboeckh/Documents/acts/run ) -find_package( Acts REQUIRED COMPONENTS Core PATHS /home/tboeckh/Documents/acts/run NO_DEFAULT_PATH ) # Component(s) in the package: atlas_add_library( FaserActsKalmanFilterLib @@ -41,7 +39,6 @@ atlas_add_component(FaserActsKalmanFilter # FaserActsKalmanFilter/TruthTrackFinderTool.h FaserActsKalmanFilter/Measurement.h # FaserActsKalmanFilter/MultiTrackFinderTool.h - FaserActsKalmanFilter/MyAmbiguitySolver.h FaserActsKalmanFilter/PerformanceWriterTool.h FaserActsKalmanFilter/PlotHelpers.h FaserActsKalmanFilter/ResPlotTool.h diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CombinatorialKalmanFilterAlg.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CombinatorialKalmanFilterAlg.h index a5cdb8048c98db22d685681c48fad156566e6b7d..3044a2fcc24068da110975f06415608752d38cdc 100644 --- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CombinatorialKalmanFilterAlg.h +++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/CombinatorialKalmanFilterAlg.h @@ -20,6 +20,10 @@ #include "FaserActsKalmanFilter/RootTrajectorySummaryWriterTool.h" #include "FaserActsKalmanFilter/PerformanceWriterTool.h" #include "FaserActsKalmanFilter/TrackSeedWriterTool.h" +#include "FaserActsKalmanFilter/FaserActsRecMultiTrajectory.h" +#include <boost/dynamic_bitset.hpp> +using ConstTrackStateProxy = Acts::detail_lt::TrackStateProxy<IndexSourceLink, 6, true>; +using ClusterSet = boost::dynamic_bitset<>; class FaserSCT_ID; @@ -45,8 +49,8 @@ public: Acts::CombinatorialKalmanFilterOptions<IndexSourceLinkAccessor, MeasurementCalibrator, Acts::MeasurementSelector>; - using TrackFitterResult = - Acts::Result<Acts::CombinatorialKalmanFilterResult<IndexSourceLink>>; + using CKFResult = Acts::CombinatorialKalmanFilterResult<IndexSourceLink>; + using TrackFitterResult = Acts::Result<CKFResult>; using TrackFinderResult = std::vector<TrackFitterResult>; using TrackParameters = Acts::CurvilinearTrackParameters; using TrackParametersContainer = std::vector<TrackParameters>; @@ -59,6 +63,29 @@ public: const TrackFinderOptions&) const = 0; }; + struct TrajectoryInfo { + TrajectoryInfo(const FaserActsRecMultiTrajectory &traj) : + trajectory{traj}, clusterSet{nClusters} { + auto state = Acts::MultiTrajectoryHelpers::trajectoryState(traj.multiTrajectory(), traj.tips().front()); + traj.multiTrajectory().visitBackwards(traj.tips().front(), [&](const ConstTrackStateProxy& state) { + auto typeFlags = state.typeFlags(); + if (not typeFlags.test(Acts::TrackStateFlag::MeasurementFlag)) { + return true; + } + clusterSet.set(state.uncalibrated().index()); + return true; + }); + nMeasurements = state.nMeasurements; + chi2 = state.chi2Sum; + } + + static size_t nClusters; + FaserActsRecMultiTrajectory trajectory; + ClusterSet clusterSet; + size_t nMeasurements; + double chi2; + }; + static std::shared_ptr<TrackFinderFunction> makeTrackFinderFunction( std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry, bool resolvePassive, bool resolveMaterial, bool resolveSensitive); @@ -88,7 +115,6 @@ public: Gaudi::Property<double> m_maxSteps {this, "maxSteps", 10000}; Gaudi::Property<double> m_chi2Max {this, "chi2Max", 15}; Gaudi::Property<unsigned long> m_nMax {this, "nMax", 10}; - Gaudi::Property<size_t> m_nTrajectories {this, "nTrajectories", 2}; SG::ReadCondHandleKey<FaserFieldCacheCondObj> m_fieldCondObjInputKey {this, "FaserFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"}; ToolHandle<ITrackSeedTool> m_trackSeedTool {this, "TrackSeed", "ClusterTrackSeedTool"}; ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool {this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"}; @@ -97,9 +123,10 @@ public: ToolHandle<RootTrajectorySummaryWriterTool> m_trajectorySummaryWriterTool {this, "RootTrajectorySummaryWriterTool", "RootTrajectorySummaryWriterTool"}; ToolHandle<TrackSeedWriterTool> m_trackSeedWriterTool {this, "TrackSeedWriterTool", "TrackSeedWriterTool"}; - std::unique_ptr<Trk::Track> makeTrack(Acts::GeometryContext& tgContext, TrackFitterResult& fitResult) const; - const Trk::TrackParameters* ConvertActsTrackParameterToATLAS(const Acts::BoundTrackParameters &actsParameter, const Acts::GeometryContext& gctx) const; - SG::WriteHandleKey<TrackCollection> m_trackCollection { this, "CKFTrackCollection", "CKFTrackCollection" }; + std::unique_ptr<Trk::Track> makeTrack(Acts::GeometryContext& tgContext, TrackFitterResult& fitResult) const; + std::unique_ptr<Trk::Track> makeTrack(const Acts::GeometryContext &geoCtx, const FaserActsRecMultiTrajectory &traj) const; + const Trk::TrackParameters* ConvertActsTrackParameterToATLAS(const Acts::BoundTrackParameters &actsParameter, const Acts::GeometryContext& gctx) const; + SG::WriteHandleKey<TrackCollection> m_trackCollection { this, "CKFTrackCollection", "CKFTrackCollection" }; }; #endif // COMBINATORIALKALMANFILTERALG_H diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/MyAmbiguitySolver.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/MyAmbiguitySolver.h deleted file mode 100644 index 65c8dc6795a96397fa2e6a77e9a14b690b4d58ca..0000000000000000000000000000000000000000 --- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/MyAmbiguitySolver.h +++ /dev/null @@ -1,82 +0,0 @@ -#ifndef FASERACTSKALMANFILTER_AMBIGUITYSOLVER_H -#define FASERACTSKALMANFILTER_AMBIGUITYSOLVER_H - -/* -#include "Acts/TrackFinding/CombinatorialKalmanFilter.hpp" -#include "FaserActsKalmanFilter/FaserActsRecMultiTrajectory.h" - -using CombinatorialKalmanFilterResult = Acts::CombinatorialKalmanFilterResult<IndexSourceLink>; -using TrackFitterResult = Acts::Result<CombinatorialKalmanFilterResult>; -using TrackFinderResult = std::vector<TrackFitterResult>; - - -size_t numberMeasurements(const CombinatorialKalmanFilterResult& ckfResult) { - auto traj = FaserActsRecMultiTrajectory(ckfResult.fittedStates, ckfResult.lastMeasurementIndices, ckfResult.fittedParameters); - const auto& mj = traj.multiTrajectory(); - const auto& trackTips = traj.tips(); - size_t maxMeasurements = 0; - for (const auto& trackTip : trackTips) { - auto trajState = Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip); - size_t nMeasurements = trajState.nMeasurements; - if (nMeasurements > maxMeasurements) { - maxMeasurements = nMeasurements; - } - std::cout << "# measurements: " << trajState.nMeasurements << std::endl; - } - return maxMeasurements; -} - -int countSharedHits(const CombinatorialKalmanFilterResult& result1, const CombinatorialKalmanFilterResult& result2) { - int count = 0; - std::vector<size_t> hitIndices {}; - - for (auto measIndex : result1.lastMeasurementIndices) { - result1.fittedStates.visitBackwards(measIndex, [&](const auto& state) { - if (not state.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag)) - return; - size_t hitIndex = state.uncalibrated().index(); - hitIndices.emplace_back(hitIndex); - }); - } - - for (auto measIndex : result2.lastMeasurementIndices) { - result2.fittedStates.visitBackwards(measIndex, [&](const auto& state) { - if (not state.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag)) - return; - size_t hitIndex = state.uncalibrated().index(); - if (std::find(hitIndices.begin(), hitIndices.end(), hitIndex) != hitIndices.end()) { - count += 1; - } - }); - } - return count; -} - - -std::pair<int, int> solveAmbiguity(TrackFinderResult& results, size_t minMeasurements = 13) { - std::map<std::pair<size_t, size_t>, size_t> trackPairs {}; - for (size_t i = 0; i < results.size(); ++i) { - // if (not results.at(i).ok()) continue; - // if (numberMeasurements(results.at(i).value()) < minMeasurements) continue; - for (size_t j = i+1; j < results.size(); ++j) { - // if (not results.at(j).ok()) continue; - // if (numberMeasurements(results.at(j).value()) < minMeasurements) continue; - int n = countSharedHits(results.at(i).value(), results.at(j).value()); - trackPairs[std::make_pair(i, j)] = n; - } - } - - std::pair<size_t, size_t> bestTrackPair; - size_t minSharedHits = std::numeric_limits<size_t>::max(); - for (const auto& trackPair : trackPairs) { - if (trackPair.second < minSharedHits) { - minSharedHits = trackPair.second; - bestTrackPair = trackPair.first; - } - } - - return bestTrackPair; -} -*/ - -#endif //FASERACTSKALMANFILTER_AMBIGUITYSOLVER_H diff --git a/Tracking/Acts/FaserActsKalmanFilter/python/CombinatorialKalmanFilterConfig.py b/Tracking/Acts/FaserActsKalmanFilter/python/CombinatorialKalmanFilterConfig.py index 7f163f34e8c3fb14af4ddbea595539e295fdcf8a..317e432194b9c200d4e68d7986ae12ab858f130e 100644 --- a/Tracking/Acts/FaserActsKalmanFilter/python/CombinatorialKalmanFilterConfig.py +++ b/Tracking/Acts/FaserActsKalmanFilter/python/CombinatorialKalmanFilterConfig.py @@ -73,8 +73,6 @@ def CombinatorialKalmanFilterCfg(flags, **kwargs): ckf.RootTrajectoryStatesWriterTool = trajectory_states_writer_tool ckf.RootTrajectorySummaryWriterTool = trajectory_summary_writer_tool ckf.PerformanceWriterTool = performance_writer_tool - ckf.nTrajectories = 2 - ckf.PerformanceWriter = False ckf.nMax = 10 ckf.chi2Max = 15 diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx index fba67f83598daab4021aa38b1f16d80ce6f27969..06d203b4df754b49cf53fefe9f0935e3001300bf 100644 --- a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx +++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx @@ -19,7 +19,6 @@ #include "Acts/Surfaces/PerigeeSurface.hpp" #include "Acts/MagneticField/MagneticFieldContext.hpp" #include "FaserActsKalmanFilter/TrackSelection.h" -#include "FaserActsKalmanFilter/MyAmbiguitySolver.h" #include <algorithm> #include "FaserActsGeometry/FASERMagneticFieldWrapper.h" diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CombinatorialKalmanFilterAlg.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CombinatorialKalmanFilterAlg.cxx index be770952a50bb875bb5b4e7224f68fafab7fc92c..bea4c44f2b9ab5784bee1951fe8ebc8b76243385 100644 --- a/Tracking/Acts/FaserActsKalmanFilter/src/CombinatorialKalmanFilterAlg.cxx +++ b/Tracking/Acts/FaserActsKalmanFilter/src/CombinatorialKalmanFilterAlg.cxx @@ -19,9 +19,10 @@ #include "Acts/Surfaces/PerigeeSurface.hpp" #include "Acts/MagneticField/MagneticFieldContext.hpp" #include "FaserActsKalmanFilter/TrackSelection.h" -#include "FaserActsKalmanFilter/MyAmbiguitySolver.h" #include <algorithm> +size_t CombinatorialKalmanFilterAlg::TrajectoryInfo::nClusters {0}; + using TrajectoriesContainer = std::vector<FaserActsRecMultiTrajectory>; std::array<Acts::BoundIndices, 2> indices = {Acts::eBoundLoc0, Acts::eBoundLoc1}; @@ -53,7 +54,6 @@ StatusCode CombinatorialKalmanFilterAlg::initialize() { if (m_actsLogging == "VERBOSE") { m_logger = Acts::getDefaultLogger("KalmanFitter", Acts::Logging::VERBOSE); } else if (m_actsLogging == "DEBUG") { - m_logger = Acts::getDefaultLogger("KalmanFitter", Acts::Logging::DEBUG); } else { m_logger = Acts::getDefaultLogger("KalmanFitter", Acts::Logging::INFO); } @@ -90,6 +90,8 @@ StatusCode CombinatorialKalmanFilterAlg::execute() { std::shared_ptr<std::vector<const Tracker::FaserSCT_Cluster*>> clusters = m_trackSeedTool->clusters(); std::shared_ptr<std::vector<std::array<std::vector<const Tracker::FaserSCT_Cluster*>, 3>>> seedClusters = m_trackSeedTool->seedClusters(); + TrajectoryInfo::nClusters = sourceLinks->size(); + m_trackSeedWriterTool->writeout(geoctx, initialParameters); TrajectoriesContainer trajectories; @@ -127,111 +129,60 @@ StatusCode CombinatorialKalmanFilterAlg::execute() { } auto results = (*m_fit)(tmp, *initialParameters, options); - for (auto& result : results) { + // results contains a MultiTrajectory for each track seed with a trajectory of each branch of the CKF. + // To simplify the ambiguity solving a list of MultiTrajectories is created, each containing only a single track. + std::list<TrajectoryInfo> allTrajectories; + for (auto &result : results) { if (not result.ok()) { continue; } - Acts::CombinatorialKalmanFilterResult<IndexSourceLink> ckfResult = result.value(); - std::cout << "next seed" << std::endl; - for (const auto& index : ckfResult.lastMeasurementIndices) { - std::cout << "?? " << index << std::endl; + CKFResult ckfResult = result.value(); + for (size_t trackTip : ckfResult.lastMeasurementIndices) { + allTrajectories.emplace_back(TrajectoryInfo(FaserActsRecMultiTrajectory( + ckfResult.fittedStates, {trackTip}, {{trackTip, ckfResult.fittedParameters.at(trackTip)}}))); } } - std::vector<TrackQuality> trackQuality; - selectTracks(results, trackQuality); - - TrackFinderResult selectedResults; - for (size_t i = 0; i < std::min(trackQuality.size(), (size_t)m_nTrajectories); ++i) { - selectedResults.push_back(Acts::Result<Acts::CombinatorialKalmanFilterResult<IndexSourceLink>>(trackQuality[i].track)); - } - /* - TrackFinderResult minTrackRequirements {}; - for (auto& result : results) { - if (not result.ok()) { - continue; - } - auto& ckfResult = result.value(); - auto traj = FaserActsRecMultiTrajectory(ckfResult.fittedStates, ckfResult.lastMeasurementIndices, ckfResult.fittedParameters); - const auto& mj = traj.multiTrajectory(); - const auto& trackTips = traj.tips(); - size_t maxMeasurements = 0; - for (const auto& trackTip : trackTips) { - auto trajState = Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip); - size_t nMeasurements = trajState.nMeasurements; - if (nMeasurements > maxMeasurements) { - maxMeasurements = nMeasurements; - } - } - if (maxMeasurements >= m_minNumberMeasurements) { - minTrackRequirements.push_back(Acts::Result<Acts::CombinatorialKalmanFilterResult<IndexSourceLink>>(ckfResult)); - } + // the list of MultiTrajectories is sorted by the number of measurements using the chi2 value as a tie-breaker + allTrajectories.sort([](const TrajectoryInfo &left, const TrajectoryInfo &right) { + if (left.nMeasurements > right.nMeasurements) return true; + if (left.nMeasurements < right.nMeasurements) return false; + if (left.chi2 < right.chi2) return true; + else return false; + }); + + // select all tracks with at least 13 heats and with 6 or less shared hits, starting from the best track + // TODO use Gaudi parameters for the number of hits and shared hits + // TODO allow shared hits only in the first station? + std::vector<FaserActsRecMultiTrajectory> selectedTrajectories {}; + while (not allTrajectories.empty()) { + TrajectoryInfo selected = allTrajectories.front(); + selectedTrajectories.push_back(selected.trajectory); + allTrajectories.remove_if([&](const TrajectoryInfo &p) { + return (p.nMeasurements <= 12) || ((p.clusterSet & selected.clusterSet).count() > 6); + }); } - TrackFinderResult selectedResults {}; - if (minTrackRequirements.size() > 2) { - std::pair<std::size_t, std::size_t> trackIndices = solveAmbiguity(results); - selectedResults.push_back( - Acts::Result<Acts::CombinatorialKalmanFilterResult<IndexSourceLink>>(results.at(trackIndices.first).value())); - selectedResults.push_back( - Acts::Result<Acts::CombinatorialKalmanFilterResult<IndexSourceLink>>(results.at(trackIndices.second).value())); - } else { - for (auto& result : minTrackRequirements) { - selectedResults.push_back( - Acts::Result<Acts::CombinatorialKalmanFilterResult<IndexSourceLink>>(result.value())); - } - } - */ - computeSharedHits(sourceLinks.get(), selectedResults); - - // Loop over the track finding results for all initial parameters - for (std::size_t iseed = 0; iseed < selectedResults.size(); ++iseed) { - // The result for this seed - auto& result = selectedResults[iseed]; - if (result.ok()) { - // Get the track finding output object - const auto& trackFindingOutput = result.value(); - - std::unique_ptr<Trk::Track> track = makeTrack(geoctx, result); - if (track) { - outputTracks->push_back(std::move(track)); - } - - if (!trackFindingOutput.fittedParameters.empty()) { - const std::unordered_map<size_t, Acts::BoundTrackParameters>& params = trackFindingOutput.fittedParameters; - for (const auto& surface_params : params) { - ATH_MSG_VERBOSE("Fitted parameters for track " << iseed << ", surface " << surface_params.first); - ATH_MSG_VERBOSE(" parameters: " << surface_params.second.parameters()); - ATH_MSG_VERBOSE(" position: " << surface_params.second.position(geoctx).transpose()); - ATH_MSG_VERBOSE(" momentum: " << surface_params.second.momentum().transpose()); -// const auto& [currentTip, tipState] = trackFindingOutput.activeTips.back(); -// ATH_MSG_VERBOSE(" #measurements: " << tipState.nMeasurements); - } - } else { - ATH_MSG_DEBUG("No fitted parameters for track " << iseed); - } - // Create a Trajectories result struct - trajectories.emplace_back(trackFindingOutput.fittedStates, - trackFindingOutput.lastMeasurementIndices, - trackFindingOutput.fittedParameters); - } else { - ATH_MSG_WARNING("Track finding failed for seed " << iseed << " with error" << result.error()); - // Track finding failed. Add an empty result so the output container has - // the same number of entries as the input. - trajectories.push_back(FaserActsRecMultiTrajectory()); + // create Trk::Tracks from the trajectories + for (const FaserActsRecMultiTrajectory &traj : selectedTrajectories) { + std::unique_ptr<Trk::Track> track = makeTrack(geoctx, traj); + if (track) { + outputTracks->push_back(std::move(track)); } } + // run the performance writer if (m_statesWriter) { - ATH_CHECK(m_trajectoryStatesWriterTool->write(geoctx, trajectories)); + ATH_CHECK(m_trajectoryStatesWriterTool->write(geoctx, selectedTrajectories)); } if (m_summaryWriter) { - ATH_CHECK(m_trajectorySummaryWriterTool->write(geoctx, trajectories)); + ATH_CHECK(m_trajectorySummaryWriterTool->write(geoctx, selectedTrajectories)); } - if (m_performanceWriter) { - ATH_CHECK(m_performanceWriterTool->write(geoctx, trajectories)); + if (m_performanceWriter) { + ATH_MSG_DEBUG("?? performance writer tool"); + ATH_CHECK(m_performanceWriterTool->write(geoctx, selectedTrajectories)); } -// ATH_CHECK(trackContainer.record(std::move(outputTracks))); + ATH_CHECK(trackContainer.record(std::move(outputTracks))); return StatusCode::SUCCESS; } @@ -254,83 +205,80 @@ Acts::MagneticFieldContext CombinatorialKalmanFilterAlg::getMagneticFieldContext } std::unique_ptr<Trk::Track> -CombinatorialKalmanFilterAlg::makeTrack(Acts::GeometryContext& geoCtx, TrackFitterResult& fitResult) const { +CombinatorialKalmanFilterAlg::makeTrack(const Acts::GeometryContext &geoCtx, const FaserActsRecMultiTrajectory &traj) const { using ConstTrackStateProxy = - Acts::detail_lt::TrackStateProxy<IndexSourceLink, 6, true>; + Acts::detail_lt::TrackStateProxy<IndexSourceLink, 6, true>; std::unique_ptr<Trk::Track> newtrack = nullptr; //Get the fit output object - const auto& fitOutput = fitResult.value(); - if (fitOutput.fittedParameters.size() > 0) { - DataVector<const Trk::TrackStateOnSurface>* finalTrajectory = new DataVector<const Trk::TrackStateOnSurface>{}; - std::vector<std::unique_ptr<const Acts::BoundTrackParameters>> actsSmoothedParam; - // Loop over all the output state to create track state - fitOutput.fittedStates.visitBackwards(fitOutput.lastMeasurementIndices.front(), [&](const ConstTrackStateProxy& state) { - auto flag = state.typeFlags(); - if (state.referenceSurface().associatedDetectorElement() != nullptr) { - // We need to determine the type of state - std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> typePattern; - const Trk::TrackParameters *parm; - - // State is a hole (no associated measurement), use predicted para meters - if (flag[Acts::TrackStateFlag::HoleFlag] == true) { - const Acts::BoundTrackParameters actsParam(state.referenceSurface().getSharedPtr(), - state.predicted(), - state.predictedCovariance()); - parm = ConvertActsTrackParameterToATLAS(actsParam, geoCtx); - // auto boundaryCheck = m_boundaryCheckTool->boundaryCheck(*p arm); - typePattern.set(Trk::TrackStateOnSurface::Hole); - } - // The state was tagged as an outlier, use filtered parameters - else if (flag[Acts::TrackStateFlag::OutlierFlag] == true) { - const Acts::BoundTrackParameters actsParam(state.referenceSurface().getSharedPtr(), - state.filtered(), state.filteredCovariance()); - parm = ConvertActsTrackParameterToATLAS(actsParam, geoCtx); - typePattern.set(Trk::TrackStateOnSurface::Outlier); - } - // The state is a measurement state, use smoothed parameters - else { - const Acts::BoundTrackParameters actsParam(state.referenceSurface().getSharedPtr(), - state.smoothed(), state.smoothedCovariance()); - actsSmoothedParam.push_back(std::make_unique<const Acts::BoundTrackParameters>(Acts::BoundTrackParameters(actsParam))); - // const auto& psurface=actsParam.referenceSurface(); - Acts::Vector2 local(actsParam.parameters()[Acts::eBoundLoc0], actsParam.parameters()[Acts::eBoundLoc1]); - // const Acts::Vector3 dir = Acts::makeDirectionUnitFromPhiTheta(actsParam.parameters()[Acts::eBoundPhi], actsParam.parameters()[Acts::eBoundTheta]); - // auto pos=actsParam.position(tgContext); - parm = ConvertActsTrackParameterToATLAS(actsParam, geoCtx); - typePattern.set(Trk::TrackStateOnSurface::Measurement); - } - Tracker::FaserSCT_ClusterOnTrack* measState = nullptr; - if (state.hasUncalibrated()) { - const Tracker::FaserSCT_Cluster* fitCluster = state.uncalibrated().hit(); - if (fitCluster->detectorElement() != nullptr) { - measState = new Tracker::FaserSCT_ClusterOnTrack{ - fitCluster, - Trk::LocalParameters{ - Trk::DefinedParameter{fitCluster->localPosition()[0], Trk::loc1}, - Trk::DefinedParameter{fitCluster->localPosition()[1], Trk::loc2} - }, - fitCluster->localCovariance(), - m_idHelper->wafer_hash(fitCluster->detectorElement()->identify()) - }; - } - } - double nDoF = state.calibratedSize(); - const Trk::FitQualityOnSurface *quality = new Trk::FitQualityOnSurface(state.chi2(), nDoF); - const Trk::TrackStateOnSurface *perState = new Trk::TrackStateOnSurface(measState, parm, quality, nullptr, typePattern); - // If a state was succesfully created add it to the trajectory - if (perState) { - finalTrajectory->insert(finalTrajectory->begin(), perState); + DataVector<const Trk::TrackStateOnSurface>* finalTrajectory = new DataVector<const Trk::TrackStateOnSurface>{}; + std::vector<std::unique_ptr<const Acts::BoundTrackParameters>> actsSmoothedParam; + // Loop over all the output state to create track state + traj.multiTrajectory().visitBackwards(traj.tips().front(), [&](const ConstTrackStateProxy& state) { + auto flag = state.typeFlags(); + if (state.referenceSurface().associatedDetectorElement() != nullptr) { + // We need to determine the type of state + std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> typePattern; + const Trk::TrackParameters *parm; + + // State is a hole (no associated measurement), use predicted para meters + if (flag[Acts::TrackStateFlag::HoleFlag] == true) { + const Acts::BoundTrackParameters actsParam(state.referenceSurface().getSharedPtr(), + state.predicted(), + state.predictedCovariance()); + parm = ConvertActsTrackParameterToATLAS(actsParam, geoCtx); + // auto boundaryCheck = m_boundaryCheckTool->boundaryCheck(*p arm); + typePattern.set(Trk::TrackStateOnSurface::Hole); + } + // The state was tagged as an outlier, use filtered parameters + else if (flag[Acts::TrackStateFlag::OutlierFlag] == true) { + const Acts::BoundTrackParameters actsParam(state.referenceSurface().getSharedPtr(), + state.filtered(), state.filteredCovariance()); + parm = ConvertActsTrackParameterToATLAS(actsParam, geoCtx); + typePattern.set(Trk::TrackStateOnSurface::Outlier); + } + // The state is a measurement state, use smoothed parameters + else { + const Acts::BoundTrackParameters actsParam(state.referenceSurface().getSharedPtr(), + state.smoothed(), state.smoothedCovariance()); + actsSmoothedParam.push_back(std::make_unique<const Acts::BoundTrackParameters>(Acts::BoundTrackParameters(actsParam))); + // const auto& psurface=actsParam.referenceSurface(); + Acts::Vector2 local(actsParam.parameters()[Acts::eBoundLoc0], actsParam.parameters()[Acts::eBoundLoc1]); + // const Acts::Vector3 dir = Acts::makeDirectionUnitFromPhiTheta(actsParam.parameters()[Acts::eBoundPhi], actsParam.parameters()[Acts::eBoundTheta]); + // auto pos=actsParam.position(tgContext); + parm = ConvertActsTrackParameterToATLAS(actsParam, geoCtx); + typePattern.set(Trk::TrackStateOnSurface::Measurement); + } + Tracker::FaserSCT_ClusterOnTrack* measState = nullptr; + if (state.hasUncalibrated()) { + const Tracker::FaserSCT_Cluster* fitCluster = state.uncalibrated().hit(); + if (fitCluster->detectorElement() != nullptr) { + measState = new Tracker::FaserSCT_ClusterOnTrack{ + fitCluster, + Trk::LocalParameters{ + Trk::DefinedParameter{fitCluster->localPosition()[0], Trk::loc1}, + Trk::DefinedParameter{fitCluster->localPosition()[1], Trk::loc2} + }, + fitCluster->localCovariance(), + m_idHelper->wafer_hash(fitCluster->detectorElement()->identify()) + }; } } - return; - }); - - // Create the track using the states - const Trk::TrackInfo newInfo(Trk::TrackInfo::TrackFitter::KalmanFitter, Trk::ParticleHypothesis::muon); - // Trk::FitQuality* q = nullptr; - // newInfo.setTrackFitter(Trk::TrackInfo::TrackFitter::KalmanFitter ); //Mark the fitter as KalmanFitter - newtrack = std::make_unique<Trk::Track>(newInfo, std::move(*finalTrajectory), nullptr); - } + double nDoF = state.calibratedSize(); + const Trk::FitQualityOnSurface *quality = new Trk::FitQualityOnSurface(state.chi2(), nDoF); + const Trk::TrackStateOnSurface *perState = new Trk::TrackStateOnSurface(measState, parm, quality, nullptr, typePattern); + // If a state was succesfully created add it to the trajectory + if (perState) { + finalTrajectory->insert(finalTrajectory->begin(), perState); + } + } + return; + }); + + // Create the track using the states + const Trk::TrackInfo newInfo(Trk::TrackInfo::TrackFitter::KalmanFitter, Trk::ParticleHypothesis::muon); + // Trk::FitQuality* q = nullptr; + // newInfo.setTrackFitter(Trk::TrackInfo::TrackFitter::KalmanFitter ); //Mark the fitter as KalmanFitter + newtrack = std::make_unique<Trk::Track>(newInfo, std::move(*finalTrajectory), nullptr); return newtrack; } diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx index 1c5bde3eae4439da2d3201ebdf9163017f5aa165..d0b40a58a723009b2c9d602e26db19631e79ac8b 100644 --- a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx +++ b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx @@ -205,10 +205,17 @@ KalmanFitterTool::getParametersFromTrack(const Trk::TrackParameters *inputParame auto atlasParam = inputParameters->parameters(); auto center = inputParameters->associatedSurface().center(); - Acts::BoundVector params {center[Trk::locY], center[Trk::locX], - atlasParam[Trk::phi0], atlasParam[Trk::theta], - charge / (inputParameters->momentum().norm() * 1_MeV), - 0.}; + // Acts::BoundVector params {center[Trk::locY], center[Trk::locX], + // atlasParam[Trk::phi0], atlasParam[Trk::theta], + // charge / (inputParameters->momentum().norm() * 1_MeV), + // 0.}; + Acts::BoundVector params; + params(0.0) = center[Trk::locY]; + params(1.0) = center[Trk::locX]; + params(2.0) = atlasParam[Trk::phi0]; + params(3.0) = atlasParam[Trk::theta]; + params(4.0) = charge / (inputParameters->momentum().norm() * 1_MeV); + params(5.0) = 0.; Acts::BoundSymMatrix cov = Acts::BoundSymMatrix::Identity(); cov.topLeftCorner(5, 5) = *inputParameters->covariance(); diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/MyAmbiguitySolver.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/MyAmbiguitySolver.cxx deleted file mode 100644 index 998269a83cf6228092dd3b31987696f929820e4e..0000000000000000000000000000000000000000 --- a/Tracking/Acts/FaserActsKalmanFilter/src/MyAmbiguitySolver.cxx +++ /dev/null @@ -1,3 +0,0 @@ -#include "FaserActsKalmanFilter/MyAmbiguitySolver.h" - - diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/ThreeStationTrackSeedTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/ThreeStationTrackSeedTool.cxx index 0203835ffc4fb4c717990e31cc621ae4cdcf3eea..8c6eb8c0328b0ca194cb06083cc3c019cdbc0fd7 100644 --- a/Tracking/Acts/FaserActsKalmanFilter/src/ThreeStationTrackSeedTool.cxx +++ b/Tracking/Acts/FaserActsKalmanFilter/src/ThreeStationTrackSeedTool.cxx @@ -126,9 +126,7 @@ Acts::CurvilinearTrackParameters ThreeStationTrackSeedTool::get_params( Acts::Vector3 dir = position_st2 - position_st1; Acts::Vector3 pos = position_st1 - (position_st1.z() - origin)/dir.z() * dir; Acts::Vector4 pos4 {pos.x(), pos.y(), pos.z(), 0}; - std::cout << "calc momentum\n"; auto [abs_momentum, charge] = momentum({{1, position_st1}, {2, position_st2}, {3, position_st3}}); - std::cout << abs_momentum << std::endl; return Acts::CurvilinearTrackParameters(pos4, dir, abs_momentum, charge, cov); } std::pair<double, double> ThreeStationTrackSeedTool::momentum(const std::map<int, Amg::Vector3D>& pos, double B) { diff --git a/Tracking/Acts/FaserActsKalmanFilter/test/CombinatorialKalmanFilterAlg.py b/Tracking/Acts/FaserActsKalmanFilter/test/CombinatorialKalmanFilterAlg.py index c8f99e85b9c263b9aadc8c884a0ddf322ab6da86..1724e04ede08bc93e4c775810f3984818c82ada0 100644 --- a/Tracking/Acts/FaserActsKalmanFilter/test/CombinatorialKalmanFilterAlg.py +++ b/Tracking/Acts/FaserActsKalmanFilter/test/CombinatorialKalmanFilterAlg.py @@ -15,7 +15,7 @@ from FaserActsKalmanFilter.CombinatorialKalmanFilterConfig import CombinatorialK log.setLevel(DEBUG) Configurable.configurableRun3Behavior = True -ConfigFlags.Input.Files = ['../single_muon_1000GeV/my.RDO.pool.root'] +ConfigFlags.Input.Files = ['my.RDO.pool.root'] ConfigFlags.Output.ESDFileName = "CKF.ESD.pool.root" ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01" ConfigFlags.GeoModel.FaserVersion = "FASER-01" @@ -29,7 +29,7 @@ acc.merge(PoolReadCfg(ConfigFlags)) # acc.merge(PoolWriteCfg(ConfigFlags)) acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags)) -acc.merge(SegmentFitAlgCfg(ConfigFlags)) +acc.merge(SegmentFitAlgCfg(ConfigFlags, SharedHitFraction=0.5, MinClustersPerFit=5)) acc.merge(CombinatorialKalmanFilterCfg(ConfigFlags)) acc.getEventAlgo("CombinatorialKalmanFilterAlg").OutputLevel = VERBOSE