diff --git a/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt b/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt index a23682b085311..ffdbdf1990a32 100644 --- a/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt +++ b/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt @@ -128,3 +128,9 @@ o2_add_test_root_macro(CheckStaggering.C O2::DetectorsVertexing O2::ReconstructionDataFormats LABELS its COMPILE_ONLY) + +o2_add_test_root_macro(CheckSeeding.C + PUBLIC_LINK_LIBRARIES O2::DataFormatsITS + O2::SimulationDataFormat + O2::Steer + LABELS its COMPILE_ONLY) diff --git a/Detectors/ITSMFT/ITS/macros/test/CheckSeeding.C b/Detectors/ITSMFT/ITS/macros/test/CheckSeeding.C new file mode 100644 index 0000000000000..915f2dda75032 --- /dev/null +++ b/Detectors/ITSMFT/ITS/macros/test/CheckSeeding.C @@ -0,0 +1,706 @@ +// 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. + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "DataFormatsITS/Vertex.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTrack.h" +#include "SimulationDataFormat/O2DatabasePDG.h" +#include "Steer/MCKinematicsReader.h" +#endif + +constexpr const char* tracFile = "o2trac_its.root"; +constexpr const char* collContextFile = "collisioncontext.root"; + +namespace +{ +namespace fs = std::filesystem; + +constexpr float MinPt = 0.05f; +constexpr float MaxEta = 1.1f; +constexpr int NMultiplicityBins = 11; +constexpr std::array MultiplicityLabels{{"2", "3", "4-5", "6-8", "9-13", "14-21", "22-33", "34-52", "53-83", "84-128", "128+"}}; + +struct TruthInfo { + int multiplicity = 0; + float x = 0.f; + float y = 0.f; + float z = 0.f; +}; + +struct BestRecoInfo { + o2::its::Vertex vertex; + float purity = -1.f; +}; + +struct GaussianSummary { + bool valid = false; + double mean = 0.; + double sigma = 0.; +}; + +bool isTrueVertexLabel(const o2::MCCompLabel& label) +{ + return label.isValid() && !label.isFake() && label.getSourceID() == 0; +} + +bool isChargedPrimary(const o2::MCTrack& track) +{ + if (!track.isPrimary() || track.GetPt() < MinPt || std::abs(track.GetEta()) > MaxEta) { + return false; + } + auto* pdg = o2::O2DatabasePDG::Instance()->GetParticle(track.GetPdgCode()); + return pdg != nullptr && pdg->Charge() != 0.; +} + +bool isBetterReco(const o2::its::Vertex& candidate, float candidatePurity, const o2::its::Vertex& current, float currentPurity) +{ + if (candidatePurity != currentPurity) { + return candidatePurity > currentPurity; + } + if (candidate.getNContributors() != current.getNContributors()) { + return candidate.getNContributors() > current.getNContributors(); + } + return candidate.getChi2() < current.getChi2(); +} + +int getMultiplicityCategory(int multiplicity) +{ + if (multiplicity <= 2) { + return 1; + } + if (multiplicity <= 3) { + return 2; + } + if (multiplicity <= 5) { + return 3; + } + if (multiplicity <= 8) { + return 4; + } + if (multiplicity <= 13) { + return 5; + } + if (multiplicity <= 21) { + return 6; + } + if (multiplicity <= 33) { + return 7; + } + if (multiplicity <= 52) { + return 8; + } + if (multiplicity <= 83) { + return 9; + } + if (multiplicity <= 128) { + return 10; + } + return 11; +} + +void fillMultiplicityHistogram(TH1* hist, int multiplicity) +{ + hist->Fill(getMultiplicityCategory(multiplicity)); +} + +GaussianSummary fitGaussianCore(TH1* hist, const char* funcName) +{ + if (hist == nullptr || hist->GetEntries() < 20) { + return {}; + } + const auto rms = hist->GetRMS(); + if (!(rms > 0.)) { + return {}; + } + + TF1 fit(funcName, "gaus", hist->GetMean() - 2. * rms, hist->GetMean() + 2. * rms); + fit.SetParameters(hist->GetMaximum(), hist->GetMean(), rms); + hist->Fit(&fit, "Q0R"); + + const auto mean = fit.GetParameter(1); + const auto sigma = std::abs(fit.GetParameter(2)); + if (!(sigma > 0.)) { + return {}; + } + + fit.SetRange(mean - 2. * sigma, mean + 2. * sigma); + hist->Fit(&fit, "Q0R"); + return {true, fit.GetParameter(1), std::abs(fit.GetParameter(2))}; +} + +TH1D* makeNormalizedCopy(const TH1D* source, const char* name, const char* title) +{ + auto* copy = static_cast(source->Clone(name)); + copy->SetTitle(title); + const auto integral = copy->Integral("width"); + if (integral > 0.) { + copy->Scale(1. / integral); + } + return copy; +} + +void setMultiplicityBinLabels(TH1* hist) +{ + for (int i = 0; i < NMultiplicityBins; ++i) { + hist->GetXaxis()->SetBinLabel(i + 1, MultiplicityLabels[i]); + } +} + +void setHistogramStyle(TH1* hist, int color, int marker) +{ + hist->SetLineColor(color); + hist->SetMarkerColor(color); + hist->SetMarkerStyle(marker); + hist->SetLineWidth(2); +} + +void printGaussianByMultiplicity(const std::array& summaries, const char* title) +{ + std::printf("%s:\n", title); + for (int i = 0; i < NMultiplicityBins; ++i) { + if (summaries[i].valid) { + std::printf(" %-4s : mean=%.6g sigma=%.6g\n", MultiplicityLabels[i], summaries[i].mean, summaries[i].sigma); + } else { + std::printf(" %-4s : n/a\n", MultiplicityLabels[i]); + } + } +} + +fs::path resolveContextFile(const fs::path& dir) +{ + const std::array candidates{ + dir / collContextFile, + dir.parent_path() / collContextFile, + fs::current_path() / collContextFile}; + for (const auto& candidate : candidates) { + if (!candidate.empty() && fs::exists(candidate) && fs::is_regular_file(candidate)) { + return candidate; + } + } + return {}; +} + +void printBinnedFractions(const TH1* numerator, const TH1* denominator, const char* title) +{ + if (numerator == nullptr || denominator == nullptr) { + return; + } + std::printf("%s:\n", title); + for (int iBin = 1; iBin <= denominator->GetNbinsX(); ++iBin) { + const auto den = denominator->GetBinContent(iBin); + const auto num = numerator->GetBinContent(iBin); + const auto value = den > 0. ? num / den : 0.; + std::printf(" %-4s : %.4f (%g / %g)\n", denominator->GetXaxis()->GetBinLabel(iBin), value, num, den); + } +} + +std::vector findDirs(const std::string& roots) +{ + fs::path root = roots.empty() ? fs::current_path() : fs::path{roots}; + std::vector result; + const auto hasFiles = [](const fs::path& dir) { + const auto tracPath = dir / tracFile; + return fs::exists(tracPath) && fs::is_regular_file(tracPath); + }; + + if (fs::is_directory(root) && hasFiles(root)) { + result.push_back(root); + return result; + } + + for (const auto& entry : fs::recursive_directory_iterator(root)) { + if (entry.is_directory() && hasFiles(entry.path())) { + result.push_back(entry.path()); + } + } + std::sort(result.begin(), result.end()); + return result; +} + +} // namespace + +void CheckSeeding(const std::string& dir = "") +{ + using Vertex = o2::its::Vertex; + const auto cwd = fs::current_path(); + gStyle->SetOptStat(0); + TH1::AddDirectory(kFALSE); + + auto dirs = findDirs(dir); + std::printf("Will iterate over %zu input dirs\n", dirs.size()); + if (dirs.empty()) { + std::printf("No input directories containing %s were found.\n", tracFile); + return; + } + + auto* hTruthMultiplicityFindable = new TH1D("hTruthMultiplicityFindable", + "Findable truth vertices;truth multiplicity bin;vertices", + NMultiplicityBins, 0.5, NMultiplicityBins + 0.5); + auto* hTruthMultiplicityFound = new TH1D("hTruthMultiplicityFound", + "Found truth vertices;truth multiplicity bin;vertices", + NMultiplicityBins, 0.5, NMultiplicityBins + 0.5); + auto* hRecoMultiplicityTrue = new TH1D("hRecoMultiplicityTrue", + "True reconstructed vertices;reco multiplicity bin;vertices", + NMultiplicityBins, 0.5, NMultiplicityBins + 0.5); + auto* hRecoMultiplicityFake = new TH1D("hRecoMultiplicityFake", + "Fake reconstructed vertices;reco multiplicity bin;vertices", + NMultiplicityBins, 0.5, NMultiplicityBins + 0.5); + auto* hDx = new TH1D("hDx", "Matched vertex residuals;x_{reco}-x_{MC} (cm);vertices", 400, -0.02, 0.02); + auto* hDy = new TH1D("hDy", "Matched vertex residuals;y_{reco}-y_{MC} (cm);vertices", 400, -0.02, 0.02); + auto* hDz = new TH1D("hDz", "Matched vertex residuals;z_{reco}-z_{MC} (cm);vertices", 400, -0.02, 0.02); + auto* hPullX = new TH1D("hPullX", "Matched vertex pulls;x pull;vertices", 600, -30., 30.); + auto* hPullY = new TH1D("hPullY", "Matched vertex pulls;y pull;vertices", 600, -30., 30.); + auto* hPullZ = new TH1D("hPullZ", "Matched vertex pulls;z pull;vertices", 600, -30., 30.); + std::array hPullXByMult{}; + std::array hPullYByMult{}; + std::array hPullZByMult{}; + for (int i = 0; i < NMultiplicityBins; ++i) { + const auto nameX = std::string("hPullX_") + std::to_string(i + 1); + const auto nameY = std::string("hPullY_") + std::to_string(i + 1); + const auto nameZ = std::string("hPullZ_") + std::to_string(i + 1); + const auto titleX = std::string("x pull ") + MultiplicityLabels[i] + ";x pull;vertices"; + const auto titleY = std::string("y pull ") + MultiplicityLabels[i] + ";y pull;vertices"; + const auto titleZ = std::string("z pull ") + MultiplicityLabels[i] + ";z pull;vertices"; + hPullXByMult[i] = new TH1D(nameX.c_str(), titleX.c_str(), 600, -30., 30.); + hPullYByMult[i] = new TH1D(nameY.c_str(), titleY.c_str(), 600, -30., 30.); + hPullZByMult[i] = new TH1D(nameZ.c_str(), titleZ.c_str(), 600, -30., 30.); + } + + setMultiplicityBinLabels(hTruthMultiplicityFindable); + setMultiplicityBinLabels(hTruthMultiplicityFound); + setMultiplicityBinLabels(hRecoMultiplicityTrue); + setMultiplicityBinLabels(hRecoMultiplicityFake); + setHistogramStyle(hTruthMultiplicityFindable, kGray + 2, 20); + setHistogramStyle(hTruthMultiplicityFound, kAzure + 2, 20); + setHistogramStyle(hRecoMultiplicityTrue, kAzure + 2, 20); + setHistogramStyle(hRecoMultiplicityFake, kOrange + 7, 24); + setHistogramStyle(hDx, kAzure + 2, 20); + setHistogramStyle(hDy, kGreen + 2, 21); + setHistogramStyle(hDz, kRed + 1, 24); + setHistogramStyle(hPullX, kAzure + 2, 20); + setHistogramStyle(hPullY, kGreen + 2, 21); + setHistogramStyle(hPullZ, kRed + 1, 24); + + size_t findable = 0; + size_t totalFound = 0; + size_t trueFound = 0; + size_t fakeFound = 0; + size_t uniqueTrueReco = 0; + size_t uniqueFindableFound = 0; + size_t sigmaXCount = 0; + size_t sigmaYCount = 0; + size_t sigmaZCount = 0; + double sumSigmaX = 0.; + double sumSigmaY = 0.; + double sumSigmaZ = 0.; + + for (const auto& inputDir : dirs) { + fs::current_path(inputDir); + std::printf("Working on %s\n", inputDir.c_str()); + const auto contextPath = resolveContextFile(inputDir); + if (contextPath.empty()) { + std::printf("Skipping %s: could not locate %s\n", inputDir.c_str(), collContextFile); + continue; + } + + o2::steer::MCKinematicsReader mcReader(contextPath.string()); + if (!mcReader.isInitialized()) { + std::printf("Skipping %s: failed to initialize MCKinematicsReader from %s\n", inputDir.c_str(), contextPath.c_str()); + continue; + } + + std::unordered_map findableTruths; + std::unordered_set uniqueTrueLabelsReco; + std::unordered_set uniqueFindableTruthFound; + std::unordered_map bestRecoByTruth; + + const int iSrc = 0; + const auto nEvents = static_cast(mcReader.getNEvents(iSrc)); + for (int iEve = 0; iEve < nEvents; ++iEve) { + const auto& tracks = mcReader.getTracks(iSrc, iEve); + const auto contributors = static_cast(std::count_if(tracks.begin(), tracks.end(), isChargedPrimary)); + if (contributors >= 2) { + const auto& header = mcReader.getMCEventHeader(iSrc, iEve); + findableTruths.emplace(iEve, TruthInfo{contributors, (float)header.GetX(), (float)header.GetY(), (float)header.GetZ()}); + fillMultiplicityHistogram(hTruthMultiplicityFindable, contributors); + } + mcReader.releaseTracksForSourceAndEvent(iSrc, iEve); + } + + auto* tracFileHandle = TFile::Open((inputDir / tracFile).c_str()); + if (tracFileHandle == nullptr || tracFileHandle->IsZombie()) { + std::printf("Skipping %s: failed to open %s\n", inputDir.c_str(), tracFile); + delete tracFileHandle; + continue; + } + + auto* tracTree = tracFileHandle->Get("o2sim"); + if (tracTree == nullptr) { + std::printf("Skipping %s: missing o2sim tree in %s\n", inputDir.c_str(), tracFile); + tracFileHandle->Close(); + delete tracFileHandle; + continue; + } + + if (tracTree->GetBranch("Vertices") == nullptr || tracTree->GetBranch("ITSVertexMCTruth") == nullptr) { + std::printf("Skipping %s: missing vertex branches in %s\n", inputDir.c_str(), tracFile); + tracFileHandle->Close(); + delete tracFileHandle; + continue; + } + + std::vector* vertices = nullptr; + std::vector* labels = nullptr; + std::vector* purities = nullptr; + const bool hasPurityBranch = tracTree->GetBranch("ITSVertexMCPurity") != nullptr; + + tracTree->SetBranchAddress("Vertices", &vertices); + tracTree->SetBranchAddress("ITSVertexMCTruth", &labels); + if (hasPurityBranch) { + tracTree->SetBranchAddress("ITSVertexMCPurity", &purities); + } + + const auto nEntries = tracTree->GetEntriesFast(); + for (Long64_t iEntry = 0; iEntry < nEntries; ++iEntry) { + tracTree->GetEntry(iEntry); + if (vertices == nullptr || labels == nullptr) { + continue; + } + auto nVertices = std::min(vertices->size(), labels->size()); + if (hasPurityBranch && purities != nullptr) { + nVertices = std::min(nVertices, purities->size()); + } + + for (size_t iVtx = 0; iVtx < nVertices; ++iVtx) { + const auto& vertex = (*vertices)[iVtx]; + const auto& label = (*labels)[iVtx]; + const auto multiplicity = static_cast(vertex.getNContributors()); + ++totalFound; + + if (!isTrueVertexLabel(label)) { + ++fakeFound; + fillMultiplicityHistogram(hRecoMultiplicityFake, multiplicity); + continue; + } + + ++trueFound; + const auto eventID = label.getEventID(); + uniqueTrueLabelsReco.insert(eventID); + fillMultiplicityHistogram(hRecoMultiplicityTrue, multiplicity); + + const auto truthIt = findableTruths.find(eventID); + if (truthIt == findableTruths.end()) { + continue; + } + + uniqueFindableTruthFound.insert(eventID); + const auto purity = (hasPurityBranch && purities != nullptr) ? (*purities)[iVtx] : -1.f; + const auto bestIt = bestRecoByTruth.find(eventID); + if (bestIt == bestRecoByTruth.end() || isBetterReco(vertex, purity, bestIt->second.vertex, bestIt->second.purity)) { + bestRecoByTruth[eventID] = BestRecoInfo{vertex, purity}; + } + } + } + + tracFileHandle->Close(); + delete tracFileHandle; + + findable += findableTruths.size(); + uniqueTrueReco += uniqueTrueLabelsReco.size(); + uniqueFindableFound += uniqueFindableTruthFound.size(); + + for (const auto eventID : uniqueFindableTruthFound) { + const auto truthIt = findableTruths.find(eventID); + if (truthIt != findableTruths.end()) { + fillMultiplicityHistogram(hTruthMultiplicityFound, truthIt->second.multiplicity); + } + } + + for (const auto& [eventID, reco] : bestRecoByTruth) { + const auto truthIt = findableTruths.find(eventID); + if (truthIt == findableTruths.end()) { + continue; + } + const auto dx = reco.vertex.getX() - truthIt->second.x; + const auto dy = reco.vertex.getY() - truthIt->second.y; + const auto dz = reco.vertex.getZ() - truthIt->second.z; + hDx->Fill(dx); + hDy->Fill(dy); + hDz->Fill(dz); + if (reco.vertex.getSigmaX() > 0.f) { + const auto pullX = dx / reco.vertex.getSigmaX(); + hPullX->Fill(pullX); + hPullXByMult[getMultiplicityCategory(reco.vertex.getNContributors()) - 1]->Fill(pullX); + sumSigmaX += reco.vertex.getSigmaX(); + ++sigmaXCount; + } + if (reco.vertex.getSigmaY() > 0.f) { + const auto pullY = dy / reco.vertex.getSigmaY(); + hPullY->Fill(pullY); + hPullYByMult[getMultiplicityCategory(reco.vertex.getNContributors()) - 1]->Fill(pullY); + sumSigmaY += reco.vertex.getSigmaY(); + ++sigmaYCount; + } + if (reco.vertex.getSigmaZ() > 0.f) { + const auto pullZ = dz / reco.vertex.getSigmaZ(); + hPullZ->Fill(pullZ); + hPullZByMult[getMultiplicityCategory(reco.vertex.getNContributors()) - 1]->Fill(pullZ); + sumSigmaZ += reco.vertex.getSigmaZ(); + ++sigmaZCount; + } + } + fs::current_path(cwd); + } + + auto* hTruthMultiplicityEfficiency = static_cast(hTruthMultiplicityFound->Clone("hTruthMultiplicityEfficiency")); + hTruthMultiplicityEfficiency->SetTitle("Unique efficiency vs truth multiplicity;truth multiplicity bin;efficiency"); + hTruthMultiplicityEfficiency->Divide(hTruthMultiplicityFound, hTruthMultiplicityFindable, 1., 1., "B"); + setMultiplicityBinLabels(hTruthMultiplicityEfficiency); + hTruthMultiplicityEfficiency->SetMinimum(0.); + hTruthMultiplicityEfficiency->SetMaximum(1.05); + + auto* hRecoMultiplicityTotal = static_cast(hRecoMultiplicityTrue->Clone("hRecoMultiplicityTotal")); + hRecoMultiplicityTotal->SetTitle("All reconstructed vertices;reco multiplicity bin;vertices"); + hRecoMultiplicityTotal->Add(hRecoMultiplicityFake); + setMultiplicityBinLabels(hRecoMultiplicityTotal); + + auto* hRecoMultiplicityPurity = static_cast(hRecoMultiplicityTrue->Clone("hRecoMultiplicityPurity")); + hRecoMultiplicityPurity->SetTitle("Purity vs reconstructed multiplicity;reco multiplicity bin;purity"); + hRecoMultiplicityPurity->Divide(hRecoMultiplicityTrue, hRecoMultiplicityTotal, 1., 1., "B"); + setMultiplicityBinLabels(hRecoMultiplicityPurity); + hRecoMultiplicityPurity->SetMinimum(0.); + hRecoMultiplicityPurity->SetMaximum(1.05); + + const auto duplicates = trueFound >= uniqueTrueReco ? (trueFound - uniqueTrueReco) : 0UL; + + const double uniqueEfficiency = findable > 0 ? static_cast(uniqueFindableFound) / findable : 0.; + const double purity = totalFound > 0 ? static_cast(trueFound) / totalFound : 0.; + const double fakeRate = totalFound > 0 ? static_cast(fakeFound) / totalFound : 0.; + const double duplicateRate = trueFound > 0 ? static_cast(duplicates) / trueFound : 0.; + const double f1 = (uniqueEfficiency + purity) > 0. ? 2. * uniqueEfficiency * purity / (uniqueEfficiency + purity) : 0.; + + const auto dxFit = fitGaussianCore(hDx, "fitDx"); + const auto dyFit = fitGaussianCore(hDy, "fitDy"); + const auto dzFit = fitGaussianCore(hDz, "fitDz"); + const auto pullXFit = fitGaussianCore(hPullX, "fitPullX"); + const auto pullYFit = fitGaussianCore(hPullY, "fitPullY"); + const auto pullZFit = fitGaussianCore(hPullZ, "fitPullZ"); + std::array pullXByMultFit{}; + std::array pullYByMultFit{}; + std::array pullZByMultFit{}; + for (int i = 0; i < NMultiplicityBins; ++i) { + const auto fitX = std::string("fitPullX_") + std::to_string(i + 1); + const auto fitY = std::string("fitPullY_") + std::to_string(i + 1); + const auto fitZ = std::string("fitPullZ_") + std::to_string(i + 1); + pullXByMultFit[i] = fitGaussianCore(hPullXByMult[i], fitX.c_str()); + pullYByMultFit[i] = fitGaussianCore(hPullYByMult[i], fitY.c_str()); + pullZByMultFit[i] = fitGaussianCore(hPullZByMult[i], fitZ.c_str()); + } + + std::printf("\nVertex validation summary\n"); + std::printf(" findable truth vertices : %zu\n", findable); + std::printf(" total reconstructed vertices : %zu\n", totalFound); + std::printf(" true reconstructed vertices : %zu\n", trueFound); + std::printf(" fake reconstructed vertices : %zu\n", fakeFound); + std::printf(" unique true labels (all) : %zu\n", uniqueTrueReco); + std::printf(" unique findable truth found : %zu\n", uniqueFindableFound); + std::printf(" unique efficiency : %.5f\n", uniqueEfficiency); + std::printf(" purity : %.5f\n", purity); + std::printf(" fake rate : %.5f\n", fakeRate); + std::printf(" duplicate rate : %.5f\n", duplicateRate); + std::printf(" F1(purity,efficiency) : %.5f\n", f1); + std::printf(" mean reported sigma x/y/z : %.6g / %.6g / %.6g cm\n", + sigmaXCount > 0 ? sumSigmaX / sigmaXCount : 0., + sigmaYCount > 0 ? sumSigmaY / sigmaYCount : 0., + sigmaZCount > 0 ? sumSigmaZ / sigmaZCount : 0.); + + if (dxFit.valid) { + std::printf(" x residual Gaussian: mean=%.6g cm sigma=%.6g cm\n", dxFit.mean, dxFit.sigma); + } + if (dyFit.valid) { + std::printf(" y residual Gaussian: mean=%.6g cm sigma=%.6g cm\n", dyFit.mean, dyFit.sigma); + } + if (dzFit.valid) { + std::printf(" z residual Gaussian: mean=%.6g cm sigma=%.6g cm\n", dzFit.mean, dzFit.sigma); + } + if (pullXFit.valid) { + std::printf(" x pull Gaussian : mean=%.6g sigma=%.6g\n", pullXFit.mean, pullXFit.sigma); + } + if (pullYFit.valid) { + std::printf(" y pull Gaussian : mean=%.6g sigma=%.6g\n", pullYFit.mean, pullYFit.sigma); + } + if (pullZFit.valid) { + std::printf(" z pull Gaussian : mean=%.6g sigma=%.6g\n", pullZFit.mean, pullZFit.sigma); + } + printGaussianByMultiplicity(pullXByMultFit, "x pull Gaussian by reconstructed multiplicity"); + printGaussianByMultiplicity(pullYByMultFit, "y pull Gaussian by reconstructed multiplicity"); + printGaussianByMultiplicity(pullZByMultFit, "z pull Gaussian by reconstructed multiplicity"); + + printBinnedFractions(hTruthMultiplicityFound, hTruthMultiplicityFindable, "Efficiency vs truth multiplicity"); + printBinnedFractions(hRecoMultiplicityTrue, hRecoMultiplicityTotal, "Purity vs reconstructed multiplicity"); + + auto* cValidation = new TCanvas("cVertexValidation", "Vertex validation summary", 1800, 1000); + cValidation->Divide(3, 2); + + cValidation->cd(1); + gPad->SetMargin(0.05, 0.05, 0.05, 0.05); + auto* summary = new TPaveText(0.02, 0.02, 0.98, 0.98, "NDC"); + summary->SetBorderSize(0); + summary->SetFillColor(0); + summary->SetTextAlign(12); + summary->SetTextFont(42); + char line[256]; + std::snprintf(line, sizeof(line), "Findable truth vertices : %zu", findable); + summary->AddText(line); + std::snprintf(line, sizeof(line), "Total reconstructed : %zu", totalFound); + summary->AddText(line); + std::snprintf(line, sizeof(line), "True reconstructed : %zu", trueFound); + summary->AddText(line); + std::snprintf(line, sizeof(line), "Fake reconstructed : %zu", fakeFound); + summary->AddText(line); + summary->AddText(""); + std::snprintf(line, sizeof(line), "Unique truth found : %zu", uniqueFindableFound); + summary->AddText(line); + std::snprintf(line, sizeof(line), "Unique efficiency : %.5f", uniqueEfficiency); + summary->AddText(line); + std::snprintf(line, sizeof(line), "Purity : %.5f", purity); + summary->AddText(line); + std::snprintf(line, sizeof(line), "Fake rate : %.5f", fakeRate); + summary->AddText(line); + std::snprintf(line, sizeof(line), "Duplicate rate : %.5f", duplicateRate); + summary->AddText(line); + std::snprintf(line, sizeof(line), "F1 : %.5f", f1); + summary->AddText(line); + std::snprintf(line, sizeof(line), "mean sigma x/y/z cm : %.3g / %.3g / %.3g", + sigmaXCount > 0 ? sumSigmaX / sigmaXCount : 0., + sigmaYCount > 0 ? sumSigmaY / sigmaYCount : 0., + sigmaZCount > 0 ? sumSigmaZ / sigmaZCount : 0.); + summary->AddText(line); + summary->AddText(""); + if (dxFit.valid) { + std::snprintf(line, sizeof(line), "dx fit mean/sigma cm : %.3g / %.3g", dxFit.mean, dxFit.sigma); + summary->AddText(line); + } + if (dyFit.valid) { + std::snprintf(line, sizeof(line), "dy fit mean/sigma cm : %.3g / %.3g", dyFit.mean, dyFit.sigma); + summary->AddText(line); + } + if (dzFit.valid) { + std::snprintf(line, sizeof(line), "dz fit mean/sigma cm : %.3g / %.3g", dzFit.mean, dzFit.sigma); + summary->AddText(line); + } + if (pullXFit.valid) { + std::snprintf(line, sizeof(line), "pull x sigma : %.3g", pullXFit.sigma); + summary->AddText(line); + } + if (pullYFit.valid) { + std::snprintf(line, sizeof(line), "pull y sigma : %.3g", pullYFit.sigma); + summary->AddText(line); + } + if (pullZFit.valid) { + std::snprintf(line, sizeof(line), "pull z sigma : %.3g", pullZFit.sigma); + summary->AddText(line); + } + summary->Draw(); + + cValidation->cd(2); + gPad->SetGridy(); + hTruthMultiplicityEfficiency->Draw("hist e1"); + + cValidation->cd(3); + gPad->SetGridy(); + const auto maxReco = std::max(hRecoMultiplicityTrue->GetMaximum(), hRecoMultiplicityFake->GetMaximum()); + hRecoMultiplicityTrue->SetMaximum(1.2 * std::max(1., maxReco)); + hRecoMultiplicityTrue->Draw("hist e1"); + hRecoMultiplicityFake->Draw("hist e1 same"); + auto hRecoMultiplicitySum = (TH1D*)hRecoMultiplicityTrue->Clone("hRecoMultiplicitySum"); + hRecoMultiplicitySum->Add(hRecoMultiplicityFake); + setHistogramStyle(hRecoMultiplicitySum, kBlack, 23); + hRecoMultiplicitySum->Draw("hist e1 same"); + { + auto* legend = new TLegend(0.58, 0.75, 0.88, 0.88); + legend->SetBorderSize(0); + legend->AddEntry(hRecoMultiplicityTrue, "true", "lep"); + legend->AddEntry(hRecoMultiplicityFake, "fake", "lep"); + legend->AddEntry(hRecoMultiplicitySum, "sum", "lep"); + legend->Draw(); + } + + cValidation->cd(4); + gPad->SetGridy(); + hRecoMultiplicityPurity->Draw("hist e1"); + + cValidation->cd(5); + gPad->SetGridy(); + auto* hDxNorm = makeNormalizedCopy(hDx, "hDxNorm", "Matched vertex residuals;residual (cm);normalized entries"); + auto* hDyNorm = makeNormalizedCopy(hDy, "hDyNorm", "Matched vertex residuals;residual (cm);normalized entries"); + auto* hDzNorm = makeNormalizedCopy(hDz, "hDzNorm", "Matched vertex residuals;residual (cm);normalized entries"); + const auto maxResidual = std::max({hDxNorm->GetMaximum(), hDyNorm->GetMaximum(), hDzNorm->GetMaximum()}); + hDzNorm->SetMaximum(1.2 * std::max(1., maxResidual)); + hDzNorm->Draw("hist"); + hDxNorm->Draw("hist same"); + hDyNorm->Draw("hist same"); + { + auto* legend = new TLegend(0.62, 0.72, 0.88, 0.88); + legend->SetBorderSize(0); + legend->AddEntry(hDxNorm, "dx", "l"); + legend->AddEntry(hDyNorm, "dy", "l"); + legend->AddEntry(hDzNorm, "dz", "l"); + legend->Draw(); + } + + cValidation->cd(6); + gPad->SetGridy(); + auto* hPullXNorm = makeNormalizedCopy(hPullX, "hPullXNorm", "Matched vertex pulls;pull;normalized entries"); + auto* hPullYNorm = makeNormalizedCopy(hPullY, "hPullYNorm", "Matched vertex pulls;pull;normalized entries"); + auto* hPullZNorm = makeNormalizedCopy(hPullZ, "hPullZNorm", "Matched vertex pulls;pull;normalized entries"); + const auto maxPull = std::max({hPullXNorm->GetMaximum(), hPullYNorm->GetMaximum(), hPullZNorm->GetMaximum()}); + hPullZNorm->SetMaximum(1.2 * std::max(1., maxPull)); + hPullZNorm->Draw("hist"); + hPullXNorm->Draw("hist same"); + hPullYNorm->Draw("hist same"); + { + auto* legend = new TLegend(0.62, 0.72, 0.88, 0.88); + legend->SetBorderSize(0); + legend->AddEntry(hPullXNorm, "pull x", "l"); + legend->AddEntry(hPullYNorm, "pull y", "l"); + legend->AddEntry(hPullZNorm, "pull z", "l"); + legend->Draw(); + } + + cValidation->cd(); + cValidation->Update(); + cValidation->SaveAs("checkSeeding.pdf"); +} diff --git a/Detectors/ITSMFT/ITS/tracking/CMakeLists.txt b/Detectors/ITSMFT/ITS/tracking/CMakeLists.txt index 8d8304d16764f..c84040e87a30f 100644 --- a/Detectors/ITSMFT/ITS/tracking/CMakeLists.txt +++ b/Detectors/ITSMFT/ITS/tracking/CMakeLists.txt @@ -16,6 +16,7 @@ o2_add_library(ITStracking src/Configuration.cxx src/FastMultEstConfig.cxx src/FastMultEst.cxx + src/LineVertexerHelpers.cxx src/TimeFrame.cxx src/IOUtils.cxx src/Tracker.cxx diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ClusterLines.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ClusterLines.h index 6fbc6d7da7721..bcb8a98a62cab 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ClusterLines.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ClusterLines.h @@ -12,8 +12,9 @@ #ifndef O2_ITS_CLUSTERLINES_H #define O2_ITS_CLUSTERLINES_H -#include +#include #include +#include #include #include #include "ITStracking/Cluster.h" @@ -59,6 +60,7 @@ class ClusterLines final public: ClusterLines() = default; ClusterLines(const int firstLabel, const Line& firstLine, const int secondLabel, const Line& secondLine); + ClusterLines(gsl::span lineLabels, gsl::span lines); void add(const int lineLabel, const Line& line); void computeClusterCentroid(); void accumulate(const Line& line); @@ -67,7 +69,7 @@ class ClusterLines final const float* getRMS2() const { return mRMS2.Array(); } float getAvgDistance2() const { return mAvgDistance2; } auto getSize() const noexcept { return mLabels.size(); } - auto& getLabels() noexcept { return mLabels; } + auto& getLabels() const noexcept { return mLabels; } const auto& getTimeStamp() const noexcept { return mTime; } bool operator==(const ClusterLines& rhs) const noexcept; float getR2() const noexcept { return (mVertex[0] * mVertex[0]) + (mVertex[1] * mVertex[1]); } diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h index 02dbeb8cf3992..1f55a95ca0d65 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h @@ -18,11 +18,9 @@ #include #ifndef GPUCA_GPUCODE_DEVICE -#include #include #include #include -#include #endif #include "DetectorsBase/Propagator.h" @@ -89,21 +87,28 @@ struct VertexingParameters { int nIterations = 1; // Number of vertexing passes to perform 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 ZBins{1}; - int PhiBins{128}; - float zCut = 0.002f; - float phiCut = 0.005f; - float pairCut = 0.04f; - float clusterCut = 0.8f; - float histPairCut = 0.04f; - float tanLambdaCut = 0.002f; // tanLambda = deltaZ/deltaR - float lowMultBeamDistCut = 0.1f; // XY cut for low-multiplicity pile up - int vertNsigmaCut = 6; // N sigma cut for vertex XY - float vertRadiusSigma = 0.33f; // sigma of vertex XY - float trackletSigma = 0.01f; // tracklet to vertex sigma - float maxZPositionAllowed = 25.f; - int clusterContributorsCut = 16; - int maxTrackletsPerCluster = 2e3; + int ZBins = 1; + int PhiBins = 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 vertNsigmaCut = -1.f; + float vertRadiusSigma = -1.f; + float trackletSigma = -1.f; + 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; diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/LineVertexerHelpers.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/LineVertexerHelpers.h new file mode 100644 index 0000000000000..0e3807aba8efb --- /dev/null +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/LineVertexerHelpers.h @@ -0,0 +1,46 @@ +// Copyright 2019-2026 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 O2_ITS_TRACKING_LINE_VERTEXER_HELPERS_H_ +#define O2_ITS_TRACKING_LINE_VERTEXER_HELPERS_H_ + +#include +#include + +#include "ITStracking/BoundedAllocator.h" +#include "ITStracking/ClusterLines.h" + +namespace o2::its::line_vertexer +{ + +struct Settings { + float beamX = 0.f; + float beamY = 0.f; + float pairCut = 0.f; + float pairCut2 = 0.f; + float clusterCut = 0.f; + float coarseZWindow = 0.f; + float seedDedupZCut = 0.f; + float refitDedupZCut = 0.f; + float duplicateZCut = 0.f; + float duplicateDistance2Cut = 0.f; + float finalSelectionZCut = 0.f; + float maxZ = 0.f; + int seedMemberRadiusTime = 1; + int seedMemberRadiusZ = 2; + std::shared_ptr memoryPool; +}; + +bounded_vector buildClusters(std::span lines, const Settings& settings); + +} // namespace o2::its::line_vertexer + +#endif diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/MathUtils.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/MathUtils.h index 95e0b4554e32c..ab3c7d5d29873 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/MathUtils.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/MathUtils.h @@ -94,6 +94,16 @@ GPUhdi() constexpr float Sq(float v) return v * v; } +GPUhdi() constexpr float SqSum(float v, float w) +{ + return Sq(v) + Sq(w); +} + +GPUhdi() constexpr float SqSum(float u, float v, float w) +{ + return Sq(u) + SqSum(v, w); +} + GPUhdi() constexpr float SqDiff(float x, float y) { return Sq(x - y); diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h index e77200a1432d1..cb291b46f5e44 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h @@ -24,24 +24,30 @@ struct VertexerParamConfig : public o2::conf::ConfigurableParamHelper50% + elem.setFakeFlag(); + } return std::make_pair(elem, static_cast(maxCount) / static_cast(elements.size())); } diff --git a/Detectors/ITSMFT/ITS/tracking/src/ClusterLines.cxx b/Detectors/ITSMFT/ITS/tracking/src/ClusterLines.cxx index f561fe0436c4a..3e3e1b8b46338 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/ClusterLines.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/ClusterLines.cxx @@ -148,6 +148,38 @@ ClusterLines::ClusterLines(const int firstLabel, const Line& firstLine, const in mAvgDistance2 += (Line::getDistance2FromPoint(secondLine, mVertex) - mAvgDistance2) / (float)getSize(); } +ClusterLines::ClusterLines(gsl::span lineLabels, gsl::span lines) +{ + if (lineLabels.size() < 2) { + return; + } + + mLabels.reserve(lineLabels.size()); + mTime = lines[lineLabels[0]].mTime; + for (size_t index = 0; index < lineLabels.size(); ++index) { + const auto lineLabel = lineLabels[index]; + if (index > 0) { + mTime += lines[lineLabel].mTime; + } + mLabels.push_back(lineLabel); + accumulate(lines[lineLabel]); + } + + computeClusterCentroid(); + if (!mIsValid) { + return; + } + + mRMS2 = Line::getDCAComponents(lines[lineLabels[0]], mVertex); + mAvgDistance2 = Line::getDistance2FromPoint(lines[lineLabels[0]], mVertex); + for (size_t index = 1; index < lineLabels.size(); ++index) { + const auto lineLabel = lineLabels[index]; + const auto tmpRMS2 = Line::getDCAComponents(lines[lineLabel], mVertex); + mRMS2 += (tmpRMS2 - mRMS2) * (1.f / static_cast(index + 1)); + mAvgDistance2 += (Line::getDistance2FromPoint(lines[lineLabel], mVertex) - mAvgDistance2) / static_cast(index + 1); + } +} + void ClusterLines::add(const int lineLabel, const Line& line) { mTime += line.mTime; diff --git a/Detectors/ITSMFT/ITS/tracking/src/Configuration.cxx b/Detectors/ITSMFT/ITS/tracking/src/Configuration.cxx index c447bb6bcc880..6c88b61f2df07 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/Configuration.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/Configuration.cxx @@ -45,8 +45,8 @@ std::string TrackingParameters::asString() const } if (!AddTimeError.empty()) { str += " AddTimeError:"; - for (size_t i = 0; i < AddTimeError.size(); i++) { - str += std::format("{} ", AddTimeError[i]); + for (unsigned int i : AddTimeError) { + str += std::format("{} ", i); } } if (std::numeric_limits::max() != MaxMemory) { @@ -57,7 +57,8 @@ std::string TrackingParameters::asString() const std::string VertexingParameters::asString() const { - std::string str = std::format("NZb:{} NPhB:{} ClsCont:{} MaxTrkltCls:{} ZCut:{} PhCut:{}", ZBins, PhiBins, clusterContributorsCut, maxTrackletsPerCluster, zCut, phiCut); + std::string str = std::format("NZb:{} NPhB:{} MinVtxCont:{} SupLowMultDebris:{} MaxTrkltCls:{} ZCut:{} PhCut:{} PairCut:{} ClCut:{} SeedRad:{}x{}", + ZBins, PhiBins, clusterContributorsCut, suppressLowMultDebris, maxTrackletsPerCluster, zCut, phiCut, pairCut, clusterCut, seedMemberRadiusTime, seedMemberRadiusZ); if (std::numeric_limits::max() != MaxMemory) { str += std::format(" MemLimit {:.2f} GB", double(MaxMemory) / constants::GB); } @@ -173,8 +174,8 @@ std::vector TrackingMode::getTrackingParameters(TrackingMode LOGP(fatal, "Unsupported ITS tracking mode {} ", toString(mode)); } - float bFactor = std::abs(o2::base::Propagator::Instance()->getNominalBz()) / 5.0066791; - float bFactorTracklets = bFactor < 0.01 ? 1. : bFactor; // for tracklets only + float bFactor = std::abs(o2::base::Propagator::Instance()->getNominalBz()) / 5.0066791f; + float bFactorTracklets = bFactor < 0.01f ? 1.f : bFactor; // for tracklets only // global parameters set for every iteration for (auto& p : trackParams) { @@ -262,6 +263,9 @@ std::vector TrackingMode::getVertexingParameters(TrackingMo p.trackletSigma = vc.trackletSigma; p.maxZPositionAllowed = vc.maxZPositionAllowed; p.clusterContributorsCut = vc.clusterContributorsCut; + p.suppressLowMultDebris = vc.suppressLowMultDebris; + p.seedMemberRadiusTime = vc.seedMemberRadiusTime; + p.seedMemberRadiusZ = vc.seedMemberRadiusZ; p.phiSpan = vc.phiSpan; p.nThreads = vc.nThreads; p.ZBins = vc.ZBins; @@ -273,12 +277,16 @@ std::vector TrackingMode::getVertexingParameters(TrackingMo vertParams[0].vertNsigmaCut = vc.vertNsigmaCut; vertParams[0].vertRadiusSigma = vc.vertRadiusSigma; vertParams[0].maxTrackletsPerCluster = vc.maxTrackletsPerCluster; - vertParams[0].lowMultBeamDistCut = vc.lowMultBeamDistCut; vertParams[0].zCut = vc.zCut; vertParams[0].phiCut = vc.phiCut; vertParams[0].pairCut = vc.pairCut; vertParams[0].clusterCut = vc.clusterCut; - vertParams[0].histPairCut = vc.histPairCut; + vertParams[0].coarseZWindow = vc.coarseZWindow; + vertParams[0].seedDedupZCut = vc.seedDedupZCut; + vertParams[0].refitDedupZCut = vc.refitDedupZCut; + vertParams[0].duplicateZCut = vc.duplicateZCut; + vertParams[0].finalSelectionZCut = vc.finalSelectionZCut; + vertParams[0].duplicateDistance2Cut = vc.duplicateDistance2Cut; vertParams[0].tanLambdaCut = vc.tanLambdaCut; return vertParams; diff --git a/Detectors/ITSMFT/ITS/tracking/src/LineVertexerHelpers.cxx b/Detectors/ITSMFT/ITS/tracking/src/LineVertexerHelpers.cxx new file mode 100644 index 0000000000000..592c22dedf347 --- /dev/null +++ b/Detectors/ITSMFT/ITS/tracking/src/LineVertexerHelpers.cxx @@ -0,0 +1,1036 @@ +// Copyright 2019-2026 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 +#include + +#include +#include + +#include "ITStracking/Constants.h" +#include "ITStracking/MathUtils.h" +#include "ITStracking/LineVertexerHelpers.h" + +namespace o2::its::line_vertexer +{ +namespace +{ +using SymMatrix3 = ROOT::Math::SMatrix>; +using SVector3 = ROOT::Math::SVector; + +constexpr float TukeyC = 4.685f; +constexpr float TukeyC2 = TukeyC * TukeyC; +constexpr float InitialScale2 = 5.f; +constexpr float MinScale2 = 1.f; +constexpr float MedianToSigma = 1.4826f; +constexpr float VertexShiftZTol = 0.01f; +constexpr float VertexShiftR2Tol = 1.e-4f; +constexpr int MaxFitIterations = 10; +constexpr int MaxSeedsPerCluster = 32; +constexpr float MinRelativePeakSupport = 0.1f; +constexpr int MaxHistogramBins = 0x7fff; +constexpr float TieTolerance = 1e-5f; + +struct LineRef { + LineRef(const Line& line, const int index, const float beamX, const float beamY, const float maxZ) : lineIndex(index) + { + const auto symTime = line.mTime.makeSymmetrical(); + tCenter = symTime.getTimeStamp(); + tHalfWidth = symTime.getTimeStampError(); + const auto dx = line.originPoint(0) - beamX; + const auto dy = line.originPoint(1) - beamY; + const auto ux = line.cosinesDirector(0); + const auto uy = line.cosinesDirector(1); + const auto uz = line.cosinesDirector(2); + const auto den = math_utils::SqSum(ux, uy); + if (den <= constants::Tolerance) { + lineIndex = constants::UnusedIndex; + return; + } + const auto s0 = -((dx * ux) + (dy * uy)) / den; + const auto xb = dx + (s0 * ux); + const auto yb = dy + (s0 * uy); + zBeam = line.originPoint(2) + s0 * uz; + if (!std::isfinite(zBeam) || o2::gpu::CAMath::Abs(zBeam) > maxZ) { + lineIndex = constants::UnusedIndex; + } + } + bool isDead() const noexcept { return lineIndex == constants::UnusedIndex; } + + int lineIndex = constants::UnusedIndex; + float zBeam = 0.f; + float tCenter = 0.f; + float tHalfWidth = 0.f; +}; + +struct VertexSeed { + explicit VertexSeed(const std::shared_ptr& mr) : contributors(mr.get()), assigned(mr.get()) {} + + std::array vertex = {}; + TimeEstBC time; + float scale2 = InitialScale2; + bounded_vector contributors; + bounded_vector assigned; + bool valid = false; + bool isUsableSeed() const noexcept + { + return valid && contributors.size() >= 2; + } +}; + +void compactSeeds(bounded_vector& seeds) +{ + seeds.erase(std::remove_if(seeds.begin(), seeds.end(), [](const VertexSeed& seed) { + return !seed.isUsableSeed(); + }), + seeds.end()); +} + +struct Histogram2D { + explicit Histogram2D(const std::shared_ptr& mr) : bins(mr.get()) {} + + int nTimeBins = 0; + int nZBins = 0; + float timeMin = 0.f; + float zMin = 0.f; + float timeBinSize = 1.f; + float zBinSize = 1.f; + bounded_vector bins; + + int getIndex(const int tBin, const int zBin) const noexcept + { + return (tBin * nZBins) + zBin; + } + + std::pair decodeIndex(const int index) const noexcept + { + return {index / nZBins, index % nZBins}; + } + + int getTimeBin(const float time) const noexcept + { + if (time < timeMin) { + return -1; + } + const auto bin = static_cast((time - timeMin) / timeBinSize); + return (bin >= 0 && bin < nTimeBins) ? bin : -1; + } + + int getZBin(const float z) const noexcept + { + if (z < zMin) { + return -1; + } + const auto bin = static_cast((z - zMin) / zBinSize); + return (bin >= 0 && bin < nZBins) ? bin : -1; + } + + void fill(const float time, const float z, const float weight) noexcept + { + const auto tBin = getTimeBin(time); + const auto zBin = getZBin(z); + if (tBin < 0 || zBin < 0) { + return; + } + bins[getIndex(tBin, zBin)] += weight; + } + + int findPeakBin() const noexcept + { + float bestWeight = 0.f; + int bestIndex = -1; + for (int index = 0; index < static_cast(bins.size()); ++index) { + if (bins[index] > bestWeight) { + bestWeight = bins[index]; + bestIndex = index; + } + } + return bestIndex; + } + + void suppressBin(const int index) noexcept + { + if (index >= 0 && index < static_cast(bins.size())) { + bins[index] = -1.f; + } + } + + void suppressNeighborhood(const int index, const int radiusTime, const int radiusZ) noexcept + { + if (index < 0) { + return; + } + const auto [tBin, zBin] = decodeIndex(index); + for (int dt = -radiusTime; dt <= radiusTime; ++dt) { + const auto tt = tBin + dt; + if (tt < 0 || tt >= nTimeBins) { + continue; + } + for (int dz = -radiusZ; dz <= radiusZ; ++dz) { + const auto zz = zBin + dz; + if (zz < 0 || zz >= nZBins) { + continue; + } + bins[getIndex(tt, zz)] = -1.f; + } + } + } + + float getNeighborhoodSum(const int index, const int radiusTime, const int radiusZ) const noexcept + { + if (index < 0) { + return 0.f; + } + const auto [tBin, zBin] = decodeIndex(index); + float sum = 0.f; + for (int dt = -radiusTime; dt <= radiusTime; ++dt) { + const auto tt = tBin + dt; + if (tt < 0 || tt >= nTimeBins) { + continue; + } + for (int dz = -radiusZ; dz <= radiusZ; ++dz) { + const auto zz = zBin + dz; + if (zz < 0 || zz >= nZBins) { + continue; + } + const auto value = bins[getIndex(tt, zz)]; + if (value > 0.f) { + sum += value; + } + } + } + return sum; + } + + float getTimeBinCenter(const int tBin) const noexcept + { + return timeMin + ((static_cast(tBin) + 0.5f) * timeBinSize); + } + + float getZBinCenter(const int zBin) const noexcept + { + return zMin + ((static_cast(zBin) + 0.5f) * zBinSize); + } + + TimeEstBC getTimeInterval(const int tBin) const noexcept + { + const auto lowFloat = timeMin + (static_cast(tBin) * timeBinSize); + const auto highFloat = lowFloat + timeBinSize; + const auto low = std::max(0., std::floor(lowFloat)); + const auto high = std::max(low + 1., (double)std::ceil(highFloat)); + constexpr auto maxTS = std::numeric_limits::max(); + const auto clampedLow = std::min(low, maxTS - 1.); + const auto width = std::min(high - clampedLow, std::numeric_limits::max()); + return {static_cast(clampedLow), static_cast(std::max(1., width))}; + } + + TimeEstBC getTimeNeighborhoodInterval(const int tBin, const int radius) const noexcept + { + const auto lowBin = std::max(0, tBin - radius); + const auto highBin = std::min(nTimeBins - 1, tBin + radius); + const auto lowFloat = timeMin + (static_cast(lowBin) * timeBinSize); + const auto highFloat = timeMin + (static_cast(highBin + 1) * timeBinSize); + const auto low = std::max(0., std::floor(lowFloat)); + const auto high = std::max(low + 1., (double)std::ceil(highFloat)); + constexpr auto maxTS = std::numeric_limits::max(); + const auto clampedLow = std::min(low, maxTS - 1.); + const auto width = std::min(high - clampedLow, std::numeric_limits::max()); + return {static_cast(clampedLow), static_cast(std::max(1., width))}; + } +}; + +class SeedHistogram +{ + public: + SeedHistogram(std::span members, + std::span lineRefs, + std::span lines, + const Settings& settings) + : mMembers(members), mLineRefs(lineRefs), mSeedMemberRadiusTime(settings.seedMemberRadiusTime), mSeedMemberRadiusZ(settings.seedMemberRadiusZ), mMemoryPool(settings.memoryPool), mHistogram(mMemoryPool) + { + const auto zBinSize = 0.25f * settings.clusterCut; + const auto timeBinSize = medianTimeError(lines); + + float minZ = std::numeric_limits::max(); + float maxZ = std::numeric_limits::lowest(); + float minTime = std::numeric_limits::max(); + float maxTime = std::numeric_limits::lowest(); + for (const auto lineRefIdx : mMembers) { + minZ = std::min(minZ, mLineRefs[lineRefIdx].zBeam); + maxZ = std::max(maxZ, mLineRefs[lineRefIdx].zBeam); + minTime = std::min(minTime, mLineRefs[lineRefIdx].tCenter); + maxTime = std::max(maxTime, mLineRefs[lineRefIdx].tCenter); + } + + const auto dz = std::max(0.f, maxZ - minZ); + const auto dt = std::max(0.f, maxTime - minTime); + mHistogram.nZBins = 1 + static_cast(dz / zBinSize); + mHistogram.nTimeBins = 1 + static_cast(dt / timeBinSize); + if (mHistogram.nTimeBins * mHistogram.nZBins > MaxHistogramBins) { + if (mHistogram.nTimeBins > mHistogram.nZBins) { + mHistogram.nTimeBins = std::max(1, (MaxHistogramBins - 1) / std::max(1, mHistogram.nZBins)); + } else { + mHistogram.nZBins = std::max(1, (MaxHistogramBins - 1) / std::max(1, mHistogram.nTimeBins)); + } + } + + mHistogram.timeBinSize = std::max(timeBinSize, dt / (float)std::max(1, mHistogram.nTimeBins)); + mHistogram.zBinSize = std::max(zBinSize, dz / (float)std::max(1, mHistogram.nZBins)); + const auto paddedTime = 0.5f * ((float)mHistogram.nTimeBins * mHistogram.timeBinSize - dt); + const auto paddedZ = 0.5f * ((float)mHistogram.nZBins * mHistogram.zBinSize - dz); + mHistogram.timeMin = minTime - paddedTime; + mHistogram.zMin = minZ - paddedZ; + mHistogram.bins.assign((size_t)mHistogram.nTimeBins * (size_t)mHistogram.nZBins, 0.f); + + for (const auto lineRefIdx : mMembers) { + mHistogram.fill(mLineRefs[lineRefIdx].tCenter, mLineRefs[lineRefIdx].zBeam, 1.f); + } + } + + int findPeakBin() const noexcept + { + return mHistogram.findPeakBin(); + } + + float getPeakSupport(const int peakIndex) const noexcept + { + return mHistogram.getNeighborhoodSum(peakIndex, mSeedMemberRadiusTime, mSeedMemberRadiusZ); + } + + bounded_vector collectLocalMembers(const int peakIndex, const int radiusTime, const int radiusZ) const + { + bounded_vector localMembers(mMemoryPool.get()); + localMembers.reserve(mMembers.size()); + const auto [timeBin, zBin] = mHistogram.decodeIndex(peakIndex); + for (const auto lineRefIdx : mMembers) { + const auto memberTimeBin = mHistogram.getTimeBin(mLineRefs[lineRefIdx].tCenter); + const auto memberZBin = mHistogram.getZBin(mLineRefs[lineRefIdx].zBeam); + if (memberTimeBin < 0 || memberZBin < 0) { + continue; + } + if (o2::gpu::GPUCommonMath::Abs(memberTimeBin - timeBin) > radiusTime) { + continue; + } + if (o2::gpu::GPUCommonMath::Abs(memberZBin - zBin) > radiusZ) { + continue; + } + localMembers.push_back(lineRefIdx); + } + return localMembers; + } + + TimeEstBC getPeakTimeInterval(const int peakIndex, const int radius = 0) const noexcept + { + return mHistogram.getTimeNeighborhoodInterval(mHistogram.decodeIndex(peakIndex).first, radius); + } + + float getPeakZCenter(const int peakIndex) const noexcept + { + return mHistogram.getZBinCenter(mHistogram.decodeIndex(peakIndex).second); + } + + void suppressPeak(const int peakIndex) noexcept + { + mHistogram.suppressBin(peakIndex); + } + + void suppressPeakNeighborhood(const int peakIndex) noexcept + { + mHistogram.suppressNeighborhood(peakIndex, mSeedMemberRadiusTime, mSeedMemberRadiusZ); + } + + private: + float medianTimeError(std::span lines) const + { + bounded_vector errors(mMemoryPool.get()); + errors.reserve(mMembers.size()); + for (const auto lineRefIdx : mMembers) { + errors.push_back(static_cast(lines[mLineRefs[lineRefIdx].lineIndex].mTime.getTimeStampError())); + } + std::sort(errors.begin(), errors.end()); + return errors.empty() ? 1.f : std::max(1.f, errors[errors.size() / 2]); + } + + std::span mMembers; + std::span mLineRefs; + int mSeedMemberRadiusTime = 1; + int mSeedMemberRadiusZ = 2; + std::shared_ptr mMemoryPool; + Histogram2D mHistogram; +}; + +float updateScale2(const std::span chi2s, const std::shared_ptr& mr) noexcept +{ + if (chi2s.empty()) { + return MinScale2; + } + + bounded_vector sorted(chi2s.begin(), chi2s.end(), mr.get()); + std::sort(sorted.begin(), sorted.end()); + const auto median = sorted[sorted.size() / 2]; + + for (auto& value : sorted) { + value = o2::gpu::GPUCommonMath::Abs(value - median); + } + std::sort(sorted.begin(), sorted.end()); + const auto mad = sorted[sorted.size() / 2]; + if (!std::isfinite(mad) || mad <= constants::Tolerance) { + return MinScale2; + } + return std::max(MinScale2, MedianToSigma * mad); +} + +class VertexFit +{ + public: + void add(const Line& line, const float weight) noexcept + { + const auto& direction = line.cosinesDirector; + const auto& origin = line.originPoint; + const auto det = ROOT::Math::Dot(direction, direction); + if (det <= constants::Tolerance) { + return; + } + + for (int i = 0; i < 3; ++i) { + for (int j = i; j < 3; ++j) { + mMatrix(i, j) += weight * (((i == j ? det : 0.f) - direction(i) * direction(j)) / det); + } + } + + const auto dDotO = ROOT::Math::Dot(direction, origin); + for (int i = 0; i < 3; ++i) { + mRhs(i) += weight * ((direction(i) * dDotO - det * origin(i)) / det); + } + } + + bool solve(std::array& vertexOut) const noexcept + { + SymMatrix3 inv{mMatrix}; + if (!inv.InvertFast()) { + return false; + } + const auto solution = inv * mRhs; + vertexOut[0] = static_cast(-solution(0)); + vertexOut[1] = static_cast(-solution(1)); + vertexOut[2] = static_cast(-solution(2)); + return std::isfinite(vertexOut[0]) && std::isfinite(vertexOut[1]) && std::isfinite(vertexOut[2]); + } + + private: + SymMatrix3 mMatrix; + SVector3 mRhs; +}; + +VertexSeed fitSeed(const VertexSeed& initialSeed, + std::span members, + std::span lineRefs, + std::span lines, + const std::shared_ptr& mr, + const float pairCut2) +{ + VertexSeed seed{mr}; + seed.vertex = initialSeed.vertex; + seed.time = initialSeed.time; + seed.scale2 = initialSeed.scale2; + seed.valid = false; + seed.contributors.clear(); + seed.assigned.clear(); + if (members.size() < 2) { + return seed; + } + + for (int iteration = 0; iteration < MaxFitIterations; ++iteration) { + VertexFit vertexFit; + TimeEstBC commonTime{}; + bool hasCommonTime = false; + bounded_vector contributors{mr.get()}; + const auto scale2 = std::max(seed.scale2, MinScale2); + const auto tukeyFactor = 1.f / (scale2 * TukeyC2); + + for (const auto lineRefIdx : members) { + const auto lineIdx = lineRefs[lineRefIdx].lineIndex; + const auto& line = lines[lineIdx]; + if (!line.mTime.isCompatible(seed.time)) { + continue; + } + if (hasCommonTime && !line.mTime.isCompatible(commonTime)) { + continue; + } + + const auto chi2 = Line::getDistance2FromPoint(line, seed.vertex) / pairCut2; + auto weight = 1.f - (chi2 * tukeyFactor); + if (weight <= 0.f) { + continue; + } + weight *= weight; + + if (!hasCommonTime) { + commonTime = line.mTime; + hasCommonTime = true; + } else { + commonTime += line.mTime; + } + + contributors.push_back(lineRefIdx); + vertexFit.add(line, weight); + } + + if (!hasCommonTime || contributors.size() < 2) { + return seed; + } + + std::sort(contributors.begin(), contributors.end()); + + std::array updatedVertex{}; + if (!vertexFit.solve(updatedVertex)) { + return seed; + } + + const auto sameContributors = contributors == seed.contributors; + const auto dz = o2::gpu::GPUCommonMath::Abs(updatedVertex[2] - seed.vertex[2]); + const auto oldR2 = (seed.vertex[0] * seed.vertex[0]) + (seed.vertex[1] * seed.vertex[1]); + const auto newR2 = (updatedVertex[0] * updatedVertex[0]) + (updatedVertex[1] * updatedVertex[1]); + const auto dr2 = o2::gpu::GPUCommonMath::Abs(newR2 - oldR2); + + seed.vertex = updatedVertex; + seed.time = commonTime; + bounded_vector updatedChi2s{mr.get()}; + updatedChi2s.reserve(contributors.size()); + for (const auto lineRefIx : contributors) { + updatedChi2s.push_back(Line::getDistance2FromPoint(lines[lineRefs[lineRefIx].lineIndex], seed.vertex) / pairCut2); + } + seed.scale2 = updateScale2(updatedChi2s, mr); + seed.contributors = std::move(contributors); + seed.valid = true; + + if (sameContributors && dz < VertexShiftZTol && dr2 < VertexShiftR2Tol) { + break; + } + } + + return seed; +} + +size_t countSharedContributors(std::span lhs, std::span rhs) noexcept +{ + size_t shared = 0; + auto lhsIt = lhs.begin(); + auto rhsIt = rhs.begin(); + while (lhsIt != lhs.end() && rhsIt != rhs.end()) { + if (*lhsIt == *rhsIt) { + ++shared; + ++lhsIt; + ++rhsIt; + } else if (*lhsIt < *rhsIt) { + ++lhsIt; + } else { + ++rhsIt; + } + } + return shared; +} + +bounded_vector collectCompatibleContributors(const VertexSeed& seed, + std::span members, + std::span lineRefs, + std::span lines, + const std::shared_ptr& mr, + const float pairCut2) +{ + bounded_vector contributors{mr.get()}; + contributors.reserve(members.size()); + for (const auto lineRefIdx : members) { + const auto lineIdx = lineRefs[lineRefIdx].lineIndex; + const auto& line = lines[lineIdx]; + if (!line.mTime.isCompatible(seed.time)) { + continue; + } + if (Line::getDistance2FromPoint(line, seed.vertex) >= pairCut2) { + continue; + } + contributors.push_back(lineRefIdx); + } + std::sort(contributors.begin(), contributors.end()); + return contributors; +} + +void deduplicateSeeds(bounded_vector& seeds, const Settings& settings) +{ + if (seeds.size() < 2) { + return; + } + + std::sort(seeds.begin(), seeds.end(), [](const VertexSeed& lhs, const VertexSeed& rhs) { + if (lhs.contributors.size() != rhs.contributors.size()) { + return lhs.contributors.size() > rhs.contributors.size(); + } + if (o2::gpu::GPUCommonMath::Abs(lhs.scale2 - rhs.scale2) > constants::Tolerance) { + return lhs.scale2 < rhs.scale2; + } + return lhs.vertex[2] < rhs.vertex[2]; + }); + + const auto dedupZCut = settings.seedDedupZCut > 0.f ? settings.seedDedupZCut : 0.25f * settings.clusterCut; + for (size_t i = 0; i < seeds.size(); ++i) { + auto& candidate = seeds[i]; + if (!candidate.isUsableSeed()) { + candidate.valid = false; + continue; + } + bool duplicate = false; + for (size_t j = 0; j < i; ++j) { + const auto& kept = seeds[j]; + if (!kept.isUsableSeed()) { + continue; + } + if (!candidate.time.isCompatible(kept.time)) { + continue; + } + const auto shared = countSharedContributors(candidate.contributors, kept.contributors); + const auto minSize = std::min(candidate.contributors.size(), kept.contributors.size()); + const auto zDelta = o2::gpu::GPUCommonMath::Abs(candidate.vertex[2] - kept.vertex[2]); + const bool clearlyWorse = kept.contributors.size() > candidate.contributors.size() || + kept.scale2 + constants::Tolerance < 0.9f * candidate.scale2; + const bool overlapDuplicate = shared > 0 && shared * 2 >= minSize; + const bool nearbyDuplicate = zDelta < dedupZCut && (shared > 0 || clearlyWorse); + if (overlapDuplicate || nearbyDuplicate) { + duplicate = true; + break; + } + } + if (duplicate) { + candidate.valid = false; + } + } + compactSeeds(seeds); +} + +void deduplicateRefittedSeeds(bounded_vector& seeds, const Settings& settings) +{ + if (seeds.size() < 2) { + return; + } + + std::sort(seeds.begin(), seeds.end(), [](const VertexSeed& lhs, const VertexSeed& rhs) { + if (lhs.contributors.size() != rhs.contributors.size()) { + return lhs.contributors.size() > rhs.contributors.size(); + } + if (o2::gpu::GPUCommonMath::Abs(lhs.scale2 - rhs.scale2) > constants::Tolerance) { + return lhs.scale2 < rhs.scale2; + } + return lhs.vertex[2] < rhs.vertex[2]; + }); + + const auto zCut = settings.refitDedupZCut > 0.f ? settings.refitDedupZCut : 0.25f * settings.clusterCut; + for (size_t i = 0; i < seeds.size(); ++i) { + auto& candidate = seeds[i]; + if (!candidate.isUsableSeed()) { + candidate.valid = false; + continue; + } + bool duplicate = false; + for (size_t j = 0; j < i; ++j) { + const auto& kept = seeds[j]; + if (!kept.isUsableSeed()) { + continue; + } + if (!candidate.time.isCompatible(kept.time)) { + continue; + } + const auto shared = countSharedContributors(candidate.contributors, kept.contributors); + const auto minSize = std::min(candidate.contributors.size(), kept.contributors.size()); + const auto zDelta = o2::gpu::GPUCommonMath::Abs(candidate.vertex[2] - kept.vertex[2]); + const bool overlapDuplicate = shared > 0 && shared * 2 >= minSize; + const bool lowSupportPair = std::min(candidate.contributors.size(), kept.contributors.size()) < 4; + const bool clearlyWorse = kept.contributors.size() > candidate.contributors.size() || + kept.scale2 + constants::Tolerance < 0.9f * candidate.scale2; + const bool geometricDuplicate = zDelta < zCut && (lowSupportPair || clearlyWorse); + if (overlapDuplicate || geometricDuplicate) { + duplicate = true; + break; + } + } + if (duplicate) { + candidate.valid = false; + } + } + compactSeeds(seeds); +} + +struct OrderedComponent { + explicit OrderedComponent(const std::shared_ptr& mr) : members(mr.get()) {} + float center = 0.f; + bounded_vector members; +}; + +bounded_vector> buildCoarseClusters(std::span lineRefs, + std::span lines, + const Settings& settings) +{ + bounded_vector> clusters(settings.memoryPool.get()); + if (lineRefs.size() < 2) { + return clusters; + } + + bounded_vector sortedByLower(lineRefs.size(), settings.memoryPool.get()); + std::iota(sortedByLower.begin(), sortedByLower.end(), 0); + std::sort(sortedByLower.begin(), sortedByLower.end(), [&](const int lhs, const int rhs) { + const auto lhsLower = lines[lineRefs[lhs].lineIndex].mTime.lower(); + const auto rhsLower = lines[lineRefs[rhs].lineIndex].mTime.lower(); + if (lhsLower != rhsLower) { + return lhsLower < rhsLower; + } + return lineRefs[lhs].lineIndex < lineRefs[rhs].lineIndex; + }); + + const auto coarseZWindow = settings.coarseZWindow > 0.f ? settings.coarseZWindow : settings.clusterCut; + bounded_vector parent(lineRefs.size(), settings.memoryPool.get()); + bounded_vector componentSize(lineRefs.size(), 1, settings.memoryPool.get()); + std::iota(parent.begin(), parent.end(), 0); + float minZ = std::numeric_limits::max(); + float maxZ = std::numeric_limits::lowest(); + for (const auto& lineRef : lineRefs) { + minZ = std::min(minZ, lineRef.zBeam); + maxZ = std::max(maxZ, lineRef.zBeam); + } + const auto nZBins = std::max(1, 1 + static_cast((maxZ - minZ) / coarseZWindow)); + auto getZBin = [&](const float z) { + return std::clamp(static_cast((z - minZ) / coarseZWindow), 0, nZBins - 1); + }; + + auto findRoot = [&](int idx) { + int root = idx; + while (parent[root] != root) { + root = parent[root]; + } + while (parent[idx] != idx) { + const auto next = parent[idx]; + parent[idx] = root; + idx = next; + } + return root; + }; + + auto unite = [&](const int lhs, const int rhs) { + auto lhsRoot = findRoot(lhs); + auto rhsRoot = findRoot(rhs); + if (lhsRoot == rhsRoot) { + return; + } + if (componentSize[lhsRoot] < componentSize[rhsRoot]) { + std::swap(lhsRoot, rhsRoot); + } + parent[rhsRoot] = lhsRoot; + componentSize[lhsRoot] += componentSize[rhsRoot]; + }; + + using ActiveEntry = std::pair; + bounded_vector activeEntries(settings.memoryPool.get()); + std::priority_queue, std::greater> activeByUpper(std::greater{}, std::move(activeEntries)); + bounded_vector activeMask(lineRefs.size(), 0, settings.memoryPool.get()); + bounded_vector> activeByZBin(settings.memoryPool.get()); + activeByZBin.reserve(nZBins); + for (int iBin = 0; iBin < nZBins; ++iBin) { + activeByZBin.emplace_back(); + } + for (const auto lineRefIdx : sortedByLower) { + const auto& lineRef = lineRefs[lineRefIdx]; + const auto& line = lines[lineRef.lineIndex]; + const auto currentLower = line.mTime.lower(); + + while (!activeByUpper.empty() && activeByUpper.top().first < currentLower) { + activeMask[activeByUpper.top().second] = 0; + activeByUpper.pop(); + } + + const auto zBin = getZBin(lineRef.zBeam); + for (int neighborBin = std::max(0, zBin - 1); neighborBin <= std::min(nZBins - 1, zBin + 1); ++neighborBin) { + auto& bucket = activeByZBin[neighborBin]; + size_t writePos = 0; + for (size_t readPos = 0; readPos < bucket.size(); ++readPos) { + const auto oLineRefIdx = bucket[readPos]; + if (!activeMask[oLineRefIdx]) { + continue; + } + bucket[writePos++] = oLineRefIdx; + const auto& oLineRef = lineRefs[oLineRefIdx]; + if (o2::gpu::GPUCommonMath::Abs(lineRef.zBeam - oLineRef.zBeam) >= coarseZWindow) { + continue; + } + const auto& otherLine = lines[oLineRef.lineIndex]; + if (line.mTime.isCompatible(otherLine.mTime)) { + unite(lineRefIdx, oLineRefIdx); + } + } + bucket.resize(writePos); + } + + activeMask[lineRefIdx] = 1; + activeByUpper.emplace(line.mTime.upper(), lineRefIdx); + activeByZBin[zBin].push_back(lineRefIdx); + } + + std::unordered_map> components; + components.reserve(lineRefs.size()); + for (int lineRefIdx = 0; lineRefIdx < static_cast(lineRefs.size()); ++lineRefIdx) { + const auto root = findRoot(lineRefIdx); + auto [it, inserted] = components.try_emplace(root, std::pmr::polymorphic_allocator{settings.memoryPool.get()}); + (void)inserted; + it->second.push_back(lineRefIdx); + } + + bounded_vector orderedComponents(settings.memoryPool.get()); + orderedComponents.reserve(components.size()); + for (auto& [root, members] : components) { + (void)root; + if (members.size() < 2) { + continue; + } + std::sort(members.begin(), members.end(), [&](const int lhs, const int rhs) { + const auto lhsLower = lines[lineRefs[lhs].lineIndex].mTime.lower(); + const auto rhsLower = lines[lineRefs[rhs].lineIndex].mTime.lower(); + if (lhsLower != rhsLower) { + return lhsLower < rhsLower; + } + return lineRefs[lhs].lineIndex < lineRefs[rhs].lineIndex; + }); + orderedComponents.emplace_back(settings.memoryPool); + orderedComponents.back().center = lineRefs[members.front()].tCenter; + orderedComponents.back().members = std::move(members); + } + + std::sort(orderedComponents.begin(), orderedComponents.end(), [](const auto& lhs, const auto& rhs) { + if (o2::gpu::GPUCommonMath::Abs(lhs.center - rhs.center) > TieTolerance) { + return lhs.center < rhs.center; + } + return lhs.members.front() < rhs.members.front(); + }); + clusters.reserve(orderedComponents.size()); + for (auto& component : orderedComponents) { + clusters.push_back(std::move(component.members)); + } + return clusters; +} + +bounded_vector buildSeeds(std::span members, + std::span lineRefs, + std::span lines, + const Settings& settings) +{ + SeedHistogram histogram(members, lineRefs, lines, settings); + bounded_vector seeds(settings.memoryPool.get()); + seeds.reserve(MaxSeedsPerCluster); + float leadingPeakSupport = 0.f; + + while (static_cast(seeds.size()) < MaxSeedsPerCluster) { + const auto peak = histogram.findPeakBin(); + if (peak < 0) { + break; + } + const auto peakSupport = histogram.getPeakSupport(peak); + if (peakSupport < 2.f) { + break; + } + if (leadingPeakSupport <= 0.f) { + leadingPeakSupport = peakSupport; + } else if (peakSupport < std::max(2.f, MinRelativePeakSupport * leadingPeakSupport)) { + break; + } + auto localMembers = histogram.collectLocalMembers(peak, 0, 0); + if (localMembers.size() < 2) { + localMembers = histogram.collectLocalMembers(peak, settings.seedMemberRadiusTime, settings.seedMemberRadiusZ); + } + if (localMembers.size() < 2) { + histogram.suppressPeak(peak); + continue; + } + + VertexSeed seed(settings.memoryPool); + seed.vertex = {settings.beamX, settings.beamY, histogram.getPeakZCenter(peak)}; + seed.time = histogram.getPeakTimeInterval(peak); + seed.scale2 = InitialScale2; + + auto fitted = fitSeed(seed, localMembers, lineRefs, lines, settings.memoryPool, settings.pairCut2); + if (fitted.valid && fitted.contributors.size() >= 2) { + seeds.push_back(std::move(fitted)); + histogram.suppressPeakNeighborhood(peak); + } else { + histogram.suppressPeak(peak); + } + } + + return seeds; +} + +void assignLinesToSeeds(bounded_vector& seeds, + std::span members, + std::span lineRefs, + std::span lines, + const float pairCut2) +{ + for (auto& seed : seeds) { + seed.assigned.clear(); + } + + for (const auto lineRefIdx : members) { + const auto lineIdx = lineRefs[lineRefIdx].lineIndex; + const auto& line = lines[lineIdx]; + + int bestSeed = -1; + float bestScore = std::numeric_limits::max(); + size_t bestMult = 0; + float bestZResidual = std::numeric_limits::max(); + + for (int seedIdx = 0; seedIdx < static_cast(seeds.size()); ++seedIdx) { + const auto& seed = seeds[seedIdx]; + if (!seed.valid || seed.contributors.size() < 2) { + continue; + } + if (!line.mTime.isCompatible(seed.time)) { + continue; + } + + const auto distance2 = Line::getDistance2FromPoint(line, seed.vertex); + if (distance2 >= pairCut2) { + continue; + } + + const auto score = distance2 / std::max(seed.scale2, MinScale2); + const auto zResidual = o2::gpu::GPUCommonMath::Abs(lineRefs[lineRefIdx].zBeam - seed.vertex[2]); + const auto multiplicity = seed.contributors.size(); + + const auto betterScore = score + TieTolerance < bestScore; + const auto betterMultiplicity = o2::gpu::GPUCommonMath::Abs(score - bestScore) <= TieTolerance && multiplicity > bestMult; + const auto betterZ = o2::gpu::GPUCommonMath::Abs(score - bestScore) <= TieTolerance && + multiplicity == bestMult && zResidual + constants::Tolerance < bestZResidual; + if (betterScore || betterMultiplicity || betterZ) { + bestSeed = seedIdx; + bestScore = score; + bestMult = multiplicity; + bestZResidual = zResidual; + } + } + + if (bestSeed >= 0) { + seeds[bestSeed].assigned.push_back(lineRefIdx); + } + } +} + +ClusterLines materializeCluster(const VertexSeed& seed, + std::span lineRefs, + std::span lines, + const std::shared_ptr& mr) +{ + bounded_vector lineIndices{mr.get()}; + lineIndices.reserve(seed.contributors.size()); + for (const auto lineRefIdx : seed.contributors) { + lineIndices.push_back(lineRefs[lineRefIdx].lineIndex); + } + std::sort(lineIndices.begin(), lineIndices.end()); + lineIndices.erase(std::unique(lineIndices.begin(), lineIndices.end()), lineIndices.end()); + + if (lineIndices.size() < 2) { + return {}; + } + + return {std::span{lineIndices.data(), lineIndices.size()}, lines}; +} + +} // namespace + +bounded_vector buildClusters(std::span lines, const Settings& settings) +{ + bounded_vector clusters(settings.memoryPool.get()); + if (lines.size() < 2) { + return clusters; + } + + bounded_vector refs(settings.memoryPool.get()); + refs.reserve(lines.size()); + for (int lineIdx = 0; lineIdx < static_cast(lines.size()); ++lineIdx) { + LineRef ref(lines[lineIdx], lineIdx, settings.beamX, settings.beamY, settings.maxZ); + if (!ref.isDead()) { + refs.push_back(ref); + } + } + + if (refs.size() < 2) { + return clusters; + } + + const auto coarseClusters = buildCoarseClusters(refs, lines, settings); + + for (const auto& members : coarseClusters) { + auto seeds = buildSeeds(members, refs, lines, settings); + if (seeds.empty()) { + continue; + } + + for (auto& seed : seeds) { + if (!seed.isUsableSeed()) { + seed.valid = false; + continue; + } + auto contributors = collectCompatibleContributors(seed, members, refs, lines, settings.memoryPool, settings.pairCut2); + if (contributors.size() < 2) { + seed.valid = false; + continue; + } + seed.contributors = std::move(contributors); + } + compactSeeds(seeds); + if (seeds.empty()) { + continue; + } + deduplicateSeeds(seeds, settings); + if (seeds.empty()) { + continue; + } + assignLinesToSeeds(seeds, members, refs, lines, settings.pairCut2); + for (auto& seed : seeds) { + if (seed.assigned.size() < 2) { + seed.valid = false; + continue; + } + seed = fitSeed(seed, seed.assigned, refs, lines, settings.memoryPool, settings.pairCut2); + if (!seed.isUsableSeed()) { + seed.valid = false; + continue; + } + } + compactSeeds(seeds); + deduplicateRefittedSeeds(seeds, settings); + for (auto& refit : seeds) { + auto cluster = materializeCluster(refit, refs, lines, settings.memoryPool); + if (cluster.getSize() < 2) { + continue; + } + if (!cluster.isValid()) { + continue; + } + clusters.push_back(std::move(cluster)); + } + } + + return clusters; +} + +} // namespace o2::its::line_vertexer diff --git a/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx b/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx index 5a32b3d3b1a95..5b412ea4eea69 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx @@ -14,19 +14,16 @@ /// #include -#include #include "Framework/Logger.h" #include "ITStracking/TimeFrame.h" #include "ITStracking/MathUtils.h" -#include "DataFormatsITSMFT/Cluster.h" #include "DataFormatsITSMFT/CompCluster.h" #include "DataFormatsITSMFT/ROFRecord.h" #include "DataFormatsITSMFT/TopologyDictionary.h" #include "ITSBase/GeometryTGeo.h" #include "ITSMFTBase/SegmentationAlpide.h" #include "ITStracking/BoundedAllocator.h" -#include "ITStracking/TrackingConfigParam.h" namespace { @@ -271,7 +268,7 @@ void TimeFrame::initialise(const int iteration, const TrackingParameter 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.SystErrorZ2[iLayer] + trkParam.SystErrorY2[iLayer]) + trkParam.LayerResolution[iLayer] * trkParam.LayerResolution[iLayer]); + mPositionResolution[iLayer] = o2::gpu::CAMath::Sqrt((0.5f * (trkParam.SystErrorZ2[iLayer] + trkParam.SystErrorY2[iLayer])) + (trkParam.LayerResolution[iLayer] * trkParam.LayerResolution[iLayer])); } clearResizeBoundedVector(mLines, getNrof(1), mMemoryPool.get()); clearResizeBoundedVector(mTrackletClusters, getNrof(1), mMemoryPool.get()); @@ -312,17 +309,17 @@ void TimeFrame::initialise(const int iteration, const TrackingParameter float oneOverR{0.001f * 0.3f * std::abs(mBz) / trkParam.TrackletMinPt}; for (unsigned int iLayer{0}; iLayer < NLayers; ++iLayer) { mMSangles[iLayer] = math_utils::MSangle(0.14f, trkParam.TrackletMinPt, trkParam.LayerxX0[iLayer]); - mPositionResolution[iLayer] = o2::gpu::CAMath::Sqrt(0.5f * (trkParam.SystErrorZ2[iLayer] + trkParam.SystErrorY2[iLayer]) + trkParam.LayerResolution[iLayer] * trkParam.LayerResolution[iLayer]); + mPositionResolution[iLayer] = o2::gpu::CAMath::Sqrt((0.5f * (trkParam.SystErrorZ2[iLayer] + trkParam.SystErrorY2[iLayer])) + (trkParam.LayerResolution[iLayer] * trkParam.LayerResolution[iLayer])); if (iLayer < mClusters.size() - 1) { const float& r1 = trkParam.LayerRadii[iLayer]; const float& r2 = trkParam.LayerRadii[iLayer + 1]; - oneOverR = (0.5 * oneOverR >= 1.f / r2) ? 2.f / r2 - o2::constants::math::Almost0 : oneOverR; + 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[iLayer]); const float res2 = o2::gpu::CAMath::Hypot(trkParam.PVres, mPositionResolution[iLayer + 1]); 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))); + 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) mPhiCuts[iLayer] = std::min(o2::gpu::CAMath::ASin(0.5f * x * oneOverR) + 2.f * mMSangles[iLayer] + delta, o2::constants::math::PI * 0.5f); } diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx b/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx index a41560e2e9e9a..eb0841888b03e 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx @@ -10,6 +10,8 @@ // or submit itself to any jurisdiction. #include +#include +#include #include #include @@ -128,12 +130,12 @@ void ITSTrackingInterface::run(framework::ProcessingContext& pc) gsl::span physTriggers; std::vector fromTRD; if (mUseTriggers == 2) { // use TRD triggers - o2::InteractionRecord ir{0, tfInfo.firstTForbit}; + o2::InteractionRecord irFirstTF{0, tfInfo.firstTForbit}; auto trdTriggers = pc.inputs().get>("phystrig"); for (const auto& trig : trdTriggers) { - if (trig.getBCData() >= ir && trig.getNumberOfTracklets()) { - ir = trig.getBCData(); - fromTRD.emplace_back(o2::itsmft::PhysTrigger{.ir = ir, .data = 0}); + if (trig.getBCData() >= irFirstTF && trig.getNumberOfTracklets()) { + irFirstTF = trig.getBCData(); + fromTRD.emplace_back(o2::itsmft::PhysTrigger{.ir = irFirstTF, .data = 0}); } } physTriggers = gsl::span(fromTRD.data(), fromTRD.size()); @@ -215,7 +217,8 @@ void ITSTrackingInterface::run(framework::ProcessingContext& pc) auto vtxSpan = mTimeFrame->getPrimaryVertices(clockLayerId, iRof); if (o2::its::TrackerParamConfig::Instance().doUPCIteration) { if (!vtxSpan.empty()) { - if (vtxSpan[0].isFlagSet(Vertex::UPCMode) == 1) { // at least one vertex in this ROF and it is from second vertex iteration + bool hasUPC = std::any_of(vtxSpan.begin(), vtxSpan.end(), [](const auto& v) { return v.isFlagSet(Vertex::UPCMode); }); + if (hasUPC) { // at least one vertex in this ROF and it is from second vertex iteration LOGP(debug, "ROF {} rejected as vertices are from the UPC iteration", iRof); processUPCMask.selectROF({clockTiming.getROFStartInBC(iRof), clockTiming.getROFEndInBC(iRof)}); vtxROF.setFlag(o2::itsmft::ROFRecord::VtxUPCMode); diff --git a/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx b/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx index 5e27e20b3ddee..a22d2d6c60990 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx @@ -1,4 +1,4 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// Copyright 2019-2026 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. // @@ -13,6 +13,8 @@ #include #include #include +#include +#include #include #include @@ -22,6 +24,7 @@ #include "ITStracking/BoundedAllocator.h" #include "ITStracking/ClusterLines.h" #include "ITStracking/Definitions.h" +#include "ITStracking/LineVertexerHelpers.h" #include "ITStracking/Tracklet.h" #include "SimulationDataFormat/DigitizationContext.h" #include "SimulationDataFormat/O2DatabasePDG.h" @@ -31,13 +34,14 @@ namespace o2::its { - +namespace +{ template -static void trackleterKernelHost( +void trackleterKernelHost( const gsl::span& clustersNextLayer, // 0 2 const gsl::span& clustersCurrentLayer, // 1 1 const gsl::span& usedClustersNextLayer, // 0 2 - int* indexTableNext, + const int* indexTableNext, const float phiCut, bounded_vector& tracklets, gsl::span foundTracklets, @@ -94,7 +98,7 @@ static void trackleterKernelHost( } } -static void trackletSelectionKernelHost( +void trackletSelectionKernelHost( const Cluster* clusters0, // global layer 0 clusters const Cluster* clusters1, // global layer 1 clusters gsl::span usedClusters0, // global layer 0 used clusters @@ -145,6 +149,7 @@ static void trackletSelectionKernelHost( offset12 += foundTracklets12[iCurrentLayerClusterIndex]; } } +} // namespace template void VertexerTraits::updateVertexingParameters(const std::vector& vrtPar) @@ -255,7 +260,7 @@ void VertexerTraits::computeTracklets(const int iteration) }); }); - /// Create tracklets labels for L0-L1, information is as flat as in tracklets vector (no rofId) + /// Create flat L0-L1 tracklet labels (no rofId) if (mTimeFrame->hasMCinformation()) { for (const auto& trk : mTimeFrame->getTracklets()[0]) { o2::MCCompLabel label; @@ -309,15 +314,7 @@ void VertexerTraits::computeTrackletMatching(const int iteration) static_cast(mTimeFrame->getClustersOnLayer(pivotRofId, 1).size()), mVrtParams[iteration].tanLambdaCut, mVrtParams[iteration].phiCut); - auto& lines = mTimeFrame->getLines(pivotRofId); - totalLines.local() += lines.size(); - std::stable_sort(lines.begin(), lines.end(), [](const Line& a, const Line& b) { - // sort by lower edge and secondly prefer wider windows - if (a.mTime.lower() != b.mTime.lower()) { - return a.mTime.lower() < b.mTime.lower(); - } - return a.mTime.upper() > b.mTime.upper(); - }); + totalLines.local() += mTimeFrame->getLines(pivotRofId).size(); } }); mTimeFrame->setNLinesTotal(totalLines.combine(std::plus())); @@ -330,125 +327,214 @@ void VertexerTraits::computeTrackletMatching(const int iteration) template void VertexerTraits::computeVertices(const int iteration) { - const auto nsigmaCut{std::min(mVrtParams[iteration].vertNsigmaCut * mVrtParams[iteration].vertNsigmaCut * (mVrtParams[iteration].vertRadiusSigma * mVrtParams[iteration].vertRadiusSigma + mVrtParams[iteration].trackletSigma * mVrtParams[iteration].trackletSigma), 1.98f)}; - const auto pairCut2{mVrtParams[iteration].pairCut * mVrtParams[iteration].pairCut}; const int nRofs = mTimeFrame->getNrof(1); - const bool hasMC = mTimeFrame->hasMCinformation(); std::vector> rofVertices(nRofs); std::vector> rofLabels(nRofs); + const float nsigmaCut = std::min(mVrtParams[iteration].vertNsigmaCut * mVrtParams[iteration].vertNsigmaCut * (mVrtParams[iteration].vertRadiusSigma * mVrtParams[iteration].vertRadiusSigma + mVrtParams[iteration].trackletSigma * mVrtParams[iteration].trackletSigma), 1.98f); + const float pairCut2 = mVrtParams[iteration].pairCut * mVrtParams[iteration].pairCut; + const float duplicateZCut = mVrtParams[iteration].duplicateZCut > 0.f ? mVrtParams[iteration].duplicateZCut : std::max(4.f * mVrtParams[iteration].pairCut, 0.5f * mVrtParams[iteration].clusterCut); + const float duplicateDistance2Cut = mVrtParams[iteration].duplicateDistance2Cut > 0.f ? mVrtParams[iteration].duplicateDistance2Cut : std::max(16.f * pairCut2, 0.0625f * mVrtParams[iteration].clusterCut * mVrtParams[iteration].clusterCut); + line_vertexer::Settings settings; + settings.beamX = mTimeFrame->getBeamX(); + settings.beamY = mTimeFrame->getBeamY(); + settings.pairCut = mVrtParams[iteration].pairCut; + settings.pairCut2 = pairCut2; + settings.clusterCut = mVrtParams[iteration].clusterCut; + settings.coarseZWindow = mVrtParams[iteration].coarseZWindow; + settings.seedDedupZCut = mVrtParams[iteration].seedDedupZCut; + settings.refitDedupZCut = mVrtParams[iteration].refitDedupZCut; + settings.duplicateZCut = duplicateZCut; + settings.duplicateDistance2Cut = duplicateDistance2Cut; + settings.finalSelectionZCut = mVrtParams[iteration].finalSelectionZCut; + settings.maxZ = mVrtParams[iteration].maxZPositionAllowed; + settings.seedMemberRadiusTime = mVrtParams[iteration].seedMemberRadiusTime; + settings.seedMemberRadiusZ = mVrtParams[iteration].seedMemberRadiusZ; + settings.memoryPool = mMemoryPool; const auto processROF = [&](const int rofId) { auto& lines = mTimeFrame->getLines(rofId); - const int nLines{static_cast(lines.size())}; - bounded_vector usedTracklets(nLines, 0, mMemoryPool.get()); - auto& clusters = mTimeFrame->getTrackletClusters(rofId); - - for (int iLine1{0}; iLine1 < nLines; ++iLine1) { - if (usedTracklets[iLine1]) { - continue; + auto clusters = line_vertexer::buildClusters(std::span{lines.data(), lines.size()}, settings); + deepVectorClear(lines); // not needed after + auto clusterBeamDistance2 = [&](const ClusterLines& cluster) { + return (mTimeFrame->getBeamX() - cluster.getVertex()[0]) * (mTimeFrame->getBeamX() - cluster.getVertex()[0]) + + (mTimeFrame->getBeamY() - cluster.getVertex()[1]) * (mTimeFrame->getBeamY() - cluster.getVertex()[1]); + }; + auto clusterBetter = [&](const ClusterLines& lhs, const ClusterLines& rhs) { + if (lhs.getSize() != rhs.getSize()) { + return lhs.getSize() > rhs.getSize(); } - const auto& line1 = lines[iLine1]; - for (int iLine2{iLine1 + 1}; iLine2 < nLines; ++iLine2) { - if (usedTracklets[iLine2]) { - continue; - } - const auto& line2 = lines[iLine2]; - if (!line1.mTime.isCompatible(line2.mTime)) { + if (o2::gpu::GPUCommonMath::Abs(lhs.getAvgDistance2() - rhs.getAvgDistance2()) > constants::Tolerance) { + return lhs.getAvgDistance2() < rhs.getAvgDistance2(); + } + const auto lhsBeam = clusterBeamDistance2(lhs); + const auto rhsBeam = clusterBeamDistance2(rhs); + if (o2::gpu::GPUCommonMath::Abs(lhsBeam - rhsBeam) > constants::Tolerance) { + return lhsBeam < rhsBeam; + } + return lhs.getVertex()[2] < rhs.getVertex()[2]; + }; + + // Cluster deduplication by local non-maximum suppression in time/space + std::sort(clusters.begin(), clusters.end(), clusterBetter); + float minClusterZ = std::numeric_limits::max(); + for (const auto& cluster : clusters) { + minClusterZ = std::min(minClusterZ, cluster.getVertex()[2]); + } + bounded_vector deduplicated(mMemoryPool.get()); + deduplicated.reserve(clusters.size()); + std::unordered_map> keptByZBin; + for (auto& candidate : clusters) { + bool duplicate = false; + const auto candidateZ = candidate.getVertex()[2]; + const auto zBin = static_cast(std::floor((candidateZ - minClusterZ) / settings.duplicateZCut)); + for (int neighborBin = zBin - 1; neighborBin <= zBin + 1 && !duplicate; ++neighborBin) { + const auto found = keptByZBin.find(neighborBin); + if (found == keptByZBin.end()) { continue; } - auto dca2{Line::getDCA2(line1, line2)}; - if (dca2 < pairCut2) { - auto& cluster = clusters.emplace_back(iLine1, line1, iLine2, line2); - if (!cluster.isValid() || cluster.getR2() > 4.f) { - clusters.pop_back(); + for (const auto ownerId : found->second) { + const auto& owner = deduplicated[ownerId]; + if (!candidate.getTimeStamp().isCompatible(owner.getTimeStamp())) { continue; } - - usedTracklets[iLine1] = 1; - usedTracklets[iLine2] = 1; - for (int iLine3{0}; iLine3 < nLines; ++iLine3) { - if (usedTracklets[iLine3]) { - continue; - } - const auto& line3 = lines[iLine3]; - if (!line3.mTime.isCompatible(cluster.getTimeStamp())) { - continue; - } - const auto distance2 = Line::getDistance2FromPoint(line3, cluster.getVertex()); - if (distance2 < pairCut2) { - cluster.add(iLine3, line3); - usedTracklets[iLine3] = 1; - } + if (o2::gpu::GPUCommonMath::Abs(candidate.getVertex()[2] - owner.getVertex()[2]) >= settings.duplicateZCut) { + continue; + } + const auto dx = candidate.getVertex()[0] - owner.getVertex()[0]; + const auto dy = candidate.getVertex()[1] - owner.getVertex()[1]; + const auto dz = candidate.getVertex()[2] - owner.getVertex()[2]; + const auto distance2 = math_utils::SqSum(dx, dy, dz); + if (distance2 < settings.duplicateDistance2Cut) { + duplicate = true; + break; } - break; } } - } + if (duplicate) { + continue; + } - // Cluster merging - std::sort(clusters.begin(), clusters.end(), - [](ClusterLines& cluster1, ClusterLines& cluster2) { return cluster1.getSize() > cluster2.getSize(); }); + const auto ownerId = static_cast(deduplicated.size()); + keptByZBin[zBin].push_back(ownerId); + deduplicated.push_back(std::move(candidate)); + } + clusters = std::move(deduplicated); int nClusters = static_cast(clusters.size()); - for (int iCluster1{0}; iCluster1 < nClusters; ++iCluster1) { - std::array vertex1{clusters[iCluster1].getVertex()}; - std::array vertex2{}; - for (int iCluster2{iCluster1 + 1}; iCluster2 < nClusters; ++iCluster2) { - if (clusters[iCluster1].getTimeStamp().isCompatible(clusters[iCluster2].getTimeStamp())) { - vertex2 = clusters[iCluster2].getVertex(); - if (o2::gpu::GPUCommonMath::Abs(vertex1[2] - vertex2[2]) < mVrtParams[iteration].clusterCut) { - float distance{((vertex1[0] - vertex2[0]) * (vertex1[0] - vertex2[0])) + - ((vertex1[1] - vertex2[1]) * (vertex1[1] - vertex2[1])) + - ((vertex1[2] - vertex2[2]) * (vertex1[2] - vertex2[2]))}; - if (distance < mVrtParams[iteration].pairCut * mVrtParams[iteration].pairCut) { - for (auto label : clusters[iCluster2].getLabels()) { - clusters[iCluster1].add(label, lines[label]); - vertex1 = clusters[iCluster1].getVertex(); - } - clusters.erase(clusters.begin() + iCluster2); - --iCluster2; - --nClusters; - } - } - } + + // Vertex filtering with score-based local NMS + std::sort(clusters.begin(), clusters.end(), clusterBetter); + std::vector candidateIndices; + candidateIndices.reserve(nClusters); + for (int iCluster{0}; iCluster < nClusters; ++iCluster) { + const bool zCompatible = o2::gpu::GPUCommonMath::Abs(clusters[iCluster].getVertex()[2]) < mVrtParams[iteration].maxZPositionAllowed; + + if (zCompatible) { + candidateIndices.push_back(iCluster); } } - // Vertex filtering - std::sort(clusters.begin(), clusters.end(), - [](const ClusterLines& cluster1, const ClusterLines& cluster2) { return cluster1.getSize() > cluster2.getSize(); }); - bool atLeastOneFound{false}; - for (int iCluster{0}; iCluster < nClusters; ++iCluster) { - bool lowMultCandidate{false}; - double beamDistance2{(mTimeFrame->getBeamX() - clusters[iCluster].getVertex()[0]) * (mTimeFrame->getBeamX() - clusters[iCluster].getVertex()[0]) + - (mTimeFrame->getBeamY() - clusters[iCluster].getVertex()[1]) * (mTimeFrame->getBeamY() - clusters[iCluster].getVertex()[1])}; - if (atLeastOneFound && (lowMultCandidate = clusters[iCluster].getSize() < mVrtParams[iteration].clusterContributorsCut)) { - lowMultCandidate &= (beamDistance2 < mVrtParams[iteration].lowMultBeamDistCut * mVrtParams[iteration].lowMultBeamDistCut); - if (!lowMultCandidate) { - clusters.erase(clusters.begin() + iCluster); - nClusters--; - continue; + if (candidateIndices.empty()) { + return; + } + + auto countSharedLabels = [](const ClusterLines& lhs, const ClusterLines& rhs) { + size_t shared = 0; + auto lhsIt = lhs.getLabels().begin(); + auto rhsIt = rhs.getLabels().begin(); + while (lhsIt != lhs.getLabels().end() && rhsIt != rhs.getLabels().end()) { + if (*lhsIt == *rhsIt) { + ++shared; + ++lhsIt; + ++rhsIt; + } else if (*lhsIt < *rhsIt) { + ++lhsIt; + } else { + ++rhsIt; } } + return shared; + }; - if (beamDistance2 < nsigmaCut && o2::gpu::GPUCommonMath::Abs(clusters[iCluster].getVertex()[2]) < mVrtParams[iteration].maxZPositionAllowed) { - atLeastOneFound = true; - Vertex vertex{clusters[iCluster].getVertex().data(), - clusters[iCluster].getRMS2(), - (ushort)clusters[iCluster].getSize(), - clusters[iCluster].getAvgDistance2()}; - - if (iteration) { - vertex.setFlags(Vertex::UPCMode); + float minCandidateZ = std::numeric_limits::max(); + for (const auto clusterId : candidateIndices) { + minCandidateZ = std::min(minCandidateZ, clusters[clusterId].getVertex()[2]); + } + std::unordered_map> selectedByZBin; + std::vector selectedIndices; + selectedIndices.reserve(candidateIndices.size()); + for (const auto clusterId : candidateIndices) { + const auto& candidate = clusters[clusterId]; + const auto candidateZ = candidate.getVertex()[2]; + const auto zBin = static_cast((candidateZ - minCandidateZ) / settings.finalSelectionZCut); + bool suppressed = false; + for (int neighborBin = zBin - 1; neighborBin <= zBin + 1 && !suppressed; ++neighborBin) { + const auto found = selectedByZBin.find(neighborBin); + if (found == selectedByZBin.end()) { + continue; } - vertex.setTimeStamp(clusters[iCluster].getTimeStamp()); - rofVertices[rofId].push_back(vertex); - if (hasMC) { - bounded_vector labels(mMemoryPool.get()); - for (auto& index : clusters[iCluster].getLabels()) { - labels.push_back(mTimeFrame->getLinesLabel(rofId)[index]); + for (const auto selectedId : found->second) { + const auto& selected = clusters[selectedId]; + if (!candidate.getTimeStamp().isCompatible(selected.getTimeStamp())) { + continue; } - rofLabels[rofId].push_back(computeMain(labels)); + const auto zDelta = o2::gpu::GPUCommonMath::Abs(candidateZ - selected.getVertex()[2]); + const auto sharedLabels = countSharedLabels(candidate, selected); + const auto minSize = std::min(candidate.getSize(), selected.getSize()); + const bool overlapDuplicate = sharedLabels > 0 && sharedLabels * 4 >= minSize; + const bool strongZDuplicate = zDelta < settings.finalSelectionZCut; + const bool clearlyBetterMultiplicity = selected.getSize() >= candidate.getSize() + 3; + const bool clearlyBetterQuality = selected.getSize() > candidate.getSize() && + selected.getAvgDistance2() + constants::Tolerance < 0.8f * candidate.getAvgDistance2(); + const bool weakCandidate = clearlyBetterMultiplicity || clearlyBetterQuality; + if (overlapDuplicate || (strongZDuplicate && weakCandidate)) { + suppressed = true; + break; + } + } + } + if (suppressed) { + continue; + } + selectedByZBin[zBin].push_back(clusterId); + selectedIndices.push_back(clusterId); + } + + // sort vertices by their multiplicity to opt. suppress lower mult. debris + std::vector sortedIndices(selectedIndices.size()); + std::iota(sortedIndices.begin(), sortedIndices.end(), 0); + std::sort(sortedIndices.begin(), sortedIndices.end(), [&selectedIndices, &clusters](int i, int j) { + return clusters[selectedIndices[i]].getSize() > clusters[selectedIndices[j]].getSize(); + }); + for (const auto sortedId : sortedIndices) { + const auto& cluster = clusters[selectedIndices[sortedId]]; + const auto beamDistance2 = clusterBeamDistance2(cluster); + if (!(beamDistance2 < nsigmaCut)) { + continue; + } + if (cluster.getSize() < mVrtParams[iteration].clusterContributorsCut) { + continue; + } + if (!rofVertices[rofId].empty() && cluster.getSize() < mVrtParams[iteration].suppressLowMultDebris) { + continue; + } + + Vertex vertex{cluster.getVertex().data(), + cluster.getRMS2(), + (ushort)cluster.getSize(), + cluster.getAvgDistance2()}; + if (iteration) { + vertex.setFlags(Vertex::UPCMode); + } + vertex.setTimeStamp(cluster.getTimeStamp()); + rofVertices[rofId].push_back(vertex); + if (mTimeFrame->hasMCinformation()) { + auto& lineLabels = mTimeFrame->getLinesLabel(rofId); + bounded_vector labels(mMemoryPool.get()); + for (auto& index : cluster.getLabels()) { + labels.push_back(lineLabels[index]); } + const auto mainLabel = computeMain(labels); + rofLabels[rofId].push_back(mainLabel); } } }; @@ -469,7 +555,7 @@ void VertexerTraits::computeVertices(const int iteration) for (auto& vertex : rofVertices[rofId]) { mTimeFrame->addPrimaryVertex(vertex); } - if (hasMC) { + if (mTimeFrame->hasMCinformation()) { for (auto& label : rofLabels[rofId]) { mTimeFrame->addPrimaryVertexLabel(label); } @@ -504,7 +590,8 @@ void VertexerTraits::addTruthSeedingVertices() if (!trk.isPrimary() || trk.GetPt() < 0.05 || std::abs(trk.GetEta()) > 1.1) { return false; } - return o2::O2DatabasePDG::Instance()->GetParticle(trk.GetPdgCode())->Charge() != 0; + const auto* p = o2::O2DatabasePDG::Instance()->GetParticle(trk.GetPdgCode()); + return (!p) ? false : p->Charge() != 0; }))); vert.setXYZ((float)eve.GetX(), (float)eve.GetY(), (float)eve.GetZ()); vert.setChi2(1); // not used as constraint diff --git a/prodtests/full-system-test/dpl-workflow.sh b/prodtests/full-system-test/dpl-workflow.sh index 9f982513fdffd..d54c05ff0f20e 100755 --- a/prodtests/full-system-test/dpl-workflow.sh +++ b/prodtests/full-system-test/dpl-workflow.sh @@ -117,6 +117,11 @@ EVE_OPT=" --jsons-folder $EDJSONS_DIR" [[ "0$ITSSTAGGERED" == "01" ]] && ITS_STAGGERED=" --enable-its-staggering " || ITS_STAGGERED= [[ "0$MFTSTAGGERED" == "01" ]] && MFT_STAGGERED=" --enable-its-staggering " || MFT_STAGGERED= +# ITS vertexing settings +if [[ $BEAMTYPE == "pp" || $LIGHTNUCLEI == "1" ]]; then + ITS_CONFIG_KEY+=";ITSVertexerParam.pairCut=0.0317563;ITSVertexerParam.clusterCut=0.6640964;ITSVertexerParam.coarseZWindow=0.2049018;ITSVertexerParam.seedDedupZCut=0.0711793;ITSVertexerParam.refitDedupZCut=0.0680009;ITSVertexerParam.duplicateZCut=0.1582193;ITSVertexerParam.finalSelectionZCut=0.1081465;ITSVertexerParam.duplicateDistance2Cut=0.0117033;ITSVertexerParam.clusterContributorsCut=2;ITSVertexerParam.seedMemberRadiusZ=0;ITSVertexerParam.vertNsigmaCut=4.0;ITSVertexerParam.vertRadiusSigma=0.0452309;ITSVertexerParam.trackletSigma=0.0025941;ITSVertexerParam.suppressLowMultDebris=0;" +fi + if [[ $CTFINPUT != 1 ]]; then GPU_OUTPUT+=",tpc-triggers" fi @@ -126,10 +131,10 @@ if [[ $SYNCMODE == 1 ]]; then MFT_STF_DEC_CONFIG+="MFTClustererParam.maxBCDiffToMaskBias=-1;" [[ $BEAMTYPE == "PbPb" || $BEAMTYPE == "pp" || $LIGHTNUCLEI == "1" ]] && MFT_CONFIG_KEY+="MFTTracking.cutMultClusLow=0;MFTTracking.cutMultClusHigh=4000;" if [[ $BEAMTYPE == "PbPb" ]]; then - ITS_CONFIG_KEY+="fastMultConfig.cutMultClusLow=${CUT_MULT_MIN_ITS:-0};fastMultConfig.cutMultClusHigh=${CUT_MULT_MAX_ITS:-400};fastMultConfig.cutMultVtxHigh=${CUT_MULT_VTX_ITS:-20};" + ITS_CONFIG_KEY+=";fastMultConfig.cutMultClusLow=${CUT_MULT_MIN_ITS:-0};fastMultConfig.cutMultClusHigh=${CUT_MULT_MAX_ITS:-400};fastMultConfig.cutMultVtxHigh=${CUT_MULT_VTX_ITS:-20};" MCH_CONFIG_KEY="MCHTracking.maxCandidates=50000;MCHTracking.maxTrackingDuration=20;" elif [[ $BEAMTYPE == "pp" || $LIGHTNUCLEI == "1" ]]; then - ITS_CONFIG_KEY+="fastMultConfig.cutMultClusLow=${CUT_MULT_MIN_ITS:--1};fastMultConfig.cutMultClusHigh=${CUT_MULT_MAX_ITS:--1};fastMultConfig.cutMultVtxHigh=${CUT_MULT_VTX_ITS:--1};ITSVertexerParam.phiCut=0.5;ITSVertexerParam.clusterContributorsCut=3;ITSVertexerParam.tanLambdaCut=0.2;" + ITS_CONFIG_KEY+=";fastMultConfig.cutMultClusLow=${CUT_MULT_MIN_ITS:--1};fastMultConfig.cutMultClusHigh=${CUT_MULT_MAX_ITS:--1};fastMultConfig.cutMultVtxHigh=${CUT_MULT_VTX_ITS:--1};" MCH_CONFIG_KEY="MCHTracking.maxCandidates=20000;MCHTracking.maxTrackingDuration=10;" fi [[ -n ${CUT_RANDOM_FRACTION_ITS:-} ]] && ITS_CONFIG_KEY+="fastMultConfig.cutRandomFraction=$CUT_RANDOM_FRACTION_ITS;" @@ -157,11 +162,6 @@ if [[ $SYNCMODE == 1 ]]; then has_detector ITS && TRD_FILTER_CONFIG+=" --filter-trigrec" else has_detectors_gpu TPC ITS && ITS_CONFIG_KEY+="ITSCATrackerParam.trackingMode=1;" # sets ITS gpu reco to async - if [[ $BEAMTYPE == "pp" || $LIGHTNUCLEI == "1" ]]; then - ITS_CONFIG_KEY+="ITSVertexerParam.phiCut=0.5;ITSVertexerParam.clusterContributorsCut=3;ITSVertexerParam.tanLambdaCut=0.2;" - elif [[ $BEAMTYPE == "PbPb" ]]; then - ITS_CONFIG_KEY+="ITSVertexerParam.lowMultBeamDistCut=0;" - fi if [[ $IS_SIMULATED_DATA == 0 && $CTFINPUT == 1 ]]; then # Enable fixes to the MCH readout mapping for async processing of real data MCH_CONFIG_KEY+="MCHDigitModifier.updateST1=true;MCHDigitModifier.updateST2=true;" @@ -585,7 +585,7 @@ has_detector_gpu ITS && GPU_OUTPUT+=",its-tracks" # --------------------------------------------------------------------------------------------------------------------- # Common reconstruction workflows -(has_detector_reco ITS && ! has_detector_gpu ITS) && ! has_detector_from_global_reader ITS && add_W o2-its-reco-workflow "--trackerCA $ITS_CONFIG $ITS_STAGGERED $DISABLE_MC ${DISABLE_DIGIT_CLUSTER_INPUT:-} $DISABLE_ROOT_OUTPUT --pipeline $(get_N its-tracker ITS REST 1 ITSTRK),$(get_N its-clusterer ITS REST 1 ITSCL)" "$ITS_CONFIG_KEY;$ITSMFT_STROBES;$ITSEXTRAERR" +(has_detector_reco ITS && ! has_detector_gpu ITS) && ! has_detector_from_global_reader ITS && add_W o2-its-reco-workflow "$ITS_CONFIG $ITS_STAGGERED $DISABLE_MC ${DISABLE_DIGIT_CLUSTER_INPUT:-} $DISABLE_ROOT_OUTPUT --pipeline $(get_N its-tracker ITS REST 1 ITSTRK),$(get_N its-clusterer ITS REST 1 ITSCL)" "$ITS_CONFIG_KEY;$ITSMFT_STROBES;$ITSEXTRAERR" [[ ${DISABLE_DIGIT_CLUSTER_INPUT:-} =~ "--digits-from-upstream" ]] && has_detector_gpu ITS && ! has_detector_from_global_reader ITS && add_W o2-its-reco-workflow "--disable-tracking ${DISABLE_DIGIT_CLUSTER_INPUT:-} $ITS_STAGGERED $DISABLE_MC $DISABLE_ROOT_OUTPUT --pipeline $(get_N its-clusterer ITS REST 1 ITSCL)" "$ITS_CONFIG_KEY;$ITSMFT_STROBES;$ITSEXTRAERR" (has_detector_reco TPC || has_detector_ctf TPC) && ! has_detector_from_global_reader TPC && add_W o2-gpu-reco-workflow "--gpu-reconstruction \"$GPU_CONFIG_SELF\" --input-type=$GPU_INPUT $DISABLE_MC --output-type $GPU_OUTPUT $TPC_CORR_OPT $ITS_STAGGERED --pipeline gpu-reconstruction:${N_TPCTRK:-1},gpu-reconstruction-prepare:${N_TPCTRK:-1} $GPU_CONFIG" "GPU_global.deviceType=$GPUTYPE;GPU_proc.debugLevel=0;$GPU_CONFIG_KEY;$TRACKTUNETPCINNER;$TPC_CORR_KEY" (has_detector_reco TOF || has_detector_ctf TOF) && ! has_detector_from_global_reader TOF && add_W o2-tof-reco-workflow "$TOF_CONFIG --input-type $TOF_INPUT --output-type $TOF_OUTPUT $DISABLE_DIGIT_ROOT_INPUT $DISABLE_ROOT_OUTPUT $DISABLE_MC --pipeline $(get_N tof-compressed-decoder TOF RAW 1),$(get_N TOFClusterer TOF REST 1)" diff --git a/prodtests/full_system_test.sh b/prodtests/full_system_test.sh index 46739e76f103b..e89d8ee09dee9 100755 --- a/prodtests/full_system_test.sh +++ b/prodtests/full_system_test.sh @@ -188,7 +188,7 @@ taskwrapper digi.log o2-sim-digitizer-workflow -n $NEvents ${DIGIQED} ${NOMCLABE touch digiTRD.log_done if [[ "0$GENERATE_ITSMFT_DICTIONARIES" == "01" ]]; then - taskwrapper itsmftdict1.log o2-its-reco-workflow --trackerCA --disable-mc --configKeyValues '"fastMultConfig.cutMultClusLow=30000;fastMultConfig.cutMultClusHigh=2000000;fastMultConfig.cutMultVtxHigh=500;"' + taskwrapper itsmftdict1.log o2-its-reco-workflow --disable-mc --configKeyValues '"fastMultConfig.cutMultClusLow=30000;fastMultConfig.cutMultClusHigh=2000000;fastMultConfig.cutMultVtxHigh=500;"' cp ~/alice/O2/Detectors/ITSMFT/ITS/macros/test/CreateDictionaries.C . taskwrapper itsmftdict2.log root -b -q CreateDictionaries.C++ rm -f CreateDictionaries_C* CreateDictionaries.C diff --git a/prodtests/sim_challenge.sh b/prodtests/sim_challenge.sh index 8c7cfb1a024b0..f5bbf8ab74ff8 100755 --- a/prodtests/sim_challenge.sh +++ b/prodtests/sim_challenge.sh @@ -153,7 +153,7 @@ if [ "$doreco" == "1" ]; then echo "Return status of tpcreco: $?" echo "Running ITS reco flow" - taskwrapper itsreco.log o2-its-reco-workflow --trackerCA --tracking-mode async $gloOpt $ITSRecOpt + taskwrapper itsreco.log o2-its-reco-workflow --tracking-mode async $gloOpt $ITSRecOpt echo "Return status of itsreco: $?" # existing checks