diff --git a/Detectors/ITSMFT/MFT/workflow/CMakeLists.txt b/Detectors/ITSMFT/MFT/workflow/CMakeLists.txt index 5c78df4dd9651..cee277a16ad7f 100644 --- a/Detectors/ITSMFT/MFT/workflow/CMakeLists.txt +++ b/Detectors/ITSMFT/MFT/workflow/CMakeLists.txt @@ -22,11 +22,13 @@ o2_add_library(MFTWorkflow O2::SimConfig O2::SimulationDataFormat O2::ITSMFTReconstruction + O2::ITSMFTTracking O2::MFTTracking O2::MFTAssessment O2::DataFormatsMFT O2::ITSMFTWorkflow O2::ITStracking + TBB::tbb O2::MFTAlignment O2::GlobalTrackingWorkflowReaders) o2_add_executable(reco-workflow diff --git a/Detectors/ITSMFT/MFT/workflow/include/MFTWorkflow/CATrackerSpec.h b/Detectors/ITSMFT/MFT/workflow/include/MFTWorkflow/CATrackerSpec.h index ace31aa83af9b..c1414d2a79f8b 100644 --- a/Detectors/ITSMFT/MFT/workflow/include/MFTWorkflow/CATrackerSpec.h +++ b/Detectors/ITSMFT/MFT/workflow/include/MFTWorkflow/CATrackerSpec.h @@ -15,28 +15,26 @@ #define O2_MFT_CATRACKERSPEC_H_ #include -#include -#include - -#include "DataFormatsITSMFT/ROFRecord.h" -#include "DataFormatsITSMFT/TopologyDictionary.h" -#include "DataFormatsParameters/GRPObject.h" #include "DetectorsBase/GRPGeomHelper.h" #include "Framework/DataProcessorSpec.h" #include "Framework/Task.h" -#include "ITStracking/ROFLookupTables.h" +#include "ITSMFTTracking/Configuration.h" +#include "ITSMFTTracking/TrackingInterface.h" namespace o2::mft { +/// MFT CA tracker DPL task. Delegates reconstruction to ITSMFTTrackingInterfaceMFT. class CATrackerDPL : public o2::framework::Task { public: - static constexpr int NCALayers = 5; - using ROFOverlapTable = o2::its::ROFOverlapTable; - - CATrackerDPL(std::shared_ptr gr, bool useMC) : mGGCCDBRequest(std::move(gr)), mUseMC(useMC) {} + CATrackerDPL(std::shared_ptr gr, + bool useMC, + o2::itsmft::TrackingMode::Type trMode) + : mGGCCDBRequest(std::move(gr)), mUseMC(useMC), mTracking(useMC, trMode, false) + { + } ~CATrackerDPL() override = default; void init(framework::InitContext& ic) final; @@ -45,18 +43,15 @@ class CATrackerDPL : public o2::framework::Task private: void updateTimeDependentParams(framework::ProcessingContext& pc); - void initialiseROFTable(gsl::span rofs); - bool mUseMC = false; std::shared_ptr mGGCCDBRequest; - const o2::itsmft::TopologyDictionary* mDict = nullptr; - ROFOverlapTable mROFTable; - ROFOverlapTable::View mROFTableView; + bool mUseMC = false; + bool mTrackingInitialised = false; + o2::itsmft::tracking::ITSMFTTrackingInterfaceMFT mTracking; }; -/// create a processor spec that reads all MFT tracker inputs and provides the -/// future insertion point for wiring them into a generalized ITS TimeFrame. -o2::framework::DataProcessorSpec getCATrackerSpec(bool useMC, bool useGeom, bool useIRFrames); +o2::framework::DataProcessorSpec getCATrackerSpec(bool useMC, bool useGeom, bool useIRFrames, + o2::itsmft::TrackingMode::Type trMode); } // namespace o2::mft diff --git a/Detectors/ITSMFT/MFT/workflow/include/MFTWorkflow/TrackWriterSpec.h b/Detectors/ITSMFT/MFT/workflow/include/MFTWorkflow/TrackWriterSpec.h index 5a8d50939a25a..a3cab9723315e 100644 --- a/Detectors/ITSMFT/MFT/workflow/include/MFTWorkflow/TrackWriterSpec.h +++ b/Detectors/ITSMFT/MFT/workflow/include/MFTWorkflow/TrackWriterSpec.h @@ -26,7 +26,7 @@ namespace mft /// create a processor spec /// write MFT tracks a root file -o2::framework::DataProcessorSpec getTrackWriterSpec(bool useMC); +o2::framework::DataProcessorSpec getTrackWriterSpec(bool useMC, bool useCA = false); } // namespace mft } // namespace o2 diff --git a/Detectors/ITSMFT/MFT/workflow/src/CATrackerSpec.cxx b/Detectors/ITSMFT/MFT/workflow/src/CATrackerSpec.cxx index 78ba9b9fc7a19..7587f1bff9f93 100644 --- a/Detectors/ITSMFT/MFT/workflow/src/CATrackerSpec.cxx +++ b/Detectors/ITSMFT/MFT/workflow/src/CATrackerSpec.cxx @@ -13,26 +13,24 @@ #include "MFTWorkflow/CATrackerSpec.h" +#include #include #include -#include "CommonConstants/LHCConstants.h" #include "CommonDataFormat/IRFrame.h" #include "DataFormatsITSMFT/CompCluster.h" -#include "DataFormatsITSMFT/DPLAlpideParam.h" #include "DataFormatsITSMFT/ROFRecord.h" #include "DataFormatsITSMFT/TopologyDictionary.h" -#include "DataFormatsParameters/GRPObject.h" +#include "DataFormatsMFT/TrackMFT.h" #include "DetectorsBase/GeometryManager.h" -#include "DetectorsCommonDataFormats/DetectorNameConf.h" #include "Framework/CCDBParamSpec.h" -#include "Framework/ConfigParamRegistry.h" -#include "Framework/ControlService.h" #include "Framework/DataProcessorSpec.h" -#include "Framework/DeviceSpec.h" #include "Framework/Logger.h" +#include "ITSMFTTracking/MFTCATrack.h" #include "MFTBase/GeometryTGeo.h" +#include "MFTTracking/Constants.h" +#include "MFTTracking/MFTTrackingParam.h" #include "SimulationDataFormat/MCCompLabel.h" #include "SimulationDataFormat/MCTruthContainer.h" @@ -41,6 +39,76 @@ using namespace o2::framework; namespace o2::mft { +static_assert(o2::itsmft::tracking::ITSMFTTrackingInterfaceMFT::DetId == o2::detectors::DetID::MFT); + +namespace +{ +template +void fillMFTOutputs(const o2::itsmft::tracking::TimeFrameMFT& tf, + gsl::span inputROFs, + TracksVec& tracks, + ClusterIdxVec& clusterIndices, + ROFVec& trackROFs, + LabelsVec& trackLabels, + SeedPatternVec& seedPatterns, + bool useMC) +{ + trackROFs.assign(inputROFs.begin(), inputROFs.end()); + for (auto& rof : trackROFs) { + rof.setFirstEntry(0); + rof.setNEntries(0); + } + + const auto& tracksIn = tf.getTracks(); + tracks.reserve(tracksIn.size()); + seedPatterns.reserve(tracksIn.size()); + if (useMC) { + trackLabels.reserve(tracksIn.size()); + } + + const auto& clockLayer = tf.getROFOverlapTableView().getClockLayer(); + std::vector rofEntries(trackROFs.size() + 1, 0); + + for (size_t iTrk = 0; iTrk < tracksIn.size(); ++iTrk) { + const auto& src = tracksIn[iTrk]; + auto dst = src.getTrack(); + dst.setExternalClusterIndexOffset(clusterIndices.size()); + int nPoints = 0; + for (int layer = o2::itsmft::tracking::MFTCATrack::MaxClusters; layer--;) { + if (!src.hasHitOnLayer(layer)) { + continue; + } + const int extIdx = src.getClusterIndex(layer); + if (extIdx < 0) { + continue; + } + const int clsSize = src.getClusterSize(layer); + dst.setClusterSize(layer, clsSize > 0 ? clsSize : tf.getClusterSize(0, extIdx)); + clusterIndices.push_back(extIdx); + ++nPoints; + } + dst.setNumberOfPoints(nPoints); + tracks.push_back(dst); + seedPatterns.push_back(src.getSeedPattern()); + + if (useMC && iTrk < tf.getTracksLabel().size()) { + trackLabels.push_back(tf.getTracksLabel()[iTrk]); + } + + const auto rof = clockLayer.getROF(src.getTimeStamp()); + if (rof >= 0 && rof < static_cast(trackROFs.size())) { + ++rofEntries[rof]; + } + } + + std::exclusive_scan(rofEntries.begin(), rofEntries.end(), rofEntries.begin(), 0); + for (size_t iROF = 0; iROF < trackROFs.size(); ++iROF) { + trackROFs[iROF].setFirstEntry(rofEntries[iROF]); + trackROFs[iROF].setNEntries(rofEntries[iROF + 1] - rofEntries[iROF]); + } +} +} // namespace + void CATrackerDPL::init(InitContext&) { o2::base::GRPGeomHelper::instance().setRequest(mGGCCDBRequest); @@ -50,37 +118,66 @@ void CATrackerDPL::run(ProcessingContext& pc) { updateTimeDependentParams(pc); + auto rofsinput = pc.inputs().get>("ROframes"); + auto& trackROFs = pc.outputs().make>(Output{"MFT", "MFTTrackROF", 0}, + rofsinput.begin(), rofsinput.end()); + auto& allTracksMFT = pc.outputs().make>(Output{"MFT", "TRACKS", 0}); + auto& allClusIdx = pc.outputs().make>(Output{"MFT", "TRACKCLSID", 0}); + auto& allSeedPatterns = pc.outputs().make>(Output{"MFT", "TRACKSEEDPAT", 0}); + std::vector allTrackLabels; + + if (!mTracking.isActive()) { + return; + } + auto compClusters = pc.inputs().get>("compClusters"); gsl::span patterns = pc.inputs().get>("patterns"); - auto rofs = pc.inputs().get>("ROframes"); - initialiseROFTable(gsl::span(rofs.data(), rofs.size())); const dataformats::MCTruthContainer* labels = nullptr; - if (mUseMC) { + if (mUseMC && pc.inputs().getPos("labels") >= 0) { labels = pc.inputs().get*>("labels").release(); } - LOGP(info, "MFT CA input pulled {} compressed clusters, {} pattern bytes, {} RO frames, dictionary={}, MC labels={}", - compClusters.size(), patterns.size(), rofs.size(), mDict != nullptr, labels != nullptr); - - size_t nClusterRefs = 0; - for (const auto& rof : rofs) { - nClusterRefs += rof.getNEntries(); + gsl::span irFrames; + if (pc.inputs().getPos("IRFramesITS") >= 0) { + irFrames = pc.inputs().get>("IRFramesITS"); } - LOGP(info, "MFT CA input ROF cluster references: {}", nClusterRefs); - if (pc.inputs().getPos("IRFramesITS") >= 0) { - auto irFrames = pc.inputs().get>("IRFramesITS"); - LOGP(info, "MFT CA input pulled {} ITS IR frames", irFrames.size()); + LOGP(info, "MFT CA input pulled {} compressed clusters in {} RO frames ({} pattern bytes)", + compClusters.size(), rofsinput.size(), patterns.size()); + + mTracking.processTimeFrame(gsl::span(rofsinput.data(), rofsinput.size()), + gsl::span(compClusters.data(), compClusters.size()), + patterns, + labels, + irFrames); + + fillMFTOutputs(mTracking.getTimeFrame(), + gsl::span(rofsinput.data(), rofsinput.size()), + allTracksMFT, + allClusIdx, + trackROFs, + allTrackLabels, + allSeedPatterns, + mUseMC); + + LOGP(info, "MFT CA pushed {} tracks in {} ROFs", allTracksMFT.size(), trackROFs.size()); + + if (mUseMC) { + pc.outputs().snapshot(Output{"MFT", "TRACKSMCTR", 0}, allTrackLabels); + LOGP(info, "MFT CA pushed {} track MC labels", allTrackLabels.size()); } - // Next step: replace this inspection point with a loader into a generalized - // ITS TimeFrame interface that accepts MFT layer geometry and cluster data. + mTracking.clearTimeFrame(); } void CATrackerDPL::updateTimeDependentParams(ProcessingContext& pc) { o2::base::GRPGeomHelper::instance().checkUpdates(pc); + if (!mTrackingInitialised) { + mTrackingInitialised = true; + mTracking.initialise(); + } static bool initOnceDone = false; if (!initOnceDone) { initOnceDone = true; @@ -90,43 +187,11 @@ void CATrackerDPL::updateTimeDependentParams(ProcessingContext& pc) pc.inputs().get("cldict"); o2::mft::GeometryTGeo::Instance()->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::T2GRot, - o2::math_utils::TransformType::T2G)); + o2::math_utils::TransformType::T2G, + o2::math_utils::TransformType::L2G)); } } -void CATrackerDPL::initialiseROFTable(gsl::span rofs) -{ - if (!o2::base::GRPGeomHelper::instance().getGRPECS()->isDetContinuousReadOut(o2::detectors::DetID::MFT)) { - LOGP(fatal, "MFT CA tracker currently supports only continuous readout"); - } - - const auto& par = o2::itsmft::DPLAlpideParam::Instance(); - const int nOrbitsPerTF = o2::base::GRPGeomHelper::getNHBFPerTF(); - const int timingSourceLayer = 0; // MFT CA disks share timing for now; keep per-CA-layer definition for future staggering. - const unsigned int nROFsPerOrbit = o2::constants::lhc::LHCMaxBunches / par.getROFLengthInBC(timingSourceLayer); - const auto nROFsTF = nROFsPerOrbit * nOrbitsPerTF; - - ROFOverlapTable rofTable; - for (int caLayer = 0; caLayer < NCALayers; ++caLayer) { - const o2::its::LayerTiming timing{ - .mNROFsTF = nROFsTF, - .mROFLength = static_cast(par.getROFLengthInBC(timingSourceLayer)), - .mROFDelay = static_cast(par.getROFDelayInBC(timingSourceLayer)), - .mROFBias = static_cast(par.getROFBiasInBC(timingSourceLayer)), - .mROFAddTimeErr = 0}; - rofTable.defineLayer(caLayer, timing); - } - rofTable.init(); - mROFTable = std::move(rofTable); - mROFTableView = mROFTable.getView(); - - if (rofs.size() != nROFsTF) { - LOGP(warn, "MFT CA ROF count differs from continuous timing expectation: received {} expected {}", rofs.size(), nROFsTF); - } - LOGP(info, "MFT CA ROF lookup table initialised for {} CA layers, {} ROFs/TF, ROF length {} BC, bias {} BC", - NCALayers, nROFsTF, par.getROFLengthInBC(timingSourceLayer), par.getROFBiasInBC(timingSourceLayer)); -} - void CATrackerDPL::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) { if (o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) { @@ -134,17 +199,21 @@ void CATrackerDPL::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) } if (matcher == ConcreteDataMatcher("MFT", "CLUSDICT", 0)) { LOG(info) << "MFT CA input cluster dictionary updated"; - mDict = static_cast(obj); + mTracking.setClusterDictionary(static_cast(obj)); return; } if (matcher == ConcreteDataMatcher("MFT", "GEOMTGEO", 0)) { LOG(info) << "MFT CA input GeometryTGeo loaded from CCDB"; o2::mft::GeometryTGeo::adopt(static_cast(obj)); + o2::mft::GeometryTGeo::Instance()->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, + o2::math_utils::TransformType::T2GRot, + o2::math_utils::TransformType::T2G, + o2::math_utils::TransformType::L2G)); return; } } -DataProcessorSpec getCATrackerSpec(bool useMC, bool useGeom, bool useIRFrames) +DataProcessorSpec getCATrackerSpec(bool useMC, bool useGeom, bool useIRFrames, o2::itsmft::TrackingMode::Type trMode) { std::vector inputs; inputs.emplace_back("compClusters", "MFT", "COMPCLUSTERS", 0, Lifetime::Timeframe); @@ -155,27 +224,38 @@ DataProcessorSpec getCATrackerSpec(bool useMC, bool useGeom, bool useIRFrames) if (useMC) { inputs.emplace_back("labels", "MFT", "CLUSTERSMCTR", 0, Lifetime::Timeframe); } - if (useIRFrames) { + + const auto& trackingParam = MFTTrackingParam::Instance(); + if (useIRFrames || trackingParam.irFramesOnly) { inputs.emplace_back("IRFramesITS", "ITS", "IRFRAMES", 0, Lifetime::Timeframe); } - auto ggRequest = std::make_shared(false, // orbitResetTime - true, // GRPECS=true - false, // GRPLHCIF - true, // GRPMagField - false, // askMatLUT - useGeom ? o2::base::GRPGeomRequest::Aligned : o2::base::GRPGeomRequest::None, // geometry + auto ggRequest = std::make_shared(false, + true, + false, + true, + true, + useGeom ? o2::base::GRPGeomRequest::Aligned : o2::base::GRPGeomRequest::None, inputs, true); if (!useGeom) { ggRequest->addInput({"mftTGeo", "MFT", "GEOMTGEO", 0, Lifetime::Condition, framework::ccdbParamSpec("MFT/Config/Geometry")}, inputs); } + std::vector outputs; + outputs.emplace_back("MFT", "TRACKS", 0, Lifetime::Timeframe); + outputs.emplace_back("MFT", "MFTTrackROF", 0, Lifetime::Timeframe); + outputs.emplace_back("MFT", "TRACKCLSID", 0, Lifetime::Timeframe); + outputs.emplace_back("MFT", "TRACKSEEDPAT", 0, Lifetime::Timeframe); + if (useMC) { + outputs.emplace_back("MFT", "TRACKSMCTR", 0, Lifetime::Timeframe); + } + return DataProcessorSpec{ "mft-ca-tracker", inputs, - {}, - AlgorithmSpec{adaptFromTask(ggRequest, useMC)}, + outputs, + AlgorithmSpec{adaptFromTask(ggRequest, useMC, trMode)}, Options{}}; } diff --git a/Detectors/ITSMFT/MFT/workflow/src/TrackWriterSpec.cxx b/Detectors/ITSMFT/MFT/workflow/src/TrackWriterSpec.cxx index f8a848f6fde32..d12fbd125b5b2 100644 --- a/Detectors/ITSMFT/MFT/workflow/src/TrackWriterSpec.cxx +++ b/Detectors/ITSMFT/MFT/workflow/src/TrackWriterSpec.cxx @@ -34,7 +34,7 @@ template using BranchDefinition = MakeRootTreeWriterSpec::BranchDefinition; using namespace o2::header; -DataProcessorSpec getTrackWriterSpec(bool useMC) +DataProcessorSpec getTrackWriterSpec(bool useMC, bool useCA) { // Spectators for logging // this is only to restore the original behavior @@ -53,6 +53,10 @@ DataProcessorSpec getTrackWriterSpec(bool useMC) tracksSizeGetter}, BranchDefinition>{InputSpec{"trackClIdx", "MFT", "TRACKCLSID", 0}, "MFTTrackClusIdx"}, + BranchDefinition>{InputSpec{"trackSeedPat", "MFT", "TRACKSEEDPAT", 0}, + "MFTTrackSeedPattern", + (useCA ? 1 : 0), + ""}, BranchDefinition>{InputSpec{"ROframes", "MFT", "MFTTrackROF", 0}, "MFTTracksROF", logger}, diff --git a/Detectors/ITSMFT/MFT/workflow/src/mft-ca-tracker-workflow.cxx b/Detectors/ITSMFT/MFT/workflow/src/mft-ca-tracker-workflow.cxx index c8d236ce8983e..8439580f0d59f 100644 --- a/Detectors/ITSMFT/MFT/workflow/src/mft-ca-tracker-workflow.cxx +++ b/Detectors/ITSMFT/MFT/workflow/src/mft-ca-tracker-workflow.cxx @@ -12,28 +12,63 @@ /// @file mft-ca-tracker-workflow.cxx #include +#include +#include "CommonUtils/ConfigurableParam.h" +#include "DataFormatsITSMFT/DPLAlpideParamInitializer.h" +#include "DetectorsRaw/HBFUtilsInitializer.h" +#include "Framework/CallbacksPolicy.h" +#include "Framework/CompletionPolicyHelpers.h" #include "Framework/ConfigParamSpec.h" +#include "ITSMFTTracking/Configuration.h" +#include "ITSMFTTracking/TrackingConfigParam.h" #include "MFTWorkflow/CATrackerSpec.h" +#include "MFTWorkflow/TrackWriterSpec.h" using namespace o2::framework; +void customize(std::vector& policies) +{ + o2::raw::HBFUtilsInitializer::addNewTimeSliceCallback(policies); +} + +void customize(std::vector& policies) +{ + policies.push_back(CompletionPolicyHelpers::consumeWhenAllOrdered(".*(?:MFT|mft).*[W,w]riter.*")); +} + void customize(std::vector& workflowOptions) { workflowOptions.push_back(ConfigParamSpec{"disable-mc", VariantType::Bool, false, {"disable MC labels"}}); + workflowOptions.push_back(ConfigParamSpec{"disable-root-output", VariantType::Bool, false, {"do not write output root files"}}); workflowOptions.push_back(ConfigParamSpec{"use-geom", VariantType::Bool, false, {"use geometry from the global geometry manager"}}); workflowOptions.push_back(ConfigParamSpec{"use-irframes", VariantType::Bool, false, {"consume ITS IR frames"}}); + workflowOptions.push_back(ConfigParamSpec{"tracking-mode", VariantType::String, "sync", {"sync,async,cosmics,unset,off"}}); + workflowOptions.push_back(ConfigParamSpec{"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings (e.g. MFTCATrackerParam.nIterations=4;MFTAlpideParam.roFrameLengthInBC=594)"}}); + o2::itsmft::DPLAlpideParamInitializer::addMFTConfigOption(workflowOptions); + o2::raw::HBFUtilsInitializer::addConfigOption(workflowOptions); } #include "Framework/runDataProcessing.h" WorkflowSpec defineDataProcessing(ConfigContext const& config) { + std::ignore = o2::itsmft::TrackerParamConfig::Instance(); + o2::conf::ConfigurableParam::updateFromString(config.options().get("configKeyValues")); + const bool useMC = !config.options().get("disable-mc"); + const bool disableRootOutput = config.options().get("disable-root-output"); const bool useGeom = config.options().get("use-geom"); const bool useIRFrames = config.options().get("use-irframes"); + const auto trMode = o2::itsmft::TrackingMode::fromString(config.options().get("tracking-mode")); WorkflowSpec specs; - specs.emplace_back(o2::mft::getCATrackerSpec(useMC, useGeom, useIRFrames)); + specs.emplace_back(o2::mft::getCATrackerSpec(useMC, useGeom, useIRFrames, trMode)); + if (!disableRootOutput) { + specs.emplace_back(o2::mft::getTrackWriterSpec(useMC, true)); + } + + o2::raw::HBFUtilsInitializer hbfIni(config, specs); + return specs; } diff --git a/Detectors/ITSMFT/common/CMakeLists.txt b/Detectors/ITSMFT/common/CMakeLists.txt index 3991f3e67a82b..57664249f41f7 100644 --- a/Detectors/ITSMFT/common/CMakeLists.txt +++ b/Detectors/ITSMFT/common/CMakeLists.txt @@ -12,5 +12,6 @@ add_subdirectory(base) add_subdirectory(simulation) add_subdirectory(reconstruction) +add_subdirectory(tracking) add_subdirectory(workflow) add_subdirectory(data) \ No newline at end of file diff --git a/Detectors/ITSMFT/common/tracking/CMakeLists.txt b/Detectors/ITSMFT/common/tracking/CMakeLists.txt new file mode 100644 index 0000000000000..0aaf05f8d8e3f --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/CMakeLists.txt @@ -0,0 +1,59 @@ +# Copyright 2019-2020 CERN and copyright holders of ALICE O2. +# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +# All rights not expressly granted are reserved. +# +# This software is distributed under the terms of the GNU General Public +# License v3 (GPL Version 3), copied verbatim in the file "COPYING". +# +# In applying this license CERN does not waive the privileges and immunities +# granted to it by virtue of its status as an Intergovernmental Organization +# or submit itself to any jurisdiction. + +# Lightweight param registration so MFTCATrackerParam keys in --configKeyValues +# are known to workflows that do not link O2::ITSMFTTracking (e.g. ctf-reader). +o2_add_library(ITSMFTTrackingParams + SOURCES src/TrackingConfigParam.cxx + PUBLIC_LINK_LIBRARIES O2::CommonUtils + O2::DetectorsCommonDataFormats) + +o2_target_root_dictionary(ITSMFTTrackingParams + HEADERS include/ITSMFTTracking/TrackingConfigParam.h + LINKDEF src/ITSMFTTrackingLinkDef.h) + +# Shared CA tracking infrastructure for ITS and MFT. Currently wired from the +# MFT CA tracker workflow via ITSMFTTrackingInterface; ITS production tracking +# remains on O2::ITSTrackingInterface / O2::ITStracking. +# Dependency is one-way (ITSMFTTracking -> ITStracking) to avoid a circular link. +o2_add_library(ITSMFTTracking + TARGETVARNAME targetName + SOURCES src/IOUtils.cxx + src/DetectorTraits.cxx + src/Configuration.cxx + src/TimeFrame.cxx + src/TrackerTraits.cxx + src/CATracker.cxx + src/TrackingInterface.cxx + src/MFTFwdTrackHelpers.cxx + PUBLIC_LINK_LIBRARIES + O2::ITSMFTTrackingParams + O2::ITStracking + O2::GPUCommon + O2::CommonConstants + O2::DetectorsCommonDataFormats + O2::DataFormatsITSMFT + O2::ITSMFTBase + O2::CommonUtils + O2::DetectorsBase + O2::MathUtils + Microsoft.GSL::GSL + O2::SimulationDataFormat + O2::ReconstructionDataFormats + O2::DataFormatsCalibration + O2::DataFormatsMFT + TBB::tbb + PRIVATE_LINK_LIBRARIES + O2::Framework + O2::ITSBase + O2::MFTBase + O2::MFTTracking + O2::DataFormatsITS) diff --git a/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/CATracker.h b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/CATracker.h new file mode 100644 index 0000000000000..95230baa70c90 --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/CATracker.h @@ -0,0 +1,77 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file CATracker.h +/// \brief Shared CA tracker orchestrator (same role as ITStracking/Tracker.h) +/// + +#ifndef ALICEO2_ITSMFT_TRACKING_CATRACKER_H_ +#define ALICEO2_ITSMFT_TRACKING_CATRACKER_H_ + +#include +#include + +#include + +#include "ITSMFTTracking/Configuration.h" +#include "ITSMFTTracking/TimeFrame.h" +#include "ITSMFTTracking/TrackerTraits.h" + +namespace o2::itsmft::tracking +{ + +template +class Tracker +{ + public: + using TimeFrameN = TimeFrame; + using TrackerTraitsN = TrackerTraits; + + explicit Tracker(TrackerTraitsN* traits); + + void adoptTimeFrame(TimeFrameN& tf); + void setParameters(const std::vector& p) { mTrkParams = p; } + void setMemoryPool(std::shared_ptr pool) { mMemoryPool = pool; } + void setBz(float bz) { mTraits->setBz(bz); } + void setNThreads(int n, std::shared_ptr& arena) { mTraits->setNThreads(n, arena); } + + /// Run all configured iterations; returns elapsed ms or -1 on failure. + float clustersToTracks(); + + const TimeFrameN& getTimeFrame() const { return *mTimeFrame; } + TimeFrameN& getTimeFrame() { return *mTimeFrame; } + + private: + void initialiseTimeFrame(int iteration) { mTraits->initialiseTimeFrame(iteration); } + void computeTracklets(int iteration, int iVertex) { mTraits->computeLayerTracklets(iteration, iVertex); } + void computeCells(int iteration) { mTraits->computeLayerCells(iteration); } + void findCellsNeighbours(int iteration) { mTraits->findCellsNeighbours(iteration); } + void findRoads(int iteration) { mTraits->findRoads(iteration); } + void rectifyClusterIndices(); + void sortTracks(); + + TrackerTraitsN* mTraits = nullptr; + TimeFrameN* mTimeFrame = nullptr; + std::vector mTrkParams; + std::shared_ptr mMemoryPool; +}; + +template +using CATracker = Tracker; + +using TrackerITS = Tracker; +using TrackerMFT = Tracker; +using CATrackerITS = TrackerITS; +using CATrackerMFT = TrackerMFT; + +} // namespace o2::itsmft::tracking + +#endif /* ALICEO2_ITSMFT_TRACKING_CATRACKER_H_ */ diff --git a/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/Cell.h b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/Cell.h new file mode 100644 index 0000000000000..403f71b6de49d --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/Cell.h @@ -0,0 +1,217 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file Cell.h +/// \brief CA cell/track seed types with hole-layer support (ITS PR #15390) +/// + +#ifndef ALICEO2_ITSMFT_TRACKING_INCLUDE_CACELL_H_ +#define ALICEO2_ITSMFT_TRACKING_INCLUDE_CACELL_H_ + +#include +#include +#include + +#include "DetectorsCommonDataFormats/DetID.h" +#include "DataFormatsITS/TrackITS.h" +#include "DataFormatsITS/TimeEstBC.h" +#include "ITSMFTTracking/Configuration.h" +#include "ITSMFTTracking/LayerMask.h" +#include "ITStracking/Constants.h" +#include "ReconstructionDataFormats/Track.h" +#include "ReconstructionDataFormats/TrackFwd.h" +#include "GPUCommonDef.h" + +namespace o2::itsmft::tracking +{ + +struct CellNeighbour { + int cellTopology{-1}; + int cell{-1}; + int nextCellTopology{-1}; + int nextCell{-1}; + int level{-1}; +}; + +template +class SeedBase : public TrackParT +{ + public: + GPUhd() LayerMask getHitLayerMask() const { return LayerMask{mHitLayerMask}; } + GPUhd() void setHitLayerMask(LayerMask mask) { mHitLayerMask = mask.value(); } + GPUhd() int getInnerLayer() const { return getHitLayerMask().first(); } + GPUhd() int getFirstTrackletIndex() const { return mTracklets[0]; }; + GPUhd() void setFirstTrackletIndex(int trkl) { mTracklets[0] = trkl; }; + GPUhd() int getSecondTrackletIndex() const { return mTracklets[1]; }; + GPUhd() void setSecondTrackletIndex(int trkl) { mTracklets[1] = trkl; }; + GPUhd() float getChi2() const { return mChi2; }; + GPUhd() void setChi2(float chi2) { mChi2 = chi2; }; + GPUhd() int getLevel() const { return mLevel; }; + GPUhd() void setLevel(int level) { mLevel = level; }; + GPUhd() int* getLevelPtr() { return &mLevel; } + GPUhd() auto& getTimeStamp() noexcept { return mTime; } + GPUhd() const auto& getTimeStamp() const noexcept { return mTime; } + + /// Road-length filter: barrel q/pT² for ITS, (invQPt)² for forward MFT seeds. + GPUhd() float getQ2Pt() const + { + if constexpr (std::is_same_v) { + const float invQPt = static_cast(TrackParT::getInvQPt()); + return invQPt * invQPt; + } else { + return TrackParT::getQ2Pt(); + } + } + + protected: + GPUhdDefault() SeedBase() = default; + GPUhdDefault() SeedBase(const SeedBase&) = default; + GPUhdDefault() ~SeedBase() = default; + GPUhdDefault() SeedBase(SeedBase&&) = default; + GPUhdDefault() SeedBase& operator=(const SeedBase&) = default; + GPUhdDefault() SeedBase& operator=(SeedBase&&) = default; + GPUhd() SeedBase(const TrackParT& tpc, float chi2, int level, const o2::its::TimeEstBC& time) + : TrackParT(tpc), mChi2(chi2), mLevel(level), mTime(time) + { + } + GPUhd() auto& clustersRaw() { return mClusters; } + GPUhd() const auto& clustersRaw() const { return mClusters; } + + private: + uint16_t mHitLayerMask{0}; + float mChi2{o2::its::constants::UnsetValue}; + int mLevel{o2::its::constants::UnusedIndex}; + std::array mTracklets = o2::its::constants::helpers::initArray(); + std::array mClusters = o2::its::constants::helpers::initArray(); + o2::its::TimeEstBC mTime; +}; + +template +class CellSeedTpl final : public SeedBase +{ + using Base = SeedBase; + + public: + GPUhdDefault() CellSeedTpl() = default; + GPUhd() CellSeedTpl(int innerL, int cl0, int cl1, int cl2, int trkl0, int trkl1, const TrackParT& tpc, float chi2, const o2::its::TimeEstBC& time) + : CellSeedTpl(LayerMask(innerL, innerL + 1, innerL + 2), cl0, cl1, cl2, trkl0, trkl1, tpc, chi2, time) + { + } + GPUhd() CellSeedTpl(LayerMask hitLayerMask, int cl0, int cl1, int cl2, int trkl0, int trkl1, const TrackParT& tpc, float chi2, const o2::its::TimeEstBC& time) + : Base(tpc, chi2, 1, time) + { + this->setHitLayerMask(hitLayerMask); + auto& clusters = this->clustersRaw(); + clusters[0] = cl0; + clusters[1] = cl1; + clusters[2] = cl2; + this->setFirstTrackletIndex(trkl0); + this->setSecondTrackletIndex(trkl1); + } + GPUhdDefault() CellSeedTpl(const CellSeedTpl&) = default; + GPUhdDefault() ~CellSeedTpl() = default; + GPUhdDefault() CellSeedTpl(CellSeedTpl&&) = default; + GPUhdDefault() CellSeedTpl& operator=(const CellSeedTpl&) = default; + GPUhdDefault() CellSeedTpl& operator=(CellSeedTpl&&) = default; + + GPUhd() int getFirstClusterIndex() const { return this->clustersRaw()[0]; }; + GPUhd() int getSecondClusterIndex() const { return this->clustersRaw()[1]; }; + GPUhd() int getThirdClusterIndex() const { return this->clustersRaw()[2]; }; + GPUhd() auto& getClusters() { return this->clustersRaw(); } + GPUhd() const auto& getClusters() const { return this->clustersRaw(); } + GPUhd() int getCluster(int layer) const + { + const int slot = this->getHitLayerMask().slot(layer); + return (slot >= 0 && slot < o2::its::constants::ClustersPerCell) ? this->clustersRaw()[slot] : o2::its::constants::UnusedIndex; + } +}; + +/// ITS default: barrel track parameters in cells. +using CellSeed = CellSeedTpl; + +template +class TrackSeedTpl final : public SeedBase +{ + using Base = SeedBase; + + public: + GPUhdDefault() TrackSeedTpl() = default; + GPUhd() TrackSeedTpl(const CellSeedTpl& cs) + : Base(static_cast(cs), cs.getChi2(), cs.getLevel(), cs.getTimeStamp()) + { + this->setHitLayerMask(cs.getHitLayerMask()); + this->setFirstTrackletIndex(cs.getFirstTrackletIndex()); + this->setSecondTrackletIndex(cs.getSecondTrackletIndex()); + auto& clusters = this->clustersRaw(); + int slot = 0; + const auto hitMask = cs.getHitLayerMask(); + for (int layer = 0; layer < NLayers; ++layer) { + if (hitMask.has(layer)) { + clusters[layer] = cs.getClusters()[slot++]; + } + } + } + GPUhdDefault() TrackSeedTpl(const TrackSeedTpl&) = default; + GPUhdDefault() ~TrackSeedTpl() = default; + GPUhdDefault() TrackSeedTpl(TrackSeedTpl&&) = default; + GPUhdDefault() TrackSeedTpl& operator=(const TrackSeedTpl&) = default; + GPUhdDefault() TrackSeedTpl& operator=(TrackSeedTpl&&) = default; + + GPUhd() int getFirstClusterIndex() const { return getClusterBySlot(0); } + GPUhd() int getSecondClusterIndex() const { return getClusterBySlot(1); } + GPUhd() int getThirdClusterIndex() const { return getClusterBySlot(2); } + GPUhd() auto& getClusters() { return this->clustersRaw(); } + GPUhd() const auto& getClusters() const { return this->clustersRaw(); } + GPUhd() int getCluster(int layer) const { return this->clustersRaw()[layer]; } + + private: + GPUhd() int getClusterBySlot(int requestedSlot) const + { + int slot = 0; + const auto hitMask = this->getHitLayerMask(); + for (int layer = 0; layer < NLayers; ++layer) { + if (hitMask.has(layer)) { + if (slot++ == requestedSlot) { + return this->clustersRaw()[layer]; + } + } + } + return o2::its::constants::UnusedIndex; + } +}; + +template +using TrackSeed = TrackSeedTpl; + +template +struct CATrackTypeHelper { + using type = o2::its::TrackITSExt; +}; + +template +using CATrackType = typename CATrackTypeHelper::type; + +/// Per-detector track parametrization stored in CA cells and extended seeds. +template +struct CASeedTrackPar { + static constexpr o2::detectors::DetID::ID DetId = detIdFromNLayers(); + using type = std::conditional_t; +}; + +template +using CellSeedN = CellSeedTpl::type>; + +template +using TrackSeedN = TrackSeedTpl::type>; + +} // namespace o2::itsmft::tracking + +#endif /* ALICEO2_ITSMFT_TRACKING_INCLUDE_CACELL_H_ */ diff --git a/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/Configuration.h b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/Configuration.h new file mode 100644 index 0000000000000..4275d8c3e8c0c --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/Configuration.h @@ -0,0 +1,219 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file Configuration.h +/// \brief Shared CA tracking configuration for ITS and MFT +/// + +#ifndef ALICEO2_ITSMFT_TRACKING_CONFIGURATION_H_ +#define ALICEO2_ITSMFT_TRACKING_CONFIGURATION_H_ + +#include +#ifndef GPUCA_GPUCODE_DEVICE +#include +#include +#include +#include +#endif + +#include "CommonUtils/EnumFlags.h" +#include "DetectorsBase/Propagator.h" +#include "DetectorsCommonDataFormats/DetID.h" +#include "ITSMFTTracking/LayerMask.h" +#include "ITSMFTTracking/TrackingConfigParam.h" +#include "MFTTracking/Constants.h" +#include "ITStracking/TrackingConfigParam.h" + +namespace o2::itsmft::tracking +{ + +constexpr int nLayersForDet(o2::detectors::DetID::ID detId) +{ + return detId == o2::detectors::DetID::MFT ? o2::mft::constants::mft::LayersNumber : ITSNLayers; +} + +template +constexpr o2::detectors::DetID::ID detIdFromNLayers() +{ + return NLayers == o2::mft::constants::mft::LayersNumber ? o2::detectors::DetID::MFT : o2::detectors::DetID::ITS; +} + +} // namespace o2::itsmft::tracking + +namespace o2::itsmft +{ + +inline constexpr int ClustersPerCell = 3; + +// Steering of dedicated steps in an iteration +enum class IterationStep : uint8_t { + FirstPass = 0, + RebuildClusterLUT, + UseUPCMask, + SelectUPCVertices, + ResetVertices, + SkipROFsAboveThreshold, + MarkVerticesAsUPC, +}; +using IterationSteps = o2::utils::EnumFlags; + +struct TrackingParameters { + int CellMinimumLevel() const noexcept + { + const int minClusters = MinTrackLength - (MaxHoles > 0 ? MaxHoles : 0); + const int effectiveMinClusters = minClusters > ClustersPerCell ? minClusters : ClustersPerCell; + return effectiveMinClusters - ClustersPerCell + 1; + } + int NeighboursPerRoad() const noexcept { return NLayers - 3; } + int CellsPerRoad() const noexcept { return NLayers - 2; } + int TrackletsPerRoad() const noexcept { return NLayers - 1; } + std::string asString() const; + + IterationSteps PassFlags{IterationStep::FirstPass, IterationStep::RebuildClusterLUT}; + int NLayers = tracking::ITSNLayers; + std::vector AddTimeError = {0, 0, 0, 0, 0, 0, 0}; + std::vector LayerZ = {16.333f + 1, 16.333f + 1, 16.333f + 1, 42.140f + 1, 42.140f + 1, 73.745f + 1, 73.745f + 1}; + std::vector LayerColHalfExtent{}; // index-table column half extent (ITS: z, MFT: global x); falls back to LayerZ + float IndexRowMin{0.f}; // index-table row origin (MFT: global y min; unused for ITS phi-z) + float IndexRowMax{0.f}; // index-table row span end (MFT: global y max; 0 => TwoPI for ITS) + std::vector LayerRadii = {2.33959f, 3.14076f, 3.91924f, 19.6213f, 24.5597f, 34.388f, 39.3329f}; + std::vector LayerxX0 = {5.e-3f, 5.e-3f, 5.e-3f, 1.e-2f, 1.e-2f, 1.e-2f, 1.e-2f}; + std::vector LayerResolution = {5.e-4f, 5.e-4f, 5.e-4f, 5.e-4f, 5.e-4f, 5.e-4f, 5.e-4f}; + std::vector SystError2Row = {0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f}; // systematic error^2 along local row (ALPIDE X) per layer + std::vector SystError2Col = {0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f}; // systematic error^2 along local column (ALPIDE Z) per layer + int ColBins{256}; // ITS: ZBins + int RowBins{128}; // ITS: PhiBins + bool UseDiamond = false; + float Diamond[3] = {0.f, 0.f, 0.f}; + float DiamondCov[6] = {25.e-6f, 0.f, 0.f, 25.e-6f, 0.f, 36.f}; + + /// General parameters + int MinTrackLength = 7; + int MaxHoles = 0; + tracking::LayerMask HoleLayerMask = 0; + float NSigmaCut = 5; + float PVres = 1.e-2f; + /// Trackleting cuts + float TrackletMinPt = 0.3f; + /// Cell finding cuts + float CellDeltaTanLambdaSigma = 0.007f; + float CellRoadRCut = 0.05f; // MFT: max distance to seed line (classic ROADclsRCut) + float TrackletMinAbsX = 0.f; // MFT: reject clusters/tracks with |x| below this (cm); 0 = disabled + /// Fitter parameters + o2::base::PropagatorImpl::MatCorrType CorrType = o2::base::PropagatorImpl::MatCorrType::USEMatCorrNONE; + float MaxChi2ClusterAttachment = 60.f; + float MaxChi2NDF = 30.f; + int ReseedIfShorter = 6; // reseed for the final fit track with the length shorter than this + std::vector MinPt = {0.f, 0.f, 0.f, 0.f}; + tracking::LayerMask StartLayerMask = 0x7F; + bool RepeatRefitOut = false; // repeat outward refit using inward refit as a seed + bool ShiftRefToCluster = true; // TrackFit: after update shift the linearization reference to cluster + bool PerPrimaryVertexProcessing = false; + bool SaveTimeBenchmarks = false; + bool DoUPCIteration = false; + bool FataliseUponFailure = true; + bool CreateArtefactLabels{false}; + bool PrintMemory = false; // print allocator usage in epilog report + size_t MaxMemory = std::numeric_limits::max(); + bool DropTFUponFailure = false; + + // Selections on tracks sharing clusters + bool AllowSharingFirstCluster = false; + float SharedClusterMaxDeltaPhi = 0.05f; // For tracks sharing clusters, maximum allowed delta phi at the cluster position + float SharedClusterMaxDeltaEta = 0.03f; // For tracks sharing clusters, maximum allowed delta eta at the cluster position + bool SharedClusterOppositeSign = false; // For tracks sharing clusters, require opposite sign of the tracklets + int SharedMaxClusters = 0; // Maximal allowed shared clusters (excluding first cluster) +}; + +/// Reset tracking parameters to detector geometry defaults (ITS: struct defaults; MFT: MFTTracking/Constants.h). +void resetDetectorDefaults(TrackingParameters& params, o2::detectors::DetID::ID detId); + +namespace TrackingMode +{ +enum Type : int8_t { + Unset = -1, + Sync = 0, + Async = 1, + Cosmics = 2, + Off = 3, +}; + +Type fromString(std::string_view str); +std::string toString(Type mode); +std::vector getTrackingParameters(o2::detectors::DetID::ID detId, Type mode); + +} // namespace TrackingMode + +struct VertexingParameters { + std::string asString() const; + + IterationSteps PassFlags{IterationStep::FirstPass, IterationStep::ResetVertices}; + std::vector LayerZ = {16.333f + 1, 16.333f + 1, 16.333f + 1, 42.140f + 1, 42.140f + 1, 73.745f + 1, 73.745f + 1}; + std::vector LayerRadii = {2.33959f, 3.14076f, 3.91924f, 19.6213f, 24.5597f, 34.388f, 39.3329f}; + int vertPerRofThreshold = 0; // Maximum number of vertices per ROF to trigger second a round + int ColBins = 1; + int RowBins = 128; + float zCut = -1.f; + float phiCut = -1.f; + float pairCut = -1.f; + float clusterCut = -1.f; + float coarseZWindow = -1.f; + float seedDedupZCut = -1.f; + float refitDedupZCut = -1.f; + float duplicateZCut = -1.f; + float finalSelectionZCut = -1.f; + float duplicateDistance2Cut = -1.f; + float tanLambdaCut = -1.f; + float NSigmaCut = -1; + float maxZPositionAllowed = -1.f; + int clusterContributorsCut = -1; + int suppressLowMultDebris = -1; + int seedMemberRadiusTime = -1; + int seedMemberRadiusZ = -1; + int maxTrackletsPerCluster = -1; + int phiSpan = -1; + int zSpan = -1; + bool SaveTimeBenchmarks = false; + + bool useTruthSeeding = false; // overwrite found vertices with MC events + + int nThreads = 1; + bool PrintMemory = false; // print allocator usage in epilog report + size_t MaxMemory = std::numeric_limits::max(); + bool DropTFUponFailure = false; +}; + +} // namespace o2::itsmft + +namespace o2::itsmft::tracking +{ + +/// MFT uses o2::itsmft::TrackerParamConfig; ITS production params stay in O2::ITStracking. +template +struct TrackerParamRef; + +template <> +struct TrackerParamRef { + using Type = o2::itsmft::TrackerParamConfig; + static const Type& get() { return Type::Instance(); } + static constexpr int nLayers() { return Type::getNLayers(); } +}; + +template <> +struct TrackerParamRef { + using Type = o2::its::TrackerParamConfig; + static const Type& get() { return Type::Instance(); } + static constexpr int nLayers() { return ITSNLayers; } +}; + +} // namespace o2::itsmft::tracking + +#endif /* ALICEO2_ITSMFT_TRACKING_CONFIGURATION_H_ */ diff --git a/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/DetectorTraits.h b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/DetectorTraits.h new file mode 100644 index 0000000000000..7f7911b80690e --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/DetectorTraits.h @@ -0,0 +1,75 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file DetectorTraits.h +/// \brief Minimal detector-specific hooks for shared CA tracking (ITS vs MFT) +/// + +#ifndef ALICEO2_ITSMFT_TRACKING_DETECTOR_TRAITS_H_ +#define ALICEO2_ITSMFT_TRACKING_DETECTOR_TRAITS_H_ + +#include "DetectorsBase/Propagator.h" +#include "DetectorsCommonDataFormats/DetID.h" +#include "DataFormatsCalibration/MeanVertexObject.h" +#include "ITSMFTTracking/Configuration.h" +#include "ITSMFTTracking/IndexTableUtils.h" +#include "ITSMFTTracking/TimeFrame.h" +#include "ITSMFTTracking/Cell.h" +#include "ITStracking/BoundedAllocator.h" +#include "ITStracking/Cluster.h" +#include "ITStracking/Tracklet.h" + +namespace o2::itsmft::tracking +{ + +/// Per-detector differences in refit, track acceptance, and index-table setup. +/// Everything else stays in TrackerTraits and matches ITS line-for-line. +template +struct DetectorTraits { + using TrackType = CATrackType; + using TrackSeedN = o2::itsmft::tracking::TrackSeedN; + using CellSeedN = o2::itsmft::tracking::CellSeedN; + using TimeFrameN = TimeFrame; + static constexpr o2::detectors::DetID::ID DetId = detIdFromNLayers(); + + static bool refitSeed(const TrackSeedN& seed, + TrackType& track, + const TrackingParameters& params, + float bz, + TimeFrameN& tf, + const o2::its::TrackingFrameInfo* const tfInfos[NLayers], + const o2::its::Cluster* const unsortedClusters[NLayers], + const o2::base::PropagatorImpl* propagator); + + static void configureIndexTableUtils(IndexTableUtils& utils, const TrackingParameters& params); + + /// ITS: barrel χ² between connected cells; MFT: forward χ² between connected cells. + static bool cellsAreCompatible(const CellSeedN& currentCell, + const CellSeedN& nextCell, + const TimeFrameN& tf, + const TrackingParameters& params, + float bz); +}; + +template +struct TrackingLoadPolicy { + static void configureBeamPosition(TimeFrame& tf, + const TrackingParameters& p, + const o2::dataformats::MeanVertexObject* meanVertex, + bool overrideBeamEstimation); +}; + +template +using TrackingLoadPolicyN = TrackingLoadPolicy(), NLayers>; + +} // namespace o2::itsmft::tracking + +#endif /* ALICEO2_ITSMFT_TRACKING_DETECTOR_TRAITS_H_ */ diff --git a/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/IOUtils.h b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/IOUtils.h new file mode 100644 index 0000000000000..348f74ff78214 --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/IOUtils.h @@ -0,0 +1,116 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file IOUtils.h +/// \brief Shared cluster I/O utilities for ITS and MFT (based on ITStracking/IOUtils.h) +/// + +#ifndef ALICEO2_ITSMFT_TRACKING_IOUTILS_H_ +#define ALICEO2_ITSMFT_TRACKING_IOUTILS_H_ + +#include +#include + +#include + +#include "DetectorsCommonDataFormats/DetID.h" +#include "ITSMFTBase/SegmentationAlpide.h" +#include "ReconstructionDataFormats/BaseCluster.h" +#include "DataFormatsITSMFT/ClusterPattern.h" +#include "DataFormatsITSMFT/CompCluster.h" +#include "DataFormatsITSMFT/TopologyDictionary.h" +#include "MathUtils/Cartesian.h" + +namespace o2::its +{ +struct TrackingFrameInfo; +} + +namespace o2::itsmft::ioutils +{ + +constexpr float DefClusErrorRow = o2::itsmft::SegmentationAlpide::PitchRow * 0.5f; +constexpr float DefClusErrorCol = o2::itsmft::SegmentationAlpide::PitchCol * 0.5f; +constexpr float DefClusError2Row = DefClusErrorRow * DefClusErrorRow; +constexpr float DefClusError2Col = DefClusErrorCol * DefClusErrorCol; + +void fillMatrixCache(o2::detectors::DetID::ID detId); +int getClusterLayer(o2::detectors::DetID::ID detId, const CompClusterExt& cluster); + +/// Decode a compact cluster into layer, size, and a TrackingFrameInfo (global + local frame). +template +void loadClusterTrackingFrameInfo(const CompClusterExt& c, + gsl::span::iterator& pattIt, + const TopologyDictionary* dict, + int& layer, + unsigned int& clusterSize, + o2::its::TrackingFrameInfo& tfInfo, + bool applySysErrors = true); + +/// Convert compact clusters to 3D spacepoints. +/// \tparam DetId o2::detectors::DetID::ITS or DetID::MFT +template +void convertCompactClusters(gsl::span clusters, + gsl::span::iterator& pattIt, + std::vector>& output, + const TopologyDictionary* dict); + +template +o2::math_utils::Point3D extractClusterData(const CompClusterExt& c, iterator& iter, const TopologyDictionary* dict, T& sig2Row, T& sig2Col, unsigned int* clusterSize = nullptr) +{ + auto pattID = c.getPatternID(); + sig2Row = DefClusError2Row; + sig2Col = DefClusError2Col; // Dummy COG errors (about half pixel size) + if (pattID != CompCluster::InvalidPatternID) { + sig2Row = dict->getErr2X(pattID); + sig2Col = dict->getErr2Z(pattID); + if (!dict->isGroup(pattID)) { + if (clusterSize != nullptr) { + *clusterSize = dict->getNpixels(pattID); + } + return dict->getClusterCoordinates(c); + } + ClusterPattern patt(iter); + if (clusterSize != nullptr) { + *clusterSize = patt.getNPixels(); + } + return dict->getClusterCoordinates(c, patt); + } + ClusterPattern patt(iter); + if (clusterSize != nullptr) { + *clusterSize = patt.getNPixels(); + } + return dict->getClusterCoordinates(c, patt, false); +} + +// same method returning coordinates as an array (suitable for the TGeoMatrix) +template +std::array extractClusterDataA(const CompClusterExt& c, iterator& iter, const TopologyDictionary* dict, T& sig2Row, T& sig2Col) +{ + auto pattID = c.getPatternID(); + sig2Row = DefClusError2Row; + sig2Col = DefClusError2Col; // Dummy COG errors (about half pixel size) + if (pattID != CompCluster::InvalidPatternID) { + sig2Row = dict->getErr2X(pattID); + sig2Col = dict->getErr2Z(pattID); + if (!dict->isGroup(pattID)) { + return dict->getClusterCoordinatesA(c); + } + ClusterPattern patt(iter); + return dict->getClusterCoordinatesA(c, patt); + } + ClusterPattern patt(iter); + return dict->getClusterCoordinatesA(c, patt, false); +} + +} // namespace o2::itsmft::ioutils + +#endif /* ALICEO2_ITSMFT_TRACKING_IOUTILS_H_ */ diff --git a/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/IndexTableUtils.h b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/IndexTableUtils.h new file mode 100644 index 0000000000000..3bf00c727d272 --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/IndexTableUtils.h @@ -0,0 +1,262 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file IndexTableUtils.h +/// \brief Shared index-table utilities for ITS (phi-z) and MFT (x-y) +/// + +#ifndef ALICEO2_ITSMFT_TRACKING_INDEXTABLEUTILS_H_ +#define ALICEO2_ITSMFT_TRACKING_INDEXTABLEUTILS_H_ + +#include + +#include "CommonConstants/MathConstants.h" +#include "GPUCommonMath.h" +#include "GPUCommonDef.h" +#include "ITSMFTTracking/Configuration.h" +#include "ITStracking/Cluster.h" +#include "MFTTracking/Constants.h" + +namespace o2::itsmft +{ + +enum class IndexTableCoordType : uint8_t { PhiZ, XY }; + +namespace index_table_utils +{ +GPUhdi() float getNormalizedPhi(float phi) +{ + phi -= o2::constants::math::TwoPI * o2::gpu::GPUCommonMath::Floor(phi * (1.f / o2::constants::math::TwoPI)); + return phi; +} +} // namespace index_table_utils + +/// Row/column LUT helper (ITS: row=phi, col=z; MFT: row=y, col=x). +template +class IndexTableUtils +{ + public: + /// Configure LUT geometry. ITS (PhiZ): row = phi [0, TwoPI), col = z; MFT (XY): row = y, col = x. + void setIndexTableParams(IndexTableCoordType coordType, int nRowBins, int nColBins, + float rowMin, float rowMax, + const std::array& layerColHalfExtent) + { + mCoordType = coordType; + mRowOrigin = (coordType == IndexTableCoordType::PhiZ) ? 0.f : rowMin; + mRowCoordinateSpan = rowMax - rowMin; + mInverseRowBinSize = (mRowCoordinateSpan > 0.f) ? static_cast(nRowBins) / mRowCoordinateSpan : 0.f; + mNcolBins = nColBins; + mNrowBins = nRowBins; + for (int iLayer{0}; iLayer < nLayers; ++iLayer) { + mLayerColHalfExtent[iLayer] = layerColHalfExtent[iLayer]; + mInverseColBinSize[iLayer] = 0.5f * nColBins / layerColHalfExtent[iLayer]; + } + } + + /// Fill LUT geometry from any struct exposing RowBins, ColBins and LayerZ (ITS phi-z). + template + void setTrackingParameters(const T& params) + { + setIndexTableParams(IndexTableCoordType::PhiZ, params.RowBins, params.ColBins, + 0.f, o2::constants::math::TwoPI, layerColHalfExtentFrom(params)); + } + + /// Fill LUT geometry for MFT (row = global y, col = global x). + template + void setTrackingParametersXY(const T& params, float rowMin, float rowMax) + { + setIndexTableParams(IndexTableCoordType::XY, params.RowBins, params.ColBins, + rowMin, rowMax, layerColHalfExtentFrom(params)); + } + + GPUhdi() float getInverseColCoordinate(const int layerIndex) const + { + return 0.5f * mNcolBins / mLayerColHalfExtent[layerIndex]; + } + + GPUhdi() int getColBinIndex(const int layerIndex, const float colCoordinate) const + { + return (colCoordinate + mLayerColHalfExtent[layerIndex]) * mInverseColBinSize[layerIndex]; + } + + GPUhdi() int getRowBinIndex(const float rowCoordinate) const + { + if (mCoordType == IndexTableCoordType::PhiZ) { + return rowCoordinate * mInverseRowBinSize; + } + return (rowCoordinate - mRowOrigin) * mInverseRowBinSize; + } + + GPUhdi() int getBinIndex(const int colIndex, const int rowIndex) const + { + return o2::gpu::GPUCommonMath::Min(rowIndex * mNcolBins + colIndex, (mNcolBins * mNrowBins) - 1); + } + + GPUhdi() int countRowSelectedBins(const int* indexTable, const int rowBinIndex, + const int minColBinIndex, const int maxColBinIndex) const + { + const int firstBinIndex{getBinIndex(minColBinIndex, rowBinIndex)}; + const int maxBinIndex{firstBinIndex + maxColBinIndex - minColBinIndex + 1}; + + return indexTable[maxBinIndex] - indexTable[firstBinIndex]; + } + + void print() const; + + GPUhdi() int getNcolBins() const { return mNcolBins; } + GPUhdi() int getNrowBins() const { return mNrowBins; } + GPUhdi() float getLayerColHalfExtent(int i) const { return mLayerColHalfExtent[i]; } + GPUhdi() void setNcolBins(const int colBins) { mNcolBins = colBins; } + GPUhdi() void setNrowBins(const int rowBins) { mNrowBins = rowBins; } + GPUhdi() IndexTableCoordType getCoordType() const { return mCoordType; } + + private: + template + static std::array layerColHalfExtentFrom(const T& params) + { + std::array extents{}; + if constexpr (requires { params.LayerColHalfExtent; }) { + const auto& colExtents = params.LayerColHalfExtent.empty() ? params.LayerZ : params.LayerColHalfExtent; + for (int iLayer{0}; iLayer < nLayers && iLayer < static_cast(colExtents.size()); ++iLayer) { + extents[iLayer] = colExtents[iLayer]; + } + } else { + for (int iLayer{0}; iLayer < nLayers && iLayer < static_cast(params.LayerZ.size()); ++iLayer) { + extents[iLayer] = params.LayerZ[iLayer]; + } + } + return extents; + } + + int mNcolBins = 0; + int mNrowBins = 0; + float mInverseRowBinSize = 0.f; + float mRowOrigin = 0.f; + float mRowCoordinateSpan = o2::constants::math::TwoPI; + IndexTableCoordType mCoordType{IndexTableCoordType::PhiZ}; + std::array mLayerColHalfExtent{}; + std::array mInverseColBinSize{}; +}; + +template +void IndexTableUtils::print() const +{ + printf("NcolBins: %d, NrowBins: %d, InverseRowBinSize: %f\n", mNcolBins, mNrowBins, mInverseRowBinSize); + for (int iLayer{0}; iLayer < nLayers; ++iLayer) { + printf("Layer %d: ColHalfExtent: %f, InverseColBinSize: %f\n", iLayer, mLayerColHalfExtent[iLayer], mInverseColBinSize[iLayer]); + } +} + +/// ITS: row = phi, col = z. +template +GPUhdi() int4 getBinsPhiZ(float phi, const int layerIndex, + float z1, float z2, float maxDeltaCol, float maxDeltaRow, + const IndexTableUtils& utils) +{ + const float colRangeMin = o2::gpu::GPUCommonMath::Min(z1, z2) - maxDeltaCol; + const float rowRangeMin = (maxDeltaRow > o2::constants::math::PI) ? 0.f : phi - maxDeltaRow; + const float colRangeMax = o2::gpu::GPUCommonMath::Max(z1, z2) + maxDeltaCol; + const float rowRangeMax = (maxDeltaRow > o2::constants::math::PI) ? o2::constants::math::TwoPI : phi + maxDeltaRow; + + if (colRangeMax < -utils.getLayerColHalfExtent(layerIndex) || + colRangeMin > utils.getLayerColHalfExtent(layerIndex) || colRangeMin > colRangeMax) { + return int4{-1, -1, -1, -1}; + } + + return int4{o2::gpu::GPUCommonMath::Max(0, utils.getColBinIndex(layerIndex, colRangeMin)), + utils.getRowBinIndex(index_table_utils::getNormalizedPhi(rowRangeMin)), + o2::gpu::GPUCommonMath::Min(utils.getNcolBins() - 1, utils.getColBinIndex(layerIndex, colRangeMax)), + utils.getRowBinIndex(index_table_utils::getNormalizedPhi(rowRangeMax))}; +} + +/// MFT: row = y, col = x. +template +GPUhdi() int4 getBinsXY(float x, float y, const int layerIndex, + float x1, float x2, float y1, float y2, + float maxDeltaCol, float maxDeltaRow, + const IndexTableUtils& utils) +{ + const float colRangeMin = o2::gpu::GPUCommonMath::Min(x1, x2) - maxDeltaCol; + const float rowRangeMin = o2::gpu::GPUCommonMath::Min(y1, y2) - maxDeltaRow; + const float colRangeMax = o2::gpu::GPUCommonMath::Max(x1, x2) + maxDeltaCol; + const float rowRangeMax = o2::gpu::GPUCommonMath::Max(y1, y2) + maxDeltaRow; + + const float colHalf = utils.getLayerColHalfExtent(layerIndex); + if (colRangeMax < -colHalf || colRangeMin > colHalf || colRangeMin > colRangeMax) { + return int4{-1, -1, -1, -1}; + } + + return int4{o2::gpu::GPUCommonMath::Max(0, utils.getColBinIndex(layerIndex, colRangeMin)), + o2::gpu::GPUCommonMath::Max(0, utils.getRowBinIndex(rowRangeMin)), + o2::gpu::GPUCommonMath::Min(utils.getNcolBins() - 1, utils.getColBinIndex(layerIndex, colRangeMax)), + utils.getRowBinIndex(rowRangeMax)}; +} + +/// MFT: extrapolate cluster to toLayer on the line through the primary vertex. +template +GPUhdi() void mftConeProject(const o2::its::Cluster& cluster, int fromLayer, int toLayer, + float pvX, float pvY, float pvZ, float& xProj, float& yProj) +{ + const auto& layerZ = o2::mft::constants::mft::LayerZCoordinate(); + const float zFrom = layerZ[fromLayer]; + const float zTo = layerZ[toLayer]; + const float dz0 = zFrom - pvZ; + if (o2::gpu::GPUCommonMath::Abs(dz0) < 1e-6f) { + xProj = cluster.xCoordinate; + yProj = cluster.yCoordinate; + return; + } + const float w = (zTo - pvZ) / dz0; + xProj = pvX + w * (cluster.xCoordinate - pvX); + yProj = pvY + w * (cluster.yCoordinate - pvY); +} + +/// MFT LUT window around a precomputed (x, y) projection on toLayer. +template +GPUhdi() int4 getBinsRectClusterAtProj(float xProj, float yProj, int toLayer, + float colRangeMin, float colRangeMax, float maxDeltaCol, float maxDeltaRow, + const IndexTableUtils& utils) +{ + const float rProj = o2::gpu::GPUCommonMath::Hypot(xProj, yProj); + float x1 = xProj; + float x2 = xProj; + float y1 = yProj; + float y2 = yProj; + if (rProj > 0.f) { + const float invRProj = 1.f / rProj; + x1 = colRangeMin * xProj * invRProj; + x2 = colRangeMax * xProj * invRProj; + y1 = colRangeMin * yProj * invRProj; + y2 = colRangeMax * yProj * invRProj; + } + return getBinsXY(xProj, yProj, toLayer, x1, x2, y1, y2, maxDeltaCol, maxDeltaRow, utils); +} + +/// Cluster-driven LUT window: phi-z for ITS, projected x-y for MFT. +/// ITS: colRangeMin/Max = z window; MFT: colRangeMin/Max = rMin/rMax at toLayer from diamond z spread. +template +GPUhdi() int4 getBinsRectCluster(const o2::its::Cluster& cluster, int fromLayer, int toLayer, + float colRangeMin, float colRangeMax, float maxDeltaCol, float maxDeltaRow, + const IndexTableUtils& utils, + float pvX = 0.f, float pvY = 0.f, float pvZ = 0.f) +{ + if (utils.getCoordType() == IndexTableCoordType::XY) { + float xProj = 0.f; + float yProj = 0.f; + mftConeProject(cluster, fromLayer, toLayer, pvX, pvY, pvZ, xProj, yProj); + return getBinsRectClusterAtProj(xProj, yProj, toLayer, colRangeMin, colRangeMax, maxDeltaCol, maxDeltaRow, utils); + } + return getBinsPhiZ(cluster.phi, toLayer, colRangeMin, colRangeMax, maxDeltaCol, maxDeltaRow, utils); +} + +} // namespace o2::itsmft + +#endif /* ALICEO2_ITSMFT_TRACKING_INDEXTABLEUTILS_H_ */ diff --git a/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/LayerMask.h b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/LayerMask.h new file mode 100644 index 0000000000000..e9afc26a916bd --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/LayerMask.h @@ -0,0 +1,115 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef ALICEO2_ITSMFT_TRACKING_INCLUDE_LAYERMASK_H_ +#define ALICEO2_ITSMFT_TRACKING_INCLUDE_LAYERMASK_H_ + +#include +#include + +#ifndef GPUCA_GPUCODE +#include +#include +#endif + +#include "GPUCommonDef.h" +#include "GPUCommonMath.h" +#include "ITStracking/Constants.h" + +namespace o2::itsmft::tracking +{ + +struct LayerMask { + GPUhdDefault() constexpr LayerMask() noexcept = default; + GPUhdDefault() constexpr LayerMask(uint16_t mask) noexcept : mBits{mask} {} + GPUhdDefault() constexpr LayerMask(int layer0, int layer1, int layer2) noexcept + : mBits{static_cast((uint16_t(1) << layer0) | (uint16_t(1) << layer1) | (uint16_t(1) << layer2))} + { + } + GPUhdi() constexpr operator uint16_t() const noexcept { return mBits; } + GPUhdi() constexpr uint16_t value() const noexcept { return mBits; } + GPUhdi() constexpr void set(int layer) noexcept { mBits |= (uint16_t(1) << layer); } + + GPUhdi() LayerMask operator~() const noexcept { return LayerMask{static_cast(~mBits)}; } + GPUhdi() LayerMask operator&(LayerMask other) const noexcept { return LayerMask{static_cast(mBits & other.mBits)}; } + GPUhdi() LayerMask operator|(LayerMask other) const noexcept { return LayerMask{static_cast(mBits | other.mBits)}; } + GPUhdi() LayerMask& operator&=(LayerMask other) noexcept + { + mBits &= other.mBits; + return *this; + } + GPUhdi() LayerMask& operator|=(LayerMask other) noexcept + { + mBits |= other.mBits; + return *this; + } + + GPUhdi() bool empty() const noexcept { return mBits == 0; } + GPUhdi() bool has(int layer) const noexcept { return mBits & (uint16_t(1) << layer); } + GPUhdi() bool isSubsetOf(LayerMask allowed) const noexcept { return (*this & ~allowed).empty(); } + GPUhdi() bool isAllowedHoleMask(int maxHoles, LayerMask allowedHoleMask) const noexcept + { + const int allowedHoles = maxHoles > 0 ? maxHoles : 0; + return count() <= allowedHoles && isSubsetOf(allowedHoleMask); + } + GPUhdi() bool isAllowed(int maxHoles, LayerMask allowedHoleMask) const noexcept + { + return holeMask().isAllowedHoleMask(maxHoles, allowedHoleMask); + } + GPUhdi() int length() const noexcept { return empty() ? 0 : last() - first() + 1; } + GPUhdi() int count() const noexcept { return static_cast(o2::gpu::GPUCommonMath::Popcount(mBits)); } + GPUhdi() int first() const noexcept { return mBits ? static_cast(o2::gpu::GPUCommonMath::Ctz(mBits)) : o2::its::constants::UnusedIndex; } + GPUhdi() int last() const noexcept { return mBits ? 31 - static_cast(o2::gpu::GPUCommonMath::Clz(mBits)) : o2::its::constants::UnusedIndex; } + GPUhdi() LayerMask holeMask() const noexcept + { + return empty() ? LayerMask{0} : (span(first(), last()) & ~(*this)); + } + + GPUhdi() int slot(int layer) const noexcept + { + if (!has(layer)) { + return o2::its::constants::UnusedIndex; + } + const uint32_t lowerLayers = (uint32_t(1) << layer) - 1; + return static_cast(o2::gpu::GPUCommonMath::Popcount(static_cast(mBits) & lowerLayers)); + } + + static GPUhdi() LayerMask span(int fromLayer, int toLayer) noexcept + { + if (fromLayer > toLayer) { + return 0; + } + const uint32_t upper = (uint32_t(1) << (toLayer + 1)) - 1; + const uint32_t lower = (uint32_t(1) << fromLayer) - 1; + return static_cast(upper & ~lower); + } + + static GPUhdi() LayerMask skipped(int fromLayer, int toLayer) noexcept + { + return (toLayer - fromLayer <= 1) ? LayerMask{0} : span(fromLayer + 1, toLayer - 1); + } + +#ifndef GPUCA_GPUCODE + std::string asString() const { return fmt::format("{:016b}", mBits); } +#endif + + private: + uint16_t mBits{0}; +}; + +static_assert(std::is_standard_layout_v); +static_assert(std::is_trivially_copyable_v); +static_assert(sizeof(LayerMask) == sizeof(uint16_t)); +static_assert(alignof(LayerMask) == alignof(uint16_t)); + +} // namespace o2::itsmft::tracking + +#endif /* ALICEO2_ITSMFT_TRACKING_INCLUDE_LAYERMASK_H_ */ diff --git a/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/MFTCATrack.h b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/MFTCATrack.h new file mode 100644 index 0000000000000..54b88c948cac6 --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/MFTCATrack.h @@ -0,0 +1,157 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file MFTCATrack.h +/// \brief MFT CA track with native forward parameters and layer-indexed cluster refs +/// + +#ifndef ALICEO2_ITSMFT_TRACKING_MFTCATRACK_H_ +#define ALICEO2_ITSMFT_TRACKING_MFTCATRACK_H_ + +#include +#include + +#include "DataFormatsITS/TimeEstBC.h" +#include "DataFormatsITS/TrackITS.h" +#include "DataFormatsMFT/TrackMFT.h" +#include "ITSMFTTracking/Cell.h" +#include "ITStracking/Constants.h" +#include "MFTTracking/Constants.h" + +namespace o2::itsmft::tracking +{ + +/// MFT CA output track: forward fit parameters plus ITS-style cluster indexing for CA bookkeeping. +class MFTCATrack +{ + enum UserBits { + kSharedClusters = 1 << 28 + }; + + public: + static constexpr int MaxClusters = o2::its::TrackITSExt::MaxClusters; + + MFTCATrack() + { + mIndex.fill(o2::its::constants::UnusedIndex); + } + + const o2::mft::TrackMFT& getTrack() const { return mTrack; } + o2::mft::TrackMFT& getTrack() { return mTrack; } + + float getChi2() const { return static_cast(mTrack.getTrackChi2()); } + float getPhi() const { return static_cast(mTrack.getPhi()); } + float getEta() const { return static_cast(mTrack.getEta()); } + double getCharge() const { return mTrack.getCharge(); } + + int getNumberOfClusters() const + { + int n{0}; + for (int layer = 0; layer < MaxClusters; ++layer) { + if (hasHitOnLayer(layer)) { + ++n; + } + } + return n; + } + + uint32_t getPattern() const { return mPattern; } + void setPattern(uint32_t p) { mPattern = p; } + bool hasHitOnLayer(uint32_t layer) const { return mPattern & (0x1u << layer); } + + uint32_t getFirstClusterLayer() const + { + uint32_t s{0}; + while (!(mPattern & (1u << s)) && s < MaxClusters) { + ++s; + } + return s; + } + + uint32_t getLastClusterLayer() const + { + uint32_t r{0}, v{mPattern & ((1u << MaxClusters) - 1)}; + while (v >>= 1) { + ++r; + } + return r; + } + + int getClusterIndex(int layer) const { return mIndex[layer]; } + + int getFirstLayerClusterIndex() const { return getClusterIndex(getFirstClusterLayer()); } + + void setClusterIndex(int layer, int idx) + { + mIndex[layer] = idx; + if (idx != o2::its::constants::UnusedIndex) { + mPattern |= 0x1u << layer; + } + } + + void setExternalClusterIndex(int layer, int idx, bool newCluster = false) + { + if (newCluster) { + mPattern |= 0x1u << layer; + } + mIndex[layer] = idx; + } + + void setClusterSize(int layer, int size) + { + if (layer >= 10) { + return; + } + if (size > 63) { + size = 63; + } + mClusterSizes &= ~(0x3fULL << (layer * 6)); + mClusterSizes |= (static_cast(size) << (layer * 6)); + } + + int getClusterSize(int layer) const + { + if (layer >= 10) { + return 0; + } + return (mClusterSizes >> (layer * 6)) & 0x3f; + } + + auto& getTimeStamp() { return mTime; } + const auto& getTimeStamp() const { return mTime; } + + uint16_t getSeedPattern() const { return mSeedPattern; } + void setSeedPattern(uint16_t pattern) { mSeedPattern = pattern; } + + void setSharedClusters(bool toggle = true) + { + mClusterSizes = toggle ? (mClusterSizes | kSharedClusters) : (mClusterSizes & ~kSharedClusters); + } + + bool hasSharedClusters() const { return mClusterSizes & kSharedClusters; } + + private: + o2::mft::TrackMFT mTrack; + std::array mIndex = {}; + uint32_t mPattern = 0; + o2::its::TimeStamp mTime; + uint16_t mSeedPattern{0}; + uint64_t mClusterSizes = 0; +}; + +template <> +struct CATrackTypeHelper { + using type = MFTCATrack; +}; + +} // namespace o2::itsmft::tracking + +#endif /* ALICEO2_ITSMFT_TRACKING_MFTCATRACK_H_ */ diff --git a/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/MFTFwdTrackHelpers.h b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/MFTFwdTrackHelpers.h new file mode 100644 index 0000000000000..7302fb004109a --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/MFTFwdTrackHelpers.h @@ -0,0 +1,331 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file MFTFwdTrackHelpers.h +/// \brief Forward-track helpers for MFT CA cell fitting and final track refit +/// + +#ifndef ALICEO2_ITSMFT_TRACKING_MFTFWDTRACKHELPERS_H_ +#define ALICEO2_ITSMFT_TRACKING_MFTFWDTRACKHELPERS_H_ + +#include +#include +#include + +#include "CommonConstants/MathConstants.h" +#include "ITSMFTTracking/Cell.h" +#include "ITSMFTTracking/Configuration.h" +#include "ITSMFTTracking/MFTCATrack.h" +#include "ITSMFTTracking/TimeFrame.h" +#include "ITStracking/Cluster.h" +#include "ITStracking/Constants.h" +#include "MFTTracking/Constants.h" +#include "ReconstructionDataFormats/TrackFwd.h" + +namespace o2::itsmft::tracking::detail +{ + +/// MFT CA uses o2::mft::constants::mft::LayersNumber half-disk layers (same index as GeometryTGeo::getLayer). +/// Physical disk index is halfLayer / 2; ROFOverlapTable stores one LayerTiming per half-layer. + +inline float mftLayerZ(int layer) +{ + return o2::mft::constants::mft::LayerZCoordinate()[layer]; +} + +inline float mftLayerMSAngle(int layer, const TrackingParameters& params) +{ + const float invP = 1.f / params.TrackletMinPt; + const float zLayer = mftLayerZ(layer); + const float rRef = params.LayerRadii[layer]; + const float tanlRef = (std::abs(rRef) > 1e-6f) ? zLayer / rRef : 0.f; + const float absTanl = std::abs(tanlRef); + const float cscLambda = (absTanl > 1e-6f) ? std::sqrt(1.f + tanlRef * tanlRef) / absTanl : 1e6f; + return 0.0136f * invP * std::sqrt(params.LayerxX0[layer] * cscLambda); +} + +inline void mftTrackletProject(float xCl, float yCl, float zCl, float pvX, float pvY, float pvZ, + int fromLayer, int toLayer, float bz, float minPt, + float& xProj, float& yProj) +{ + if (std::abs(bz) > 0.01f && minPt > 0.f) { + const float zTo = mftLayerZ(toLayer); + const float dxTan = xCl - pvX; + const float dyTan = yCl - pvY; + const float dzTan = zCl - pvZ; + const float drTan = std::sqrt(dxTan * dxTan + dyTan * dyTan); + float invQPt = 1.f / minPt; + float tanl = (drTan > 1e-6f) ? -std::abs(dzTan) / drTan : -1.f; + float phi = (drTan > 1e-6f) ? std::atan2(dyTan, dxTan) : 0.f; + if (std::abs(tanl) > 1e-6f) { + const float k = std::abs(o2::constants::math::B2C * bz); + const float hz = (bz > 0.f) ? 1.f : -1.f; + phi -= 0.5f * hz * invQPt * dzTan * k / tanl; + } + ROOT::Math::SVector params{xCl, yCl, phi, tanl, invQPt}; + ROOT::Math::SMatrix> cov{}; + cov(0, 0) = cov(1, 1) = cov(2, 2) = cov(3, 3) = 1.; + const double qptSigma = std::clamp(static_cast(std::abs(invQPt)), 1., 10.); + cov(4, 4) = qptSigma * qptSigma; + o2::track::TrackParCovFwd track{zCl, params, cov, 0.}; + track.propagateToZhelix(zTo, bz); + xProj = static_cast(track.getX()); + yProj = static_cast(track.getY()); + } else { + const float dz0 = mftLayerZ(fromLayer) - pvZ; + if (std::abs(dz0) < 1e-6f) { + xProj = xCl; + yProj = yCl; + return; + } + const float w = (mftLayerZ(toLayer) - pvZ) / dz0; + xProj = pvX + w * (xCl - pvX); + yProj = pvY + w * (yCl - pvY); + } +} + +inline void mftTrackletSigmaXY(float x0, float y0, float pvX, float pvY, float pvZ, + float sigma2X0, float sigma2Y0, float sigma2PvX, float sigma2PvY, float sigma2PvZ, + int fromLayer, int toLayer, float rLayerFrom, float meanDeltaZ, float msAngle, + float bendingAngle, float xProj, float yProj, float& sigmaX, float& sigmaY) +{ + const float zFrom = mftLayerZ(fromLayer); + const float zTo = mftLayerZ(toLayer); + const float dz0 = zFrom - pvZ; + const float tanlRef = (std::abs(rLayerFrom) > 1e-6f) ? zFrom / rLayerFrom : 0.f; + const float sigma2MS = meanDeltaZ * meanDeltaZ * msAngle * msAngle * (tanlRef * tanlRef + 1.f); + if (std::abs(dz0) < o2::its::constants::Tolerance) { + sigmaX = std::sqrt(sigma2X0 + sigma2PvX + sigma2MS); + sigmaY = std::sqrt(sigma2Y0 + sigma2PvY + sigma2MS); + } else { + const float w = (zTo - pvZ) / dz0; + const float invDz0 = w / dz0; + const float sigma2W = invDz0 * invDz0 * sigma2PvZ; + const float dx0 = x0 - pvX; + const float dy0 = y0 - pvY; + const float oneMinusW = 1.f - w; + sigmaX = std::sqrt(oneMinusW * oneMinusW * sigma2PvX + w * w * sigma2X0 + dx0 * dx0 * sigma2W + sigma2MS); + sigmaY = std::sqrt(oneMinusW * oneMinusW * sigma2PvY + w * w * sigma2Y0 + dy0 * dy0 * sigma2W + sigma2MS); + } + const float rProj = std::hypot(xProj, yProj); + if (rProj > 1e-6f && bendingAngle > 0.f) { + const float dr = rProj * bendingAngle; + const float invR = 1.f / rProj; + const float sinPhi = yProj * invR; + const float cosPhi = xProj * invR; + sigmaX = std::sqrt(sigmaX * sigmaX + dr * dr * sinPhi * sinPhi); + sigmaY = std::sqrt(sigmaY * sigmaY + dr * dr * cosPhi * cosPhi); + } +} + +inline void mftFwdPropagateToZ(o2::track::TrackParCovFwd& track, float z, float bz) +{ + if (std::abs(bz) > 0.01f) { + track.propagateToZhelix(z, bz); + } else { + track.propagateToZlinear(z); + } +} + +inline float mftFwdPredictedChi2(const o2::track::TrackParCovFwd& track, float x, float y, float sigma2X, float sigma2Y) +{ + const float dx = x - static_cast(track.getX()); + const float dy = y - static_cast(track.getY()); + const float vx = static_cast(track.getSigma2X()) + sigma2X; + const float vy = static_cast(track.getSigma2Y()) + sigma2Y; + if (vx <= 0.f || vy <= 0.f) { + return o2::constants::math::VeryBig; + } + return dx * dx / vx + dy * dy / vy; +} + +inline float mftFwdStateChi2(const o2::track::TrackParCovFwd& current, const o2::track::TrackParCovFwd& rhs) +{ + ROOT::Math::SVector diff{ + rhs.getX() - current.getX(), + rhs.getY() - current.getY(), + rhs.getPhi() - current.getPhi(), + rhs.getTanl() - current.getTanl(), + rhs.getInvQPt() - current.getInvQPt()}; + auto cov = current.getCovariances(); + cov += rhs.getCovariances(); + if (!cov.Invert()) { + return o2::constants::math::VeryBig; + } + return static_cast(ROOT::Math::Similarity(cov, diff)); +} + +inline bool mftFwdAttachCluster(o2::track::TrackParCovFwd& track, float z, float x, float y, + float sigma2X, float sigma2Y, float xOverX0, float bz, float maxChi2, + float& chi2, bool checkChi2OnLast = false) +{ + mftFwdPropagateToZ(track, z, bz); + if (xOverX0 > 0.f) { + track.addMCSEffect(xOverX0); + } + const float predChi2 = mftFwdPredictedChi2(track, x, y, sigma2X, sigma2Y); + if (checkChi2OnLast && predChi2 > maxChi2) { + return false; + } + const std::array p{x, y}; + const std::array cov{sigma2X, sigma2Y}; + if (!track.update(p, cov)) { + return false; + } + chi2 += predChi2; + return true; +} + +/// Squared transverse distance from cluster c to the seed line c1→c2 (legacy MFT getDistanceToSeed). +inline float mftDistanceToSeedSquared(const o2::its::Cluster& c1, const o2::its::Cluster& c2, const o2::its::Cluster& c) +{ + const float dxSeed = c2.xCoordinate - c1.xCoordinate; + const float dySeed = c2.yCoordinate - c1.yCoordinate; + const float dzSeed = c2.zCoordinate - c1.zCoordinate; + if (std::abs(dzSeed) < 1e-9f) { + return std::numeric_limits::max(); + } + const float invdzSeed = (c.zCoordinate - c1.zCoordinate) / dzSeed; + const float xSeed = c1.xCoordinate + dxSeed * invdzSeed; + const float ySeed = c1.yCoordinate + dySeed * invdzSeed; + const float dx = c.xCoordinate - xSeed; + const float dy = c.yCoordinate - ySeed; + return dx * dx + dy * dy; +} + +/// Conical road scale (1 + dz/z_from)^2 between half-layers (legacy ROADclsRCut behaviour). +inline float mftConicalRoadR2Scale(int layerFrom, int layerTo) +{ + const float zFrom = mftLayerZ(layerFrom); + if (std::abs(zFrom) < 1e-6f) { + return 1.f; + } + const float dCone = 1.f + (mftLayerZ(layerTo) - zFrom) / zFrom; + return dCone * dCone; +} + +/// Cheap geometric pre-cut before forward cell fit (CellRoadRCut / ROADclsRCut). +inline bool validateMFTCellClusters(const o2::its::Cluster& c0, int layer0, + const o2::its::Cluster& c1, int layer1, + const o2::its::Cluster& c2, int layer2, + float r2Cut) +{ + return mftDistanceToSeedSquared(c0, c2, c1) < r2Cut * mftConicalRoadR2Scale(layer0, layer1) && + mftDistanceToSeedSquared(c0, c1, c2) < r2Cut * mftConicalRoadR2Scale(layer0, layer2) && + mftDistanceToSeedSquared(c1, c2, c0) < r2Cut * mftConicalRoadR2Scale(layer1, layer0); +} + +/// Build inward forward seed at the outer cluster and Kalman-fit the three cell clusters. +template +inline bool mftFwdFitCellClusters(const int hitLayers[3], const int clusIds[3], + const TimeFrame& tf, const TrackingParameters& params, + float bz, o2::track::TrackParCovFwd& track, float& chi2) +{ + const auto& cInner = tf.getUnsortedClusters()[hitLayers[0]][clusIds[0]]; + const auto& cMid = tf.getUnsortedClusters()[hitLayers[1]][clusIds[1]]; + const auto& cOuter = tf.getUnsortedClusters()[hitLayers[2]][clusIds[2]]; + if (cInner.zCoordinate <= cOuter.zCoordinate + 1.e-6f) { + return false; + } + + const float dxTan = cMid.xCoordinate - cInner.xCoordinate; + const float dyTan = cMid.yCoordinate - cInner.yCoordinate; + const float dzTan = cMid.zCoordinate - cInner.zCoordinate; + const float drTan = std::sqrt(dxTan * dxTan + dyTan * dyTan); + const float dxPhi = cOuter.xCoordinate - cInner.xCoordinate; + const float dyPhi = cOuter.yCoordinate - cInner.yCoordinate; + const float dzPhi = cOuter.zCoordinate - cInner.zCoordinate; + const float drPhi = std::sqrt(dxPhi * dxPhi + dyPhi * dyPhi); + if (drTan < 1.e-6f || std::abs(dzTan) < 1.e-6f || drPhi < 1.e-6f || std::abs(dzPhi) < 1.e-6f) { + return false; + } + + const float minPt = params.TrackletMinPt; + const float invQPt = (minPt > 0.f) ? 1.f / minPt : 0.f; + float tanl{0.f}; + float phi{0.f}; + if (std::abs(bz) > 0.01f) { + tanl = -std::abs(dzTan) / drTan; + phi = std::atan2(dyPhi, dxPhi); + if (std::abs(tanl) > 1.e-6f) { + const float k = std::abs(o2::constants::math::B2C * bz); + const float hz = (bz > 0.f) ? 1.f : -1.f; + phi -= 0.5f * hz * invQPt * dzPhi * k / tanl; + } + } else { + tanl = -std::abs(dzPhi) / drPhi; + phi = std::atan2(dyPhi, dxPhi); + } + + const auto& tfOuter = tf.getTrackingFrameInfoOnLayer(hitLayers[2])[clusIds[2]]; + ROOT::Math::SVector seedParams{cOuter.xCoordinate, cOuter.yCoordinate, phi, tanl, invQPt}; + ROOT::Math::SMatrix> seedCov{}; + seedCov(0, 0) = tfOuter.covarianceTrackingFrame[0] > 0.f ? tfOuter.covarianceTrackingFrame[0] : 1.f; + seedCov(1, 1) = tfOuter.covarianceTrackingFrame[2] > 0.f ? tfOuter.covarianceTrackingFrame[2] : 1.f; + seedCov(2, 2) = seedCov(3, 3) = 1.; + const double qptSigma = std::clamp(static_cast(std::abs(invQPt)), 1., 10.); + seedCov(4, 4) = qptSigma * qptSigma; + track = {cOuter.zCoordinate, seedParams, seedCov, 0.}; + + chi2 = 0.f; + for (int iC{2}; iC >= 0; --iC) { + const int layer = hitLayers[iC]; + const int clIdx = clusIds[iC]; + const auto& tfInfo = tf.getTrackingFrameInfoOnLayer(layer)[clIdx]; + if (!mftFwdAttachCluster(track, tfInfo.zCoordinate, tfInfo.xCoordinate, tfInfo.yCoordinate, + tfInfo.covarianceTrackingFrame[0], tfInfo.covarianceTrackingFrame[2], + params.LayerxX0[layer], bz, params.MaxChi2ClusterAttachment, chi2, iC == 0)) { + return false; + } + } + return true; +} + +template +inline bool mftFwdAttachNeighbourToSeed(const SeedT& seed, int neighbourLayer, int neighbourClIdx, + const TimeFrame& tf, const TrackingParameters& params, + float bz, o2::track::TrackParCovFwd& fwdTrack, float& chi2) +{ + fwdTrack = static_cast(seed); + chi2 = seed.getChi2(); + const auto& trHit = tf.getTrackingFrameInfoOnLayer(neighbourLayer)[neighbourClIdx]; + return mftFwdAttachCluster(fwdTrack, trHit.zCoordinate, trHit.xCoordinate, trHit.yCoordinate, + trHit.covarianceTrackingFrame[0], trHit.covarianceTrackingFrame[2], + params.LayerxX0[neighbourLayer], bz, params.MaxChi2ClusterAttachment, chi2, true); +} + +inline bool mftFwdCellsAreCompatible(const CellSeedTpl& currentCell, + const CellSeedTpl& nextCell, + float bz, float maxChi2) +{ + if (currentCell.getSecondClusterIndex() != nextCell.getFirstClusterIndex()) { + return false; + } + o2::track::TrackParCovFwd nextFwd = static_cast(nextCell); + const auto& currentFwd = static_cast(currentCell); + mftFwdPropagateToZ(nextFwd, static_cast(currentFwd.getZ()), bz); + return mftFwdStateChi2(currentFwd, nextFwd) <= maxChi2; +} + +} // namespace o2::itsmft::tracking::detail + +namespace o2::itsmft::tracking +{ + +bool refitTrackFwd(const TrackSeedN& seed, + MFTCATrack& track, + const TimeFrame& tf, + const TrackingParameters& params, + float bz); + +} // namespace o2::itsmft::tracking + +#endif /* ALICEO2_ITSMFT_TRACKING_MFTFWDTRACKHELPERS_H_ */ diff --git a/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/TimeFrame.h b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/TimeFrame.h new file mode 100644 index 0000000000000..dd976bde96185 --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/TimeFrame.h @@ -0,0 +1,648 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// + +#ifndef ALICEO2_ITSMFT_TRACKING_TIMEFRAME_H_ +#define ALICEO2_ITSMFT_TRACKING_TIMEFRAME_H_ + +#include +#include +#include +#include +#include +#include + +#include "DataFormatsITS/TrackITS.h" +#include "DataFormatsITS/Vertex.h" + +#include "ITStracking/BoundedAllocator.h" +#include "ITSMFTTracking/Cell.h" +#include "ITStracking/Cluster.h" +#include "ITStracking/ClusterLines.h" +#include "ITStracking/ExternalAllocator.h" +#include "ITStracking/ROFLookupTables.h" +#include "ITStracking/Tracklet.h" + +#include "ITSMFTTracking/MFTCATrack.h" +#include "ITSMFTTracking/Configuration.h" +#include "ITSMFTTracking/IndexTableUtils.h" +#include "ITSMFTTracking/LayerMask.h" +#include "ITSMFTTracking/TrackingTopology.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTruthContainer.h" + +#include "DetectorsBase/Propagator.h" +#include "DetectorsCommonDataFormats/DetID.h" + +namespace o2 +{ +namespace gpu +{ +class GPUChainITS; +} + +namespace itsmft +{ +class CompClusterExt; +class TopologyDictionary; +class ROFRecord; +} // namespace itsmft + +namespace itsmft::tracking +{ + +// Re-use ITS CA tracking data structures; only TimeFrame and index-table I/O are detector-aware. +using Cluster = o2::its::Cluster; +using TrackingFrameInfo = o2::its::TrackingFrameInfo; +using Tracklet = o2::its::Tracklet; +using LayerMask = o2::itsmft::tracking::LayerMask; +using BoundedMemoryResource = o2::its::BoundedMemoryResource; +template +using bounded_vector = o2::its::bounded_vector; +using ExternalAllocator = o2::its::ExternalAllocator; +using Line = o2::its::Line; +using ClusterLines = o2::its::ClusterLines; +using Vertex = o2::its::Vertex; +using VertexLabel = o2::its::VertexLabel; +using TrackITSExt = o2::its::TrackITSExt; + +namespace constants +{ +using namespace o2::its::constants; +using o2::itsmft::tracking::ITSNLayers; +using o2::itsmft::tracking::nLayersForDet; +} // namespace constants + +template +struct TimeFrame { + using IndexTableUtilsN = o2::itsmft::IndexTableUtils; + using ROFOverlapTableN = o2::its::ROFOverlapTable; + using ROFVertexLookupTableN = o2::its::ROFVertexLookupTable; + using ROFMaskTableN = o2::its::ROFMaskTable; + using TrackingTopologyN = o2::itsmft::tracking::TrackingTopology; + using CellSeedN = o2::itsmft::tracking::CellSeedN; + using TrackSeedN = o2::itsmft::tracking::TrackSeedN; + + TimeFrame() = default; + virtual ~TimeFrame() = default; + + const Vertex& getPrimaryVertex(const int ivtx) const { return mPrimaryVertices[ivtx]; } + auto& getPrimaryVertices() { return mPrimaryVertices; }; + auto getPrimaryVerticesNum() { return mPrimaryVertices.size(); }; + const auto& getPrimaryVertices() const { return mPrimaryVertices; }; + auto& getPrimaryVerticesLabels() { return mPrimaryVerticesLabels; }; + gsl::span getPrimaryVertices(int layer, int rofId) const; + void addPrimaryVertex(const Vertex& vertex); + void addPrimaryVertexLabel(const VertexLabel& label) { mPrimaryVerticesLabels.push_back(label); } + + // read-in data + void loadROFrameData(gsl::span rofs, + gsl::span clusters, + gsl::span::iterator& pattIt, + const itsmft::TopologyDictionary* dict, + int layer, + const dataformats::MCTruthContainer* mcLabels = nullptr, + o2::detectors::DetID::ID detId = o2::detectors::DetID::ITS); + o2::detectors::DetID::ID getDetId() const { return mDetId; } + void setDetId(o2::detectors::DetID::ID detId) { mDetId = detId; } + void resetROFrameData(int iLayer); + void prepareROFrameData(gsl::span clusters, int layer); + + int getTotalClusters() const; + bool empty() const { return getTotalClusters() == 0; } + int getSortedIndex(int rofId, int layer, int idx) const { return mROFramesClusters[layer][rofId] + idx; } + int getSortedStartIndex(const int rofId, const int layer) const { return mROFramesClusters[layer][rofId]; } + int getNrof(int layer) const { return mROFramesClusters[layer].size() - 1; } + + void resetBeamXY(const float x, const float y, const float w = 0); + void setBeamPosition(const float x, const float y, const float s2, const float base = 50.f, const float systematic = 0.f) + { + isBeamPositionOverridden = true; + resetBeamXY(x, y, s2 / o2::gpu::CAMath::Sqrt((base * base) + systematic)); + } + + float getBeamX() const { return mBeamPos[0]; } + float getBeamY() const { return mBeamPos[1]; } + std::array& getBeamXY() { return mBeamPos; } + + auto& getMinRs() { return mMinR; } + auto& getMaxRs() { return mMaxR; } + float getMinR(int layer) const { return mMinR[layer]; } + float getMaxR(int layer) const { return mMaxR[layer]; } + float getTransitionPhiCut(int transitionId) const { return mTransitionPhiCuts[transitionId]; } + float getTransitionMSAngle(int transitionId) const { return mTransitionMSAngles[transitionId]; } + auto& getTransitionPhiCuts() { return mTransitionPhiCuts; } + auto& getTransitionMSAngles() { return mTransitionMSAngles; } + float getPositionResolution(int layer) const { return mPositionResolution[layer]; } + auto& getPositionResolutions() { return mPositionResolution; } + + gsl::span getClustersOnLayer(int rofId, int layerId); + gsl::span getClustersOnLayer(int rofId, int layerId) const; + gsl::span getClustersPerROFrange(int rofMin, int range, int layerId) const; + gsl::span getUnsortedClustersOnLayer(int rofId, int layerId) const; + gsl::span getUsedClustersROF(int rofId, int layerId); + gsl::span getUsedClustersROF(int rofId, int layerId) const; + gsl::span getROFramesClustersPerROFrange(int rofMin, int range, int layerId) const; + gsl::span getROFrameClusters(int layerId) const; + gsl::span getNClustersROFrange(int rofMin, int range, int layerId) const; + gsl::span getIndexTable(int rofId, int layerId); + const auto& getTrackingFrameInfoOnLayer(int layerId) const { return mTrackingFrameInfo[layerId]; } + + // navigation tables + const auto& getIndexTableUtils() const { return mIndexTableUtils; } + const auto& getROFOverlapTable() const { return mROFOverlapTable; } + const auto& getROFOverlapTableView() const { return mROFOverlapTableView; } + const auto& getTrackerTopologies() const { return mTrackerTopologies; } + const auto& getTrackingTopologyView() const { return mTrackingTopologyView; } + void setROFOverlapTable(ROFOverlapTableN table) + { + mROFOverlapTable = std::move(table); + mROFOverlapTableView = mROFOverlapTable.getView(); + } + const auto& getROFVertexLookupTable() const { return mROFVertexLookupTable; } + const auto& getROFVertexLookupTableView() const { return mROFVertexLookupTableView; } + void setROFVertexLookupTable(ROFVertexLookupTableN table) + { + mROFVertexLookupTable = std::move(table); + mROFVertexLookupTableView = mROFVertexLookupTable.getView(); + } + void updateROFVertexLookupTable() { mROFVertexLookupTable.update(mPrimaryVertices.data(), mPrimaryVertices.size()); } + void setMultiplicityCutMask(ROFMaskTableN cutMask) + { + mMultiplicityCutMask = std::move(cutMask); + mROFMaskView = mROFMask->getView(); + } + void useMultiplictyMask() noexcept + { + mROFMask = &mMultiplicityCutMask; + mROFMaskView = mROFMask->getView(); + } + void setUPCCutMask(ROFMaskTableN cutMask) { mUPCCutMask = std::move(cutMask); } + void useUPCMask() noexcept + { + mROFMask = &mUPCCutMask; + mROFMaskView = mROFMask->getView(); + } + const auto& getROFMaskView() const { return mROFMaskView; } + + const TrackingFrameInfo& getClusterTrackingFrameInfo(int layerId, const Cluster& cl) const; + gsl::span getClusterLabels(int layerId, const Cluster& cl) const { return getClusterLabels(layerId, cl.clusterId); } + gsl::span getClusterLabels(int layerId, const int clId) const { return mClusterLabels[((mIsStaggered) ? layerId : 0)]->getLabels(mClusterExternalIndices[layerId][clId]); } + int getClusterExternalIndex(int layerId, const int clId) const { return mClusterExternalIndices[layerId][clId]; } + int getClusterSize(int layer, int clusterId) const { return mClusterSize[layer][clusterId]; } + void setClusterSize(int layer, bounded_vector& v) { mClusterSize[layer] = std::move(v); } + + auto& getTrackletsLabel(int layer) { return mTrackletLabels[layer]; } + auto& getCellsLabel(int layer) { return mCellLabels[layer]; } + + bool hasMCinformation() const { return mClusterLabels[0] != nullptr; } + void initVertexingTopology(const TrackingParameters& trkParam); + void initDefaultTrackingTopology(const TrackingParameters& trkParam, const int maxLayers = NLayers); + void initTrackerTopologies(gsl::span trkParams, const int maxLayers = NLayers); + void initialise(const TrackingParameters& trkParam, const int maxLayers = NLayers, const int iteration = constants::UnusedIndex); + + bool isClusterUsed(int layer, int clusterId) const { return mUsedClusters[layer][clusterId]; } + void markUsedCluster(int layer, int clusterId) { mUsedClusters[layer][clusterId] = true; } + gsl::span getUsedClusters(const int layer); + + auto& getTracklets() { return mTracklets; } + auto& getTrackletsLookupTable() { return mTrackletsLookupTable; } + + auto& getClusters() { return mClusters; } + auto& getUnsortedClusters() { return mUnsortedClusters; } + const auto& getUnsortedClusters() const { return mUnsortedClusters; } + int getClusterROF(int iLayer, int iCluster); + auto& getCells() { return mCells; } + + auto& getCellsLookupTable() { return mCellsLookupTable; } + auto& getCellsNeighbours() { return mCellsNeighbours; } + auto& getCellsNeighboursTopology() { return mCellsNeighboursTopology; } + auto& getCellsNeighboursLUT() { return mCellsNeighboursLUT; } + auto& getTracks() { return mTracks; } + const auto& getTracks() const { return mTracks; } + auto& getTracksLabel() { return mTracksLabel; } + const auto& getTracksLabel() const { return mTracksLabel; } + auto& getLinesLabel(const int rofId) { return mLinesLabels[rofId]; } + + size_t getNumberOfClusters() const; + virtual size_t getNumberOfCells() const; + virtual size_t getNumberOfTracklets() const; + virtual size_t getNumberOfNeighbours() const; + size_t getNumberOfTracks() const; + size_t getNumberOfUsedClusters() const; + + /// memory management + void setMemoryPool(std::shared_ptr pool); + auto& getMemoryPool() const noexcept { return mMemoryPool; } + bool checkMemory(unsigned long max) { return getArtefactsMemory() < max; } + unsigned long getArtefactsMemory() const; + void printArtefactsMemory() const; + + /// staggering + void setIsStaggered(bool b) noexcept { mIsStaggered = b; } + + // Vertexer + void computeTrackletsPerROFScans(); + void computeTracletsPerClusterScans(); + int& getNTrackletsROF(int rofId, int combId) { return mNTrackletsPerROF[combId][rofId]; } + auto& getLines(int rofId) { return mLines[rofId]; } + int getNLinesTotal() const noexcept { return mTotalLines; } + void setNLinesTotal(uint32_t a) noexcept { mTotalLines = a; } + auto& getTrackletClusters(int rofId) { return mTrackletClusters[rofId]; } + gsl::span getFoundTracklets(int rofId, int combId) const; + gsl::span getFoundTracklets(int rofId, int combId); + gsl::span getLabelsFoundTracklets(int rofId, int combId) const; + gsl::span getNTrackletsCluster(int rofId, int combId); + gsl::span getExclusiveNTrackletsCluster(int rofId, int combId); + uint32_t getTotalTrackletsTF(const int iLayer) { return mTotalTracklets[iLayer]; } + int getTotalClustersPerROFrange(int rofMin, int range, int layerId) const; + // \Vertexer + + int hasBogusClusters() const { return std::accumulate(mBogusClusters.begin(), mBogusClusters.end(), 0); } + + void setBz(float bz) { mBz = bz; } + float getBz() const { return mBz; } + + /// State if memory will be externally managed by the GPU framework + ExternalAllocator* mExternalAllocator{nullptr}; + std::shared_ptr mExtMemoryPool; // host memory pool managed by the framework + auto getFrameworkAllocator() { return mExternalAllocator; }; + void setFrameworkAllocator(ExternalAllocator* ext); + bool hasFrameworkAllocator() const noexcept { return mExternalAllocator != nullptr; } + std::pmr::memory_resource* getMaybeFrameworkHostResource(bool forceHost = false) { return (hasFrameworkAllocator() && !forceHost) ? mExtMemoryPool.get() : mMemoryPool.get(); } + + // Propagator + const o2::base::PropagatorImpl* getDevicePropagator() const { return mPropagatorDevice; } + virtual void setDevicePropagator(const o2::base::PropagatorImpl* /*unused*/) {}; + + template + void addClusterToLayer(int layer, T&&... args); + template + void addTrackingFrameInfoToLayer(int layer, T&&... args); + void addTrackingFrameInfoToLayer(int layer, const TrackingFrameInfo& tfInfo) { mTrackingFrameInfo[layer].push_back(tfInfo); } + void addClusterExternalIndexToLayer(int layer, const int idx) { mClusterExternalIndices[layer].push_back(idx); } + + std::array, NLayers> mClusters; + std::array, NLayers> mTrackingFrameInfo; + std::array, NLayers> mClusterExternalIndices; + std::array, NLayers> mROFramesClusters; + std::array*, NLayers> mClusterLabels{nullptr}; + std::array, 2> mNTrackletsPerCluster; + std::array, 2> mNTrackletsPerClusterSum; + std::array, NLayers> mNClustersPerROF; + std::array, NLayers> mIndexTables; + std::vector> mTrackletsLookupTable; + std::array, NLayers> mUsedClusters; + + std::array, NLayers> mUnsortedClusters; + std::vector> mTracklets; + std::vector> mCells; + bounded_vector> mTracks; + bounded_vector mTracksLabel; + std::vector> mCellsNeighbours; + std::vector> mCellsNeighboursTopology; + std::vector> mCellsLookupTable; + + const o2::base::PropagatorImpl* mPropagatorDevice = nullptr; // Needed only for GPU + o2::detectors::DetID::ID mDetId{o2::detectors::DetID::ITS}; + + virtual void wipe(); + + // interface + virtual bool isGPU() const noexcept { return false; } + virtual const char* getName() const noexcept { return "CPU"; } + + protected: + void prepareClusters(const TrackingParameters& trkParam, const int maxLayers = NLayers); + float mBz = 5.; + unsigned int mNTotalLowPtVertices = 0; + int mBeamPosWeight = 0; + std::array mBeamPos = {0.f, 0.f}; + bool isBeamPositionOverridden = false; + std::array mMinR; + std::array mMaxR; + bounded_vector mTransitionPhiCuts; + bounded_vector mTransitionMSAngles; + bounded_vector mPositionResolution; + std::array, NLayers> mClusterSize; + + bounded_vector> mPValphaX; /// PV x and alpha for track propagation + std::vector> mTrackletLabels; + std::vector> mCellLabels; + std::vector> mCellsNeighboursLUT; + bounded_vector mBogusClusters; /// keep track of clusters with wild coordinates + + // Vertexer + bounded_vector mPrimaryVertices; + bounded_vector mPrimaryVerticesLabels; + std::vector> mNTrackletsPerROF; + std::vector> mLines; + std::vector> mTrackletClusters; + std::array, 2> mTrackletsIndexROF; + std::vector> mLinesLabels; + std::array mTotalTracklets = {0, 0}; + uint32_t mTotalLines = 0; + // \Vertexer + + // lookup tables + IndexTableUtilsN mIndexTableUtils; + ROFOverlapTableN mROFOverlapTable; + ROFOverlapTableN::View mROFOverlapTableView; + TrackingTopologyN mVertexingTopology; + TrackingTopologyN mDefaultTrackingTopology; + std::vector mTrackerTopologies; + typename TrackingTopologyN::View mTrackingTopologyView; + ROFVertexLookupTableN mROFVertexLookupTable; + ROFVertexLookupTableN::View mROFVertexLookupTableView; + ROFMaskTableN mMultiplicityCutMask; + ROFMaskTableN mUPCCutMask; + ROFMaskTableN* mROFMask = &mMultiplicityCutMask; + ROFMaskTableN::View mROFMaskView; + + bool mIsStaggered{false}; + + std::shared_ptr mMemoryPool; +}; + +template +gsl::span TimeFrame::getPrimaryVertices(int layer, int rofId) const +{ + if (rofId < 0 || rofId >= getNrof(layer)) { + return {}; + } + const auto& entry = mROFVertexLookupTableView.getVertices(layer, rofId); + return {&mPrimaryVertices[entry.getFirstEntry()], static_cast::size_type>(entry.getEntries())}; +} + +template +inline void TimeFrame::resetBeamXY(const float x, const float y, const float w) +{ + mBeamPos[0] = x; + mBeamPos[1] = y; + mBeamPosWeight = w; +} + +template +inline gsl::span TimeFrame::getROFrameClusters(int layerId) const +{ + return {&mROFramesClusters[layerId][0], static_cast::size_type>(mROFramesClusters[layerId].size())}; +} + +template +inline gsl::span TimeFrame::getClustersOnLayer(int rofId, int layerId) +{ + if (rofId < 0 || rofId >= getNrof(layerId)) { + return {}; + } + int startIdx{mROFramesClusters[layerId][rofId]}; + return {&mClusters[layerId][startIdx], static_cast::size_type>(mROFramesClusters[layerId][rofId + 1] - startIdx)}; +} + +template +inline gsl::span TimeFrame::getClustersOnLayer(int rofId, int layerId) const +{ + if (rofId < 0 || rofId >= getNrof(layerId)) { + return {}; + } + int startIdx{mROFramesClusters[layerId][rofId]}; + return {&mClusters[layerId][startIdx], static_cast::size_type>(mROFramesClusters[layerId][rofId + 1] - startIdx)}; +} + +template +inline gsl::span TimeFrame::getUsedClustersROF(int rofId, int layerId) +{ + if (rofId < 0 || rofId >= getNrof(layerId)) { + return {}; + } + int startIdx{mROFramesClusters[layerId][rofId]}; + return {&mUsedClusters[layerId][startIdx], static_cast::size_type>(mROFramesClusters[layerId][rofId + 1] - startIdx)}; +} + +template +inline gsl::span TimeFrame::getUsedClustersROF(int rofId, int layerId) const +{ + if (rofId < 0 || rofId >= getNrof(layerId)) { + return {}; + } + int startIdx{mROFramesClusters[layerId][rofId]}; + return {&mUsedClusters[layerId][startIdx], static_cast::size_type>(mROFramesClusters[layerId][rofId + 1] - startIdx)}; +} + +template +inline gsl::span TimeFrame::getClustersPerROFrange(int rofMin, int range, int layerId) const +{ + if (rofMin < 0 || rofMin >= getNrof(layerId)) { + return {}; + } + int startIdx{mROFramesClusters[layerId][rofMin]}; // First cluster of rofMin + int endIdx{mROFramesClusters[layerId][o2::gpu::CAMath::Min(rofMin + range, getNrof(layerId))]}; + return {&mClusters[layerId][startIdx], static_cast::size_type>(endIdx - startIdx)}; +} + +template +inline gsl::span TimeFrame::getROFramesClustersPerROFrange(int rofMin, int range, int layerId) const +{ + int chkdRange{o2::gpu::CAMath::Min(range, getNrof(layerId) - rofMin)}; + return {&mROFramesClusters[layerId][rofMin], static_cast::size_type>(chkdRange)}; +} + +template +inline gsl::span TimeFrame::getNClustersROFrange(int rofMin, int range, int layerId) const +{ + int chkdRange{o2::gpu::CAMath::Min(range, getNrof(layerId) - rofMin)}; + return {&mNClustersPerROF[layerId][rofMin], static_cast::size_type>(chkdRange)}; +} + +template +inline int TimeFrame::getTotalClustersPerROFrange(int rofMin, int range, int layerId) const +{ + int startIdx{rofMin}; // First cluster of rofMin + int endIdx{o2::gpu::CAMath::Min(rofMin + range, getNrof(layerId))}; + return mROFramesClusters[layerId][endIdx] - mROFramesClusters[layerId][startIdx]; +} + +template +inline int TimeFrame::getClusterROF(int iLayer, int iCluster) +{ + return std::lower_bound(mROFramesClusters[iLayer].begin(), mROFramesClusters[iLayer].end(), iCluster + 1) - mROFramesClusters[iLayer].begin() - 1; +} + +template +inline gsl::span TimeFrame::getUnsortedClustersOnLayer(int rofId, int layerId) const +{ + if (rofId < 0 || rofId >= getNrof(layerId)) { + return {}; + } + int startIdx{mROFramesClusters[layerId][rofId]}; + return {&mUnsortedClusters[layerId][startIdx], static_cast::size_type>(mROFramesClusters[layerId][rofId + 1] - startIdx)}; +} + +template +inline gsl::span TimeFrame::getIndexTable(int rofId, int layer) +{ + if (rofId < 0 || rofId >= getNrof(layer)) { + return {}; + } + const int tableSize = mIndexTableUtils.getNrowBins() * mIndexTableUtils.getNcolBins() + 1; + return {&mIndexTables[layer][rofId * tableSize], static_cast::size_type>(tableSize)}; +} + +template +template +void TimeFrame::addClusterToLayer(int layer, T&&... values) +{ + mUnsortedClusters[layer].emplace_back(std::forward(values)...); +} + +template +template +void TimeFrame::addTrackingFrameInfoToLayer(int layer, T&&... values) +{ + mTrackingFrameInfo[layer].emplace_back(std::forward(values)...); +} + +template +inline gsl::span TimeFrame::getUsedClusters(const int layer) +{ + return {&mUsedClusters[layer][0], static_cast::size_type>(mUsedClusters[layer].size())}; +} + +template +inline gsl::span TimeFrame::getNTrackletsCluster(int rofId, int combId) +{ + if (rofId < 0 || rofId >= getNrof(1)) { + return {}; + } + auto startIdx{mROFramesClusters[1][rofId]}; + return {&mNTrackletsPerCluster[combId][startIdx], static_cast::size_type>(mROFramesClusters[1][rofId + 1] - startIdx)}; +} + +template +inline gsl::span TimeFrame::getExclusiveNTrackletsCluster(int rofId, int combId) +{ + if (rofId < 0 || rofId >= getNrof(1)) { + return {}; + } + auto clusStartIdx{mROFramesClusters[1][rofId]}; + + return {&mNTrackletsPerClusterSum[combId][clusStartIdx], static_cast::size_type>(mROFramesClusters[1][rofId + 1] - clusStartIdx)}; +} + +template +inline gsl::span TimeFrame::getFoundTracklets(int rofId, int combId) +{ + if (rofId < 0 || rofId >= getNrof(1) || mTracklets[combId].empty()) { + return {}; + } + auto startIdx{mNTrackletsPerROF[combId][rofId]}; + return {&mTracklets[combId][startIdx], static_cast::size_type>(mNTrackletsPerROF[combId][rofId + 1] - startIdx)}; +} + +template +inline gsl::span TimeFrame::getFoundTracklets(int rofId, int combId) const +{ + if (rofId < 0 || rofId >= getNrof(1)) { + return {}; + } + auto startIdx{mNTrackletsPerROF[combId][rofId]}; + return {&mTracklets[combId][startIdx], static_cast::size_type>(mNTrackletsPerROF[combId][rofId + 1] - startIdx)}; +} + +template +inline gsl::span TimeFrame::getLabelsFoundTracklets(int rofId, int combId) const +{ + if (rofId < 0 || rofId >= getNrof(1) || !hasMCinformation()) { + return {}; + } + auto startIdx{mNTrackletsPerROF[combId][rofId]}; + return {&mTrackletLabels[combId][startIdx], static_cast::size_type>(mNTrackletsPerROF[combId][rofId + 1] - startIdx)}; +} + +template +inline int TimeFrame::getTotalClusters() const +{ + size_t totalClusters{0}; + for (const auto& clusters : mUnsortedClusters) { + totalClusters += clusters.size(); + } + return int(totalClusters); +} + +template +inline size_t TimeFrame::getNumberOfClusters() const +{ + size_t nClusters{0}; + for (const auto& layer : mClusters) { + nClusters += layer.size(); + } + return nClusters; +} + +template +inline size_t TimeFrame::getNumberOfCells() const +{ + size_t nCells{0}; + for (const auto& layer : mCells) { + nCells += layer.size(); + } + return nCells; +} + +template +inline size_t TimeFrame::getNumberOfTracklets() const +{ + size_t nTracklets{0}; + for (const auto& layer : mTracklets) { + nTracklets += layer.size(); + } + return nTracklets; +} + +template +inline size_t TimeFrame::getNumberOfNeighbours() const +{ + size_t neigh{0}; + for (const auto& l : mCellsNeighbours) { + neigh += l.size(); + } + return neigh; +} + +template +inline size_t TimeFrame::getNumberOfTracks() const +{ + return mTracks.size(); +} + +template +inline size_t TimeFrame::getNumberOfUsedClusters() const +{ + size_t nClusters = 0; + for (const auto& layer : mUsedClusters) { + nClusters += std::count(layer.begin(), layer.end(), true); + } + return nClusters; +} + +template +inline const TrackingFrameInfo& TimeFrame::getClusterTrackingFrameInfo(int layerId, const Cluster& cl) const +{ + return mTrackingFrameInfo[layerId][cl.clusterId]; +} + +using TimeFrameITS = TimeFrame; +/// MFT CA TimeFrame: NLayers = half-disk CA layers (see MFTFwdTrackHelpers.h). +using TimeFrameMFT = TimeFrame; + +} // namespace itsmft::tracking +} // namespace o2 + +#endif /* ALICEO2_ITSMFT_TRACKING_TIMEFRAME_H_ */ diff --git a/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/Tracker.h b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/Tracker.h new file mode 100644 index 0000000000000..ff1b3035050e5 --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/Tracker.h @@ -0,0 +1,21 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file Tracker.h +/// \brief ITS-style alias for the shared CA tracker orchestrator +/// + +#ifndef ALICEO2_ITSMFT_TRACKING_TRACKER_H_ +#define ALICEO2_ITSMFT_TRACKING_TRACKER_H_ + +#include "ITSMFTTracking/CATracker.h" + +#endif /* ALICEO2_ITSMFT_TRACKING_TRACKER_H_ */ diff --git a/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/TrackerTraits.h b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/TrackerTraits.h new file mode 100644 index 0000000000000..5cbe6553f9a5e --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/TrackerTraits.h @@ -0,0 +1,87 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file TrackerTraits.h +/// \brief Shared CA tracker traits: same ITS-style tracklet/cell/road logic; MFT uses x-y LUT and forward refit +/// + +#ifndef ALICEO2_ITSMFT_TRACKING_TRACKERTRAITS_H_ +#define ALICEO2_ITSMFT_TRACKING_TRACKERTRAITS_H_ + +#include + +#include "ITSMFTTracking/Configuration.h" +#include "ITSMFTTracking/TimeFrame.h" +#include "ITStracking/BoundedAllocator.h" + +namespace o2::itsmft::tracking +{ + +template +class TrackerTraits +{ + public: + using TimeFrameN = TimeFrame; + using IndexTableUtilsN = o2::itsmft::IndexTableUtils; + using CellSeedN = typename TimeFrameN::CellSeedN; + using TrackSeedN = typename TimeFrameN::TrackSeedN; + + virtual ~TrackerTraits() = default; + virtual void adoptTimeFrame(TimeFrameN* tf) { mTimeFrame = tf; } + virtual void initialiseTimeFrame(const int iteration) + { + mTimeFrame->initialise(mTrkParams[iteration], mTrkParams[iteration].NLayers, iteration); + } + + virtual void computeLayerTracklets(const int iteration, int iVertex); + virtual void computeLayerCells(const int iteration); + virtual void findCellsNeighbours(const int iteration); + virtual void findRoads(const int iteration); + + template + void processNeighbours(int iteration, int defaultCellTopologyId, int iLevel, const bounded_vector& currentCellSeed, const bounded_vector& currentCellId, const bounded_vector& currentCellTopologyId, bounded_vector& updatedCellSeed, bounded_vector& updatedCellId, bounded_vector& updatedCellTopologyId); + + void acceptTracks(int iteration, bounded_vector>& tracks, bounded_vector>& firstClusters); + void markTracks(int iteration); + + void updateTrackingParameters(const std::vector& trkPars) { mTrkParams = trkPars; } + TimeFrameN* getTimeFrame() { return mTimeFrame; } + + virtual void setBz(float bz); + float getBz() const { return mBz; } + virtual const char* getName() const noexcept { return "CPU"; } + virtual bool isGPU() const noexcept { return false; } + void setMemoryPool(std::shared_ptr pool) noexcept { mMemoryPool = pool; } + auto getMemoryPool() const noexcept { return mMemoryPool; } + + void setNThreads(int n, std::shared_ptr& arena); + int getNThreads() { return mTaskArena->max_concurrency(); } + + int getTFNumberOfClusters() const { return mTimeFrame->getNumberOfClusters(); } + int getTFNumberOfTracklets() const { return mTimeFrame->getNumberOfTracklets(); } + int getTFNumberOfCells() const { return mTimeFrame->getNumberOfCells(); } + + private: + std::shared_ptr mMemoryPool; + std::shared_ptr mTaskArena; + + protected: + TimeFrameN* mTimeFrame = nullptr; + std::vector mTrkParams; + float mBz{-999.f}; +}; + +using TrackerTraitsITS = TrackerTraits; +using TrackerTraitsMFT = TrackerTraits; + +} // namespace o2::itsmft::tracking + +#endif /* ALICEO2_ITSMFT_TRACKING_TRACKERTRAITS_H_ */ diff --git a/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/TrackingConfigParam.h b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/TrackingConfigParam.h new file mode 100644 index 0000000000000..8793b15c2285f --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/TrackingConfigParam.h @@ -0,0 +1,168 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef ALICEO2_ITSMFT_TRACKING_CONFIG_PARAM_H_ +#define ALICEO2_ITSMFT_TRACKING_CONFIG_PARAM_H_ + +#include +#include + +#include "CommonUtils/ConfigurableParam.h" +#include "CommonUtils/ConfigurableParamHelper.h" +#include "DetectorsCommonDataFormats/DetID.h" + +namespace o2::itsmft::tracking +{ +/// ITS CA layer count (matches legacy o2::its::TrackingParameters::NLayers default). +constexpr int ITSNLayers = 7; +/// MFT CA half-disk layer count (same as o2::mft::constants::mft::LayersNumber). +constexpr int MFTNLayers = 10; +/// Max CA iterations (same as o2::its::constants::MaxIter). +constexpr int MaxIter = 4; +} // namespace o2::itsmft::tracking + +namespace o2::itsmft +{ + +/// ITS vertexer settings (header only for now; not registered or used while ITS +/// tracking stays on O2::ITStracking) +struct VertexerParamConfig : public o2::conf::ConfigurableParamHelper { + bool saveTimeBenchmarks = false; // dump metrics on file + + int nIterations = 1; // Number of vertexing passes to perform. + int vertPerRofThreshold = 0; // Maximum number of vertices per ROF to trigger second a iteration. + + // geometrical cuts for tracklet selection for Pb-Pb + float zCut = 0.002f; + float phiCut = 0.005f; + float pairCut = 0.017321f; + float clusterCut = 0.170048f; + float coarseZWindow = 0.055458f; + float seedDedupZCut = 0.116685f; + float refitDedupZCut = 0.039855f; + float duplicateZCut = 0.200097f; + float finalSelectionZCut = 0.034535f; + float duplicateDistance2Cut = 0.005117f; + float tanLambdaCut = 0.002f; // tanLambda = deltaZ/deltaR + float nSigmaCut = 0.0164651f; + float maxZPositionAllowed = 25.f; // 4x sZ of the beam + + // Artefacts selections + int clusterContributorsCut = 3; // minimum number of contributors for an accepted final vertex + int suppressLowMultDebris = 16; // suppress all vertices below this threshold if a vertex was already found in a rof + int seedMemberRadiusTime = 0; + int seedMemberRadiusZ = 2; + int maxTrackletsPerCluster = 100; + int phiSpan = -1; + int zSpan = -1; + int ZBins = 1; // z-phi index table configutation: number of z bins + int PhiBins = 128; // z-phi index table configutation: number of phi bins + + bool useTruthSeeding{false}; // overwrite seeding vertices with MC truth + + int nThreads = 1; + bool printMemory = false; + size_t maxMemory = std::numeric_limits::max(); + bool dropTFUponFailure = false; + + O2ParamDef(VertexerParamConfig, "ITSVertexerParam"); +}; + +template +struct TrackerParamConfig : public o2::conf::ConfigurableParamHelper> { + static constexpr std::string_view getParamName() + { + return N == o2::detectors::DetID::ITS ? TrackerParamName[0] : TrackerParamName[1]; + } + + static constexpr int MinTrackLength = N == o2::detectors::DetID::ITS ? 4 : 5; + static constexpr int MaxTrackLength = N == o2::detectors::DetID::ITS ? o2::itsmft::tracking::ITSNLayers + : o2::itsmft::tracking::MFTNLayers; + + static constexpr int getNLayers() + { + return N == o2::detectors::DetID::ITS ? o2::itsmft::tracking::ITSNLayers : o2::itsmft::tracking::MFTNLayers; + } + + bool useMatCorrTGeo = false; // use full geometry to corect for material budget accounting in the fits. Default is to use the material budget LUT. + bool useFastMaterial = false; // use faster material approximation for material budget accounting in the fits. + int addTimeError[getNLayers()] = {0}; // configure the width of the window in BC to be considered for the tracking. + int minTrackLgtIter[o2::itsmft::tracking::MaxIter] = {}; // minimum track length at each iteration, used only if >0, otherwise use code defaults + uint8_t startLayerMask[o2::itsmft::tracking::MaxIter] = {}; // mask of start layer for this iteration (if >0) + int maxHolesIter[o2::itsmft::tracking::MaxIter] = {}; // maximum number of missing internal layers allowed in the CA topology for each iteration + uint16_t holeLayerMaskIter[o2::itsmft::tracking::MaxIter] = {}; // layers that may be skipped by the CA topology for each iteration + float minPtIterLgt[o2::itsmft::tracking::MaxIter * (MaxTrackLength - MinTrackLength + 1)] = {}; // min.pT for given track length at this iteration, used only if >0, otherwise use code defaults + float sysErr2Row[getNLayers()] = {0}; // systematic error^2 along ALPIDE rows (local X) per layer + float sysErr2Col[getNLayers()] = {0}; // systematic error^2 along ALPIDE columns (local Z) per layer + float maxChi2ClusterAttachment = -1.f; + float maxChi2NDF = -1.f; + float nSigmaCut = -1.f; + float deltaTanLres = -1.f; + float minPt = -1.f; + float pvRes = -1.f; + int LUTbinsU = N == o2::detectors::DetID::MFT ? 64 : -1; // number of LUT bins along the first coordinate (ITS: phi, MFT: global x) + int LUTbinsV = N == o2::detectors::DetID::MFT ? 128 : -1; // number of LUT bins along the second coordinate (ITS: z, MFT: global y) + float diamondPos[3] = {0.f, 0.f, 0.f}; // diamond vertex position (MFT CA tracklet seed; ITS when useDiamond) + bool useDiamond = N == o2::detectors::DetID::MFT; // MFT CA: always diamond; ITS: opt-in via param + bool perPrimaryVertexProcessing = false; // perform the full tracking considering the vertex hypotheses one at the time. + bool saveTimeBenchmarks = false; // dump metrics on file + bool overrideBeamEstimation = false; // ITS only: meanVertex CCDB beam seed (MFT CA always uses diamond) + int trackingMode = -1; // -1: unset (use --tracking-mode), 0=sync, 1=async, 2=cosmics + bool doUPCIteration = false; // Perform an additional iteration for UPC events on tagged vertices. You want to combine this config with VertexerParamConfig.nIterations=2 + int nIterations = N == o2::detectors::DetID::MFT ? 1 : o2::itsmft::tracking::MaxIter; // overwrite the number of iterations + int reseedIfShorter = 6; // for the final refit reseed the track with circle if they are shorter than this value + bool shiftRefToCluster{true}; // TrackFit: after update shift the linearization reference to cluster + bool repeatRefitOut{false}; // repeat outward refit using inward refit as a seed + bool createArtefactLabels{false}; // create on-the-fly labels for the artefacts + + int nThreads = 1; + bool printMemory = false; + size_t maxMemory = std::numeric_limits::max(); + bool dropTFUponFailure = false; + bool fataliseUponFailure = true; // granular management of the fatalisation in async mode + + // Selections on tracks sharing clusters + bool allowSharingFirstCluster = false; // allow first cluster sharing among tracks + float sharedClusterMaxDeltaPhi = 0.05f; // Maximum allowed delta phi at the cluster position + float sharedClusterMaxDeltaEta = 0.03f; // Maximum allowed delta eta at the cluster position + bool sharedClusterOppositeSign = false; // Require opposite sign of the tracklets + float cellRoadRCut = -1.f; // MFT: max distance to seed line; <=0 uses default (0.05 cm) + float trackletMinAbsX = -1.f; // MFT: min |x| (cm) for tracklet seeds and accepted tracks; <=0 uses code default + + O2ParamDef(TrackerParamConfig, getParamName().data()); + + private: + static constexpr std::string_view TrackerParamName[2] = {"ITSCATrackerParam", "MFTCATrackerParam"}; + + static_assert(N == o2::detectors::DetID::ITS || N == o2::detectors::DetID::MFT, "only DetID::ITS or DetID::MFT are allowed"); +}; + +template +TrackerParamConfig TrackerParamConfig::sInstance; + +} // namespace o2::itsmft + +namespace framework +{ +template +struct is_messageable; +template <> +struct is_messageable : std::true_type { +}; +template <> +struct is_messageable> : std::true_type { +}; +template <> +struct is_messageable> : std::true_type { +}; +} // namespace framework + +#endif /* ALICEO2_ITSMFT_TRACKING_CONFIG_PARAM_H_ */ diff --git a/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/TrackingInterface.h b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/TrackingInterface.h new file mode 100644 index 0000000000000..790cb5b991ccc --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/TrackingInterface.h @@ -0,0 +1,123 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file TrackingInterface.h +/// \brief Shared DPL-facing CA tracking interface for ITS and MFT +/// + +#ifndef ALICEO2_ITSMFT_TRACKING_TRACKINGINTERFACE_H_ +#define ALICEO2_ITSMFT_TRACKING_TRACKINGINTERFACE_H_ + +#include +#include +#include + +#include + +#include "CommonDataFormat/IRFrame.h" +#include "DataFormatsCalibration/MeanVertexObject.h" +#include "DataFormatsITSMFT/CompCluster.h" +#include "DataFormatsITSMFT/ROFRecord.h" +#include "DataFormatsITSMFT/TopologyDictionary.h" +#include "DetectorsCommonDataFormats/DetID.h" +#include "ITSMFTTracking/DetectorTraits.h" +#include "ITSMFTTracking/Tracker.h" +#include "ITSMFTTracking/Configuration.h" +#include "ITSMFTTracking/TimeFrame.h" +#include "ITSMFTTracking/TrackerTraits.h" +#include "ITStracking/BoundedAllocator.h" +#include "ITStracking/ROFLookupTables.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTruthContainer.h" + +namespace o2::itsmft::tracking +{ + +template +class ITSMFTTrackingInterface +{ + public: + static_assert(NLayers == ITSNLayers || NLayers == o2::mft::constants::mft::LayersNumber, + "ITSMFTTrackingInterface supports ITS (7) and MFT (10) layer counts only"); + static constexpr o2::detectors::DetID::ID DetId = detIdFromNLayers(); + + using TimeFrameN = TimeFrame; + using TrackerN = Tracker; + using TrackerTraitsN = TrackerTraits; + using ROFOverlapTableN = o2::its::ROFOverlapTable; + using ROFVertexLookupTableN = o2::its::ROFVertexLookupTable; + using ROFMaskTableN = o2::its::ROFMaskTable; + using BoundedMemoryResourceN = BoundedMemoryResource; + + ITSMFTTrackingInterface(bool useMC, o2::itsmft::TrackingMode::Type mode, bool overrideBeamEst); + + void setTrackingMode(o2::itsmft::TrackingMode::Type mode) { mTrackingMode = mode; } + void setClusterDictionary(const o2::itsmft::TopologyDictionary* dict) { mDict = dict; } + void setMeanVertex(const o2::dataformats::MeanVertexObject* v) { mMeanVertex = v; } + + void initialise(); + + /// Phase 2 (TimeFrame load) + Phase 3 (CA tracking). Returns tracker elapsed ms or -1 on failure. + /// Does not wipe the TimeFrame; call clearTimeFrame() after outputs are extracted. + float processTimeFrame(gsl::span rofs, + gsl::span clusters, + gsl::span patterns, + const o2::dataformats::MCTruthContainer* labels, + gsl::span irFrames = {}); + + void clearTimeFrame() { mTimeFrame.wipe(); } + + TimeFrameN& getTimeFrame() { return mTimeFrame; } + const TimeFrameN& getTimeFrame() const { return mTimeFrame; } + const std::vector& getTrackingParameters() const { return mTrackParams; } + bool isActive() const { return !mTrackParams.empty(); } + + protected: + virtual void onTimeFrameLoaded() {} + virtual void onTrackingFinished(float elapsedMs) {} + + private: + void resolveTrackingParameters(); + void initialiseMemoryPool(); + void initialiseTracker(); + void loadTimeFrame(gsl::span rofs, + gsl::span clusters, + gsl::span patterns, + const o2::dataformats::MCTruthContainer* labels, + gsl::span irFrames); + float runTracking(); + void configureROFLookupTables(); + void configureBeamPosition(); + void configureTrackingTopology(); + void configureROFMask(gsl::span rofs, + gsl::span irFrames); + void validateROFInput(gsl::span rofs) const; + + bool mUseMC = false; + bool mOverrideBeamEstimation = false; + o2::itsmft::TrackingMode::Type mTrackingMode = o2::itsmft::TrackingMode::Unset; + std::vector mTrackParams; + std::shared_ptr mMemoryPool; + std::unique_ptr mTrackerTraits; + std::unique_ptr mTracker; + const o2::itsmft::TopologyDictionary* mDict = nullptr; + const o2::dataformats::MeanVertexObject* mMeanVertex = nullptr; + TimeFrameN mTimeFrame; + int mMFTROFrameLengthInBC = 0; + bool mMFTTriggered = false; +}; + +using ITSMFTTrackingInterfaceITS = ITSMFTTrackingInterface; +using ITSMFTTrackingInterfaceMFT = ITSMFTTrackingInterface; + +} // namespace o2::itsmft::tracking + +#endif /* ALICEO2_ITSMFT_TRACKING_TRACKINGINTERFACE_H_ */ diff --git a/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/TrackingTopology.h b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/TrackingTopology.h new file mode 100644 index 0000000000000..822369508af31 --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/TrackingTopology.h @@ -0,0 +1,215 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef ALICEO2_ITSMFT_TRACKING_INCLUDE_TRACKINGTOPOLOGY_H_ +#define ALICEO2_ITSMFT_TRACKING_INCLUDE_TRACKINGTOPOLOGY_H_ + +#include +#include +#include +#include + +#ifndef GPUCA_GPUCODE +#include +#include +#include "Framework/Logger.h" +#endif + +#include "CommonDataFormat/RangeReference.h" +#include "GPUCommonDef.h" +#include "GPUCommonMath.h" +#include "ITSMFTTracking/LayerMask.h" + +namespace o2::itsmft::tracking +{ + +template +class TrackingTopology +{ + public: + using Id = uint8_t; + using Mask = LayerMask; + using Range = o2::dataformats::RangeReference; + static constexpr int MaxTransitions = (NLayers * (NLayers - 1)) / 2; + static constexpr int MaxCells = (NLayers * (NLayers - 1) * (NLayers - 2)) / 6; + static_assert(NLayers < std::numeric_limits::max()); + static_assert(MaxTransitions <= std::numeric_limits::max()); + static_assert(MaxCells <= std::numeric_limits::max()); + + struct LayerTransition { + Id fromLayer{0}; + Id toLayer{0}; + }; + static_assert(std::is_standard_layout_v); + static_assert(std::is_trivially_copyable_v); + static_assert(sizeof(LayerTransition) == (2 * sizeof(Id))); + + struct CellTopology { + Id firstTransition{0}; + Id secondTransition{0}; + Mask hitLayerMask{0}; + }; + static_assert(std::is_standard_layout_v); + static_assert(std::is_trivially_copyable_v); + static_assert(sizeof(CellTopology) == (2 * sizeof(Id)) + sizeof(Mask)); + + struct View { + const LayerTransition* transitions{nullptr}; + const CellTopology* cells{nullptr}; + const Range* cellsByFirstTransitionIndex{nullptr}; + const Id* cellsByFirstTransition{nullptr}; + Id nTransitions{0}; + Id nCells{0}; + Id nCellsByFirstTransition{0}; + + GPUhdi() const LayerTransition& getTransition(Id id) const { return transitions[id]; } + GPUhdi() const CellTopology& getCell(Id id) const { return cells[id]; } + GPUhdi() Range getCellsStartingWithTransition(Id transitionId) const { return cellsByFirstTransitionIndex[transitionId]; } + +#ifndef GPUCA_GPUCODE + std::string asString() const + { + std::string out = fmt::format("TrackingTopology: transitions={} cells={}", nTransitions, nCells); + out += "\n transitions:"; + for (Id transitionId = 0; transitionId < nTransitions; ++transitionId) { + const auto& t = transitions[transitionId]; + out += fmt::format("\n {}: {} -> {}", transitionId, t.fromLayer, t.toLayer); + } + out += "\n cells:"; + for (Id cellId = 0; cellId < nCells; ++cellId) { + const auto& c = cells[cellId]; + const auto& first = transitions[c.firstTransition]; + const auto& second = transitions[c.secondTransition]; + out += fmt::format("\n {}: {} -> {} -> {} hitMask={} transitions=({}, {})", cellId, first.fromLayer, first.toLayer, second.toLayer, c.hitLayerMask.asString(), c.firstTransition, c.secondTransition); + } + return out; + } + + void print() const + { + LOGP(info, "{}", asString()); + } +#endif + }; + + void init(int maxLayers, int maxHoles, Mask holeLayerMask) + { + clear(); + mMaxLayers = o2::gpu::CAMath::Max(0, o2::gpu::CAMath::Min(maxLayers, NLayers)); + mMaxHoles = o2::gpu::CAMath::Max(maxHoles, 0); + mHoleLayerMask = holeLayerMask; + for (int fromLayer = 0; fromLayer < mMaxLayers; ++fromLayer) { + for (int toLayer = fromLayer + 1; toLayer < mMaxLayers; ++toLayer) { + if (Mask::skipped(fromLayer, toLayer).isAllowedHoleMask(mMaxHoles, mHoleLayerMask)) { + mTransitions[mNTransitions++] = LayerTransition{static_cast(fromLayer), static_cast(toLayer)}; + } + } + } + + for (Id firstId = 0; firstId < mNTransitions; ++firstId) { + const auto& first = mTransitions[firstId]; + for (Id secondId = 0; secondId < mNTransitions; ++secondId) { + const auto& second = mTransitions[secondId]; + if (first.toLayer != second.fromLayer) { + continue; + } + const Mask hitMask{first.fromLayer, first.toLayer, second.toLayer}; + if (hitMask.isAllowed(mMaxHoles, mHoleLayerMask)) { + mCells[mNCells++] = CellTopology{firstId, secondId, hitMask}; + } + } + } + + fillCellsByTransition(); + } + + View getView() const + { + return View{mTransitions.data(), + mCells.data(), + mCellsByFirstTransitionIndex.data(), + mCellsByFirstTransition.data(), + mNTransitions, + mNCells, + mNCellsByFirstTransition}; + } + + View getDeviceView(const LayerTransition* deviceTransitions, + const CellTopology* deviceCells, + const Range* deviceCellsByFirstTransitionIndex, + const Id* deviceCellsByFirstTransition) const + { + return View{deviceTransitions, + deviceCells, + deviceCellsByFirstTransitionIndex, + deviceCellsByFirstTransition, + mNTransitions, + mNCells, + mNCellsByFirstTransition}; + } + + const auto& getTransitions() const noexcept { return mTransitions; } + const auto& getCells() const noexcept { return mCells; } + const auto& getCellsByFirstTransitionIndex() const noexcept { return mCellsByFirstTransitionIndex; } + const auto& getCellsByFirstTransition() const noexcept { return mCellsByFirstTransition; } + Id getNTransitions() const noexcept { return mNTransitions; } + Id getNCells() const noexcept { return mNCells; } + Id getNCellsByFirstTransition() const noexcept { return mNCellsByFirstTransition; } + + private: + void clear() + { + mNTransitions = 0; + mNCells = 0; + mNCellsByFirstTransition = 0; + mTransitions.fill({}); + mCells.fill({}); + mCellsByFirstTransitionIndex.fill(Range{0, 0}); + mCellsByFirstTransition.fill(0); + } + + void fillCellsByTransition() + { + std::array counts{}; + for (Id cellId = 0; cellId < mNCells; ++cellId) { + ++counts[mCells[cellId].firstTransition]; + } + + Id offset = 0; + for (Id transitionId = 0; transitionId < mNTransitions; ++transitionId) { + mCellsByFirstTransitionIndex[transitionId].setFirstEntry(offset); + mCellsByFirstTransitionIndex[transitionId].setEntries(counts[transitionId]); + offset += counts[transitionId]; + } + + std::array cursor{}; + for (Id cellId = 0; cellId < mNCells; ++cellId) { + const Id transitionId = mCells[cellId].firstTransition; + mCellsByFirstTransition[mCellsByFirstTransitionIndex[transitionId].getFirstEntry() + cursor[transitionId]++] = cellId; + } + mNCellsByFirstTransition = offset; + } + + int mMaxLayers{0}; + int mMaxHoles{0}; + Mask mHoleLayerMask{0}; + Id mNTransitions{0}; + Id mNCells{0}; + Id mNCellsByFirstTransition{0}; + std::array mTransitions{}; + std::array mCells{}; + std::array mCellsByFirstTransitionIndex{}; + std::array mCellsByFirstTransition{}; +}; + +} // namespace o2::itsmft::tracking + +#endif /* ALICEO2_ITSMFT_TRACKING_INCLUDE_TRACKINGTOPOLOGY_H_ */ diff --git a/Detectors/ITSMFT/common/tracking/src/CATracker.cxx b/Detectors/ITSMFT/common/tracking/src/CATracker.cxx new file mode 100644 index 0000000000000..2a967dddd74b3 --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/src/CATracker.cxx @@ -0,0 +1,198 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file CATracker.cxx +/// \brief +/// + +#include "ITSMFTTracking/CATracker.h" + +#include +#include +#include + +#include "Framework/Logger.h" +#include "ITStracking/Constants.h" +#include "ITStracking/BoundedAllocator.h" +#include "MFTTracking/MFTTrackingParam.h" +#include "SimulationDataFormat/MCCompLabel.h" + +namespace o2::itsmft::tracking +{ + +namespace +{ +template +void computeTracksMClabels(TimeFrame& tf) +{ + auto& trackLabels = tf.getTracksLabel(); + trackLabels.clear(); + trackLabels.reserve(tf.getNumberOfTracks()); + + for (auto& track : tf.getTracks()) { + std::vector> occurrences; + for (int iLayer = 0; iLayer < NLayers; ++iLayer) { + const int index = track.getClusterIndex(iLayer); + if (index == constants::UnusedIndex) { + continue; + } + const auto labels = tf.getClusterLabels(iLayer, index); + bool found{false}; + for (auto& occurrence : occurrences) { + for (const auto& label : labels) { + if (label == occurrence.first) { + ++occurrence.second; + found = true; + } + } + } + if (!found) { + for (const auto& label : labels) { + occurrences.emplace_back(label, 1); + } + } + } + + MCCompLabel maxOccurrencesValue; + if (!occurrences.empty()) { + std::sort(occurrences.begin(), occurrences.end(), [](const auto& e1, const auto& e2) { + return e1.second > e2.second; + }); + maxOccurrencesValue = occurrences[0].first; + if constexpr (NLayers == o2::mft::constants::mft::LayersNumber) { + const float threshold = o2::mft::MFTTrackingParam::Instance().TrueTrackMCThreshold; + if (static_cast(occurrences[0].second) / track.getNumberOfClusters() < threshold) { + maxOccurrencesValue.setFakeFlag(); + } + } else if (occurrences[0].second < static_cast(track.getNumberOfClusters())) { + maxOccurrencesValue.setFakeFlag(); + } + } else { + maxOccurrencesValue.setFakeFlag(); + } + trackLabels.emplace_back(maxOccurrencesValue); + } +} +} // namespace + +template +Tracker::Tracker(TrackerTraits* traits) : mTraits(traits) +{ +} + +template +void Tracker::adoptTimeFrame(TimeFrame& tf) +{ + mTimeFrame = &tf; + mTraits->adoptTimeFrame(&tf); +} + +template +float Tracker::clustersToTracks() +{ + mTraits->updateTrackingParameters(mTrkParams); + + int maxNvertices{-1}; + if (mTrkParams[0].PerPrimaryVertexProcessing) { + maxNvertices = mTimeFrame->getROFVertexLookupTableView().getMaxVerticesPerROF(); + } + + float total{0.f}; + try { + for (int iteration = 0; iteration < static_cast(mTrkParams.size()); ++iteration) { + mMemoryPool->setMaxMemory(mTrkParams[iteration].MaxMemory); + if (mTrkParams[iteration].PassFlags[IterationStep::UseUPCMask]) { + mTimeFrame->useUPCMask(); + } + + int iVertex = std::min(maxNvertices, 0); + initialiseTimeFrame(iteration); + do { + computeTracklets(iteration, iVertex); + computeCells(iteration); + findCellsNeighbours(iteration); + findRoads(iteration); + } while (++iVertex < maxNvertices); + } + } catch (const BoundedMemoryResource::MemoryLimitExceeded& err) { + LOGP(error, "CA tracker exceeded memory limit: {}", err.what()); + if (mTrkParams[0].DropTFUponFailure) { + mTimeFrame->wipe(); + return -1.f; + } + throw; + } catch (const std::exception& err) { + LOGP(error, "CA tracker failed: {}", err.what()); + mTimeFrame->getTracks().clear(); + return -1.f; + } + + if (mTimeFrame->hasMCinformation()) { + computeTracksMClabels(*mTimeFrame); + } + rectifyClusterIndices(); + sortTracks(); + return total; +} + +template +void Tracker::rectifyClusterIndices() +{ + for (auto& track : mTimeFrame->getTracks()) { + for (int iCluster = 0; iCluster < CATrackType::MaxClusters; ++iCluster) { + const int index = track.getClusterIndex(iCluster); + if (index == constants::UnusedIndex) { + continue; + } + track.setExternalClusterIndex(iCluster, mTimeFrame->getClusterExternalIndex(iCluster, index)); + } + } +} + +template +void Tracker::sortTracks() +{ + auto& tracks = mTimeFrame->getTracks(); + bounded_vector indices(tracks.size(), mMemoryPool.get()); + std::iota(indices.begin(), indices.end(), 0); + std::sort(indices.begin(), indices.end(), [&tracks](size_t i, size_t j) { + const auto& a = tracks[i]; + const auto& b = tracks[j]; + const auto aLower = a.getTimeStamp().getTimeStamp() - a.getTimeStamp().getTimeStampError(); + const auto bLower = b.getTimeStamp().getTimeStamp() - b.getTimeStamp().getTimeStampError(); + if (aLower != bLower) { + return aLower < bLower; + } + return a.getChi2() < b.getChi2(); + }); + + bounded_vector> sortedTracks(mMemoryPool.get()); + sortedTracks.reserve(tracks.size()); + for (size_t idx : indices) { + sortedTracks.push_back(tracks[idx]); + } + tracks.swap(sortedTracks); + + if (mTimeFrame->hasMCinformation()) { + auto& trackLabels = mTimeFrame->getTracksLabel(); + bounded_vector sortedLabels(mMemoryPool.get()); + sortedLabels.reserve(trackLabels.size()); + for (size_t idx : indices) { + sortedLabels.push_back(trackLabels[idx]); + } + trackLabels.swap(sortedLabels); + } +} + +template class Tracker<7>; +template class Tracker<10>; + +} // namespace o2::itsmft::tracking diff --git a/Detectors/ITSMFT/common/tracking/src/Configuration.cxx b/Detectors/ITSMFT/common/tracking/src/Configuration.cxx new file mode 100644 index 0000000000000..bd178dd52f0d1 --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/src/Configuration.cxx @@ -0,0 +1,341 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "DetectorsBase/Propagator.h" +#include "Framework/Logger.h" +#include "ITSMFTTracking/Configuration.h" +#include "ITSMFTTracking/TrackingConfigParam.h" +#include "ITStracking/Constants.h" +#include "MFTTracking/Constants.h" + +namespace +{ +constexpr bool iequals(std::string_view a, std::string_view b) +{ + return std::equal(a.begin(), a.end(), b.begin(), b.end(), + [](char x, char y) { return std::tolower(x) == std::tolower(y); }); +} +} // namespace + +namespace o2::itsmft +{ + +std::string TrackingParameters::asString() const +{ + std::string str = std::format("NColB:{} NRowB:{} PerVtx:{} DropFail:{} TtklMinPt:{:.2f} MinCl:{}", ColBins, RowBins, PerPrimaryVertexProcessing, DropTFUponFailure, TrackletMinPt, MinTrackLength); + auto isSet = [](auto e) { return e >= 0; }; + auto isAnySet = [&isSet](auto v) { return !v.empty() && std::any_of(v.begin(), v.end(), isSet); }; + bool first = true; + for (int il = NLayers; il >= MinTrackLength; il--) { + int slot = NLayers - il; + if (slot < (int)MinPt.size() && MinPt[slot] > 0) { + if (first) { + first = false; + str += " MinPt: "; + } + str += std::format("L{}:{:.2f} ", il, MinPt[slot]); + } + } + if (isAnySet(SystError2Row) || isAnySet(SystError2Col)) { + str += " SystErrRow/Col:"; + for (size_t i = 0; i < SystError2Row.size(); i++) { + str += std::format("{:.2e}/{:.2e} ", SystError2Row[i], SystError2Col[i]); + } + } + if (isAnySet(AddTimeError)) { + str += " AddTimeError:"; + for (unsigned int i : AddTimeError) { + str += std::format("{} ", i); + } + } + if (SharedMaxClusters) { + str += std::format(" ShaMaxCls:{} ", SharedMaxClusters); + } + if (AllowSharingFirstCluster) { + str += std::format(" ShaClsDPhi:{} ShaClsDEta:{} ShaClsSign:{}", SharedClusterMaxDeltaPhi, SharedClusterMaxDeltaEta, SharedClusterOppositeSign); + } + if (MaxHoles) { + str += std::format(" MaxHoles:{} HoleMask:{}", MaxHoles, HoleLayerMask.asString()); + } + if (std::numeric_limits::max() != MaxMemory) { + str += std::format(" MemLimit {:.2f} GB", double(MaxMemory) / (1024.f * 1024.f * 1024.f)); + } + return str; +} + +std::string VertexingParameters::asString() const +{ + std::string str = std::format("NColB:{} NRowB:{} MinVtxCont:{} SupLowMultDebris:{} MaxTrkltCls:{} ZCut:{} PhCut:{} PairCut:{} ClCut:{} SeedRad:{}x{}", + ColBins, RowBins, clusterContributorsCut, suppressLowMultDebris, maxTrackletsPerCluster, zCut, phiCut, pairCut, clusterCut, seedMemberRadiusTime, seedMemberRadiusZ); + if (std::numeric_limits::max() != MaxMemory) { + str += std::format(" MemLimit {:.2f} GB", double(MaxMemory) / (1024.f * 1024.f * 1024.f)); + } + return str; +} + +void resetDetectorDefaults(TrackingParameters& p, detectors::DetID::ID detId) +{ + if (detId == detectors::DetID::ITS) { + p = TrackingParameters{}; + p.MinPt.assign(TrackerParamConfig::MaxTrackLength - TrackerParamConfig::MinTrackLength + 1, 0.f); + return; + } + + if (detId == detectors::DetID::MFT) { + namespace mftc = o2::mft::constants; + namespace mft = mftc::mft; + constexpr int nLayers = o2::mft::constants::mft::LayersNumber; + + p = TrackingParameters{}; + p.NLayers = nLayers; + p.LayerZ.clear(); + p.LayerZ.reserve(nLayers); + for (float z : mft::LayerZCoordinate()) { + p.LayerZ.push_back(std::abs(z)); + } + p.LayerColHalfExtent.assign(mftc::index_table::RMax.begin(), mftc::index_table::RMax.end()); + p.IndexRowMin = -20.f; + p.IndexRowMax = 20.f; + p.LayerRadii.resize(nLayers); + for (int i{0}; i < nLayers; ++i) { + p.LayerRadii[i] = 0.5f * (mftc::index_table::RMin[i] + mftc::index_table::RMax[i]); + } + p.LayerxX0.assign(nLayers, 0.042f / mft::DisksNumber); // MFTRadLength / disks + p.LayerResolution.assign(nLayers, mft::Resolution); + p.SystError2Row.assign(nLayers, 0.f); + p.SystError2Col.assign(nLayers, 0.f); + p.AddTimeError.assign(nLayers, 0u); + p.ColBins = 64; + p.RowBins = 128; + p.UseDiamond = true; + p.PerPrimaryVertexProcessing = false; + p.StartLayerMask = (1u << nLayers) - 1u; + p.TrackletMinAbsX = 0.05f; + p.MinPt.assign(TrackerParamConfig::MaxTrackLength - TrackerParamConfig::MinTrackLength + 1, 0.f); + return; + } + + LOGP(fatal, "Unsupported detector id {} in resetDetectorDefaults", static_cast(detId)); +} + +namespace TrackingMode +{ + +Type fromString(std::string_view str) +{ + constexpr std::array smodes = { + std::pair{"sync", Sync}, + std::pair{"async", Async}, + std::pair{"cosmics", Cosmics}, + std::pair{"unset", Unset}, + std::pair{"off", Off}}; + + const auto it = std::find_if(smodes.begin(), smodes.end(), [&str](const auto& pair) { + return iequals(str, pair.first); + }); + if (it == smodes.end()) { + LOGP(fatal, "Unrecognized CA tracking mode '{}'", str); + } + return it->second; +} + +std::string toString(Type mode) +{ + switch (mode) { + case Sync: + return "sync"; + case Async: + return "async"; + case Cosmics: + return "cosmics"; + case Unset: + return "unset"; + case Off: + return "off"; + } + LOGP(fatal, "Unrecognized CA tracking mode {}", static_cast(mode)); + return ""; +} + +std::vector getTrackingParameters(detectors::DetID::ID detId, Type mode) +{ + if (detId == detectors::DetID::ITS) { + LOGP(fatal, "ITS CA tracking via O2::ITSMFTTracking is not enabled yet; use O2::ITStracking"); + } + if (detId != detectors::DetID::MFT) { + LOGP(fatal, "Unsupported detector id {} in getTrackingParameters", static_cast(detId)); + } + + const auto& tc = TrackerParamConfig::Instance(); + std::vector trackParams; + + if (mode == Off) { + return trackParams; + } + if (mode == Unset) { + LOGP(fatal, "CA tracking mode is unset; set --tracking-mode or {}.trackingMode", TrackerParamConfig::getParamName()); + } + + if (mode == Async) { + trackParams.resize(3); + for (auto& param : trackParams) { + resetDetectorDefaults(param, detId); + } + + trackParams[1].TrackletMinPt = 0.15f; + trackParams[1].CellDeltaTanLambdaSigma *= 2.f; + trackParams[2].TrackletMinPt = 0.08f; + trackParams[2].CellDeltaTanLambdaSigma *= 4.f; + + trackParams[0].MinPt[0] = 1.f / 12.f; // 10 clusters + trackParams[1].MinPt[0] = 1.f / 12.f; + + trackParams[2].MinTrackLength = TrackerParamConfig::MinTrackLength; + trackParams[2].MinPt[0] = 1.f / 12.f; // 10 cl + trackParams[2].MinPt[1] = 1.f / 8.f; // 9 cl + trackParams[2].MinPt[2] = 1.f / 5.f; // 8 cl + trackParams[2].MinPt[3] = 1.f / 3.f; // 7 cl + trackParams[2].MinPt[4] = 1.f / 2.f; // 6 cl + trackParams[2].MinPt[5] = 1.f / 1.f; // 5 cl + + for (int ip = 0; ip < static_cast(trackParams.size()); ip++) { + auto& param = trackParams[ip]; + if (ip < o2::its::constants::MaxIter) { + if (tc.startLayerMask[ip] > 0) { + param.StartLayerMask = tc.startLayerMask[ip]; + } + if (tc.minTrackLgtIter[ip] > 0) { + param.MinTrackLength = tc.minTrackLgtIter[ip]; + } + for (int ilg = tc.MaxTrackLength; ilg >= tc.MinTrackLength; ilg--) { + const int lslot0 = tc.MaxTrackLength - ilg; + const int lslot = lslot0 + ip * (tc.MaxTrackLength - tc.MinTrackLength + 1); + if (tc.minPtIterLgt[lslot] > 0.f) { + param.MinPt[lslot0] = tc.minPtIterLgt[lslot]; + } + } + } + } + } else if (mode == Sync) { + trackParams.resize(1); + resetDetectorDefaults(trackParams[0], detId); + trackParams[0].MinTrackLength = TrackerParamConfig::MinTrackLength; + } else if (mode == Cosmics) { + trackParams.resize(1); + resetDetectorDefaults(trackParams[0], detId); + trackParams[0].MinTrackLength = TrackerParamConfig::MinTrackLength; + trackParams[0].CellDeltaTanLambdaSigma *= 10.f; + trackParams[0].ColBins = 32; + trackParams[0].RowBins = 64; + trackParams[0].PVres = 1.e5f; + trackParams[0].MaxChi2ClusterAttachment = 60.f; + trackParams[0].MaxChi2NDF = 40.f; + } else { + LOGP(fatal, "Unsupported CA tracking mode {}", toString(mode)); + } + + for (auto& param : trackParams) { + param.PassFlags.reset(); + } + if (!trackParams.empty()) { + trackParams[0].PassFlags.set(IterationStep::FirstPass, IterationStep::RebuildClusterLUT); + } + + const float bFactor = std::abs(o2::base::Propagator::Instance()->getNominalBz()) / 5.0066791f; + const float bFactorTracklets = bFactor < 0.01f ? 1.f : bFactor; + + for (auto& p : trackParams) { + p.TrackletMinPt *= bFactorTracklets; + for (int ilg = tc.MaxTrackLength; ilg >= tc.MinTrackLength; ilg--) { + const int lslot = tc.MaxTrackLength - ilg; + if (lslot < static_cast(p.MinPt.size())) { + p.MinPt[lslot] *= bFactor; + } + } + + p.ReseedIfShorter = tc.reseedIfShorter; + p.RepeatRefitOut = tc.repeatRefitOut; + p.ShiftRefToCluster = tc.shiftRefToCluster; + p.CreateArtefactLabels = tc.createArtefactLabels; + p.PrintMemory = tc.printMemory; + p.MaxMemory = tc.maxMemory; + p.DropTFUponFailure = tc.dropTFUponFailure; + p.SaveTimeBenchmarks = tc.saveTimeBenchmarks; + p.FataliseUponFailure = tc.fataliseUponFailure; + p.AllowSharingFirstCluster = tc.allowSharingFirstCluster; + p.SharedClusterMaxDeltaPhi = tc.sharedClusterMaxDeltaPhi; + p.SharedClusterMaxDeltaEta = tc.sharedClusterMaxDeltaEta; + p.SharedClusterOppositeSign = tc.sharedClusterOppositeSign; + p.PerPrimaryVertexProcessing = tc.perPrimaryVertexProcessing; + + const auto iter = &p - trackParams.data(); + if (iter < o2::its::constants::MaxIter) { + p.MaxHoles = tc.maxHolesIter[iter]; + p.HoleLayerMask = tc.holeLayerMaskIter[iter]; + } + + if (tc.useMatCorrTGeo) { + p.CorrType = o2::base::PropagatorImpl::MatCorrType::USEMatCorrTGeo; + } else if (tc.useFastMaterial) { + p.CorrType = o2::base::PropagatorImpl::MatCorrType::USEMatCorrNONE; + } else { + p.CorrType = o2::base::PropagatorImpl::MatCorrType::USEMatCorrLUT; + } + + for (int i{0}; i < TrackerParamConfig::getNLayers(); ++i) { + p.SystError2Row[i] = tc.sysErr2Row[i] > 0 ? tc.sysErr2Row[i] : p.SystError2Row[i]; + p.SystError2Col[i] = tc.sysErr2Col[i] > 0 ? tc.sysErr2Col[i] : p.SystError2Col[i]; + p.AddTimeError[i] = tc.addTimeError[i]; + } + + p.MaxChi2ClusterAttachment = tc.maxChi2ClusterAttachment > 0 ? tc.maxChi2ClusterAttachment : p.MaxChi2ClusterAttachment; + p.MaxChi2NDF = tc.maxChi2NDF > 0 ? tc.maxChi2NDF : p.MaxChi2NDF; + p.ColBins = tc.LUTbinsU > 0 ? tc.LUTbinsU : p.ColBins; + p.RowBins = tc.LUTbinsV > 0 ? tc.LUTbinsV : p.RowBins; + p.PVres = tc.pvRes > 0 ? tc.pvRes : p.PVres; + p.NSigmaCut *= tc.nSigmaCut > 0 ? tc.nSigmaCut : 1.f; + p.CellDeltaTanLambdaSigma *= tc.deltaTanLres > 0 ? tc.deltaTanLres : 1.f; + p.TrackletMinPt *= tc.minPt > 0 ? tc.minPt : 1.f; + if (tc.cellRoadRCut > 0.f) { + p.CellRoadRCut = tc.cellRoadRCut; + } + if (tc.trackletMinAbsX >= 0.f) { + p.TrackletMinAbsX = tc.trackletMinAbsX; + } + for (int iD{0}; iD < 3; ++iD) { + p.Diamond[iD] = tc.diamondPos[iD]; + } + if (detId == detectors::DetID::MFT) { + p.UseDiamond = true; + p.PerPrimaryVertexProcessing = false; + } else { + p.UseDiamond = tc.useDiamond; + } + } + + if (trackParams.size() > static_cast(tc.nIterations)) { + trackParams.resize(tc.nIterations); + } + + return trackParams; +} + +} // namespace TrackingMode +} // namespace o2::itsmft diff --git a/Detectors/ITSMFT/common/tracking/src/DetectorTraits.cxx b/Detectors/ITSMFT/common/tracking/src/DetectorTraits.cxx new file mode 100644 index 0000000000000..56dbd86e78e3b --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/src/DetectorTraits.cxx @@ -0,0 +1,136 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file DetectorTraits.cxx +/// \brief Detector-specific refit and load hooks for shared CA tracking +/// + +#include "ITSMFTTracking/DetectorTraits.h" + +#include "Framework/Logger.h" +#include "ITSMFTTracking/MFTFwdTrackHelpers.h" +#include "ITStracking/TrackHelpers.h" + +namespace o2::itsmft::tracking +{ + +namespace +{ +template +bool refitSeedITS(const typename DetectorTraits::TrackSeedN& seed, + o2::its::TrackITSExt& track, + const TrackingParameters& params, + float bz, + const o2::its::TrackingFrameInfo* const tfInfos[NLayers], + const o2::its::Cluster* const unsortedClusters[NLayers], + const o2::base::PropagatorImpl* propagator) +{ + // Common TrackSeed matches o2::its::TrackSeed layout (hole-layer mask aside); refit uses ITS helper. + const auto& itsSeed = reinterpret_cast&>(seed); + return o2::its::track::refitTrack(itsSeed, + track, + params.MaxChi2ClusterAttachment, + params.MaxChi2NDF, + bz, + tfInfos, + unsortedClusters, + params.LayerxX0.data(), + params.LayerRadii.data(), + params.MinPt.data(), + propagator, + params.CorrType, + params.ReseedIfShorter, + params.ShiftRefToCluster, + params.RepeatRefitOut); +} +} // namespace + +template +bool DetectorTraits::cellsAreCompatible(const CellSeedN& currentCell, + const CellSeedN& nextCell, + const TimeFrameN& tf, + const TrackingParameters& params, + float bz) +{ + if constexpr (DetId == o2::detectors::DetID::MFT) { + return detail::mftFwdCellsAreCompatible(currentCell, nextCell, bz, params.MaxChi2ClusterAttachment); + } else { + auto nextCellSeed = nextCell; + if (!nextCellSeed.rotate(currentCell.getAlpha()) || + !nextCellSeed.propagateTo(currentCell.getX(), bz)) { + return false; + } + return currentCell.getPredictedChi2(nextCellSeed) <= params.MaxChi2ClusterAttachment; + } +} + +template +bool DetectorTraits::refitSeed(const TrackSeedN& seed, + TrackType& track, + const TrackingParameters& params, + float bz, + TimeFrameN& tf, + const o2::its::TrackingFrameInfo* const tfInfos[NLayers], + const o2::its::Cluster* const unsortedClusters[NLayers], + const o2::base::PropagatorImpl* propagator) +{ + if constexpr (DetId == o2::detectors::DetID::MFT) { + return refitTrackFwd(seed, track, tf, params, bz); + } else { + return refitSeedITS(seed, track, params, bz, tfInfos, unsortedClusters, propagator); + } +} + +template +void DetectorTraits::configureIndexTableUtils(IndexTableUtils& utils, const TrackingParameters& params) +{ + if constexpr (DetId == o2::detectors::DetID::MFT) { + constexpr float defaultYMin{-20.f}; + constexpr float defaultYMax{20.f}; + const bool hasRowRange = params.IndexRowMax != 0.f; + const float rowMin = hasRowRange ? params.IndexRowMin : defaultYMin; + const float rowMax = hasRowRange ? params.IndexRowMax : defaultYMax; + utils.setTrackingParametersXY(params, rowMin, rowMax); + } else { + utils.setTrackingParameters(params); + } +} + +template +void TrackingLoadPolicy::configureBeamPosition(TimeFrame& tf, + const TrackingParameters& p, + const o2::dataformats::MeanVertexObject* meanVertex, + bool overrideBeamEstimation) +{ + const auto& tc = TrackerParamRef::get(); + const float systErrY2 = p.SystError2Row.empty() ? 0.f : p.SystError2Row[0]; + const float layerRes = p.LayerResolution.empty() ? 0.f : p.LayerResolution[0]; + + if constexpr (DetId == o2::detectors::DetID::MFT) { + tf.setBeamPosition(p.Diamond[0], p.Diamond[1], p.DiamondCov[3], layerRes, systErrY2); + LOGP(info, "MFT CA vertex seed from diamond: x={:.4f} y={:.4f} z={:.4f}", + p.Diamond[0], p.Diamond[1], p.Diamond[2]); + } else if ((overrideBeamEstimation || tc.overrideBeamEstimation) && meanVertex != nullptr) { + tf.setBeamPosition(meanVertex->getX(), meanVertex->getY(), meanVertex->getSigmaY2(), layerRes, systErrY2); + LOGP(info, "ITS CA beam position from MeanVertex: x={:.4f} y={:.4f}", meanVertex->getX(), meanVertex->getY()); + } else if (p.UseDiamond) { + tf.setBeamPosition(p.Diamond[0], p.Diamond[1], p.DiamondCov[3], layerRes, systErrY2); + LOGP(info, "ITS CA beam position from diamond: x={:.4f} y={:.4f} z={:.4f}", + p.Diamond[0], p.Diamond[1], p.Diamond[2]); + } +} + +template struct DetectorTraits<7>; +template struct DetectorTraits<10>; +template struct TrackingLoadPolicy; +template struct TrackingLoadPolicy; + +} // namespace o2::itsmft::tracking diff --git a/Detectors/ITSMFT/common/tracking/src/IOUtils.cxx b/Detectors/ITSMFT/common/tracking/src/IOUtils.cxx new file mode 100644 index 0000000000000..2f12420730956 --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/src/IOUtils.cxx @@ -0,0 +1,233 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file IOUtils.cxx +/// \brief Shared cluster I/O utilities for ITS and MFT (based on ITStracking/IOUtils.cxx) +/// + +#include "ITSMFTTracking/IOUtils.h" +#include "ITSMFTTracking/Configuration.h" +#include "ITStracking/Cluster.h" + +#include "Framework/Logger.h" +#include "ITSBase/GeometryTGeo.h" +#include "MFTBase/GeometryTGeo.h" +#include "MathUtils/Utils.h" + +namespace +{ +constexpr int PrimaryVertexLayerId{-1}; +constexpr int EventLabelsSeparator{-1}; + +template +bool shouldApplySysErrors() +{ + const auto& conf = o2::itsmft::tracking::TrackerParamRef::get(); + for (int il = 0; il < o2::itsmft::tracking::TrackerParamRef::nLayers(); il++) { + if constexpr (DetId == o2::detectors::DetID::MFT) { + if (conf.sysErr2Row[il] > 0.f || conf.sysErr2Col[il] > 0.f) { + return true; + } + } else { + if (conf.sysErrY2[il] > 0.f || conf.sysErrZ2[il] > 0.f) { + return true; + } + } + } + return false; +} + +template +void addSysErrors(int layerId, float& sigma2Row, float& sigma2Col) +{ + const auto& conf = o2::itsmft::tracking::TrackerParamRef::get(); + if constexpr (DetId == o2::detectors::DetID::MFT) { + sigma2Row += conf.sysErr2Row[layerId]; + sigma2Col += conf.sysErr2Col[layerId]; + } else { + sigma2Row += conf.sysErrY2[layerId]; + sigma2Col += conf.sysErrZ2[layerId]; + } +} + +template +void loadClusterTrackingFrameInfoImpl(GeomT* geom, + const o2::itsmft::CompClusterExt& c, + gsl::span::iterator& pattIt, + const o2::itsmft::TopologyDictionary* dict, + int& layer, + unsigned int& clusterSize, + o2::its::TrackingFrameInfo& tfInfo, + bool applySysErrors) +{ + const auto sensorID = c.getSensorID(); + if (sensorID < 0 || sensorID >= geom->getSize()) { + LOGP(fatal, "Cluster sensorID {} is out of geometry range [0, {})", sensorID, geom->getSize()); + } + layer = geom->getLayer(sensorID); + if (layer < 0 || layer >= o2::itsmft::tracking::TrackerParamRef::nLayers()) { + LOGP(fatal, "Cluster sensorID {} maps to invalid layer {} (expected [0, {}))", + sensorID, layer, o2::itsmft::tracking::TrackerParamRef::nLayers()); + } + const auto pattID = c.getPatternID(); + if (pattID != o2::itsmft::CompCluster::InvalidPatternID && dict != nullptr && + (pattID < 0 || pattID >= dict->getSize())) { + LOGP(fatal, "Cluster patternID {} is out of dictionary range [0, {})", pattID, dict->getSize()); + } + + float sigma2Row{0.f}; + float sigma2Col{0.f}; + const auto locXYZ = o2::itsmft::ioutils::extractClusterData(c, pattIt, dict, sigma2Row, sigma2Col, &clusterSize); + + if (applySysErrors && shouldApplySysErrors()) { + addSysErrors(layer, sigma2Row, sigma2Col); + } + + if constexpr (DetId == o2::detectors::DetID::ITS) { + const auto trkXYZ = geom->getMatrixT2L(sensorID) ^ locXYZ; + const auto gloXYZ = geom->getMatrixL2G(sensorID) * locXYZ; + tfInfo = o2::its::TrackingFrameInfo{ + gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), trkXYZ.x(), geom->getSensorRefAlpha(sensorID), + std::array{trkXYZ.y(), trkXYZ.z()}, + std::array{sigma2Row, 0.f, sigma2Col}}; + } else { + if (!geom->getCacheL2G().isFilled() || geom->getCacheL2G().getSize() <= sensorID) { + LOGP(fatal, "MFT L2G matrix cache unavailable for sensorID {} (filled={}, size={})", + sensorID, geom->getCacheL2G().isFilled(), geom->getCacheL2G().getSize()); + } + const auto gloXYZ = geom->getMatrixL2G(sensorID) * locXYZ; + // ALPIDE row (local X) -> global X, column (local Z) -> global Y + tfInfo = o2::its::TrackingFrameInfo{ + gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), gloXYZ.x(), 0.f, + std::array{gloXYZ.y(), gloXYZ.z()}, + std::array{sigma2Row, 0.f, sigma2Col}}; + } +} + +template +void fillOutputClusters(GeomT* geom, + gsl::span clusters, + gsl::span::iterator& pattIt, + std::vector>& output, + const o2::itsmft::TopologyDictionary* dict, + bool applyMisalignment) +{ + for (const auto& c : clusters) { + float sigma2Row{0.f}, sigma2Col{0.f}; + const float sigmaRowCol{0.f}; // row-column covariance (unused) + const auto locXYZ = o2::itsmft::ioutils::extractClusterData(c, pattIt, dict, sigma2Row, sigma2Col); + if (applyMisalignment) { + addSysErrors(geom->getLayer(c.getSensorID()), sigma2Row, sigma2Col); + } + o2::math_utils::Point3D outXYZ{}; + if constexpr (DetId == o2::detectors::DetID::ITS) { + outXYZ = geom->getMatrixT2L(c.getSensorID()) ^ locXYZ; + } else { + outXYZ = geom->getMatrixL2G(c.getSensorID()) * locXYZ; + } + auto& cl3d = output.emplace_back(c.getSensorID(), outXYZ); + cl3d.setErrors(sigma2Row, sigma2Col, sigmaRowCol); + } +} + +} // namespace + +namespace o2::itsmft::ioutils +{ + +void fillMatrixCache(o2::detectors::DetID::ID detId) +{ + const auto mask = o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::L2G); + if (detId == o2::detectors::DetID::ITS) { + o2::its::GeometryTGeo::Instance()->fillMatrixCache(mask); + } else if (detId == o2::detectors::DetID::MFT) { + o2::mft::GeometryTGeo::Instance()->fillMatrixCache(mask); + } else { + LOGP(fatal, "Unsupported detector id {} in fillMatrixCache", static_cast(detId)); + } +} + +int getClusterLayer(o2::detectors::DetID::ID detId, const CompClusterExt& cluster) +{ + if (detId == o2::detectors::DetID::ITS) { + return o2::its::GeometryTGeo::Instance()->getLayer(cluster.getSensorID()); + } + if (detId == o2::detectors::DetID::MFT) { + return o2::mft::GeometryTGeo::Instance()->getLayer(cluster.getSensorID()); + } + LOGP(fatal, "Unsupported detector id {} in getClusterLayer", static_cast(detId)); + return -1; +} + +template +void loadClusterTrackingFrameInfo(const CompClusterExt& c, + gsl::span::iterator& pattIt, + const TopologyDictionary* dict, + int& layer, + unsigned int& clusterSize, + o2::its::TrackingFrameInfo& tfInfo, + bool applySysErrors) +{ + if constexpr (DetId == o2::detectors::DetID::ITS) { + loadClusterTrackingFrameInfoImpl(o2::its::GeometryTGeo::Instance(), c, pattIt, dict, layer, clusterSize, tfInfo, applySysErrors); + } else { + loadClusterTrackingFrameInfoImpl(o2::mft::GeometryTGeo::Instance(), c, pattIt, dict, layer, clusterSize, tfInfo, applySysErrors); + } +} + +template void loadClusterTrackingFrameInfo(const CompClusterExt& c, + gsl::span::iterator& pattIt, + const TopologyDictionary* dict, + int& layer, + unsigned int& clusterSize, + o2::its::TrackingFrameInfo& tfInfo, + bool applySysErrors); + +template void loadClusterTrackingFrameInfo(const CompClusterExt& c, + gsl::span::iterator& pattIt, + const TopologyDictionary* dict, + int& layer, + unsigned int& clusterSize, + o2::its::TrackingFrameInfo& tfInfo, + bool applySysErrors); + +template +void convertCompactClusters(gsl::span clusters, + gsl::span::iterator& pattIt, + std::vector>& output, + const TopologyDictionary* dict) +{ + const auto mask = o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::L2G); + + const bool applyMisalignment = shouldApplySysErrors(); + + if constexpr (DetId == o2::detectors::DetID::ITS) { + auto* geom = o2::its::GeometryTGeo::Instance(); + geom->fillMatrixCache(mask); + fillOutputClusters(geom, clusters, pattIt, output, dict, applyMisalignment); + } else { + auto* geom = o2::mft::GeometryTGeo::Instance(); + geom->fillMatrixCache(mask); + fillOutputClusters(geom, clusters, pattIt, output, dict, applyMisalignment); + } +} + +template void convertCompactClusters(gsl::span clusters, + gsl::span::iterator& pattIt, + std::vector>& output, + const TopologyDictionary* dict); + +template void convertCompactClusters(gsl::span clusters, + gsl::span::iterator& pattIt, + std::vector>& output, + const TopologyDictionary* dict); + +} // namespace o2::itsmft::ioutils diff --git a/Detectors/ITSMFT/common/tracking/src/ITSMFTTrackingLinkDef.h b/Detectors/ITSMFT/common/tracking/src/ITSMFTTrackingLinkDef.h new file mode 100644 index 0000000000000..b63c36e55b80a --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/src/ITSMFTTrackingLinkDef.h @@ -0,0 +1,21 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifdef __CLING__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class o2::itsmft::TrackerParamConfig < o2::detectors::DetID::MFT> + ; +#pragma link C++ class o2::conf::ConfigurableParamHelper < o2::itsmft::TrackerParamConfig < o2::detectors::DetID::MFT>> + ; + +#endif diff --git a/Detectors/ITSMFT/common/tracking/src/MFTFwdTrackHelpers.cxx b/Detectors/ITSMFT/common/tracking/src/MFTFwdTrackHelpers.cxx new file mode 100644 index 0000000000000..3d786bf83f9db --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/src/MFTFwdTrackHelpers.cxx @@ -0,0 +1,143 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file MFTFwdTrackHelpers.cxx +/// \brief MFT CA final forward Kalman refit (legacy TrackFitter) +/// + +#include "ITSMFTTracking/MFTFwdTrackHelpers.h" + +#include + +#include "Framework/Logger.h" +#include "ITStracking/Constants.h" +#include "MFTTracking/Cluster.h" +#include "MFTTracking/Constants.h" +#include "MFTTracking/MFTTrackingParam.h" +#include "MFTTracking/TrackFitter.h" + +namespace o2::itsmft::tracking +{ + +namespace +{ +template +bool refitTrackFwdImpl(const TrackSeedN& seed, + MFTCATrack& track, + const TimeFrame& tf, + const TrackingParameters& params, + float bz) +{ + TrackLTFType ltf(true); + const auto& unsorted = tf.getUnsortedClusters(); + const auto hitMask = seed.getHitLayerMask(); + + for (int layer = 0; layer < o2::mft::constants::mft::LayersNumber; ++layer) { + if (!hitMask.has(layer)) { + continue; + } + const int clIdx = seed.getCluster(layer); + if (clIdx == o2::its::constants::UnusedIndex) { + continue; + } + if (clIdx >= static_cast(unsorted[layer].size())) { + LOGP(warn, "MFT CA forward refit: invalid cluster index {} on layer {}", clIdx, layer); + return false; + } + const auto& cluster = unsorted[layer][clIdx]; + const auto& tfInfo = tf.getClusterTrackingFrameInfo(layer, cluster); + o2::mft::Cluster mftCluster{ + tfInfo.xCoordinate, tfInfo.yCoordinate, tfInfo.zCoordinate, + cluster.phi, cluster.radius, clIdx, 0, + tfInfo.covarianceTrackingFrame[0], tfInfo.covarianceTrackingFrame[2], 0}; + const int extIdx = tf.getClusterExternalIndex(layer, clIdx); + ltf.setPoint(mftCluster, layer, clIdx, {}, extIdx, tf.getClusterSize(0, extIdx)); + } + + if (ltf.getNumberOfPoints() < params.MinTrackLength) { + return false; + } + ltf.sort(); + + const auto& mftParam = o2::mft::MFTTrackingParam::Instance(); + o2::mft::TrackFitter fitter; + fitter.setBz(bz); + fitter.setMFTRadLength(mftParam.MFTRadLength); + fitter.setTrackModel(mftParam.trackmodel); + fitter.setAlignResiduals(mftParam.alignResidual); + + TrackLTFType outward = ltf; + if (!fitter.initTrack(ltf) || !fitter.fit(ltf) || + !fitter.initTrack(outward, true) || !fitter.fit(outward, true)) { + return false; + } + ltf.setOutParam(outward); + + const int nCl = ltf.getNumberOfPoints(); + if (nCl < params.MinTrackLength) { + return false; + } + if (ltf.getPt() < params.MinPt[o2::mft::constants::mft::LayersNumber - nCl]) { + return false; + } + if (static_cast(ltf.getTrackChi2() / std::max(1, 2 * nCl - 5)) > params.MaxChi2NDF) { + return false; + } + + if (params.TrackletMinAbsX > 0.f) { + if (std::abs(ltf.getX()) < params.TrackletMinAbsX) { + return false; + } + for (int layer = 0; layer < o2::mft::constants::mft::LayersNumber; ++layer) { + if (!hitMask.has(layer)) { + continue; + } + const int clIdx = seed.getCluster(layer); + if (clIdx != o2::its::constants::UnusedIndex && + std::abs(unsorted[layer][clIdx].xCoordinate) < params.TrackletMinAbsX) { + return false; + } + } + } + + auto& mftTr = track.getTrack(); + mftTr = static_cast(ltf); + mftTr.setCA(true); + + track.setPattern(0); + for (int layer = 0; layer < o2::mft::constants::mft::LayersNumber; ++layer) { + if (!hitMask.has(layer)) { + track.setClusterIndex(layer, o2::its::constants::UnusedIndex); + continue; + } + const int clIdx = seed.getCluster(layer); + track.setClusterIndex(layer, clIdx); + const int extIdx = tf.getClusterExternalIndex(layer, clIdx); + track.setClusterSize(layer, tf.getClusterSize(0, extIdx)); + } + return true; +} +} // namespace + +bool refitTrackFwd(const TrackSeedN& seed, + MFTCATrack& track, + const TimeFrame& tf, + const TrackingParameters& params, + float bz) +{ + const auto& mftParam = o2::mft::MFTTrackingParam::Instance(); + if (mftParam.forceZeroField || std::abs(bz) < 1e-6f) { + return refitTrackFwdImpl(seed, track, tf, params, 0.f); + } + return refitTrackFwdImpl(seed, track, tf, params, bz); +} + +} // namespace o2::itsmft::tracking diff --git a/Detectors/ITSMFT/common/tracking/src/TimeFrame.cxx b/Detectors/ITSMFT/common/tracking/src/TimeFrame.cxx new file mode 100644 index 0000000000000..252a78e68ff23 --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/src/TimeFrame.cxx @@ -0,0 +1,542 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file TimeFrame.cxx +/// \brief +/// + +#include + +#include "Framework/Logger.h" +#include "ITSMFTTracking/DetectorTraits.h" +#include "ITSMFTTracking/IOUtils.h" +#include "ITSMFTTracking/MFTFwdTrackHelpers.h" +#include "ITSMFTTracking/TimeFrame.h" +#include "ITStracking/MathUtils.h" +#include "DataFormatsITSMFT/CompCluster.h" +#include "DataFormatsITSMFT/ROFRecord.h" +#include "DataFormatsITSMFT/TopologyDictionary.h" +#include "DetectorsCommonDataFormats/DetID.h" + +namespace +{ +struct ClusterHelper { + float rowCoord; + float r; + int bin; + int ind; +}; +} // namespace + +namespace o2::itsmft::tracking +{ + +using o2::its::clearResizeBoundedVector; +using o2::its::deepVectorClear; +namespace math_utils = o2::its::math_utils; +using o2::itsmft::IndexTableCoordType; +using o2::itsmft::IterationStep; +using o2::itsmft::TrackingParameters; + +template +void TimeFrame::addPrimaryVertex(const Vertex& vert) +{ + mPrimaryVertices.emplace_back(vert); + if (!isBeamPositionOverridden) { + const float w = vert.getNContributors(); + mBeamPos[0] = (mBeamPos[0] * mBeamPosWeight + vert.getX() * w) / (mBeamPosWeight + w); + mBeamPos[1] = (mBeamPos[1] * mBeamPosWeight + vert.getY() * w) / (mBeamPosWeight + w); + mBeamPosWeight += w; + } +} + +template +void TimeFrame::loadROFrameData(gsl::span rofs, + gsl::span clusters, + gsl::span::iterator& pattIt, + const itsmft::TopologyDictionary* dict, + int layer, + const dataformats::MCTruthContainer* mcLabels, + o2::detectors::DetID::ID detId) +{ + mDetId = detId; + if (NLayers != constants::nLayersForDet(detId)) { + LOGP(fatal, "TimeFrame<{}> is incompatible with detector {} (expected {} layers)", + NLayers, static_cast(detId), constants::nLayersForDet(detId)); + } + ioutils::fillMatrixCache(detId); + resetROFrameData(layer); + prepareROFrameData(clusters, layer); + + // check for missing/empty/unset rofs + // the code requires consistent monotonically increasing input without gaps + const auto& timing = mROFOverlapTableView.getLayer(layer >= 0 ? layer : 0); + if (timing.mNROFsTF != rofs.size()) { + LOGP(fatal, "Received inconsistent number of rofs on layer:{} expected:{} received:{}", layer, timing.mNROFsTF, rofs.size()); + } + + for (int32_t iRof{0}; iRof < rofs.size(); ++iRof) { + const auto& rof = rofs[iRof]; + for (int clusterId{rof.getFirstEntry()}; clusterId < rof.getFirstEntry() + rof.getNEntries(); ++clusterId) { + const auto& c = clusters[clusterId]; + int lay{0}; + unsigned int clusterSize{0}; + TrackingFrameInfo tfInfo; + if (detId == o2::detectors::DetID::MFT) { + o2::itsmft::ioutils::loadClusterTrackingFrameInfo(c, pattIt, dict, lay, clusterSize, tfInfo); + } else { + o2::itsmft::ioutils::loadClusterTrackingFrameInfo(c, pattIt, dict, lay, clusterSize, tfInfo); + } + mClusterSize[layer >= 0 ? layer : 0][clusterId] = std::clamp(clusterSize, 0u, 255u); + addTrackingFrameInfoToLayer(lay, tfInfo); + addClusterToLayer(lay, tfInfo.xCoordinate, tfInfo.yCoordinate, tfInfo.zCoordinate, mUnsortedClusters[lay].size()); + addClusterExternalIndexToLayer(lay, clusterId); + } + // effectively calculating an exclusive sum + if (layer >= 0) { + mROFramesClusters[layer][iRof + 1] = mUnsortedClusters[layer].size(); + } else { + for (unsigned int iL{0}; iL < mUnsortedClusters.size(); ++iL) { + mROFramesClusters[iL][iRof + 1] = mUnsortedClusters[iL].size(); + } + } + } + + if (layer == 1 || layer == -1) { + for (auto i = 0; i < mNTrackletsPerCluster.size(); ++i) { + mNTrackletsPerCluster[i].resize(mUnsortedClusters[1].size()); + mNTrackletsPerClusterSum[i].resize(mUnsortedClusters[1].size() + 1); + } + } + + if (mcLabels != nullptr) { + mClusterLabels[layer >= 0 ? layer : 0] = mcLabels; + } else { + mClusterLabels[layer >= 0 ? layer : 0] = nullptr; + } +} + +template +void TimeFrame::resetROFrameData(int layer) +{ + if (layer >= 0) { + deepVectorClear(mUnsortedClusters[layer], getMaybeFrameworkHostResource()); + deepVectorClear(mTrackingFrameInfo[layer], getMaybeFrameworkHostResource()); + deepVectorClear(mClusterExternalIndices[layer], mMemoryPool.get()); + clearResizeBoundedVector(mROFramesClusters[layer], mROFOverlapTableView.getLayer(layer).mNROFsTF + 1, getMaybeFrameworkHostResource()); + } else { + for (int iLayer{0}; iLayer < NLayers; ++iLayer) { + deepVectorClear(mUnsortedClusters[iLayer], getMaybeFrameworkHostResource()); + deepVectorClear(mTrackingFrameInfo[iLayer], getMaybeFrameworkHostResource()); + deepVectorClear(mClusterExternalIndices[iLayer], mMemoryPool.get()); + clearResizeBoundedVector(mROFramesClusters[iLayer], mROFOverlapTableView.getLayer(iLayer).mNROFsTF + 1, getMaybeFrameworkHostResource()); + } + } +} + +template +void TimeFrame::prepareROFrameData(gsl::span clusters, int layer) +{ + if (layer >= 0) { + mUnsortedClusters[layer].reserve(clusters.size()); + mTrackingFrameInfo[layer].reserve(clusters.size()); + mClusterExternalIndices[layer].reserve(clusters.size()); + clearResizeBoundedVector(mClusterSize[layer], clusters.size(), mMemoryPool.get()); + } else { + clearResizeBoundedVector(mClusterSize[0], clusters.size(), mMemoryPool.get()); + std::array clusterCountPerLayer{0}; + for (const auto& cls : clusters) { + ++clusterCountPerLayer[ioutils::getClusterLayer(mDetId, cls)]; + } + for (int iLayer{0}; iLayer < NLayers; ++iLayer) { + mUnsortedClusters[iLayer].reserve(clusterCountPerLayer[iLayer]); + mTrackingFrameInfo[iLayer].reserve(clusterCountPerLayer[iLayer]); + mClusterExternalIndices[iLayer].reserve(clusterCountPerLayer[iLayer]); + } + } +} + +template +void TimeFrame::prepareClusters(const TrackingParameters& trkParam, const int maxLayers) +{ + const int numBins{trkParam.RowBins * trkParam.ColBins}; + const int stride{numBins + 1}; + bounded_vector cHelper(mMemoryPool.get()); + bounded_vector clsPerBin(numBins, 0, mMemoryPool.get()); + bounded_vector lutPerBin(numBins, 0, mMemoryPool.get()); + for (int iLayer{0}, stopLayer = std::min(trkParam.NLayers, maxLayers); iLayer < stopLayer; ++iLayer) { + for (int rof{0}; rof < getNrof(iLayer); ++rof) { + if (!mROFMaskView.isROFEnabled(iLayer, rof)) { + continue; + } + const auto& unsortedClusters{getUnsortedClustersOnLayer(rof, iLayer)}; + const int clustersNum{static_cast(unsortedClusters.size())}; + auto* tableBase = mIndexTables[iLayer].data() + rof * stride; + + cHelper.resize(clustersNum); + + const bool useXYBinning = mIndexTableUtils.getCoordType() == IndexTableCoordType::XY; + for (int iCluster{0}; iCluster < clustersNum; ++iCluster) { + const Cluster& c = unsortedClusters[iCluster]; + ClusterHelper& h = cHelper[iCluster]; + + const float x = c.xCoordinate - (useXYBinning ? 0.f : mBeamPos[0]); + const float y = c.yCoordinate - (useXYBinning ? 0.f : mBeamPos[1]); + const float z = c.zCoordinate; + + const float rowCoord = useXYBinning ? c.yCoordinate : math_utils::computePhi(x, y); + const float colCoord = useXYBinning ? c.xCoordinate : z; + int colBin{mIndexTableUtils.getColBinIndex(iLayer, colCoord)}; + if (colBin < 0 || colBin >= trkParam.ColBins) { + colBin = std::clamp(colBin, 0, trkParam.ColBins - 1); + mBogusClusters[iLayer]++; + } + int bin = mIndexTableUtils.getBinIndex(colBin, mIndexTableUtils.getRowBinIndex(rowCoord)); + h.rowCoord = rowCoord; + h.r = math_utils::hypot(x, y); + mMinR[iLayer] = o2::gpu::GPUCommonMath::Min(h.r, mMinR[iLayer]); + mMaxR[iLayer] = o2::gpu::GPUCommonMath::Max(h.r, mMaxR[iLayer]); + h.bin = bin; + h.ind = clsPerBin[bin]++; + } + std::exclusive_scan(clsPerBin.begin(), clsPerBin.end(), lutPerBin.begin(), 0); + + auto clusters2beSorted{getClustersOnLayer(rof, iLayer)}; + for (int iCluster{0}; iCluster < clustersNum; ++iCluster) { + const ClusterHelper& h = cHelper[iCluster]; + Cluster& c = clusters2beSorted[lutPerBin[h.bin] + h.ind]; + + c = unsortedClusters[iCluster]; + c.phi = useXYBinning ? math_utils::computePhi(c.xCoordinate, c.yCoordinate) : h.rowCoord; + c.radius = h.r; + c.indexTableBinIndex = h.bin; + } + std::copy_n(lutPerBin.data(), clsPerBin.size(), tableBase); + std::fill_n(tableBase + clsPerBin.size(), stride - clsPerBin.size(), clustersNum); + + std::fill(clsPerBin.begin(), clsPerBin.end(), 0); + cHelper.clear(); + } + } +} + +template +void TimeFrame::initVertexingTopology(const TrackingParameters& trkParam) +{ + mVertexingTopology.init(3, trkParam.MaxHoles, LayerMask{trkParam.HoleLayerMask}); +} + +template +void TimeFrame::initDefaultTrackingTopology(const TrackingParameters& trkParam, const int maxLayers) +{ + mDefaultTrackingTopology.init(maxLayers, trkParam.MaxHoles, LayerMask{trkParam.HoleLayerMask}); +} + +template +void TimeFrame::initTrackerTopologies(gsl::span trkParams, const int maxLayers) +{ + mTrackerTopologies.resize(trkParams.size()); + for (size_t iteration = 0; iteration < trkParams.size(); ++iteration) { + const int iterationMaxLayers = std::min(maxLayers, trkParams[iteration].NLayers); + mTrackerTopologies[iteration].init(iterationMaxLayers, trkParams[iteration].MaxHoles, LayerMask{trkParams[iteration].HoleLayerMask}); + } +} + +template +void TimeFrame::initialise(const TrackingParameters& trkParam, const int maxLayers, const int iteration) +{ + mTrackingTopologyView = iteration != constants::UnusedIndex ? mTrackerTopologies[iteration].getView() : (maxLayers == 3 ? mVertexingTopology.getView() : mDefaultTrackingTopology.getView()); + + if (trkParam.PassFlags[IterationStep::FirstPass]) { + deepVectorClear(mTracks); + deepVectorClear(mTracksLabel); + deepVectorClear(mLines); + deepVectorClear(mLinesLabels); + if (trkParam.PassFlags[IterationStep::ResetVertices]) { + deepVectorClear(mPrimaryVertices); + deepVectorClear(mPrimaryVerticesLabels); + } + clearResizeBoundedVector(mLinesLabels, getNrof(1), mMemoryPool.get()); + DetectorTraits::configureIndexTableUtils(mIndexTableUtils, trkParam); + clearResizeBoundedVector(mPositionResolution, trkParam.NLayers, mMemoryPool.get()); + clearResizeBoundedVector(mBogusClusters, trkParam.NLayers, mMemoryPool.get()); + deepVectorClear(mTrackletClusters); + for (unsigned int iLayer{0}; iLayer < std::min((int)mClusters.size(), maxLayers); ++iLayer) { + clearResizeBoundedVector(mClusters[iLayer], mUnsortedClusters[iLayer].size(), getMaybeFrameworkHostResource(maxLayers != NLayers)); + clearResizeBoundedVector(mUsedClusters[iLayer], mUnsortedClusters[iLayer].size(), getMaybeFrameworkHostResource(maxLayers != NLayers)); + mPositionResolution[iLayer] = o2::gpu::CAMath::Sqrt((0.5f * (trkParam.SystError2Col[iLayer] + trkParam.SystError2Row[iLayer])) + (trkParam.LayerResolution[iLayer] * trkParam.LayerResolution[iLayer])); + } + clearResizeBoundedVector(mLines, getNrof(1), mMemoryPool.get()); + clearResizeBoundedVector(mTrackletClusters, getNrof(1), mMemoryPool.get()); + + for (int iLayer{0}; iLayer < NLayers; ++iLayer) { + clearResizeBoundedVector(mIndexTables[iLayer], getNrof(iLayer) * ((trkParam.ColBins * trkParam.RowBins) + 1), getMaybeFrameworkHostResource()); + } + for (int iLayer{0}; iLayer < trkParam.NLayers; ++iLayer) { + if (trkParam.SystError2Row[iLayer] > 0.f || trkParam.SystError2Col[iLayer] > 0.f) { + for (auto& tfInfo : mTrackingFrameInfo[iLayer]) { + /// Account for alignment systematics in the cluster covariance matrix + tfInfo.covarianceTrackingFrame[0] += trkParam.SystError2Row[iLayer]; + tfInfo.covarianceTrackingFrame[2] += trkParam.SystError2Col[iLayer]; + } + } + } + + mMinR.fill(std::numeric_limits::max()); + mMaxR.fill(std::numeric_limits::min()); + } + clearResizeBoundedVector(mCells, mTrackingTopologyView.nCells, mMemoryPool.get()); + clearResizeBoundedVector(mCellsLookupTable, mTrackingTopologyView.nCells, mMemoryPool.get()); + clearResizeBoundedVector(mCellsNeighbours, mTrackingTopologyView.nCells, mMemoryPool.get()); + clearResizeBoundedVector(mCellsNeighboursTopology, mTrackingTopologyView.nCells, mMemoryPool.get()); + clearResizeBoundedVector(mCellsNeighboursLUT, mTrackingTopologyView.nCells, mMemoryPool.get()); + clearResizeBoundedVector(mCellLabels, mTrackingTopologyView.nCells, mMemoryPool.get()); + clearResizeBoundedVector(mTracklets, mTrackingTopologyView.nTransitions, mMemoryPool.get()); + clearResizeBoundedVector(mTrackletLabels, mTrackingTopologyView.nTransitions, mMemoryPool.get()); + clearResizeBoundedVector(mTrackletsLookupTable, mTrackingTopologyView.nTransitions, mMemoryPool.get()); + clearResizeBoundedVector(mTransitionPhiCuts, mTrackingTopologyView.nTransitions, mMemoryPool.get()); + clearResizeBoundedVector(mTransitionMSAngles, mTrackingTopologyView.nTransitions, mMemoryPool.get()); + mNTrackletsPerROF.resize(2); + for (auto& v : mNTrackletsPerROF) { + v = bounded_vector(getNrof(1) + 1, 0, mMemoryPool.get()); + } + if (trkParam.PassFlags[IterationStep::RebuildClusterLUT]) { + prepareClusters(trkParam, maxLayers); + } + mTotalTracklets = {0, 0}; + if (maxLayers < trkParam.NLayers) { // Vertexer only, but in both iterations + for (size_t iLayer{0}; iLayer < maxLayers; ++iLayer) { + deepVectorClear(mUsedClusters[iLayer]); + clearResizeBoundedVector(mUsedClusters[iLayer], mUnsortedClusters[iLayer].size(), mMemoryPool.get()); + } + } + + // estimate MS per layer + std::array msAngles{}; + for (unsigned int iLayer{0}; iLayer < NLayers; ++iLayer) { + if constexpr (NLayers == o2::mft::constants::mft::LayersNumber) { + msAngles[iLayer] = detail::mftLayerMSAngle(iLayer, trkParam); + } else { + msAngles[iLayer] = math_utils::MSangle(0.14f, trkParam.TrackletMinPt, trkParam.LayerxX0[iLayer]); + } + mPositionResolution[iLayer] = o2::gpu::CAMath::Sqrt((0.5f * (trkParam.SystError2Col[iLayer] + trkParam.SystError2Row[iLayer])) + (trkParam.LayerResolution[iLayer] * trkParam.LayerResolution[iLayer])); + } + + // for each transition calculate the phi-cuts + integrated MS + if constexpr (NLayers == o2::mft::constants::mft::LayersNumber) { + float oneOverR{0.001f * 0.3f * std::abs(mBz) / trkParam.TrackletMinPt}; + for (int transitionId{0}; transitionId < (int)mTracklets.size(); ++transitionId) { + const auto& transition = mTrackingTopologyView.getTransition(transitionId); + float ms2 = 0.f; + for (int layer = transition.fromLayer; layer < transition.toLayer; ++layer) { + ms2 += math_utils::Sq(msAngles[layer]); + } + mTransitionMSAngles[transitionId] = o2::gpu::CAMath::Sqrt(ms2); + const float& r1 = trkParam.LayerRadii[transition.fromLayer]; + const float& r2 = trkParam.LayerRadii[transition.toLayer]; + oneOverR = (0.5f * oneOverR >= 1.f / r2) ? (2.f / r2) - o2::constants::math::Almost0 : oneOverR; + const float res1 = o2::gpu::CAMath::Hypot(trkParam.PVres, mPositionResolution[transition.fromLayer]); + const float res2 = o2::gpu::CAMath::Hypot(trkParam.PVres, mPositionResolution[transition.toLayer]); + const float cosTheta1half = o2::gpu::CAMath::Sqrt(1.f - math_utils::Sq(0.5f * r1 * oneOverR)); + const float cosTheta2half = o2::gpu::CAMath::Sqrt(1.f - math_utils::Sq(0.5f * r2 * oneOverR)); + const float x = (r2 * cosTheta1half) - (r1 * cosTheta2half); + const float delta = o2::gpu::CAMath::Sqrt(1.f / (1.f - 0.25f * math_utils::Sq(x * oneOverR)) * + (math_utils::Sq((0.25f * r1 * r2 * math_utils::Sq(oneOverR) / cosTheta2half) + cosTheta1half) * math_utils::Sq(res1) + + math_utils::Sq((0.25f * r1 * r2 * math_utils::Sq(oneOverR) / cosTheta1half) + cosTheta2half) * math_utils::Sq(res2))); + // For MFT: bending angle (rad) at TrackletMinPt, used to widen x-y tracklet windows. + mTransitionPhiCuts[transitionId] = o2::gpu::CAMath::Min(o2::gpu::CAMath::ASin(0.5f * x * oneOverR) + 2.f * mTransitionMSAngles[transitionId] + delta, o2::constants::math::PI * 0.5f); + + deepVectorClear(mTracklets[transitionId]); + deepVectorClear(mTrackletLabels[transitionId]); + deepVectorClear(mTrackletsLookupTable[transitionId]); + mTrackletsLookupTable[transitionId].resize(mClusters[transition.fromLayer].size() + 1, 0); + } + } else { + float oneOverR{0.001f * 0.3f * std::abs(mBz) / trkParam.TrackletMinPt}; + for (int transitionId{0}; transitionId < (int)mTracklets.size(); ++transitionId) { + const auto& transition = mTrackingTopologyView.getTransition(transitionId); + float ms2 = 0.f; + for (int layer = transition.fromLayer; layer < transition.toLayer; ++layer) { + ms2 += math_utils::Sq(msAngles[layer]); + } + mTransitionMSAngles[transitionId] = o2::gpu::CAMath::Sqrt(ms2); + const float& r1 = trkParam.LayerRadii[transition.fromLayer]; + const float& r2 = trkParam.LayerRadii[transition.toLayer]; + oneOverR = (0.5 * oneOverR >= 1.f / r2) ? (2.f / r2) - o2::constants::math::Almost0 : oneOverR; + const float res1 = o2::gpu::CAMath::Hypot(trkParam.PVres, mPositionResolution[transition.fromLayer]); + const float res2 = o2::gpu::CAMath::Hypot(trkParam.PVres, mPositionResolution[transition.toLayer]); + const float cosTheta1half = o2::gpu::CAMath::Sqrt(1.f - math_utils::Sq(0.5f * r1 * oneOverR)); + const float cosTheta2half = o2::gpu::CAMath::Sqrt(1.f - math_utils::Sq(0.5f * r2 * oneOverR)); + float x = (r2 * cosTheta1half) - (r1 * cosTheta2half); + float delta = o2::gpu::CAMath::Sqrt(1.f / (1.f - 0.25f * math_utils::Sq(x * oneOverR)) * (math_utils::Sq((0.25f * r1 * r2 * math_utils::Sq(oneOverR) / cosTheta2half) + cosTheta1half) * math_utils::Sq(res1) + math_utils::Sq((0.25f * r1 * r2 * math_utils::Sq(oneOverR) / cosTheta1half) + cosTheta2half) * math_utils::Sq(res2))); + /// the expression std::asin(0.5f * x * oneOverR) is equivalent to std::aCos(0.5f * r1 * oneOverR) - std::acos(0.5 * r2 * oneOverR) + mTransitionPhiCuts[transitionId] = o2::gpu::CAMath::Min(o2::gpu::CAMath::ASin(0.5f * x * oneOverR) + 2.f * mTransitionMSAngles[transitionId] + delta, o2::constants::math::PI * 0.5f); + + deepVectorClear(mTracklets[transitionId]); + deepVectorClear(mTrackletLabels[transitionId]); + deepVectorClear(mTrackletsLookupTable[transitionId]); + mTrackletsLookupTable[transitionId].resize(mClusters[transition.fromLayer].size() + 1, 0); + } + } + + for (int cellId{0}; cellId < (int)mCells.size(); ++cellId) { + deepVectorClear(mCells[cellId]); + deepVectorClear(mCellsLookupTable[cellId]); + deepVectorClear(mCellsNeighbours[cellId]); + deepVectorClear(mCellsNeighboursTopology[cellId]); + deepVectorClear(mCellsNeighboursLUT[cellId]); + deepVectorClear(mCellLabels[cellId]); + } +} + +template +unsigned long TimeFrame::getArtefactsMemory() const +{ + unsigned long size{0}; + for (const auto& trkl : mTracklets) { + size += sizeof(Tracklet) * trkl.size(); + } + for (const auto& cells : mCells) { + size += sizeof(CellSeedN) * cells.size(); + } + for (const auto& cellsN : mCellsNeighbours) { + size += sizeof(int) * cellsN.size(); + } + for (const auto& cellsN : mCellsNeighboursTopology) { + size += sizeof(int) * cellsN.size(); + } + return size; +} + +template +void TimeFrame::printArtefactsMemory() const +{ + LOGP(info, "TimeFrame: Artefacts occupy {:.2f} MB", getArtefactsMemory() / constants::MB); +} + +template +void TimeFrame::computeTrackletsPerROFScans() +{ + for (ushort iLayer = 0; iLayer < 2; ++iLayer) { + for (unsigned int iRof{0}; iRof < getNrof(1); ++iRof) { + if (mROFMaskView.isROFEnabled(1, iRof)) { + mTotalTracklets[iLayer] += mNTrackletsPerROF[iLayer][iRof]; + } + } + std::exclusive_scan(mNTrackletsPerROF[iLayer].begin(), mNTrackletsPerROF[iLayer].end(), mNTrackletsPerROF[iLayer].begin(), 0); + std::exclusive_scan(mNTrackletsPerCluster[iLayer].begin(), mNTrackletsPerCluster[iLayer].end(), mNTrackletsPerClusterSum[iLayer].begin(), 0); + } +} + +template +void TimeFrame::setMemoryPool(std::shared_ptr pool) +{ + mMemoryPool = pool; + + auto initVector = [&](bounded_vector& vec, bool useExternal = false) { + std::pmr::memory_resource* mr = (useExternal) ? mExtMemoryPool.get() : mMemoryPool.get(); + deepVectorClear(vec, mr); + }; + + auto initContainers = [&](Container& container, bool useExternal = false) { + for (auto& v : container) { + initVector(v, useExternal); + } + }; + + // these will only reside on the host for the cpu part + initContainers(mClusterExternalIndices); + initContainers(mNTrackletsPerCluster); + initContainers(mNTrackletsPerClusterSum); + initContainers(mNClustersPerROF); + initVector(mPrimaryVertices); + initVector(mTransitionPhiCuts); + initVector(mTransitionMSAngles); + initVector(mPositionResolution); + initContainers(mClusterSize); + initVector(mPValphaX); + initVector(mBogusClusters); + initContainers(mTrackletsIndexROF); + initVector(mTracks); + initContainers(mTracklets); + initContainers(mCells); + initContainers(mCellsNeighbours); + initContainers(mCellsLookupTable); + // MC info (we don't know if we have MC) + initVector(mPrimaryVerticesLabels); + initContainers(mLinesLabels); + initContainers(mTrackletLabels); + initContainers(mCellLabels); + initVector(mTracksLabel); + // these will use possibly an externally provided allocator + initContainers(mClusters, hasFrameworkAllocator()); + initContainers(mUsedClusters, hasFrameworkAllocator()); + initContainers(mUnsortedClusters, hasFrameworkAllocator()); + initContainers(mIndexTables, hasFrameworkAllocator()); + initContainers(mTrackingFrameInfo, hasFrameworkAllocator()); + initContainers(mROFramesClusters, hasFrameworkAllocator()); +} + +template +void TimeFrame::setFrameworkAllocator(ExternalAllocator* ext) +{ + mExternalAllocator = ext; + mExtMemoryPool = std::make_shared(mExternalAllocator); +} + +template +void TimeFrame::wipe() +{ + deepVectorClear(mTracks); + deepVectorClear(mTracklets); + deepVectorClear(mCells); + deepVectorClear(mCellsNeighbours); + deepVectorClear(mCellsNeighboursTopology); + deepVectorClear(mCellsLookupTable); + deepVectorClear(mPrimaryVertices); + deepVectorClear(mTrackletsLookupTable); + deepVectorClear(mClusterExternalIndices); + deepVectorClear(mNTrackletsPerCluster); + deepVectorClear(mNTrackletsPerClusterSum); + deepVectorClear(mNClustersPerROF); + deepVectorClear(mTransitionPhiCuts); + deepVectorClear(mTransitionMSAngles); + deepVectorClear(mPositionResolution); + deepVectorClear(mClusterSize); + deepVectorClear(mPValphaX); + deepVectorClear(mBogusClusters); + deepVectorClear(mTrackletsIndexROF); + deepVectorClear(mTrackletClusters); + deepVectorClear(mLines); + // if we use the external host allocator then the assumption is that we + // don't clear the memory ourself + if (!hasFrameworkAllocator()) { + deepVectorClear(mClusters); + deepVectorClear(mUsedClusters); + deepVectorClear(mUnsortedClusters); + deepVectorClear(mIndexTables); + deepVectorClear(mTrackingFrameInfo); + deepVectorClear(mROFramesClusters); + } + // only needed to clear if we have MC info + if (hasMCinformation()) { + deepVectorClear(mLinesLabels); + deepVectorClear(mPrimaryVerticesLabels); + deepVectorClear(mTrackletLabels); + deepVectorClear(mCellLabels); + deepVectorClear(mTracksLabel); + } +} + +template class TimeFrame; +template class TimeFrame; + +} // namespace o2::itsmft::tracking diff --git a/Detectors/ITSMFT/common/tracking/src/TrackerTraits.cxx b/Detectors/ITSMFT/common/tracking/src/TrackerTraits.cxx new file mode 100644 index 0000000000000..166e9f504ed40 --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/src/TrackerTraits.cxx @@ -0,0 +1,1084 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file TrackerTraits.cxx +/// \brief +/// + +#include +#include +#include +#include + +#include +#include + +#include "DetectorsBase/Propagator.h" +#include "Framework/Logger.h" +#include "GPUCommonMath.h" +#include "ITStracking/BoundedAllocator.h" +#include "ITSMFTTracking/Cell.h" +#include "ITStracking/Constants.h" +#include "ITSMFTTracking/Configuration.h" +#include "ITSMFTTracking/DetectorTraits.h" +#include "ITSMFTTracking/MFTFwdTrackHelpers.h" +#include "ITSMFTTracking/IndexTableUtils.h" +#include "ITSMFTTracking/LayerMask.h" +#include "ITStracking/ROFLookupTables.h" +#include "ITSMFTTracking/TrackerTraits.h" +#include "ITStracking/TrackHelpers.h" +#include "ITStracking/Tracklet.h" +#include "SimulationDataFormat/MCCompLabel.h" + +namespace o2::itsmft::tracking +{ + +namespace math_utils = o2::its::math_utils; +using o2::its::deepVectorClear; +using o2::its::TimeEstBC; + +struct PassMode { + using OnePass = std::integral_constant; + using TwoPassCount = std::integral_constant; + using TwoPassInsert = std::integral_constant; +}; + +template +void TrackerTraits::computeLayerTracklets(const int iteration, int iVertex) +{ + const auto topology = mTimeFrame->getTrackingTopologyView(); + for (int transitionId = 0; transitionId < topology.nTransitions; ++transitionId) { + mTimeFrame->getTracklets()[transitionId].clear(); + mTimeFrame->getTrackletsLabel(transitionId).clear(); + std::fill(mTimeFrame->getTrackletsLookupTable()[transitionId].begin(), mTimeFrame->getTrackletsLookupTable()[transitionId].end(), 0); + } + + const Vertex diamondVert(mTrkParams[iteration].Diamond, mTrkParams[iteration].DiamondCov, 1, 1.f); + gsl::span diamondSpan(&diamondVert, 1); + + mTaskArena->execute([&] { + auto forTracklets = [&](auto Tag, int transitionId, int pivotROF, int base, int& offset) -> int { + const auto& transition = topology.getTransition(transitionId); + if (!mTimeFrame->getROFMaskView().isROFEnabled(transition.fromLayer, pivotROF)) { + return 0; + } + gsl::span primaryVertices = mTrkParams[iteration].UseDiamond ? diamondSpan : mTimeFrame->getPrimaryVertices(transition.fromLayer, pivotROF); + if (primaryVertices.empty()) { + return 0; + } + const int startVtx = iVertex >= 0 ? iVertex : 0; + const int endVtx = iVertex >= 0 ? o2::gpu::CAMath::Min(iVertex + 1, int(primaryVertices.size())) : int(primaryVertices.size()); + if (endVtx <= startVtx || (iVertex + 1) > primaryVertices.size()) { + return 0; + } + + const auto& rofOverlap = mTimeFrame->getROFOverlapTableView().getOverlap(transition.fromLayer, transition.toLayer, pivotROF); + if (!rofOverlap.getEntries()) { + return 0; + } + + int localCount = 0; + auto& tracklets = mTimeFrame->getTracklets()[transitionId]; + auto layer0 = mTimeFrame->getClustersOnLayer(pivotROF, transition.fromLayer); + if (layer0.empty()) { + return 0; + } + + const float meanDeltaR = mTrkParams[iteration].LayerRadii[transition.toLayer] - mTrkParams[iteration].LayerRadii[transition.fromLayer]; + const float phiCut = mTimeFrame->getTransitionPhiCut(transitionId); + const float msAngle = mTimeFrame->getTransitionMSAngle(transitionId); + const bool useDiamond = mTrkParams[iteration].UseDiamond; + const bool isMFT = DetectorTraits::DetId == o2::detectors::DetID::MFT; + const float meanDeltaZ = isMFT ? detail::mftLayerZ(transition.toLayer) - detail::mftLayerZ(transition.fromLayer) : 0.f; + + for (int iCluster = 0; iCluster < int(layer0.size()); ++iCluster) { + const Cluster& currentCluster = layer0[iCluster]; + const int currentSortedIndex = mTimeFrame->getSortedIndex(pivotROF, transition.fromLayer, iCluster); + if (mTimeFrame->isClusterUsed(transition.fromLayer, currentCluster.clusterId)) { + continue; + } + const float inverseR0 = 1.f / currentCluster.radius; + + for (int iV = startVtx; iV < endVtx; ++iV) { + const auto& pv = primaryVertices[iV]; + if (!useDiamond && !mTimeFrame->getROFVertexLookupTableView().isVertexCompatible(transition.fromLayer, pivotROF, pv)) { + continue; + } + if (pv.isFlagSet(Vertex::Flags::UPCMode) != mTrkParams[iteration].PassFlags[IterationStep::SelectUPCVertices]) { + continue; + } + float colWindow = 0.f; + float rowWindow = 0.f; + float sigmaX = 0.f; + float sigmaY = 0.f; + float sigmaZ = 0.f; + float lutRangeMin = 0.f; + float lutRangeMax = 0.f; + float xProj = 0.f; + float yProj = 0.f; + if (isMFT) { + const auto& tfInfo = mTimeFrame->getClusterTrackingFrameInfo(transition.fromLayer, currentCluster); + detail::mftTrackletProject(currentCluster.xCoordinate, currentCluster.yCoordinate, currentCluster.zCoordinate, + pv.getX(), pv.getY(), pv.getZ(), + transition.fromLayer, transition.toLayer, getBz(), mTrkParams[iteration].TrackletMinPt, + xProj, yProj); + detail::mftTrackletSigmaXY(currentCluster.xCoordinate, currentCluster.yCoordinate, + pv.getX(), pv.getY(), pv.getZ(), + tfInfo.covarianceTrackingFrame[0], tfInfo.covarianceTrackingFrame[2], + pv.getSigmaX2(), pv.getSigmaY2(), pv.getSigmaZ2(), + transition.fromLayer, transition.toLayer, + mTrkParams[iteration].LayerRadii[transition.fromLayer], + meanDeltaZ, msAngle, phiCut, xProj, yProj, sigmaX, sigmaY); + const float zSpread = mTrkParams[iteration].NSigmaCut * pv.getSigmaZ(); + const float zVtxMin = pv.getZ() - zSpread; + const float zVtxMax = pv.getZ() + zSpread; + const float zLayerFrom = detail::mftLayerZ(transition.fromLayer); + const float zLayerTo = detail::mftLayerZ(transition.toLayer); + const float absZFrom = std::abs(zLayerFrom); + const float absZTo = std::abs(zLayerTo); + const float denomMin = zVtxMax + absZFrom; + const float denomMax = absZFrom + zVtxMin; + lutRangeMin = (std::abs(denomMin) > 1.e-6f) ? currentCluster.radius * (zVtxMax + absZTo) / denomMin : currentCluster.radius; + lutRangeMax = (std::abs(denomMax) > 1.e-6f) ? currentCluster.radius * (absZTo + zVtxMin) / denomMax : currentCluster.radius; + if (lutRangeMin > lutRangeMax) { + const float tmp = lutRangeMin; + lutRangeMin = lutRangeMax; + lutRangeMax = tmp; + } + colWindow = sigmaX * mTrkParams[iteration].NSigmaCut; + rowWindow = sigmaY * mTrkParams[iteration].NSigmaCut; + } else { + const float resolution = o2::gpu::CAMath::Sqrt(math_utils::Sq(mTimeFrame->getPositionResolution(transition.fromLayer)) + math_utils::Sq(mTrkParams[iteration].PVres) / float(pv.getNContributors())); + const float tanLambda = (currentCluster.zCoordinate - pv.getZ()) * inverseR0; + lutRangeMin = tanLambda * (mTimeFrame->getMinR(transition.toLayer) - currentCluster.radius) + currentCluster.zCoordinate; + lutRangeMax = tanLambda * (mTimeFrame->getMaxR(transition.toLayer) - currentCluster.radius) + currentCluster.zCoordinate; + const float sqInvDeltaZ0 = 1.f / (math_utils::Sq(currentCluster.zCoordinate - pv.getZ()) + constants::Tolerance); + sigmaZ = o2::gpu::CAMath::Sqrt((math_utils::Sq(resolution) * math_utils::Sq(tanLambda) * ((math_utils::Sq(inverseR0) + sqInvDeltaZ0) * math_utils::Sq(meanDeltaR) + 1.f)) + math_utils::Sq(meanDeltaR * msAngle)); + colWindow = sigmaZ * mTrkParams[iteration].NSigmaCut; + rowWindow = phiCut; + } + const auto bins = isMFT + ? o2::itsmft::getBinsRectClusterAtProj(xProj, yProj, transition.toLayer, + lutRangeMin, lutRangeMax, colWindow, rowWindow, + mTimeFrame->getIndexTableUtils()) + : o2::itsmft::getBinsRectCluster(currentCluster, transition.fromLayer, transition.toLayer, + lutRangeMin, lutRangeMax, colWindow, rowWindow, + mTimeFrame->getIndexTableUtils()); + if (bins.x < 0) { + continue; + } + int rowBinsNum = bins.w - bins.y + 1; + if (rowBinsNum < 0) { + rowBinsNum += mTrkParams[iteration].RowBins; + } + + for (int targetROF = rofOverlap.getFirstEntry(); targetROF < rofOverlap.getEntriesBound(); ++targetROF) { + if (!mTimeFrame->getROFMaskView().isROFEnabled(transition.toLayer, targetROF)) { + continue; + } + auto layer1 = mTimeFrame->getClustersOnLayer(targetROF, transition.toLayer); + if (layer1.empty()) { + continue; + } + const auto ts = mTimeFrame->getROFOverlapTableView().getTimeStamp(transition.fromLayer, pivotROF, transition.toLayer, targetROF); + if (!useDiamond && !ts.isCompatible(pv.getTimeStamp())) { + continue; + } + const auto& targetIndexTable = mTimeFrame->getIndexTable(targetROF, transition.toLayer); + const int colBinRange = (bins.z - bins.x) + 1; + for (int iRow = 0; iRow < rowBinsNum; ++iRow) { + const int iRowBin = isMFT ? (bins.y + iRow) : ((bins.y + iRow) % mTrkParams[iteration].RowBins); + if (isMFT && iRowBin >= mTrkParams[iteration].RowBins) { + break; + } + const int firstBinIdx = mTimeFrame->getIndexTableUtils().getBinIndex(bins.x, iRowBin); + const int maxBinIdx = firstBinIdx + colBinRange; + const int firstRow = targetIndexTable[firstBinIdx]; + const int lastRow = targetIndexTable[maxBinIdx]; + for (int iNext = firstRow; iNext < lastRow; ++iNext) { + if (iNext >= int(layer1.size())) { + break; + } + const Cluster& nextCluster = layer1[iNext]; + if (mTimeFrame->isClusterUsed(transition.toLayer, nextCluster.clusterId)) { + continue; + } + + bool acceptTracklet = false; + float tanL = 0.f; + if (isMFT) { + const float dx = nextCluster.xCoordinate - xProj; + const float dy = nextCluster.yCoordinate - yProj; + const float invSigmaX2 = (sigmaX > 0.f) ? 1.f / (sigmaX * sigmaX) : 0.f; + const float invSigmaY2 = (sigmaY > 0.f) ? 1.f / (sigmaY * sigmaY) : 0.f; + const float transChi2 = dx * dx * invSigmaX2 + dy * dy * invSigmaY2; + const float nSigmaCut2 = math_utils::Sq(mTrkParams[iteration].NSigmaCut); + if (transChi2 < nSigmaCut2) { + acceptTracklet = std::abs(meanDeltaZ) > 1.e-6f; + tanL = (currentCluster.zCoordinate - nextCluster.zCoordinate) / meanDeltaZ; + } + } else { + const float tanLambda = (currentCluster.zCoordinate - pv.getZ()) * inverseR0; + const float deltaZ = o2::gpu::CAMath::Abs((tanLambda * (nextCluster.radius - currentCluster.radius)) + currentCluster.zCoordinate - nextCluster.zCoordinate); + if (deltaZ / sigmaZ < mTrkParams[iteration].NSigmaCut && + math_utils::isPhiDifferenceBelow(currentCluster.phi, nextCluster.phi, phiCut)) { + acceptTracklet = true; + tanL = (currentCluster.zCoordinate - nextCluster.zCoordinate) / (currentCluster.radius - nextCluster.radius); + } + } + + if (acceptTracklet) { + const float phi{o2::gpu::CAMath::ATan2(currentCluster.yCoordinate - nextCluster.yCoordinate, currentCluster.xCoordinate - nextCluster.xCoordinate)}; + if constexpr (decltype(Tag)::value == PassMode::OnePass::value) { + tracklets.emplace_back(currentSortedIndex, mTimeFrame->getSortedIndex(targetROF, transition.toLayer, iNext), tanL, phi, ts); + } else if constexpr (decltype(Tag)::value == PassMode::TwoPassCount::value) { + ++localCount; + } else if constexpr (decltype(Tag)::value == PassMode::TwoPassInsert::value) { + const int idx = base + offset++; + tracklets[idx] = Tracklet(currentSortedIndex, mTimeFrame->getSortedIndex(targetROF, transition.toLayer, iNext), tanL, phi, ts); + } + } + } + } + } + } + } + return localCount; + }; + + int dummy{0}; + if (mTaskArena->max_concurrency() <= 1) { + for (int transitionId{0}; transitionId < topology.nTransitions; ++transitionId) { + const int fromLayer = topology.getTransition(transitionId).fromLayer; + const int startROF = 0, endROF = mTimeFrame->getROFOverlapTableView().getLayer(fromLayer).mNROFsTF; + for (int pivotROF{startROF}; pivotROF < endROF; ++pivotROF) { + forTracklets(PassMode::OnePass{}, transitionId, pivotROF, 0, dummy); + } + } + } else { + tbb::parallel_for(0, static_cast(topology.nTransitions), [&](const int transitionId) { + const int fromLayer = topology.getTransition(transitionId).fromLayer; + const int startROF = 0, endROF = mTimeFrame->getROFOverlapTableView().getLayer(fromLayer).mNROFsTF; + bounded_vector perROFCount((endROF - startROF) + 1, mMemoryPool.get()); + tbb::parallel_for(startROF, endROF, [&](const int pivotROF) { + perROFCount[pivotROF - startROF] = forTracklets(PassMode::TwoPassCount{}, transitionId, pivotROF, 0, dummy); + }); + std::exclusive_scan(perROFCount.begin(), perROFCount.end(), perROFCount.begin(), 0); + const int nTracklets = perROFCount.back(); + mTimeFrame->getTracklets()[transitionId].resize(nTracklets); + if (nTracklets == 0) { + return; + } + tbb::parallel_for(startROF, endROF, [&](const int pivotROF) { + int baseIdx = perROFCount[pivotROF - startROF]; + if (baseIdx == perROFCount[pivotROF + 1 - startROF]) { + return; + } + int localIdx = 0; + forTracklets(PassMode::TwoPassInsert{}, transitionId, pivotROF, baseIdx, localIdx); + }); + }); + } + + tbb::parallel_for(0, static_cast(topology.nTransitions), [&](const int transitionId) { + /// Sort tracklets & remove duplicates + // duplicates can exist simply since we evaluate per vertex + auto& trkl{mTimeFrame->getTracklets()[transitionId]}; + std::sort(trkl.begin(), trkl.end()); + trkl.erase(std::unique(trkl.begin(), trkl.end()), trkl.end()); + trkl.shrink_to_fit(); + auto& lut{mTimeFrame->getTrackletsLookupTable()[transitionId]}; + if (!trkl.empty()) { + for (const auto& tkl : trkl) { + lut[tkl.firstClusterIndex + 1]++; + } + std::inclusive_scan(lut.begin(), lut.end(), lut.begin()); + } + }); + + /// Create tracklets labels + if (mTimeFrame->hasMCinformation() && mTrkParams[iteration].CreateArtefactLabels) { + tbb::parallel_for(0, static_cast(topology.nTransitions), [&](const int transitionId) { + const auto& transition = topology.getTransition(transitionId); + for (auto& trk : mTimeFrame->getTracklets()[transitionId]) { + MCCompLabel label; + int currentId{mTimeFrame->getClusters()[transition.fromLayer][trk.firstClusterIndex].clusterId}; + int nextId{mTimeFrame->getClusters()[transition.toLayer][trk.secondClusterIndex].clusterId}; + for (const auto& lab1 : mTimeFrame->getClusterLabels(transition.fromLayer, currentId)) { + for (const auto& lab2 : mTimeFrame->getClusterLabels(transition.toLayer, nextId)) { + if (lab1 == lab2 && lab1.isValid()) { + label = lab1; + break; + } + } + if (label.isValid()) { + break; + } + } + mTimeFrame->getTrackletsLabel(transitionId).emplace_back(label); + } + }); + } + }); +} + +template +void TrackerTraits::computeLayerCells(const int iteration) +{ + const auto topology = mTimeFrame->getTrackingTopologyView(); + for (int cellTopologyId = 0; cellTopologyId < topology.nCells; ++cellTopologyId) { + deepVectorClear(mTimeFrame->getCells()[cellTopologyId]); + deepVectorClear(mTimeFrame->getCellsLookupTable()[cellTopologyId]); + if (mTimeFrame->hasMCinformation() && mTrkParams[iteration].CreateArtefactLabels) { + deepVectorClear(mTimeFrame->getCellsLabel(cellTopologyId)); + } + } + + mTaskArena->execute([&] { + auto forTrackletCells = [&](auto Tag, int cellTopologyId, bounded_vector& layerCells, int iTracklet, int offset = 0) -> int { + const auto& cellTopology = topology.getCell(cellTopologyId); + const auto& firstTransition = topology.getTransition(cellTopology.firstTransition); + const auto& secondTransition = topology.getTransition(cellTopology.secondTransition); + const Tracklet& currentTracklet{mTimeFrame->getTracklets()[cellTopology.firstTransition][iTracklet]}; + const int nextLayerClusterIndex{currentTracklet.secondClusterIndex}; + const int nextLayerFirstTrackletIndex{mTimeFrame->getTrackletsLookupTable()[cellTopology.secondTransition][nextLayerClusterIndex]}; + const int nextLayerLastTrackletIndex{mTimeFrame->getTrackletsLookupTable()[cellTopology.secondTransition][nextLayerClusterIndex + 1]}; + int foundCells{0}; + for (int iNextTracklet{nextLayerFirstTrackletIndex}; iNextTracklet < nextLayerLastTrackletIndex; ++iNextTracklet) { + const Tracklet& nextTracklet{mTimeFrame->getTracklets()[cellTopology.secondTransition][iNextTracklet]}; + if (nextTracklet.firstClusterIndex != nextLayerClusterIndex) { + break; + } + if (!currentTracklet.getTimeStamp().isCompatible(nextTracklet.getTimeStamp())) { + continue; + } + + const float deltaTanLambdaSigma = std::abs(currentTracklet.tanLambda - nextTracklet.tanLambda) / mTrkParams[iteration].CellDeltaTanLambdaSigma; + if (deltaTanLambdaSigma < mTrkParams[iteration].NSigmaCut) { + + /// Track seed preparation. Clusters are numbered progressively from the innermost going outward. + const int clusId[3]{ + mTimeFrame->getClusters()[firstTransition.fromLayer][currentTracklet.firstClusterIndex].clusterId, + mTimeFrame->getClusters()[firstTransition.toLayer][nextTracklet.firstClusterIndex].clusterId, + mTimeFrame->getClusters()[secondTransition.toLayer][nextTracklet.secondClusterIndex].clusterId}; + const int hitLayers[3]{firstTransition.fromLayer, firstTransition.toLayer, secondTransition.toLayer}; + + float chi2{0.f}; + bool good{false}; + if constexpr (DetectorTraits::DetId == o2::detectors::DetID::MFT) { + const auto& cluster1_glo = mTimeFrame->getUnsortedClusters()[hitLayers[0]][clusId[0]]; + const auto& cluster2_glo = mTimeFrame->getUnsortedClusters()[hitLayers[1]][clusId[1]]; + const auto& cluster3_glo = mTimeFrame->getUnsortedClusters()[hitLayers[2]][clusId[2]]; + const float r2Cut = mTrkParams[iteration].CellRoadRCut * mTrkParams[iteration].CellRoadRCut; + if (!detail::validateMFTCellClusters(cluster1_glo, hitLayers[0], + cluster2_glo, hitLayers[1], + cluster3_glo, hitLayers[2], + r2Cut)) { + continue; + } + o2::track::TrackParCovFwd fwdTrack; + good = detail::mftFwdFitCellClusters(hitLayers, clusId, *mTimeFrame, mTrkParams[iteration], getBz(), fwdTrack, chi2); + if (good) { + TimeEstBC ts = currentTracklet.getTimeStamp(); + ts += nextTracklet.getTimeStamp(); + if constexpr (decltype(Tag)::value == PassMode::OnePass::value) { + layerCells.emplace_back(cellTopology.hitLayerMask, clusId[0], clusId[1], clusId[2], iTracklet, iNextTracklet, fwdTrack, chi2, ts); + ++foundCells; + } else if constexpr (decltype(Tag)::value == PassMode::TwoPassCount::value) { + ++foundCells; + } else if constexpr (decltype(Tag)::value == PassMode::TwoPassInsert::value) { + layerCells[offset++] = CellSeedN(cellTopology.hitLayerMask, clusId[0], clusId[1], clusId[2], iTracklet, iNextTracklet, fwdTrack, chi2, ts); + ++foundCells; + } else { + static_assert(false, "Unknown mode!"); + } + } + } else { + const auto& cluster1_glo = mTimeFrame->getUnsortedClusters()[firstTransition.fromLayer][clusId[0]]; + const auto& cluster2_glo = mTimeFrame->getUnsortedClusters()[firstTransition.toLayer][clusId[1]]; + const auto& cluster3_tf = mTimeFrame->getTrackingFrameInfoOnLayer(secondTransition.toLayer)[clusId[2]]; + auto track{o2::its::track::buildTrackSeed(cluster1_glo, cluster2_glo, cluster3_tf, mBz)}; + + for (int iC{2}; iC--;) { + const int hitLayer = hitLayers[iC]; + const TrackingFrameInfo& trackingHit = mTimeFrame->getTrackingFrameInfoOnLayer(hitLayer)[clusId[iC]]; + + if (!track.rotate(trackingHit.alphaTrackingFrame)) { + break; + } + + if (!track.propagateTo(trackingHit.xTrackingFrame, getBz())) { + break; + } + + if (!track.correctForMaterial(mTrkParams[iteration].LayerxX0[hitLayer], mTrkParams[iteration].LayerxX0[hitLayer] * constants::Radl * constants::Rho, true)) { + break; + } + + const auto predChi2{track.getPredictedChi2Quiet(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)}; + if (!iC && predChi2 > mTrkParams[iteration].MaxChi2ClusterAttachment) { + break; + } + + if (!track.o2::track::TrackParCov::update(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)) { + break; + } + + good = !iC; + chi2 += predChi2; + } + if (good) { + TimeEstBC ts = currentTracklet.getTimeStamp(); + ts += nextTracklet.getTimeStamp(); + if constexpr (decltype(Tag)::value == PassMode::OnePass::value) { + layerCells.emplace_back(cellTopology.hitLayerMask, clusId[0], clusId[1], clusId[2], iTracklet, iNextTracklet, track, chi2, ts); + ++foundCells; + } else if constexpr (decltype(Tag)::value == PassMode::TwoPassCount::value) { + ++foundCells; + } else if constexpr (decltype(Tag)::value == PassMode::TwoPassInsert::value) { + layerCells[offset++] = CellSeedN(cellTopology.hitLayerMask, clusId[0], clusId[1], clusId[2], iTracklet, iNextTracklet, track, chi2, ts); + ++foundCells; + } else { + static_assert(false, "Unknown mode!"); + } + } + } + } + } + return foundCells; + }; + + for (int cellTopologyId = 0; cellTopologyId < topology.nCells; ++cellTopologyId) { + const auto& cellTopology = topology.getCell(cellTopologyId); + if (mTimeFrame->getTracklets()[cellTopology.firstTransition].empty() || + mTimeFrame->getTracklets()[cellTopology.secondTransition].empty()) { + continue; + } + + auto& layerCells = mTimeFrame->getCells()[cellTopologyId]; + const int currentLayerTrackletsNum{static_cast(mTimeFrame->getTracklets()[cellTopology.firstTransition].size())}; + bounded_vector perTrackletCount(currentLayerTrackletsNum + 1, 0, mMemoryPool.get()); + if (mTaskArena->max_concurrency() <= 1) { + for (int iTracklet{0}; iTracklet < currentLayerTrackletsNum; ++iTracklet) { + perTrackletCount[iTracklet] = forTrackletCells(PassMode::OnePass{}, cellTopologyId, layerCells, iTracklet); + } + std::exclusive_scan(perTrackletCount.begin(), perTrackletCount.end(), perTrackletCount.begin(), 0); + } else { + tbb::parallel_for(0, currentLayerTrackletsNum, [&](const int iTracklet) { + perTrackletCount[iTracklet] = forTrackletCells(PassMode::TwoPassCount{}, cellTopologyId, layerCells, iTracklet); + }); + + std::exclusive_scan(perTrackletCount.begin(), perTrackletCount.end(), perTrackletCount.begin(), 0); + auto totalCells{perTrackletCount.back()}; + if (totalCells == 0) { + auto& lut = mTimeFrame->getCellsLookupTable()[cellTopologyId]; + lut.resize(currentLayerTrackletsNum + 1); + std::fill(lut.begin(), lut.end(), 0); + continue; + } + layerCells.resize(totalCells); + + tbb::parallel_for(0, currentLayerTrackletsNum, [&](const int iTracklet) { + int offset = perTrackletCount[iTracklet]; + if (offset == perTrackletCount[iTracklet + 1]) { + return; + } + forTrackletCells(PassMode::TwoPassInsert{}, cellTopologyId, layerCells, iTracklet, offset); + }); + } + + auto& lut = mTimeFrame->getCellsLookupTable()[cellTopologyId]; + lut.resize(currentLayerTrackletsNum + 1); + std::copy_n(perTrackletCount.begin(), currentLayerTrackletsNum + 1, lut.begin()); + + if (mTimeFrame->hasMCinformation() && mTrkParams[iteration].CreateArtefactLabels) { + auto& labels = mTimeFrame->getCellsLabel(cellTopologyId); + labels.reserve(layerCells.size()); + for (const auto& cell : layerCells) { + MCCompLabel currentLab{mTimeFrame->getTrackletsLabel(cellTopology.firstTransition)[cell.getFirstTrackletIndex()]}; + MCCompLabel nextLab{mTimeFrame->getTrackletsLabel(cellTopology.secondTransition)[cell.getSecondTrackletIndex()]}; + labels.emplace_back(currentLab == nextLab ? currentLab : MCCompLabel()); + } + } + } + }); + + for (int transitionId = 0; transitionId < topology.nTransitions; ++transitionId) { + deepVectorClear(mTimeFrame->getTracklets()[transitionId]); + deepVectorClear(mTimeFrame->getTrackletsLabel(transitionId)); + } +} + +template +void TrackerTraits::findCellsNeighbours(const int iteration) +{ + const auto topology = mTimeFrame->getTrackingTopologyView(); + mTaskArena->execute([&] { + std::vector> cellsNeighboursByTarget; + cellsNeighboursByTarget.reserve(topology.nCells); + for (int cellTopologyId{0}; cellTopologyId < topology.nCells; ++cellTopologyId) { + deepVectorClear(mTimeFrame->getCellsNeighbours()[cellTopologyId]); + deepVectorClear(mTimeFrame->getCellsNeighboursTopology()[cellTopologyId]); + deepVectorClear(mTimeFrame->getCellsNeighboursLUT()[cellTopologyId]); + cellsNeighboursByTarget.emplace_back(mMemoryPool.get()); + } + + for (int outerLayer{0}; outerLayer < NLayers; ++outerLayer) { + for (int cellTopologyId{0}; cellTopologyId < topology.nCells; ++cellTopologyId) { + const auto& cellTopology = topology.getCell(cellTopologyId); + if (cellTopology.hitLayerMask.last() != outerLayer || + mTimeFrame->getCells()[cellTopologyId].empty()) { + continue; + } + const auto successors = topology.getCellsStartingWithTransition(cellTopology.secondTransition); + if (!successors.getEntries()) { + continue; + } + + tbb::enumerable_thread_specific> sourceNeighbours([&]() { return bounded_vector{mMemoryPool.get()}; }); + tbb::parallel_for(0, static_cast(mTimeFrame->getCells()[cellTopologyId].size()), [&](const int iCell) { + auto& localNeighbours = sourceNeighbours.local(); + const auto& currentCellSeed{mTimeFrame->getCells()[cellTopologyId][iCell]}; + const int nextLayerTrackletIndex{currentCellSeed.getSecondTrackletIndex()}; + for (int iSuccessor{0}; iSuccessor < successors.getEntries(); ++iSuccessor) { + const int nextCellTopologyId = topology.cellsByFirstTransition[successors.getFirstEntry() + iSuccessor]; + if (mTimeFrame->getCells()[nextCellTopologyId].empty() || + mTimeFrame->getCellsLookupTable()[nextCellTopologyId].empty()) { + continue; + } + const auto& nextCellLUT = mTimeFrame->getCellsLookupTable()[nextCellTopologyId]; + if (nextLayerTrackletIndex + 1 >= static_cast(nextCellLUT.size())) { + continue; + } + const int nextLayerFirstCellIndex{nextCellLUT[nextLayerTrackletIndex]}; + const int nextLayerLastCellIndex{nextCellLUT[nextLayerTrackletIndex + 1]}; + for (int iNextCell{nextLayerFirstCellIndex}; iNextCell < nextLayerLastCellIndex; ++iNextCell) { + const auto& nextCellSeedRef{mTimeFrame->getCells()[nextCellTopologyId][iNextCell]}; + if (nextCellSeedRef.getFirstTrackletIndex() != nextLayerTrackletIndex || !currentCellSeed.getTimeStamp().isCompatible(nextCellSeedRef.getTimeStamp())) { + break; + } + + bool neighbourAccepted{false}; + auto nextCellSeed{mTimeFrame->getCells()[nextCellTopologyId][iNextCell]}; /// copy + if constexpr (DetectorTraits::DetId == o2::detectors::DetID::MFT) { + neighbourAccepted = DetectorTraits::cellsAreCompatible(currentCellSeed, + nextCellSeed, + *mTimeFrame, + mTrkParams[iteration], + getBz()); + } else { + if (!nextCellSeed.rotate(currentCellSeed.getAlpha()) || + !nextCellSeed.propagateTo(currentCellSeed.getX(), getBz())) { + continue; + } + + const float chi2 = currentCellSeed.getPredictedChi2(nextCellSeed); + if (chi2 > mTrkParams[iteration].MaxChi2ClusterAttachment) { + continue; + } + neighbourAccepted = true; + } + if (!neighbourAccepted) { + continue; + } + + const int nextLevel = currentCellSeed.getLevel() + 1; + localNeighbours.emplace_back(cellTopologyId, iCell, nextCellTopologyId, iNextCell, nextLevel); + } + } + }); + + bounded_vector count(topology.nCells, 0, mMemoryPool.get()); + for (const auto& localNeighbours : sourceNeighbours) { + for (const auto& neigh : localNeighbours) { + ++count[neigh.nextCellTopology]; + } + } + for (size_t i{0}; i < topology.nCells; ++i) { + cellsNeighboursByTarget[i].reserve(count[i]); + } + for (const auto& localNeighbours : sourceNeighbours) { + for (const auto& neigh : localNeighbours) { + cellsNeighboursByTarget[neigh.nextCellTopology].emplace_back(neigh); + if (neigh.level > mTimeFrame->getCells()[neigh.nextCellTopology][neigh.nextCell].getLevel()) { + mTimeFrame->getCells()[neigh.nextCellTopology][neigh.nextCell].setLevel(neigh.level); + } + } + } + } + } + + for (int cellTopologyId{0}; cellTopologyId < topology.nCells; ++cellTopologyId) { + auto& cellsNeighbours = cellsNeighboursByTarget[cellTopologyId]; + if (cellsNeighbours.empty()) { + continue; + } + + std::sort(cellsNeighbours.begin(), cellsNeighbours.end(), [](const auto& a, const auto& b) { + return a.nextCell < b.nextCell; + }); + + auto& cellsNeighbourLUT = mTimeFrame->getCellsNeighboursLUT()[cellTopologyId]; + cellsNeighbourLUT.assign(mTimeFrame->getCells()[cellTopologyId].size(), 0); + for (const auto& neigh : cellsNeighbours) { + ++cellsNeighbourLUT[neigh.nextCell]; + } + std::inclusive_scan(cellsNeighbourLUT.begin(), cellsNeighbourLUT.end(), cellsNeighbourLUT.begin()); + + mTimeFrame->getCellsNeighbours()[cellTopologyId].reserve(cellsNeighbours.size()); + mTimeFrame->getCellsNeighboursTopology()[cellTopologyId].reserve(cellsNeighbours.size()); + std::ranges::transform(cellsNeighbours, std::back_inserter(mTimeFrame->getCellsNeighbours()[cellTopologyId]), [](const auto& neigh) { return neigh.cell; }); + std::ranges::transform(cellsNeighbours, std::back_inserter(mTimeFrame->getCellsNeighboursTopology()[cellTopologyId]), [](const auto& neigh) { return neigh.cellTopology; }); + } + + // clean up LUTs + for (auto& cellLUT : mTimeFrame->getCellsLookupTable()) { + deepVectorClear(cellLUT); + } + }); +} + +template +template +void TrackerTraits::processNeighbours(int iteration, int defaultCellTopologyId, int iLevel, const bounded_vector& currentCellSeed, const bounded_vector& currentCellId, const bounded_vector& currentCellTopologyId, bounded_vector& updatedCellSeeds, bounded_vector& updatedCellsIds, bounded_vector& updatedCellsTopologyIds) +{ + auto propagator = o2::base::Propagator::Instance(); + + mTaskArena->execute([&] { + auto forCellNeighbours = [&](auto Tag, int iCell, int offset = 0) -> int { + const auto& currentCell{currentCellSeed[iCell]}; + const int cellTopologyId = currentCellTopologyId.empty() ? defaultCellTopologyId : currentCellTopologyId[iCell]; + + if constexpr (decltype(Tag)::value != PassMode::TwoPassInsert::value) { + if (currentCell.getLevel() != iLevel) { + return 0; + } + if (currentCellId.empty()) { + for (int layer = 0; layer < NLayers; ++layer) { + const int clusterIndex = currentCell.getCluster(layer); + if (clusterIndex != constants::UnusedIndex && mTimeFrame->isClusterUsed(layer, clusterIndex)) { + return 0; /// this we do only on the first iteration, hence the check on currentCellId + } + } + } + } + + const int cellId = currentCellId.empty() ? iCell : currentCellId[iCell]; + if (cellTopologyId < 0 || mTimeFrame->getCellsNeighboursLUT()[cellTopologyId].empty()) { + return 0; + } + const int startNeighbourId{cellId ? mTimeFrame->getCellsNeighboursLUT()[cellTopologyId][cellId - 1] : 0}; + const int endNeighbourId{mTimeFrame->getCellsNeighboursLUT()[cellTopologyId][cellId]}; + int foundSeeds{0}; + for (int iNeighbourCell{startNeighbourId}; iNeighbourCell < endNeighbourId; ++iNeighbourCell) { + const int neighbourCellTopologyId = mTimeFrame->getCellsNeighboursTopology()[cellTopologyId][iNeighbourCell]; + const int neighbourCellId = mTimeFrame->getCellsNeighbours()[cellTopologyId][iNeighbourCell]; + const auto& neighbourCell = mTimeFrame->getCells()[neighbourCellTopologyId][neighbourCellId]; + if (neighbourCell.getSecondTrackletIndex() != currentCell.getFirstTrackletIndex()) { + continue; + } + if (!currentCell.getTimeStamp().isCompatible(neighbourCell.getTimeStamp())) { + continue; + } + if (currentCell.getLevel() - 1 != neighbourCell.getLevel()) { + continue; + } + const int neighbourLayer = neighbourCell.getInnerLayer(); + const int neighbourCluster = neighbourCell.getFirstClusterIndex(); + if (mTimeFrame->isClusterUsed(neighbourLayer, neighbourCluster)) { + continue; + } + + /// Let's start the fitting procedure + TrackSeedN seed{currentCell}; + seed.getTimeStamp() = currentCell.getTimeStamp(); + seed.getTimeStamp() += neighbourCell.getTimeStamp(); + + if constexpr (DetectorTraits::DetId == o2::detectors::DetID::MFT) { + o2::track::TrackParCovFwd fwdTrack; + float chi2 = seed.getChi2(); + if (!detail::mftFwdAttachNeighbourToSeed(seed, neighbourLayer, neighbourCluster, *mTimeFrame, mTrkParams[iteration], getBz(), fwdTrack, chi2)) { + continue; + } + static_cast(seed) = fwdTrack; + seed.setChi2(chi2); + } else { + const auto& trHit = mTimeFrame->getTrackingFrameInfoOnLayer(neighbourLayer)[neighbourCluster]; + + if (!seed.rotate(trHit.alphaTrackingFrame)) { + continue; + } + + if (!propagator->propagateToX(seed, trHit.xTrackingFrame, getBz(), o2::base::PropagatorImpl::MAX_SIN_PHI, o2::base::PropagatorImpl::MAX_STEP, mTrkParams[iteration].CorrType)) { + continue; + } + + if (mTrkParams[iteration].CorrType == o2::base::PropagatorF::MatCorrType::USEMatCorrNONE) { + if (!seed.correctForMaterial(mTrkParams[iteration].LayerxX0[neighbourLayer], mTrkParams[iteration].LayerxX0[neighbourLayer] * o2::its::constants::Radl * o2::its::constants::Rho, true)) { + continue; + } + } + + auto predChi2{seed.getPredictedChi2Quiet(trHit.positionTrackingFrame, trHit.covarianceTrackingFrame)}; + if ((predChi2 > mTrkParams[iteration].MaxChi2ClusterAttachment) || predChi2 < 0.f) { + continue; + } + seed.setChi2(seed.getChi2() + predChi2); + if (!seed.o2::track::TrackParCov::update(trHit.positionTrackingFrame, trHit.covarianceTrackingFrame)) { + continue; + } + } + + if constexpr (decltype(Tag)::value != PassMode::TwoPassCount::value) { + seed.getClusters()[neighbourLayer] = neighbourCluster; + auto mask = seed.getHitLayerMask(); + mask.set(neighbourLayer); + seed.setHitLayerMask(mask); + seed.setLevel(neighbourCell.getLevel()); + seed.setFirstTrackletIndex(neighbourCell.getFirstTrackletIndex()); + seed.setSecondTrackletIndex(neighbourCell.getSecondTrackletIndex()); + } + + if constexpr (decltype(Tag)::value == PassMode::OnePass::value) { + updatedCellSeeds.push_back(seed); + updatedCellsIds.push_back(neighbourCellId); + updatedCellsTopologyIds.push_back(neighbourCellTopologyId); + } else if constexpr (decltype(Tag)::value == PassMode::TwoPassCount::value) { + ++foundSeeds; + } else if constexpr (decltype(Tag)::value == PassMode::TwoPassInsert::value) { + updatedCellSeeds[offset] = seed; + updatedCellsIds[offset] = neighbourCellId; + updatedCellsTopologyIds[offset++] = neighbourCellTopologyId; + } else { + static_assert(false, "Unknown mode!"); + } + } + return foundSeeds; + }; + + const int nCells = static_cast(currentCellSeed.size()); + if (mTaskArena->max_concurrency() <= 1) { + for (int iCell{0}; iCell < nCells; ++iCell) { + forCellNeighbours(PassMode::OnePass{}, iCell); + } + } else { + bounded_vector perCellCount(nCells + 1, 0, mMemoryPool.get()); + tbb::parallel_for(0, nCells, [&](const int iCell) { + perCellCount[iCell] = forCellNeighbours(PassMode::TwoPassCount{}, iCell); + }); + + std::exclusive_scan(perCellCount.begin(), perCellCount.end(), perCellCount.begin(), 0); + auto totalNeighbours{perCellCount.back()}; + if (totalNeighbours == 0) { + return; + } + updatedCellSeeds.resize(totalNeighbours); + updatedCellsIds.resize(totalNeighbours); + updatedCellsTopologyIds.resize(totalNeighbours); + + tbb::parallel_for(0, nCells, [&](const int iCell) { + int offset = perCellCount[iCell]; + if (offset == perCellCount[iCell + 1]) { + return; + } + forCellNeighbours(PassMode::TwoPassInsert{}, iCell, offset); + }); + } + }); +} + +template +void TrackerTraits::findRoads(const int iteration) +{ + bounded_vector> firstClusters(mTrkParams[iteration].NLayers, bounded_vector(mMemoryPool.get()), mMemoryPool.get()); + firstClusters.resize(mTrkParams[iteration].NLayers); + const auto propagator = o2::base::Propagator::Instance(); + const TrackingFrameInfo* tfInfos[NLayers]{}; + const Cluster* unsortedClusters[NLayers]{}; + for (int iLayer = 0; iLayer < NLayers; ++iLayer) { + tfInfos[iLayer] = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer).data(); + unsortedClusters[iLayer] = mTimeFrame->getUnsortedClusters()[iLayer].data(); + } + const auto topology = mTimeFrame->getTrackingTopologyView(); + for (int startLevel{mTrkParams[iteration].CellsPerRoad()}; startLevel >= mTrkParams[iteration].CellMinimumLevel(); --startLevel) { + + auto seedFilter = [&](const auto& seed) { + return seed.getHitLayerMask().isAllowed(mTrkParams[iteration].MaxHoles, mTrkParams[iteration].HoleLayerMask) && + seed.getHitLayerMask().length() >= mTrkParams[iteration].MinTrackLength && + seed.getQ2Pt() <= 1.e3 && seed.getChi2() <= mTrkParams[iteration].MaxChi2NDF * ((startLevel + 2) * 2 - 5); + }; + + bounded_vector trackSeeds(mMemoryPool.get()); + for (int startCellTopologyId{0}; startCellTopologyId < topology.nCells; ++startCellTopologyId) { + const int startLayer = topology.getCell(startCellTopologyId).hitLayerMask.last(); + if (!(mTrkParams[iteration].StartLayerMask.has(startLayer)) || mTimeFrame->getCells()[startCellTopologyId].empty()) { + continue; + } + + bounded_vector lastCellId(mMemoryPool.get()), updatedCellId(mMemoryPool.get()); + bounded_vector lastCellTopologyId(mMemoryPool.get()), updatedCellTopologyId(mMemoryPool.get()); + bounded_vector lastCellSeed(mMemoryPool.get()), updatedCellSeed(mMemoryPool.get()); + + processNeighbours(iteration, startCellTopologyId, startLevel, mTimeFrame->getCells()[startCellTopologyId], lastCellId, lastCellTopologyId, updatedCellSeed, updatedCellId, updatedCellTopologyId); + + int level = startLevel; + while (level > 2 && !updatedCellSeed.empty()) { + lastCellSeed.swap(updatedCellSeed); + lastCellId.swap(updatedCellId); + lastCellTopologyId.swap(updatedCellTopologyId); + deepVectorClear(updatedCellSeed); /// tame the memory peaks + deepVectorClear(updatedCellId); /// tame the memory peaks + deepVectorClear(updatedCellTopologyId); + processNeighbours(iteration, constants::UnusedIndex, --level, lastCellSeed, lastCellId, lastCellTopologyId, updatedCellSeed, updatedCellId, updatedCellTopologyId); + } + deepVectorClear(lastCellId); /// tame the memory peaks + deepVectorClear(lastCellTopologyId); /// tame the memory peaks + deepVectorClear(lastCellSeed); /// tame the memory peaks + + if (!updatedCellSeed.empty()) { + trackSeeds.reserve(trackSeeds.size() + std::count_if(updatedCellSeed.begin(), updatedCellSeed.end(), seedFilter)); + std::copy_if(updatedCellSeed.begin(), updatedCellSeed.end(), std::back_inserter(trackSeeds), seedFilter); + } + } + + if (trackSeeds.empty()) { + continue; + } + + using TrackT = typename DetectorTraits::TrackType; + bounded_vector tracks(mMemoryPool.get()); + mTaskArena->execute([&] { + auto forSeed = [&](auto Tag, int iSeed, int offset = 0) { + TrackT temporaryTrack; + const bool refitSuccess = DetectorTraits::refitSeed(trackSeeds[iSeed], + temporaryTrack, + mTrkParams[iteration], + mBz, + *mTimeFrame, + tfInfos, + unsortedClusters, + propagator); + if (refitSuccess) { + if constexpr (DetectorTraits::DetId == o2::detectors::DetID::MFT) { + temporaryTrack.setSeedPattern(trackSeeds[iSeed].getHitLayerMask().value()); + } + if constexpr (decltype(Tag)::value == PassMode::OnePass::value) { + tracks.push_back(temporaryTrack); + } else if constexpr (decltype(Tag)::value == PassMode::TwoPassCount::value) { + // nothing to do + } else if constexpr (decltype(Tag)::value == PassMode::TwoPassInsert::value) { + tracks[offset] = temporaryTrack; + } else { + static_assert(false, "Unknown mode!"); + } + return 1; + } + return 0; + }; + + const int nSeeds = static_cast(trackSeeds.size()); + if (mTaskArena->max_concurrency() <= 1) { + for (int iSeed{0}; iSeed < nSeeds; ++iSeed) { + forSeed(PassMode::OnePass{}, iSeed); + } + } else { + bounded_vector perSeedCount(nSeeds + 1, 0, mMemoryPool.get()); + tbb::parallel_for(0, nSeeds, [&](const int iSeed) { + perSeedCount[iSeed] = forSeed(PassMode::TwoPassCount{}, iSeed); + }); + + std::exclusive_scan(perSeedCount.begin(), perSeedCount.end(), perSeedCount.begin(), 0); + auto totalTracks{perSeedCount.back()}; + if (totalTracks == 0) { + return; + } + tracks.resize(totalTracks); + + tbb::parallel_for(0, nSeeds, [&](const int iSeed) { + if (perSeedCount[iSeed] == perSeedCount[iSeed + 1]) { + return; + } + forSeed(PassMode::TwoPassInsert{}, iSeed, perSeedCount[iSeed]); + }); + } + + deepVectorClear(trackSeeds); + }); + + // Same ordering as o2::its::track::isBetter (longer track, then lower chi2). + std::sort(tracks.begin(), tracks.end(), [](const TrackT& a, const TrackT& b) { + const auto ncla = a.getNumberOfClusters(); + const auto nclb = b.getNumberOfClusters(); + return (ncla == nclb) ? (a.getChi2() < b.getChi2()) : ncla > nclb; + }); + acceptTracks(iteration, tracks, firstClusters); + } + markTracks(iteration); +} + +template +void TrackerTraits::acceptTracks(int iteration, bounded_vector>& tracks, bounded_vector>& firstClusters) +{ + auto& trks = mTimeFrame->getTracks(); + trks.reserve(trks.size() + tracks.size()); + const float smallestROFHalf = mTimeFrame->getROFOverlapTableView().getClockLayer().mROFLength * 0.5f; + for (auto& track : tracks) { + int nShared = 0; + bool isFirstShared{false}; + int firstLayer{-1}, firstCluster{-1}; + for (int iLayer{0}; iLayer < mTrkParams[iteration].NLayers; ++iLayer) { + if (track.getClusterIndex(iLayer) == constants::UnusedIndex) { + continue; + } + bool isShared = mTimeFrame->isClusterUsed(iLayer, track.getClusterIndex(iLayer)); + nShared += int(isShared); + if (firstLayer < 0) { + firstCluster = track.getClusterIndex(iLayer); + isFirstShared = isShared && mTrkParams[iteration].AllowSharingFirstCluster && std::find(firstClusters[iLayer].begin(), firstClusters[iLayer].end(), firstCluster) != firstClusters[iLayer].end(); + firstLayer = iLayer; + } + } + + /// do not account for the first cluster in the shared clusters number if it is allowed + if (nShared - int(isFirstShared && mTrkParams[iteration].AllowSharingFirstCluster) > mTrkParams[iteration].SharedMaxClusters) { + continue; + } + + bool firstCls{true}, nominalCompatible{true}; + TimeEstBC nominalTS, expandedTS; + for (int iLayer{0}; iLayer < mTrkParams[iteration].NLayers; ++iLayer) { + if (track.getClusterIndex(iLayer) == constants::UnusedIndex) { + continue; + } + mTimeFrame->markUsedCluster(iLayer, track.getClusterIndex(iLayer)); + int currentROF = mTimeFrame->getClusterROF(iLayer, track.getClusterIndex(iLayer)); + const auto nominalROFTS = mTimeFrame->getROFOverlapTableView().getLayer(iLayer).getROFTimeBounds(currentROF); + const auto expandedROFTS = mTimeFrame->getROFOverlapTableView().getLayer(iLayer).getROFTimeBounds(currentROF, true); + if (firstCls) { + firstCls = false; + nominalTS = nominalROFTS; + expandedTS = expandedROFTS; + } else { + if (nominalCompatible) { + if (nominalTS.isCompatible(nominalROFTS)) { + nominalTS += nominalROFTS; + } else { + nominalCompatible = false; + } + } + if (!expandedTS.isCompatible(expandedROFTS)) { + LOGP(fatal, "TS {}+/-{} are incompatible with {}+/-{}, this should not happen!", expandedROFTS.getTimeStamp(), expandedROFTS.getTimeStampError(), expandedTS.getTimeStamp(), expandedTS.getTimeStampError()); + } + expandedTS += expandedROFTS; + } + } + track.getTimeStamp() = (nominalCompatible ? nominalTS : expandedTS).makeSymmetrical(); + // this is a sanity clamp + // we cannot be worse than the clock so we clamp to this + if (track.getTimeStamp().getTimeStampError() > smallestROFHalf) { + track.getTimeStamp().setTimeStampError(smallestROFHalf); + } + if constexpr (DetectorTraits::DetId == o2::detectors::DetID::ITS) { + track.setUserField(0); + track.getParamOut().setUserField(0); + } + trks.emplace_back(track); + + if (mTrkParams[iteration].AllowSharingFirstCluster) { + firstClusters[firstLayer].push_back(firstCluster); + } + } +} + +template +void TrackerTraits::markTracks(int iteration) +{ + if (mTrkParams[iteration].AllowSharingFirstCluster) { + /// Now we have to set the shared cluster flag + auto& tracks = mTimeFrame->getTracks(); + + bounded_vector fclusSort(tracks.size(), mMemoryPool.get()); + std::iota(fclusSort.begin(), fclusSort.end(), 0); + std::sort(fclusSort.begin(), fclusSort.end(), [&tracks](int a, int b) { + return tracks[a].getClusterIndex(tracks[a].getFirstClusterLayer()) < tracks[b].getClusterIndex(tracks[b].getFirstClusterLayer()); + }); + + auto areTracksSelected = [this, iteration](const auto& t1, const auto& t2) { + const auto t1FirstLayer{t1.getFirstClusterLayer()}, t2FirstLayer{t2.getFirstClusterLayer()}; + if (t1FirstLayer != t2FirstLayer) { + return false; + } + if (mTimeFrame->getClusterROF(t1FirstLayer, t1.getClusterIndex(t1FirstLayer)) != mTimeFrame->getClusterROF(t2FirstLayer, t2.getClusterIndex(t2FirstLayer))) { + return false; + } + if (!math_utils::isPhiDifferenceBelow(t1.getPhi(), t2.getPhi(), mTrkParams[iteration].SharedClusterMaxDeltaPhi)) { + return false; + } + if (std::abs(t1.getEta() - t2.getEta()) > mTrkParams[iteration].SharedClusterMaxDeltaEta) { + return false; + } + if (mTrkParams[iteration].SharedClusterOppositeSign) { + if constexpr (DetectorTraits::DetId == o2::detectors::DetID::MFT) { + if (t1.getCharge() == t2.getCharge()) { + return false; + } + } else { + if (t1.getSign() == t2.getSign()) { + return false; + } + } + } + return true; + }; + + for (int i{0}; i < static_cast(fclusSort.size()); ++i) { + auto& track = tracks[fclusSort[i]]; + for (int j{i + 1}; j < static_cast(fclusSort.size()) && tracks[fclusSort[j]].getClusterIndex(tracks[fclusSort[j]].getFirstClusterLayer()) == track.getClusterIndex(track.getFirstClusterLayer()); ++j) { + auto& track2 = tracks[fclusSort[j]]; + if (areTracksSelected(track, track2)) { + track.setSharedClusters(); + track2.setSharedClusters(); + } + } + } + } +} + +template +void TrackerTraits::setBz(float bz) +{ + mBz = bz; + mTimeFrame->setBz(bz); +} + +template +void TrackerTraits::setNThreads(int n, std::shared_ptr& arena) +{ +#if defined(OPTIMISATION_OUTPUT) + mTaskArena = std::make_shared(1); +#else + if (arena == nullptr) { + mTaskArena = std::make_shared(std::abs(n)); + LOGP(info, "Setting tracker with {} threads.", n); + } else { + mTaskArena = arena; + } +#endif +} + +template class TrackerTraits<7>; +template class TrackerTraits<10>; +template void TrackerTraits<7>::processNeighbours::CellSeedN>(int, int, int, const bounded_vector::CellSeedN>&, const bounded_vector&, const bounded_vector&, bounded_vector::TrackSeedN>&, bounded_vector&, bounded_vector&); +template void TrackerTraits<7>::processNeighbours::TrackSeedN>(int, int, int, const bounded_vector::TrackSeedN>&, const bounded_vector&, const bounded_vector&, bounded_vector::TrackSeedN>&, bounded_vector&, bounded_vector&); +template void TrackerTraits<10>::processNeighbours::CellSeedN>(int, int, int, const bounded_vector::CellSeedN>&, const bounded_vector&, const bounded_vector&, bounded_vector::TrackSeedN>&, bounded_vector&, bounded_vector&); +template void TrackerTraits<10>::processNeighbours::TrackSeedN>(int, int, int, const bounded_vector::TrackSeedN>&, const bounded_vector&, const bounded_vector&, bounded_vector::TrackSeedN>&, bounded_vector&, bounded_vector&); + +} // namespace o2::itsmft::tracking diff --git a/Detectors/ITSMFT/common/tracking/src/TrackingConfigParam.cxx b/Detectors/ITSMFT/common/tracking/src/TrackingConfigParam.cxx new file mode 100644 index 0000000000000..9d68c2719a439 --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/src/TrackingConfigParam.cxx @@ -0,0 +1,19 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "ITSMFTTracking/TrackingConfigParam.h" + +namespace o2::itsmft +{ +// Ensure MFT CA params are registered in the global param database. +// ITS production tracking uses o2::its::TrackerParamConfig in O2::ITStracking. +static auto& sMFTCATrackerParam = TrackerParamConfig::Instance(); +} // namespace o2::itsmft diff --git a/Detectors/ITSMFT/common/tracking/src/TrackingInterface.cxx b/Detectors/ITSMFT/common/tracking/src/TrackingInterface.cxx new file mode 100644 index 0000000000000..77dbbe052350d --- /dev/null +++ b/Detectors/ITSMFT/common/tracking/src/TrackingInterface.cxx @@ -0,0 +1,330 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file TrackingInterface.cxx +/// \brief +/// + +#include "ITSMFTTracking/DetectorTraits.h" +#include "ITSMFTTracking/TrackingInterface.h" + +#include + +#include "ITStracking/Constants.h" + +#include + +#include "CommonConstants/LHCConstants.h" +#include "CommonDataFormat/IRFrame.h" +#include "CommonDataFormat/InteractionRecord.h" +#include "DataFormatsITSMFT/CompCluster.h" +#include "DataFormatsITSMFT/DPLAlpideParam.h" +#include "DataFormatsITSMFT/ROFRecord.h" +#include "DetectorsBase/GRPGeomHelper.h" +#include "DetectorsBase/Propagator.h" +#include "Framework/Logger.h" +#include "MFTTracking/MFTTrackingParam.h" + +namespace o2::itsmft::tracking +{ + +namespace +{ +template +consteval const char* detName() +{ + return detId == o2::detectors::DetID::MFT ? "MFT" : "ITS"; +} + +bool rofOverlapsIRFrames(const o2::itsmft::ROFRecord& rof, int rofLengthInBC, gsl::span irFrames) +{ + o2::InteractionRecord rofStart{rof.getBCData()}; + o2::InteractionRecord rofEnd = rofStart + rofLengthInBC - 1; + o2::dataformats::IRFrame ref(rofStart, rofEnd); + for (const auto& ir : irFrames) { + if (ir.info > 0) { + const auto overlap = ref.getOverlap(ir); + if (overlap.isValid()) { + return true; + } + } + } + return false; +} +} // namespace + +template +ITSMFTTrackingInterface::ITSMFTTrackingInterface(bool useMC, + o2::itsmft::TrackingMode::Type mode, + bool overrideBeamEst) + : mUseMC(useMC), mTrackingMode(mode), mOverrideBeamEstimation(overrideBeamEst) +{ +} + +template +void ITSMFTTrackingInterface::initialise() +{ + resolveTrackingParameters(); + initialiseMemoryPool(); + initialiseTracker(); +} + +template +void ITSMFTTrackingInterface::resolveTrackingParameters() +{ + const auto& tc = o2::itsmft::tracking::TrackerParamRef::get(); + auto mode = mTrackingMode; + if (auto parmode = static_cast(tc.trackingMode); + mode == o2::itsmft::TrackingMode::Unset || + (parmode != o2::itsmft::TrackingMode::Unset && mode != parmode)) { + if (parmode != o2::itsmft::TrackingMode::Unset) { + LOGP(info, "{} CA tracking mode overwritten by configurable params from {} to {}", + detName(), o2::itsmft::TrackingMode::toString(mode), o2::itsmft::TrackingMode::toString(parmode)); + mode = parmode; + } + } + if (mode == o2::itsmft::TrackingMode::Unset) { + mode = o2::itsmft::TrackingMode::Sync; + } + mTrackingMode = mode; + mTrackParams = o2::itsmft::TrackingMode::getTrackingParameters(DetId, mode); + LOGP(info, "{} CA tracker initialized in {} mode with {} iteration(s)", + detName(), o2::itsmft::TrackingMode::toString(mode), mTrackParams.size()); + for (size_t i = 0; i < mTrackParams.size(); ++i) { + LOGP(info, " iteration {}: {}", i, mTrackParams[i].asString()); + } +} + +template +void ITSMFTTrackingInterface::initialiseMemoryPool() +{ + size_t maxMemory = std::numeric_limits::max(); + if (!mTrackParams.empty() && mTrackParams[0].MaxMemory != maxMemory) { + maxMemory = mTrackParams[0].MaxMemory; + } + mMemoryPool = std::make_shared(maxMemory); + mTimeFrame.setMemoryPool(mMemoryPool); +} + +template +void ITSMFTTrackingInterface::initialiseTracker() +{ + if (mTrackParams.empty()) { + return; + } + const auto& tc = o2::itsmft::tracking::TrackerParamRef::get(); + mTrackerTraits = std::make_unique(); + mTrackerTraits->setMemoryPool(mMemoryPool); + std::shared_ptr taskArena; + mTrackerTraits->setNThreads(tc.nThreads, taskArena); + mTracker = std::make_unique(mTrackerTraits.get()); +} + +template +float ITSMFTTrackingInterface::processTimeFrame(gsl::span rofs, + gsl::span clusters, + gsl::span patterns, + const o2::dataformats::MCTruthContainer* labels, + gsl::span irFrames) +{ + if (mTrackParams.empty()) { + LOGP(info, "{} CA tracking mode is off, skipping TimeFrame processing", detName()); + return 0.f; + } + if (mDict == nullptr) { + LOGP(fatal, "{} CA tracker cluster dictionary is not available", detName()); + } + + loadTimeFrame(rofs, clusters, patterns, labels, irFrames); + onTimeFrameLoaded(); + + const float elapsedMs = runTracking(); + onTrackingFinished(elapsedMs); + + return elapsedMs; +} + +template +void ITSMFTTrackingInterface::loadTimeFrame(gsl::span rofs, + gsl::span clusters, + gsl::span patterns, + const o2::dataformats::MCTruthContainer* labels, + gsl::span irFrames) +{ + configureROFLookupTables(); + validateROFInput(rofs); + + mTimeFrame.setBz(o2::base::Propagator::Instance()->getNominalBz()); + configureBeamPosition(); + + configureROFMask(rofs, irFrames); + + auto pattIt = patterns.begin(); + mTimeFrame.loadROFrameData(rofs, clusters, pattIt, mDict, -1, labels, DetId); + + configureTrackingTopology(); + + LOGP(info, "{} CA loaded {} clusters from {} ROFs into TimeFrame ({} pattern bytes, MC={})", + detName(), mTimeFrame.getTotalClusters(), rofs.size(), patterns.size(), labels != nullptr); + + for (int iLayer = 0; iLayer < NLayers; ++iLayer) { + LOGP(info, " layer {}: {} ROF slots", iLayer, mTimeFrame.getNrof(iLayer)); + } +} + +template +float ITSMFTTrackingInterface::runTracking() +{ + if (!mTracker || mTrackParams.empty()) { + return 0.f; + } + mTracker->adoptTimeFrame(mTimeFrame); + mTracker->setParameters(mTrackParams); + mTracker->setMemoryPool(mMemoryPool); + mTracker->setBz(mTimeFrame.getBz()); + const float elapsedMs = mTracker->clustersToTracks(); + if (elapsedMs < 0.f) { + LOGP(warn, "{} CA tracking failed for this TF", detName()); + return elapsedMs; + } + LOGP(info, "{} CA tracking produced {} tracks in {:.2f} ms", + detName(), mTimeFrame.getNumberOfTracks(), elapsedMs); + return elapsedMs; +} + +template +void ITSMFTTrackingInterface::configureTrackingTopology() +{ + if (mTrackParams.empty()) { + return; + } + mTimeFrame.initDefaultTrackingTopology(mTrackParams[0], NLayers); + mTimeFrame.initTrackerTopologies(mTrackParams); +} + +template +void ITSMFTTrackingInterface::configureBeamPosition() +{ + if (mTrackParams.empty()) { + return; + } + const auto& p = mTrackParams[0]; + TrackingLoadPolicyN::configureBeamPosition(mTimeFrame, p, mMeanVertex, mOverrideBeamEstimation); +} + +template +void ITSMFTTrackingInterface::configureROFLookupTables() +{ + if constexpr (DetId == o2::detectors::DetID::MFT) { + const bool continuous = o2::base::GRPGeomHelper::instance().getGRPECS()->isDetContinuousReadOut(DetId); + LOGP(info, "{} CA tracker RO: continuous = {}", detName(), continuous); + mMFTTriggered = !continuous; + const auto& alpParams = o2::itsmft::DPLAlpideParam::Instance(); + if (mMFTTriggered) { + mMFTROFrameLengthInBC = std::max(1, static_cast(alpParams.roFrameLengthTrig / (o2::constants::lhc::LHCBunchSpacingNS * 1e3))); + } else { + mMFTROFrameLengthInBC = alpParams.roFrameLengthInBC; + } + } + + const auto& par = o2::itsmft::DPLAlpideParam::Instance(); + const int nOrbitsPerTF = o2::base::GRPGeomHelper::getNHBFPerTF(); + + ROFOverlapTableN rofTable; + ROFVertexLookupTableN vtxTable; + for (int iLayer = 0; iLayer < NLayers; ++iLayer) { + const unsigned int nROFsPerOrbit = o2::constants::lhc::LHCMaxBunches / par.getROFLengthInBC(iLayer); + const o2::its::LayerTiming timing{ + .mNROFsTF = nROFsPerOrbit * static_cast(nOrbitsPerTF), + .mROFLength = static_cast(par.getROFLengthInBC(iLayer)), + .mROFDelay = static_cast(par.getROFDelayInBC(iLayer)), + .mROFBias = static_cast(par.getROFBiasInBC(iLayer)), + .mROFAddTimeErr = mTrackParams.empty() + ? static_cast(o2::itsmft::tracking::TrackerParamRef::get().addTimeError[iLayer]) + : mTrackParams[0].AddTimeError[iLayer]}; + rofTable.defineLayer(iLayer, timing); + vtxTable.defineLayer(iLayer, timing); + } + + const auto nROFsLayer0 = rofTable.getLayer(0).mNROFsTF; + for (int iLayer = 1; iLayer < NLayers; ++iLayer) { + if (rofTable.getLayer(iLayer).mNROFsTF != nROFsLayer0) { + LOGP(fatal, + "{} CA single CLUSTERSROF input requires identical mNROFsTF on all {} layers (layer 0: {}, layer {}: {})", + detName(), NLayers, nROFsLayer0, iLayer, rofTable.getLayer(iLayer).mNROFsTF); + } + } + + rofTable.init(); + mTimeFrame.setROFOverlapTable(std::move(rofTable)); + vtxTable.init(); + mTimeFrame.setROFVertexLookupTable(std::move(vtxTable)); +} + +template +void ITSMFTTrackingInterface::configureROFMask(gsl::span rofs, + gsl::span irFrames) +{ + ROFMaskTableN mask{mTimeFrame.getROFOverlapTable()}; + mask.resetMask(); + + if constexpr (DetId == o2::detectors::DetID::MFT) { + const auto& trackingParam = o2::mft::MFTTrackingParam::Instance(); + const bool useIrFilter = trackingParam.irFramesOnly && !irFrames.empty(); + + if (trackingParam.isMultCutRequested()) { + LOGP(info, "{} CA multiplicity filter: Min nClusters = {} ; Max nClusters = {}", + detName(), trackingParam.cutMultClusLow, trackingParam.cutMultClusHigh); + } + if (useIrFilter) { + LOGP(info, "{} CA IRFrame filter enabled with {} ITS IR frames", detName(), irFrames.size()); + } + + const auto nROFs = mTimeFrame.getROFOverlapTableView().getLayer(0).mNROFsTF; + for (int iRof = 0; iRof < static_cast(nROFs); ++iRof) { + bool accept = true; + if (iRof < static_cast(rofs.size())) { + if (useIrFilter) { + accept = rofOverlapsIRFrames(rofs[iRof], mMFTROFrameLengthInBC, irFrames); + } + if (accept && trackingParam.isMultCutRequested()) { + accept = trackingParam.isPassingMultCut(rofs[iRof].getNEntries()); + } + } + if (accept) { + for (int iLayer = 0; iLayer < NLayers; ++iLayer) { + mask.setROFEnabled(iLayer, iRof, 1); + } + } + } + } else { + for (int iLayer = 0; iLayer < NLayers; ++iLayer) { + mask.setROFsEnabled(iLayer, 0, mTimeFrame.getROFOverlapTableView().getLayer(iLayer).mNROFsTF, 1); + } + } + + mTimeFrame.setMultiplicityCutMask(std::move(mask)); +} + +template +void ITSMFTTrackingInterface::validateROFInput(gsl::span rofs) const +{ + const auto expectedROFsTF = mTimeFrame.getROFOverlapTableView().getLayer(0).mNROFsTF; + if (rofs.size() != expectedROFsTF) { + LOGP(warn, "{} CA ROF count differs from continuous timing expectation: received {} expected {}", + detName(), rofs.size(), expectedROFsTF); + } +} + +template class ITSMFTTrackingInterface<7>; +template class ITSMFTTrackingInterface<10>; + +} // namespace o2::itsmft::tracking