diff --git a/.gitignore b/.gitignore index 656d005..aefaf52 100644 --- a/.gitignore +++ b/.gitignore @@ -247,3 +247,6 @@ test/k4FWCoreTest/**/*.root test/inputFiles/*.slcio test/gaudi_opts/testConverterConstants.py +# TrackingValidation test outputs +TrackingPerformance/test/*.root +TrackingPerformance/test/*.log diff --git a/CMakeLists.txt b/CMakeLists.txt index 9d2eccf..164938b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -47,5 +47,6 @@ function(set_test_env _testname) endfunction() add_subdirectory(RecoMCTruthLinkers) +add_subdirectory(TrackingPerformance) include(cmake/CreateProjectConfig.cmake) diff --git a/TrackingPerformance/CMakeLists.txt b/TrackingPerformance/CMakeLists.txt new file mode 100644 index 0000000..ae4d405 --- /dev/null +++ b/TrackingPerformance/CMakeLists.txt @@ -0,0 +1,103 @@ +#[[ +Copyright (c) 2020-2024 Key4hep-Project. + +This file is part of Key4hep. +See https://key4hep.github.io/key4hep-doc/ for further info. + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +]] + +set(PackageName TrackingPerformance) + +project(${PackageName}) + +include(ExternalData) +list(APPEND ExternalData_URL_TEMPLATES + "https://key4hep.web.cern.ch:443/testFiles/k4RecTracker/%(hash)" +) + +file(GLOB sources + ${PROJECT_SOURCE_DIR}/src/*.cc + ${PROJECT_SOURCE_DIR}/src/*.cpp + ${PROJECT_SOURCE_DIR}/src/components/*.cpp +) + +file(GLOB headers + ${PROJECT_SOURCE_DIR}/include/*.h +) + +find_package(ROOT REQUIRED COMPONENTS Core RIO Tree MathCore MathMore Graf Graf3d Hist) + +gaudi_add_module(${PackageName} + SOURCES ${sources} + LINK + Gaudi::GaudiKernel + EDM4HEP::edm4hep + k4FWCore::k4FWCore + ROOT::MathCore + ROOT::MathMore + ROOT::Physics + ROOT::Graf + ROOT::Graf3d + ROOT::Hist +) + +target_include_directories(${PackageName} PUBLIC + $ + $ + ${MarlinUtil_INCLUDE_DIRS} + ${DELPHES_INCLUDE_DIRS} +) + +set_target_properties(${PackageName} PROPERTIES + PUBLIC_HEADER "${headers}" +) + +file(GLOB test_python_scripts + ${PROJECT_SOURCE_DIR}/test/*.py +) + +set(test_shell_scripts + ${PROJECT_SOURCE_DIR}/test/testTrackingValidation.sh +) + +set(test_scripts + ${test_python_scripts} + ${test_shell_scripts} +) + +file(COPY ${test_scripts} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/test) + +install(TARGETS ${PackageName} + EXPORT ${CMAKE_PROJECT_NAME}Targets + RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin + LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib + PUBLIC_HEADER DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}/${CMAKE_PROJECT_NAME}" COMPONENT dev +) +message(STATUS "CMAKE_SOURCE_DIR = ${CMAKE_SOURCE_DIR}") +message(STATUS "CMAKE_CURRENT_SOURCE_DIR = ${CMAKE_CURRENT_SOURCE_DIR}") +message(STATUS "PROJECT_SOURCE_DIR = ${PROJECT_SOURCE_DIR}") +install(FILES ${test_scripts} DESTINATION test) + +SET(test_name "testTrackingValidation") + +ExternalData_Add_Test(testTrackingValidation + NAME ${test_name} + COMMAND sh +x TrackingPerformance/test/testTrackingValidation.sh + DATA{${CMAKE_SOURCE_DIR}/TrackingPerformance/test/inputFiles/SimpleGatrIDEAv3o1.onnx} +) + +set_test_env(${test_name}) +set_tests_properties(${test_name} PROPERTIES WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}") +ExternalData_Add_Target(${test_name}) + diff --git a/TrackingPerformance/README.md b/TrackingPerformance/README.md new file mode 100644 index 0000000..e20af91 --- /dev/null +++ b/TrackingPerformance/README.md @@ -0,0 +1,178 @@ +# TrackingValidation + +## Overview + +`TrackingValidation` is a validation algorithm for studying the performance (efficiency, purity, residuals, resolutions) of track finding and track fitting in the tracking reconstruction. +It is designed to compare reconstructed and fitted tracks with Monte Carlo truth information and, when enabled, with tracks obtained from perfect tracking, i.e. tracks fitted using the correct simHits from the particle truth information. The algorithm writes validation information to a ROOT output file containing TTrees and summary plots that can be used later for performance studies and plotting. + +Typical use cases include: +- validation of track-finder performance, +- validation of fitted-track parameters against MC truth, +- comparison between standard reconstructed tracks and perfectly associated reference tracks. + +--- + +## Inputs + +`TrackingValidation` expects EDM4hep event content in which the relevant collections have already been produced by the preceding steps of the reconstruction chain. + +### Input collection types + +`TrackingValidation` consumes the following input collections: + +- **MC particle collection** + Type: `edm4hep::MCParticleCollection` + Used as the truth reference for particle-level validation. + +- **Planar digi-to-sim link collections** + Type: `std::vector` + Used to connect reconstructed planar hits to the originating simulated particles. + +- **Drift-chamber digi-to-sim link collections** + Type: `std::vector` + Used to connect reconstructed drift-chamber hits to the originating simulated particles. + +- **Finder track collection** + Type: `edm4hep::TrackCollection` + Collection of tracks produced by the track-finding stage. + +- **Fitted track collection** + Type: `edm4hep::TrackCollection` + Collection of tracks produced by the standard fitting stage. + +- **Perfect fitted-track collections (optional)** + Type: `std::vector` + Optional reference collections produced from perfect truth-based associations, used when perfect-fit validation is enabled. + +--- + +## Outputs + +The algorithm writes a ROOT file specified by `OutputFile`. + +The file contains validation TTrees for finder-level and fitter-level studies, together with summary performance plots produced in `finalize()`. The fitter validation trees store residuals of the reconstructed track parameters with respect to the chosen reference. + +### Output content by mode + +The exact content filled in the output depends on the validation mode selected through `Mode`: + +- **`Mode = 0` (full-pipeline mode)** + Both finder-level and fitter-level validation are performed. + The output includes the association trees and the fitter residual trees. + +- **`Mode = 1` (finder-only mode)** + Only the finder-level validation is performed. + The finder and perfect-association trees are filled, while the fitter trees are booked in the file but are not filled. + +- **`Mode = 2` (fitter-only mode)** + Only the fitter-level validation is performed. + The fitter trees are filled, while the finder and perfect-association trees are booked in the file but are not filled. + +### Effect of `DoPerfectFit` + +The flag `DoPerfectFit` controls the handling of the `fitter_vs_perfect` output: + +- if **`DoPerfectFit = true`** and perfect fitted-track collections are provided, the fitter-to-perfect comparison is filled; +- if **`DoPerfectFit = false`**, the `fitter_vs_perfect` tree is still created but its per-event content remains empty; +- if **`DoPerfectFit = true`** but no perfect fitted-track collection is provided, the tree is still written and a warning is issued. + +### Summary plots + +In `finalize()`, the algorithm also writes summary plots to the same ROOT file, including: + +- tracking efficiency vs momentum, +- `d0` resolution vs momentum, +- momentum resolution vs momentum, +- transverse-momentum resolution vs momentum. + +--- + +## Finder validation: efficiency and purity + +To evaluate finder performance, each reconstructed track is matched to the truth particle with which it shares the largest number of hits. + +For each particle-track pair, the algorithm stores two standard hit-based quantities: + +- **track hit purity**: the fraction of hits on the reconstructed track that originate from the matched truth particle; +- **track hit efficiency**: the fraction of the truth-particle hits that are recovered in the reconstructed track. + +The summary **tracking efficiency** can then be defined in more than one way. + +- **`FinderEfficiencyDefinition = 1`** + A truth particle is counted as reconstructed if it is associated to at least one finder track with + `purity >= FinderPurityThreshold`. + In the default configuration, `FinderPurityThreshold = 0.75`, following the CMS association convention in which a reconstructed track is associated to a simulated particle if more than 75% of its hits originate from that particle. The tracking efficiency is then defined as the fraction of simulated tracks associated to at least one reconstructed track. :contentReference[oaicite:0]{index=0} + +- **`FinderEfficiencyDefinition = 2`** + A truth particle is counted as reconstructed if it is associated to at least one finder track with + `purity >= 0.5` **and** `efficiency >= 0.5`. + This corresponds to the stricter two-ratio variant, where both the purity of the reconstructed track and the fraction of recovered truth hits must exceed 50%. + +In the current implementation, the denominator of the efficiency plot includes generator-level particles with status 1 and at least one truth-linked hit. +For more details on the CMS association convention and the related definitions of tracking efficiency, fake rate, and duplicate rate, see the CMS performance note [*Performance of the track selection DNN in Run 3*](https://cds.cern.ch/record/2854696/files/DP2023_009.pdf). + + + + +--- + +## How to run + +`TrackingValidation` is tested through a small end-to-end workflow driven by `ctest`. The test starts from a simulated EDM4hep file, runs the reconstruction and validation steering, and writes the final validation ROOT output. + +In the current setup: + +- the simulation step is performed with `ddsim` in the shell test, +- the reconstruction and validation steps are controlled by `runTrackingValidation.py`, +- the full test is launched through `ctest`. + +Run the test from the build directory: + +```bash +cd k4DetectorPerformance/build +ctest -V -R testTrackingValidation +``` + +The validation output is written to: + +```text +k4DetectorPerformance/TrackingPerformance/test/validation_output_test.root +``` +### Test configuration and steering options + +The current test runs the full reconstruction and validation chain with the following settings: + +- `TRACKINGPERF_RUN_SIM = 1` + simulation step in the shell test (`0` = skip simulation and use an existing EDM4hep input file via `TRACKINGPERF_INPUT_FILE_OVERRIDE`, `1` = run simulation with `ddsim`); + +- `runDigi = 1` + digitization step (`0` = skip digitization, `1` = run digitization); + +- `runFinder = 1` + track-finder step (`0` = skip track finder, `1` = run track finder); + +- `runFitter = 1` + reconstructed-track fitter (`0` = skip reco fitter, `1` = run reco fitter); + +- `runPerfectTracking = 1` + perfect-tracking and perfect-fitter chain (`0` = skip perfect tracking/perfect fitter, `1` = run them); + +- `runValidation = 1` + validation step (`0` = skip validation, `1` = run validation); + +- `useDCH = 1` + drift-chamber collections (`0` = disable DCH, `1` = use DCH); + +- `mode = 0` + validation mode (`0` = full validation, `1` = finder-only validation, `2` = fitter-only validation); + +- `doPerfectFit = 1` + fitter-versus-perfect comparison (`0` = disable fitter-vs-perfect filling, `1` = enable it); + +- `finderEfficiencyDefinition = 1` + tracking-efficiency definition (`1` = purity-based definition, `2` = purity >= 0.5 and efficiency >= 0.5); + +- `finderPurityThreshold = 0.75` + purity threshold used when `FinderEfficiencyDefinition = 1`. + +These command-line flags are defined in `runTrackingValidation.py`, which allows the same steering file to be used either for the full chain or for reduced workflows in which some reconstruction steps are skipped and only the validation is run. diff --git a/TrackingPerformance/include/TrackingValidationHelpers.h b/TrackingPerformance/include/TrackingValidationHelpers.h new file mode 100644 index 0000000..eb0ec1a --- /dev/null +++ b/TrackingPerformance/include/TrackingValidationHelpers.h @@ -0,0 +1,99 @@ +/* + * Copyright (c) 2020-2024 Key4hep-Project. + * + * This file is part of Key4hep. + * See https://key4hep.github.io/key4hep-doc/ for further info. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#ifndef TRACKINGVALIDATIONHELPERS_H +#define TRACKINGVALIDATIONHELPERS_H + +#include "edm4hep/MCParticle.h" +#include "edm4hep/Track.h" +#include "edm4hep/TrackState.h" +#include "podio/ObjectID.h" + +#include + +namespace TrackingValidationHelpers { + +/// Container for helix parameters in the fitter convention +struct HelixParams { + float D0 = 0.f; + float Z0 = 0.f; + float phi = 0.f; + float omega = 0.f; + float tanLambda = 0.f; + float p = 0.f; + float pT = 0.f; +}; + +/// Helper container for PCA position and tangent angle +struct PCAInfoHelper { + float pcaX = 0.f; + float pcaY = 0.f; + float pcaZ = 0.f; + float phi0 = 0.f; + bool ok = false; +}; + +/// Build a unique integer key from a podio object identifier +uint64_t oidKey(const podio::ObjectID& id); + +/// Safe wrapper around std::atan2 +float safeAtan2(float y, float x); + +/// Wrap a phi difference into the interval [-pi, pi] +float wrapDeltaPhi(float a, float b); + +/** + * @brief Compute the PCA position and tangent angle in mm using a GenFit-like convention. + * + * The function derives the point of closest approach to the reference point + * and the corresponding tangent direction from the input position, momentum, + * charge sign, and magnetic field. + */ + +PCAInfoHelper PCAInfo_mm(float x, float y, float z, + float px, float py, float pz, + int chargeSign, + float refX, float refY, + float Bz); + + +/** + * @brief Build truth helix parameters in the same convention used by the fitter. + * + * The returned parameters are derived from the MC particle kinematics and vertex + * using the same reference-point and helix convention adopted for fitted tracks, + * so that residuals can be computed consistently. + */ + +HelixParams truthFromMC_GenfitConvention(const edm4hep::MCParticle& mc, + float Bz, + float refX, float refY, float refZ); + +/// Retrieve the track state stored at the interaction point +bool getAtIPState(const edm4hep::Track& trk, edm4hep::TrackState& out); + +/// Compute the transverse momentum from a track state +float ptFromState(const edm4hep::TrackState& st, float Bz); + +/// Compute the total momentum from a track state +float momentumFromState(const edm4hep::TrackState& st, float Bz); + +} // namespace TrackingValidationHelpers + +#endif diff --git a/TrackingPerformance/include/TrackingValidationPlots.h b/TrackingPerformance/include/TrackingValidationPlots.h new file mode 100644 index 0000000..acec572 --- /dev/null +++ b/TrackingPerformance/include/TrackingValidationPlots.h @@ -0,0 +1,118 @@ +/* + * Copyright (c) 2020-2024 Key4hep-Project. + * + * This file is part of Key4hep. + * See https://key4hep.github.io/key4hep-doc/ for further info. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#ifndef TRACKINGVALIDATIONPLOTS_H +#define TRACKINGVALIDATIONPLOTS_H + +#include "TCanvas.h" +#include "TF1.h" +#include "TGraphErrors.h" +#include "TH1F.h" +#include "TTree.h" + +#include +#include +#include + +namespace TrackingValidationPlots { + +/// Build logarithmic momentum bins +std::vector makeLogBins(double min, double max, double step); + +/// Fit the Gaussian core of a histogram +TF1* fitGaussianCore(TH1F* h, const std::string& name); + +/** + * @brief Build the d0 resolution as a function of momentum. + * + * The function reads the fitter validation tree, fills residual histograms in + * momentum bins, fits the Gaussian core in each bin, and returns the extracted + * sigma values as a TGraphErrors. + */ + +TGraphErrors* makeD0ResolutionVsMomentum(TTree* tree, + const char* graphName = "g_d0_resolution_vs_p", + double pMin = 0.1, + double pMax = 100.0, + double logStep = 0.15); + +/// Draw the d0 resolution graph on a logarithmic momentum axis +TCanvas* drawD0ResolutionCanvas(TGraphErrors* g, + const char* canvasName = "c_d0_resolution_vs_p", + double xMin = 0.1, + double xMax = 100.0); + +/** + * @brief Build the total-momentum resolution as a function of momentum. + * + * The function reads the fitter validation tree, fills relative momentum + * residual histograms in momentum bins, fits the Gaussian core in each bin, + * and returns the extracted sigma values as a TGraphErrors. + */ +TGraphErrors* makeMomentumResolutionVsMomentum(TTree* tree, + const char* graphName = "g_p_resolution_vs_p", + double pMin = 0.1, + double pMax = 100.0, + double logStep = 0.15); + +/** + * @brief Build the transverse-momentum resolution as a function of momentum. + * + * The function reads the fitter validation tree, fills relative transverse- + * momentum residual histograms in momentum bins, fits the Gaussian core in + * each bin, and returns the extracted sigma values as a TGraphErrors. + */ +TGraphErrors* makePtResolutionVsMomentum(TTree* tree, + const char* graphName = "g_pt_resolution_vs_p", + double pMin = 0.1, + double pMax = 100.0, + double logStep = 0.15); + +/// Draw a generic resolution graph on a logarithmic momentum axis +TCanvas* drawResolutionCanvas(TGraphErrors* g, + const char* canvasName, + const char* title, + double xMin = 0.1, + double xMax = 100.0); + +/** + * @brief Build the tracking-efficiency graph as a function of momentum. + * + * The function reads the finder validation tree and computes the efficiency + * according to the selected matching definition, returning the result as a + * TGraphErrors. + */ +TGraphErrors* makeEfficiencyVsMomentum(TTree* finderTree, + const char* graphName, + int efficiencyDefinition, + double purityThreshold, + double pMin = 0.1, + double pMax = 100.0, + double logStep = 0.15); + +/// Draw the tracking-efficiency graph on a logarithmic momentum axis +TCanvas* drawEfficiencyCanvas(TGraphErrors* g, + const char* canvasName, + const char* title, + double xMin = 0.1, + double xMax = 100.0); + +} // namespace TrackingValidationPlots + +#endif diff --git a/TrackingPerformance/src/TrackingValidationHelpers.cpp b/TrackingPerformance/src/TrackingValidationHelpers.cpp new file mode 100644 index 0000000..5d89025 --- /dev/null +++ b/TrackingPerformance/src/TrackingValidationHelpers.cpp @@ -0,0 +1,181 @@ +/* + * Copyright (c) 2020-2024 Key4hep-Project. + * + * This file is part of Key4hep. + * See https://key4hep.github.io/key4hep-doc/ for further info. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "TrackingValidationHelpers.h" + +#include +#include + +namespace TrackingValidationHelpers { + +// Constants matching the fitter code +static constexpr float c_mm_s = 2.998e11f; +static constexpr float a_genfit = 1e-15f * c_mm_s; + +uint64_t oidKey(const podio::ObjectID& id) { + return (uint64_t(id.collectionID) << 32) | uint64_t(uint32_t(id.index)); +} + +float safeAtan2(float y, float x) { + return std::atan2(y, x); +} + +float wrapDeltaPhi(float a, float b) { + float d = a - b; + while (d > M_PI) d -= 2.f * M_PI; + while (d < -M_PI) d += 2.f * M_PI; + return d; +} + +PCAInfoHelper PCAInfo_mm(float x, float y, float z, + float px, float py, float pz, + int chargeSign, + float refX, float refY, + float Bz) { + PCAInfoHelper out; + + const float pt = std::sqrt(px * px + py * py); + if (pt == 0.f) return out; + if (chargeSign == 0) chargeSign = 1; + if (Bz == 0.f) return out; + + const float R = pt / (0.3f * std::abs(chargeSign) * Bz) * 1000.f; + + const float tx = px / pt; + const float ty = py / pt; + + const float nx = float(chargeSign) * ty; + const float ny = float(chargeSign) * (-tx); + + const float xc = x + R * nx; + const float yc = y + R * ny; + + const float vx = refX - xc; + const float vy = refY - yc; + const float vxy = std::sqrt(vx * vx + vy * vy); + if (vxy == 0.f) return out; + + const float ux = vx / vxy; + const float uy = vy / vxy; + + const float pcaX = xc + R * ux; + const float pcaY = yc + R * uy; + + const float rx = pcaX - xc; + const float ry = pcaY - yc; + + const int sign = (chargeSign > 0) ? 1 : -1; + float tanX = -sign * ry; + float tanY = sign * rx; + + const float tnorm = std::sqrt(tanX * tanX + tanY * tanY); + if (tnorm == 0.f) return out; + + tanX /= tnorm; + tanY /= tnorm; + + const float phi0 = std::atan2(tanY, tanX); + + const float pR = pt; + const float pZ = pz; + const float R0 = std::sqrt(x * x + y * y); + const float Z0 = z; + + const float denom = (pR * pR + pZ * pZ); + if (denom == 0.f) return out; + + const float tPCA = -(R0 * pR + Z0 * pZ) / denom; + const float ZPCA = Z0 + pZ * tPCA; + + out.pcaX = pcaX; + out.pcaY = pcaY; + out.pcaZ = ZPCA; + out.phi0 = phi0; + out.ok = true; + return out; +} + +HelixParams truthFromMC_GenfitConvention(const edm4hep::MCParticle& mc, + float Bz, + float refX, float refY, float refZ) { + HelixParams hp; + + const auto& mom = mc.getMomentum(); + const float px = float(mom.x); + const float py = float(mom.y); + const float pz = float(mom.z); + + const float pT = std::sqrt(px * px + py * py); + const float p = std::sqrt(px * px + py * py + pz * pz); + + hp.pT = pT; + hp.p = p; + + int qSign = 1; + if (mc.getCharge() < 0.f) qSign = -1; + + const auto& v = mc.getVertex(); + const float x = float(v.x); + const float y = float(v.y); + const float z = float(v.z); + + const auto info = PCAInfo_mm(x, y, z, px, py, pz, qSign, refX, refY, Bz); + if (!info.ok) { + const float NaN = std::numeric_limits::quiet_NaN(); + hp.D0 = NaN; + hp.Z0 = NaN; + hp.phi = NaN; + hp.omega = NaN; + hp.tanLambda = NaN; + return hp; + } + + hp.D0 = ((-(refX - info.pcaX)) * std::sin(info.phi0) + + (refY - info.pcaY) * std::cos(info.phi0)); + hp.Z0 = (info.pcaZ - refZ); + hp.phi = safeAtan2(py, px); + hp.tanLambda = (pT > 0.f) ? (pz / pT) : 0.f; + hp.omega = (pT > 0.f) ? (std::abs(a_genfit * Bz / pT) * float(qSign)) : 0.f; + + return hp; +} + +bool getAtIPState(const edm4hep::Track& trk, edm4hep::TrackState& out) { + for (const auto& st : trk.getTrackStates()) { + if (st.location == edm4hep::TrackState::AtIP) { + out = st; + return true; + } + } + return false; +} + +float ptFromState(const edm4hep::TrackState& st, float Bz) { + const float omega = std::abs(float(st.omega)); + if (omega == 0.f) return 0.f; + return a_genfit * std::abs(Bz) / omega; +} + +float momentumFromState(const edm4hep::TrackState& st, float Bz) { + const float pT = ptFromState(st, Bz); + const float tl = float(st.tanLambda); + return pT * std::sqrt(1.f + tl * tl); +} + +} // namespace TrackingValidationHelpers diff --git a/TrackingPerformance/src/TrackingValidationPlots.cpp b/TrackingPerformance/src/TrackingValidationPlots.cpp new file mode 100644 index 0000000..38f2ded --- /dev/null +++ b/TrackingPerformance/src/TrackingValidationPlots.cpp @@ -0,0 +1,440 @@ +#include "TrackingValidationPlots.h" +#include "TStyle.h" +#include +#include +#include +#include + +// makeLogBins +//fitGaussianCore +//makeD0ResolutionVsMomentum +//drawD0ResolutionCanvas +//makeMomentumResolutionVsMomentum +//makePtResolutionVsMomentum +//drawResolutionCanvas +//makeEfficiencyVsMomentum +//drawEfficiencyCanvas + +namespace TrackingValidationPlots { +std::vector makeLogBins(double min, double max, double step) { + std::vector bins; + for (double x = std::log10(min); x <= std::log10(max); x += step) { + bins.push_back(std::pow(10., x)); + } + if (bins.empty() || bins.back() < max) bins.push_back(max); + return bins; +} + +TF1* fitGaussianCore(TH1F* h, const std::string& name) { + if (!h || h->GetEntries() < 10) return nullptr; + + const double mean = h->GetMean(); + const double rms = h->GetRMS(); + if (rms <= 0.) return nullptr; + + TF1* g1 = new TF1((name + "_g1").c_str(), "gaus", mean - 2.0 * rms, mean + 2.0 * rms); + h->Fit(g1, "RQ0"); + + double m1 = g1->GetParameter(1); + double s1 = std::abs(g1->GetParameter(2)); + if (s1 <= 0.) s1 = rms; + + TF1* g2 = new TF1(name.c_str(), "gaus", m1 - 1.5 * s1, m1 + 1.5 * s1); + h->Fit(g2, "RQ0"); + + return g2; +} + +TGraphErrors* makeD0ResolutionVsMomentum(TTree* tree, + const char* graphName, + double pMin, + double pMax, + double logStep) { + if (!tree) return nullptr; + + std::vector bins = makeLogBins(pMin, pMax, logStep); + const int nBins = bins.size() - 1; + + std::vector> hists; + hists.reserve(nBins); + + for (int i = 0; i < nBins; ++i) { + hists.emplace_back(std::make_unique( + Form("h_d0_bin_%d", i), + Form("d0 residual bin %d;resD0 [#mum];Entries", i), + 120, -20.0, 20.0)); + } + + std::vector* resD0 = nullptr; + std::vector* p_ref_vec = nullptr; + + tree->SetBranchAddress("p_ref", &p_ref_vec); + tree->SetBranchAddress("resD0", &resD0); + + const Long64_t nEntries = tree->GetEntries(); + for (Long64_t ievt = 0; ievt < nEntries; ++ievt) { + tree->GetEntry(ievt); + + if (!p_ref_vec || !resD0) continue; + if (p_ref_vec->size() != resD0->size()) continue; + + for (size_t i = 0; i < p_ref_vec->size(); ++i) { + const double p = (*p_ref_vec)[i]; + const double d0_um = (*resD0)[i] * 1000.0; + + if (!std::isfinite(p) || !std::isfinite(d0_um)) continue; + if (p < pMin || p >= pMax) continue; + + int bin = -1; + for (int b = 0; b < nBins; ++b) { + if (p >= bins[b] && p < bins[b + 1]) { + bin = b; + break; + } + } + if (bin < 0) continue; + + hists[bin]->Fill(d0_um); + } + } + + TGraphErrors* g = new TGraphErrors(); + g->SetName(graphName); + g->SetTitle(";p_{ref} [GeV];#sigma(d_{0}) [#mum]"); + + int ip = 0; + for (int b = 0; b < nBins; ++b) { + if (hists[b]->GetEntries() < 20) continue; + + TF1* fit = fitGaussianCore(hists[b].get(), Form("fit_d0_bin_%d", b)); + if (!fit) continue; + + const double sigma = std::abs(fit->GetParameter(2)); + const double sigmaErr = fit->GetParError(2); + const double pCenter = std::sqrt(bins[b] * bins[b + 1]); + + g->SetPoint(ip, pCenter, sigma); + g->SetPointError(ip, 0.0, sigmaErr); + ++ip; + } + + return g; +} + +TCanvas* drawD0ResolutionCanvas(TGraphErrors* g, + const char* canvasName, + double xMin, + double xMax) { + if (!g) return nullptr; + + gStyle->SetOptStat(0); + + TCanvas* c = new TCanvas(canvasName, "d0 resolution vs momentum", 800, 600); + c->SetLogx(); + + g->SetMarkerStyle(20); + g->SetLineWidth(2); + g->GetXaxis()->SetLimits(xMin, xMax); + g->Draw("AP"); + + return c; +} + +TGraphErrors* makeMomentumResolutionVsMomentum(TTree* tree, + const char* graphName, + double pMin, + double pMax, + double logStep) { + if (!tree) return nullptr; + + std::vector bins = makeLogBins(pMin, pMax, logStep); + const int nBins = bins.size() - 1; + + std::vector> hists; + hists.reserve(nBins); + + for (int i = 0; i < nBins; ++i) { + hists.emplace_back(std::make_unique( + Form("h_pres_bin_%d", i), + Form("p resolution bin %d;(p_{reco}-p_{ref})/p_{ref};Entries", i), + 120, -0.2, 0.2)); + } + + std::vector* p_ref_vec = nullptr; + std::vector* p_reco_vec = nullptr; + + tree->SetBranchAddress("p_ref", &p_ref_vec); + tree->SetBranchAddress("p_reco", &p_reco_vec); + + const Long64_t nEntries = tree->GetEntries(); + for (Long64_t ievt = 0; ievt < nEntries; ++ievt) { + tree->GetEntry(ievt); + + if (!p_ref_vec || !p_reco_vec) continue; + if (p_ref_vec->size() != p_reco_vec->size()) continue; + + for (size_t i = 0; i < p_ref_vec->size(); ++i) { + const double pRef = (*p_ref_vec)[i]; + const double pReco = (*p_reco_vec)[i]; + + if (!std::isfinite(pRef) || !std::isfinite(pReco)) continue; + if (pRef <= 0.) continue; + if (pRef < pMin || pRef >= pMax) continue; + + const double res = (pReco - pRef) / pRef; + + int bin = -1; + for (int b = 0; b < nBins; ++b) { + if (pRef >= bins[b] && pRef < bins[b + 1]) { + bin = b; + break; + } + } + if (bin < 0) continue; + + hists[bin]->Fill(res); + } + } + + TGraphErrors* g = new TGraphErrors(); + g->SetName(graphName); + g->SetTitle(";p_{ref} [GeV];#sigma((p_{reco}-p_{ref})/p_{ref})"); + + int ip = 0; + for (int b = 0; b < nBins; ++b) { + if (hists[b]->GetEntries() < 20) continue; + + TF1* fit = fitGaussianCore(hists[b].get(), Form("fit_pres_bin_%d", b)); + if (!fit) continue; + + const double sigma = std::abs(fit->GetParameter(2)); + const double sigmaErr = fit->GetParError(2); + const double pCenter = std::sqrt(bins[b] * bins[b + 1]); + + g->SetPoint(ip, pCenter, sigma); + g->SetPointError(ip, 0.0, sigmaErr); + ++ip; + } + + return g; +} + +TGraphErrors* makePtResolutionVsMomentum(TTree* tree, + const char* graphName, + double pMin, + double pMax, + double logStep) { + if (!tree) return nullptr; + + std::vector bins = makeLogBins(pMin, pMax, logStep); + const int nBins = bins.size() - 1; + + std::vector> hists; + hists.reserve(nBins); + + for (int i = 0; i < nBins; ++i) { + hists.emplace_back(std::make_unique( + Form("h_ptres_bin_%d", i), + Form("pT resolution bin %d;(pT_{reco}-pT_{ref})/pT_{ref};Entries", i), + 120, -0.2, 0.2)); + } + + std::vector* p_ref_vec = nullptr; + std::vector* pt_ref_vec = nullptr; + std::vector* pt_reco_vec = nullptr; + + tree->SetBranchAddress("p_ref", &p_ref_vec); + tree->SetBranchAddress("pT_ref", &pt_ref_vec); + tree->SetBranchAddress("pT_reco", &pt_reco_vec); + + const Long64_t nEntries = tree->GetEntries(); + for (Long64_t ievt = 0; ievt < nEntries; ++ievt) { + tree->GetEntry(ievt); + + if (!p_ref_vec || !pt_ref_vec || !pt_reco_vec) continue; + if (p_ref_vec->size() != pt_ref_vec->size()) continue; + if (pt_ref_vec->size() != pt_reco_vec->size()) continue; + + for (size_t i = 0; i < p_ref_vec->size(); ++i) { + const double pRef = (*p_ref_vec)[i]; + const double ptRef = (*pt_ref_vec)[i]; + const double ptReco = (*pt_reco_vec)[i]; + + if (!std::isfinite(pRef) || !std::isfinite(ptRef) || !std::isfinite(ptReco)) continue; + if (ptRef <= 0.) continue; + if (pRef < pMin || pRef >= pMax) continue; + + const double res = (ptReco - ptRef) / ptRef; + + int bin = -1; + for (int b = 0; b < nBins; ++b) { + if (pRef >= bins[b] && pRef < bins[b + 1]) { + bin = b; + break; + } + } + if (bin < 0) continue; + + hists[bin]->Fill(res); + } + } + + TGraphErrors* g = new TGraphErrors(); + g->SetName(graphName); + g->SetTitle(";p_{ref} [GeV];#sigma((pT_{reco}-pT_{ref})/pT_{ref})"); + + int ip = 0; + for (int b = 0; b < nBins; ++b) { + if (hists[b]->GetEntries() < 20) continue; + + TF1* fit = fitGaussianCore(hists[b].get(), Form("fit_ptres_bin_%d", b)); + if (!fit) continue; + + const double sigma = std::abs(fit->GetParameter(2)); + const double sigmaErr = fit->GetParError(2); + const double pCenter = std::sqrt(bins[b] * bins[b + 1]); + + g->SetPoint(ip, pCenter, sigma); + g->SetPointError(ip, 0.0, sigmaErr); + ++ip; + } + + return g; +} + +TCanvas* drawResolutionCanvas(TGraphErrors* g, + const char* canvasName, + const char* title, + double xMin, + double xMax) { + if (!g) return nullptr; + + gStyle->SetOptStat(0); + + TCanvas* c = new TCanvas(canvasName, title, 800, 600); + c->SetLogx(); + + g->SetMarkerStyle(20); + g->SetLineWidth(2); + g->SetTitle(title); + g->GetXaxis()->SetLimits(xMin, xMax); + g->Draw("AP"); + + return c; +} + +TGraphErrors* makeEfficiencyVsMomentum(TTree* finderTree, + const char* graphName, + int efficiencyDefinition, + double purityThreshold, + double pMin, + double pMax, + double logStep) { + if (!finderTree) return nullptr; + + std::vector bins = makeLogBins(pMin, pMax, logStep); + const int nBins = bins.size() - 1; + + std::vector nDen(nBins, 0); + std::vector nNum(nBins, 0); + + std::vector* pVec = nullptr; + std::vector>* purVec = nullptr; + std::vector>* effVec = nullptr; + + finderTree->SetBranchAddress("p", &pVec); + finderTree->SetBranchAddress("matchPurity", &purVec); + finderTree->SetBranchAddress("matchEfficiency", &effVec); + + const Long64_t nEntries = finderTree->GetEntries(); + for (Long64_t ievt = 0; ievt < nEntries; ++ievt) { + finderTree->GetEntry(ievt); + + if (!pVec || !purVec || !effVec) continue; + if (pVec->size() != purVec->size()) continue; + if (pVec->size() != effVec->size()) continue; + + for (size_t i = 0; i < pVec->size(); ++i) { + const double p = (*pVec)[i]; + if (!std::isfinite(p) || p < pMin || p >= pMax) continue; + + int bin = -1; + for (int b = 0; b < nBins; ++b) { + if (p >= bins[b] && p < bins[b + 1]) { + bin = b; + break; + } + } + if (bin < 0) continue; + + nDen[bin]++; + + bool isMatched = false; + + const auto& purities = (*purVec)[i]; + const auto& efficiencies = (*effVec)[i]; + const size_t nMatches = std::min(purities.size(), efficiencies.size()); + + for (size_t j = 0; j < nMatches; ++j) { + const float purity = purities[j]; + const float efficiency = efficiencies[j]; + + if (efficiencyDefinition == 2) { + if (purity >= 0.5f && efficiency >= 0.5f) { + isMatched = true; + break; + } + } else { + if (purity >= purityThreshold) { + isMatched = true; + break; + } + } + } + + if (isMatched) nNum[bin]++; + } + } + + TGraphErrors* g = new TGraphErrors(); + g->SetName(graphName); + g->SetTitle(";p [GeV];Tracking efficiency"); + + int ip = 0; + for (int b = 0; b < nBins; ++b) { + if (nDen[b] == 0) continue; + + const double eff = double(nNum[b]) / double(nDen[b]); + const double err = std::sqrt(eff * (1.0 - eff) / double(nDen[b])); + const double pCenter = std::sqrt(bins[b] * bins[b + 1]); + + g->SetPoint(ip, pCenter, eff); + g->SetPointError(ip, 0.0, err); + ++ip; + } + + return g; +} + +TCanvas* drawEfficiencyCanvas(TGraphErrors* g, + const char* canvasName, + const char* title, + double xMin, + double xMax) { + if (!g) return nullptr; + + gStyle->SetOptStat(0); + + TCanvas* c = new TCanvas(canvasName, title, 800, 600); + c->SetLogx(); + + g->SetMarkerStyle(20); + g->SetLineWidth(2); + g->SetTitle(title); + g->GetYaxis()->SetRangeUser(0.0, 1.05); + g->GetXaxis()->SetLimits(xMin, xMax); + g->Draw("AP"); + + return c; +} +} diff --git a/TrackingPerformance/src/components/TrackingValidation.cpp b/TrackingPerformance/src/components/TrackingValidation.cpp new file mode 100644 index 0000000..1f1783a --- /dev/null +++ b/TrackingPerformance/src/components/TrackingValidation.cpp @@ -0,0 +1,740 @@ +/* + * Copyright (c) 2020-2024 Key4hep-Project. + * + * This file is part of Key4hep. + * See https://key4hep.github.io/key4hep-doc/ for further info. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. +*/ + +#include "TrackingValidationHelpers.h" +#include "TrackingValidationPlots.h" + +// k4FWCore +#include "k4FWCore/Consumer.h" + +// Gaudi +#include "Gaudi/Property.h" +#include "GaudiKernel/MsgStream.h" + +// EDM4hep +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/TrackCollection.h" +#include "edm4hep/TrackState.h" +#include "edm4hep/TrackerHitSimTrackerHitLinkCollection.h" + +// ROOT +#include "TFile.h" +#include "TTree.h" +#include "TGraphErrors.h" +#include "TCanvas.h" + + +// STL +#include +#include +#include +#include +#include +#include +#include +#include + + +/** @struct TrackingValidation + * + * Gaudi Consumer that validates the performance of track finding and track fitting + * by comparing reconstructed tracks with Monte Carlo truth information and, + * optionally, with perfectly associated fitted tracks. + * + * The consumer writes several ROOT TTrees containing finder-level associations, + * fitter residuals with respect to MC truth, and fitter residuals with respect + * to perfectly fitted reference tracks. In addition, summary performance plots + * are produced in finalize() and written to the same ROOT file. + * + * The supported validation modes are: + * - full pipeline validation, + * - finder-only validation, + * - fitter-only validation. + * + * input: + * - MC particle collection : edm4hep::MCParticleCollection + * - planar digi-to-sim link collections : std::vector + * - drift-chamber digi-to-sim link collections : std::vector + * - finder track collection : edm4hep::TrackCollection + * - fitted track collection : edm4hep::TrackCollection + * - optional perfect fitted-track collections : std::vector + * + * output: + * - ROOT file containing validation TTrees + * - summary performance plots written to the same ROOT file + * + * @author Arina Ponomareva + * @date 2026-03 + * + */ + + +// ---------- CONSUMER ---------- +struct TrackingValidation final + : k4FWCore::Consumer&, // planar links + const std::vector&, // optional DCH links + const edm4hep::TrackCollection&, // finder tracks + const edm4hep::TrackCollection&, // fitted tracks (reco) + const std::vector& // optional perfect fitted tracks + )> { + + TrackingValidation(const std::string& name, ISvcLocator* svcLoc) + : Consumer( + name, svcLoc, + { + KeyValues("MCParticles", {"MCParticles"}), + + KeyValues("PlanarLinks", + {"SiWrBSimDigiLinks", "SiWrDSimDigiLinks", "VTXBSimDigiLinks", "VTXDSimDigiLinks"}), + + + KeyValues("DCHLinks", {"DCH_DigiSimAssociationCollection"}), + + KeyValues("FinderTracks", {"GGTFTracks"}), + KeyValues("FittedTracks", {"FittedTracks"}), + + + KeyValues("PerfectFittedTracks", {"PerfectFittedTracks"}), + }) {} + + StatusCode initialize() override { + info() << "Initializing TrackingValidationConsumer" << endmsg; + + m_outFile = std::make_unique(m_outputFile.value().c_str(), "RECREATE"); + if (!m_outFile || m_outFile->IsZombie()) { + error() << "Cannot open output file: " << m_outputFile.value() << endmsg; + return StatusCode::FAILURE; + } + + bookAssocTree(m_finder_p2t, "finder_particle_to_tracks"); + bookAssocTree(m_finder_t2p, "finder_track_to_particles"); + bookAssocTree(m_perf_p2t, "perfect_particle_to_tracks"); + bookAssocTree(m_perf_t2p, "perfect_track_to_particles"); + + bookFitterTree(m_fit_vs_mc, "fitter_vs_mc"); + bookFitterTree(m_fit_vs_perfect, "fitter_vs_perfect"); + + return StatusCode::SUCCESS; + } + + void operator()(const edm4hep::MCParticleCollection& mcParts, + const std::vector& planarLinksVec, + const std::vector& dchLinksVec, + const edm4hep::TrackCollection& finderTracks, + const edm4hep::TrackCollection& fittedTracks, + const std::vector& perfectFittedTracksVec) const override { + + const int event = m_evt++; + const int mode = m_mode.value(); // 0 full, 1 finder-only, 2 fitter-only + + // ---------- Build truth maps: hit -> particle, particle -> hits ---------- + std::unordered_map> hitsPerParticle; + hitsPerParticle.reserve(mcParts.size()); + + std::unordered_map hitToParticle; + hitToParticle.reserve(200000); + + // planar links + for (const auto* links : planarLinksVec) { + if (!links) continue; + for (const auto& link : *links) { + const auto digi = link.getFrom(); + const auto sim = link.getTo(); + const auto mc = sim.getParticle(); + if (!digi.isAvailable() || !mc.isAvailable()) continue; + + const int pid = mc.getObjectID().index; + const uint64_t key = TrackingValidationHelpers::oidKey(digi.getObjectID()); + hitsPerParticle[pid].push_back(key); + hitToParticle[key] = pid; + } + } + + // DCH links + for (const auto* links : dchLinksVec) { + if (!links) continue; + for (const auto& link : *links) { + const auto digi = link.getFrom(); + const auto sim = link.getTo(); + const auto mc = sim.getParticle(); + if (!digi.isAvailable() || !mc.isAvailable()) continue; + + const int pid = mc.getObjectID().index; + const uint64_t key = TrackingValidationHelpers::oidKey(digi.getObjectID()); + hitsPerParticle[pid].push_back(key); + hitToParticle[key] = pid; + } + } + + // ---------- Finder & Perfect association trees ---------- + if (mode == 0 || mode == 1) { + fillPerfectAssoc(event, mcParts, hitsPerParticle); + fillFinderAssoc(event, mcParts, finderTracks, hitToParticle, hitsPerParticle); + } + + // ---------- Build pid -> best perfect-fitted AtIP state ---------- + struct StateWithNHits { + edm4hep::TrackState st; + int nHits = 0; + }; + + std::unordered_map perfectAtIPByPid; + + const bool wantPerfect = m_doPerfectFit.value(); + const bool havePerfectCollections = !perfectFittedTracksVec.empty(); + const bool doPerfect = wantPerfect && havePerfectCollections; + + if (wantPerfect && !havePerfectCollections && !m_warnedMissingPerfectInput) { + warning() << "DoPerfectFit=true but no PerfectFittedTracks collection was provided. " + << "fitter_vs_perfect will be filled with NaNs/empty content." << endmsg; + m_warnedMissingPerfectInput = true; + } + + if (doPerfect) { + size_t nPerfectTracks = 0; + for (const auto* coll : perfectFittedTracksVec) { + if (!coll) continue; + nPerfectTracks += coll->size(); + } + perfectAtIPByPid.reserve(nPerfectTracks); + + for (const auto* coll : perfectFittedTracksVec) { + if (!coll) continue; + + for (const auto& trk : *coll) { + edm4hep::TrackState st; + if (!TrackingValidationHelpers::getAtIPState(trk, st)) continue; + + const int pid = majorityParticleForTrack(trk, hitToParticle); + if (pid < 0 || pid >= (int)mcParts.size()) continue; + + const int nHits = (int)trk.getTrackerHits().size(); + auto it = perfectAtIPByPid.find(pid); + if (it == perfectAtIPByPid.end() || nHits > it->second.nHits) { + perfectAtIPByPid[pid] = StateWithNHits{st, nHits}; + } + } + } + } + + // ---------- Fitter trees ---------- + if (mode == 0 || mode == 2) { + fillFitterTrees(event, mcParts, fittedTracks, hitToParticle, perfectAtIPByPid, doPerfect); + } + } + + StatusCode finalize() override { + info() << "Finalizing TrackingValidation, wrote " << m_evt << " events" << endmsg; + + if (m_outFile) { + m_outFile->cd(); + + //write trees + if (m_finder_p2t.tree) m_finder_p2t.tree->Write(); + if (m_finder_t2p.tree) m_finder_t2p.tree->Write(); + if (m_perf_p2t.tree) m_perf_p2t.tree->Write(); + if (m_perf_t2p.tree) m_perf_t2p.tree->Write(); + + if (m_fit_vs_mc.tree) m_fit_vs_mc.tree->Write(); + if (m_fit_vs_perfect.tree) m_fit_vs_perfect.tree->Write(); + + //fitter summary plots + // d0 resolution vs momentum from fitter_vs_mc + TGraphErrors* g_d0_vs_p = TrackingValidationPlots::makeD0ResolutionVsMomentum(m_fit_vs_mc.tree, + "g_d0_resolution_vs_p", + 0.1, 100.0, 0.15); + if (g_d0_vs_p) { + TCanvas* c_d0_vs_p = TrackingValidationPlots::drawD0ResolutionCanvas(g_d0_vs_p, + "c_d0_resolution_vs_p", + 0.1, 100.0); + g_d0_vs_p->Write(); + if (c_d0_vs_p) c_d0_vs_p->Write(); + } + // p resolution vs momentum + TGraphErrors* g_p_vs_p = TrackingValidationPlots::makeMomentumResolutionVsMomentum(m_fit_vs_mc.tree, + "g_p_resolution_vs_p", + 0.1, 100.0, 0.15); + if (g_p_vs_p) { + TCanvas* c_p_vs_p = TrackingValidationPlots::drawResolutionCanvas(g_p_vs_p, + "c_p_resolution_vs_p", + "momentum resolution vs momentum;p_{ref} [GeV];#sigma((p_{reco}-p_{ref})/p_{ref})", + 0.1, 100.0); + g_p_vs_p->Write(); + if (c_p_vs_p) c_p_vs_p->Write(); + } + + // pT resolution vs momentum + TGraphErrors* g_pt_vs_p = TrackingValidationPlots::makePtResolutionVsMomentum(m_fit_vs_mc.tree, + "g_pt_resolution_vs_p", + 0.1, 100.0, 0.15); + if (g_pt_vs_p) { + TCanvas* c_pt_vs_p = TrackingValidationPlots::drawResolutionCanvas(g_pt_vs_p, + "c_pt_resolution_vs_p", + "pT resolution vs momentum;p_{ref} [GeV];#sigma((pT_{reco}-pT_{ref})/pT_{ref})", + 0.1, 100.0); + g_pt_vs_p->Write(); + if (c_pt_vs_p) c_pt_vs_p->Write(); + } + + // finder summary plot + TGraphErrors* g_eff_vs_p = TrackingValidationPlots::makeEfficiencyVsMomentum( + m_finder_p2t.tree, + "g_efficiency_vs_p", + m_finderEfficiencyDefinition.value(), + m_finderPurityThreshold.value(), + 0.1, 100.0, 0.15); + if (g_eff_vs_p) { + TCanvas* c_eff_vs_p = TrackingValidationPlots::drawEfficiencyCanvas( + g_eff_vs_p, + "c_efficiency_vs_p", + "tracking efficiency vs momentum;p [GeV];Efficiency", + 0.1, 100.0); + g_eff_vs_p->Write(); + if (c_eff_vs_p) c_eff_vs_p->Write(); + + } + + m_outFile->Close(); + } + return StatusCode::SUCCESS; + } + +private: + // ---------- properties ---------- + Gaudi::Property m_outputFile{this, "OutputFile", "validation.root", "Output ROOT file (TTrees)"}; + + // 0 full, 1 finder-only, 2 fitter-only + Gaudi::Property m_mode{this, "Mode", 0, "Validation mode: 0=Full, 1=FinderOnly, 2=FitterOnly"}; + + Gaudi::Property m_Bz{this, "Bz", 2.f, "Magnetic field Bz [T] used in omega convention (GenFit-style)"}; + + Gaudi::Property m_refX{this, "RefPointX", 0.f, "Reference point X [mm] (must match fitter m_VP_referencePoint)"}; + Gaudi::Property m_refY{this, "RefPointY", 0.f, "Reference point Y [mm] (must match fitter m_VP_referencePoint)"}; + Gaudi::Property m_refZ{this, "RefPointZ", 0.f, "Reference point Z [mm] (must match fitter m_VP_referencePoint)"}; + +// Definition used for the summary tracking-efficiency plot. +// Default = 1 keeps the current behaviour unchanged. + Gaudi::Property m_finderEfficiencyDefinition{ + this, "FinderEfficiencyDefinition", 1, + "Definition used for the tracking-efficiency summary plot: " + "1 = require purity >= FinderPurityThreshold; " + "2 = require purity >= 0.5 and efficiency >= 0.5" + }; + + Gaudi::Property m_finderPurityThreshold{ + this, "FinderPurityThreshold", 0.75f, + "Minimum purity for a particle-track match to count in tracking efficiency " + "when FinderEfficiencyDefinition = 1" + }; + + Gaudi::Property m_doPerfectFit{ + this, "DoPerfectFit", false, + "If true: fill fitter_vs_perfect using PerfectFitted_tracks if available. " + "If false: tree exists but is empty per event."}; + + // ---------- output structs ---------- + struct AssocTree { + TTree* tree = nullptr; + int event = 0; + std::vector index; + std::vector p; + std::vector pT; + std::vector nTrueHits; + std::vector> assoc; + + // only really used for finder_particle_to_tracks + std::vector> sharedHits; + std::vector> matchEfficiency; + std::vector> matchPurity; + + void clear() { + index.clear(); + p.clear(); + pT.clear(); + nTrueHits.clear(); + assoc.clear(); + sharedHits.clear(); + matchEfficiency.clear(); + matchPurity.clear(); + } + }; + + struct FitterTree { + TTree* tree = nullptr; + int event = 0; + + std::vector track_index; + std::vector track_location; + + std::vector resD0, resZ0, resPhi, resOmega, resTanL; + std::vector p_reco, p_ref; + std::vector pT_reco, pT_ref; + + void clear() { + track_index.clear(); + track_location.clear(); + resD0.clear(); + resZ0.clear(); + resPhi.clear(); + resOmega.clear(); + resTanL.clear(); + p_reco.clear(); + p_ref.clear(); + pT_reco.clear(); + pT_ref.clear(); + } + }; + + static void bookAssocTree(AssocTree& t, const char* name) { + t.tree = new TTree(name, name); + t.tree->Branch("event", &t.event); + t.tree->Branch("index", &t.index); + t.tree->Branch("p", &t.p); + t.tree->Branch("pT", &t.pT); + t.tree->Branch("nTrueHits", &t.nTrueHits); + t.tree->Branch("assoc", &t.assoc); + t.tree->Branch("sharedHits", &t.sharedHits); + t.tree->Branch("matchEfficiency", &t.matchEfficiency); + t.tree->Branch("matchPurity", &t.matchPurity); + } + + static void bookFitterTree(FitterTree& t, const char* name) { + t.tree = new TTree(name, name); + t.tree->Branch("event", &t.event); + t.tree->Branch("track_index", &t.track_index); + t.tree->Branch("track_location", &t.track_location); + t.tree->Branch("resD0", &t.resD0); + t.tree->Branch("resZ0", &t.resZ0); + t.tree->Branch("resPhi", &t.resPhi); + t.tree->Branch("resOmega", &t.resOmega); + t.tree->Branch("resTanLambda", &t.resTanL); + t.tree->Branch("p_reco", &t.p_reco); + t.tree->Branch("p_ref", &t.p_ref); + t.tree->Branch("pT_reco", &t.pT_reco); + t.tree->Branch("pT_ref", &t.pT_ref); + } + + // ---------- association trees ---------- + void fillPerfectAssoc(int event, const edm4hep::MCParticleCollection& mcParts, + const std::unordered_map>& hitsPerParticle) const { + + m_perf_p2t.clear(); + m_perf_t2p.clear(); + m_perf_p2t.event = event; + m_perf_t2p.event = event; + + for (int i = 0; i < (int)mcParts.size(); ++i) { + const auto& mc = mcParts[i]; + if (mc.getGeneratorStatus() != 1) continue; + + auto it = hitsPerParticle.find(i); + if (it == hitsPerParticle.end() || it->second.empty()) continue; + const auto& mom = mc.getMomentum(); + const float px = float(mom.x); + const float py = float(mom.y); + const float pz = float(mom.z); + const float p = std::sqrt(px*px + py*py + pz*pz); + const float pT = std::sqrt(px*px + py*py); + const int nHits = (int)it->second.size(); + + m_perf_p2t.index.push_back(i); + m_perf_p2t.p.push_back(p); + m_perf_p2t.pT.push_back(pT); + m_perf_p2t.nTrueHits.push_back(nHits); + m_perf_p2t.assoc.push_back({i}); + m_perf_p2t.sharedHits.push_back({}); + m_perf_p2t.matchEfficiency.push_back({}); + m_perf_p2t.matchPurity.push_back({}); + + m_perf_t2p.index.push_back(i); + m_perf_t2p.p.push_back(p); + m_perf_t2p.pT.push_back(pT); + m_perf_t2p.nTrueHits.push_back(nHits); + m_perf_t2p.assoc.push_back({i}); + m_perf_t2p.sharedHits.push_back({}); + m_perf_t2p.matchEfficiency.push_back({}); + m_perf_t2p.matchPurity.push_back({}); + } + + if (m_perf_p2t.tree) m_perf_p2t.tree->Fill(); + if (m_perf_t2p.tree) m_perf_t2p.tree->Fill(); + } + + void fillFinderAssoc(int event, const edm4hep::MCParticleCollection& mcParts, + const edm4hep::TrackCollection& finderTracks, + const std::unordered_map& hitToParticle, + const std::unordered_map>& hitsPerParticle) const { + + m_finder_p2t.clear(); + m_finder_t2p.clear(); + m_finder_p2t.event = event; + m_finder_t2p.event = event; + + std::vector> trackParticleCounts; + std::vector trackNHits; + trackParticleCounts.resize(finderTracks.size()); + trackNHits.resize(finderTracks.size(), 0); + + int tIdx = 0; + for (const auto& trk : finderTracks) { + trackNHits[tIdx] = (int)trk.getTrackerHits().size(); + for (const auto& h : trk.getTrackerHits()) { + const uint64_t hk = TrackingValidationHelpers::oidKey(h.getObjectID()); + auto it = hitToParticle.find(hk); + if (it == hitToParticle.end()) continue; + trackParticleCounts[tIdx][it->second] += 1; + } + ++tIdx; + } + + // track -> particles + for (int t = 0; t < (int)finderTracks.size(); ++t) { + m_finder_t2p.index.push_back(t); + m_finder_t2p.p.push_back(-1.f); + m_finder_t2p.pT.push_back(-1.f); + m_finder_t2p.nTrueHits.push_back(trackNHits[t]); + + std::vector parts; + std::vector sh; + std::vector effs; + std::vector purs; + + for (const auto& kv : trackParticleCounts[t]) { + const int pid = kv.first; + if (pid < 0) continue; + + const int shared = kv.second; + const int nTrackHits = trackNHits[t]; + const int nParticleHits = + hitsPerParticle.count(pid) ? (int)hitsPerParticle.at(pid).size() : 0; + + const float eff = (nParticleHits > 0) ? float(shared) / float(nParticleHits) : 0.f; + const float pur = (nTrackHits > 0) ? float(shared) / float(nTrackHits) : 0.f; + + parts.push_back(pid); + sh.push_back(shared); + effs.push_back(eff); + purs.push_back(pur); + } + + m_finder_t2p.assoc.push_back(parts); + m_finder_t2p.sharedHits.push_back(sh); + m_finder_t2p.matchEfficiency.push_back(effs); + m_finder_t2p.matchPurity.push_back(purs); + } + + // particle -> tracks + for (int p = 0; p < (int)mcParts.size(); ++p) { + const auto& mc = mcParts[p]; + if (mc.getGeneratorStatus() != 1) continue; + + auto itHits = hitsPerParticle.find(p); + if (itHits == hitsPerParticle.end() || itHits->second.empty()) continue; + + const auto& mom = mc.getMomentum(); + const float px = float(mom.x); + const float py = float(mom.y); + const float pz = float(mom.z); + const float pAbs = std::sqrt(px*px + py*py + pz*pz); + const float pT = std::sqrt(px*px + py*py); + const int nParticleHits = (int)itHits->second.size(); + + std::vector tracks; + std::vector sh; + std::vector effs; + std::vector purs; + + for (int t = 0; t < (int)finderTracks.size(); ++t) { + auto it = trackParticleCounts[t].find(p); + if (it == trackParticleCounts[t].end()) continue; + + const int shared = it->second; + const int nTrackHits = trackNHits[t]; + + const float eff = (nParticleHits > 0) ? float(shared) / float(nParticleHits) : 0.f; + const float pur = (nTrackHits > 0) ? float(shared) / float(nTrackHits) : 0.f; + + tracks.push_back(t); + sh.push_back(shared); + effs.push_back(eff); + purs.push_back(pur); + } + + m_finder_p2t.index.push_back(p); + m_finder_p2t.p.push_back(pAbs); + m_finder_p2t.pT.push_back(pT); + m_finder_p2t.nTrueHits.push_back(nParticleHits); + m_finder_p2t.assoc.push_back(tracks); + m_finder_p2t.sharedHits.push_back(sh); + m_finder_p2t.matchEfficiency.push_back(effs); + m_finder_p2t.matchPurity.push_back(purs); + } + + if (m_finder_p2t.tree) m_finder_p2t.tree->Fill(); + if (m_finder_t2p.tree) m_finder_t2p.tree->Fill(); +} + + // ---------- matching helper ---------- + int majorityParticleForTrack(const edm4hep::Track& trk, + const std::unordered_map& hitToParticle) const { + std::unordered_map counts; + for (const auto& h : trk.getTrackerHits()) { + const uint64_t hk = TrackingValidationHelpers::oidKey(h.getObjectID()); + auto it = hitToParticle.find(hk); + if (it == hitToParticle.end()) continue; + counts[it->second] += 1; + } + if (counts.empty()) return -1; + + int bestP = -1; + int bestN = -1; + for (const auto& kv : counts) { + if (kv.second > bestN) { + bestN = kv.second; + bestP = kv.first; + } + } + return bestP; + } + + // ---------- fitter trees ---------- + template + void fillFitterTrees(int event, + const edm4hep::MCParticleCollection& mcParts, + const edm4hep::TrackCollection& fittedTracks, + const std::unordered_map& hitToParticle, + const PerfectMapT& perfectAtIPByPid, + bool doPerfect) const { + + m_fit_vs_mc.clear(); + m_fit_vs_perfect.clear(); + m_fit_vs_mc.event = event; + m_fit_vs_perfect.event = event; + + const float NaN = std::numeric_limits::quiet_NaN(); + + int tIdx = 0; + for (const auto& trk : fittedTracks) { + edm4hep::TrackState stReco; + if (!TrackingValidationHelpers::getAtIPState(trk, stReco)) { + ++tIdx; + continue; + } + + const int pid = majorityParticleForTrack(trk, hitToParticle); + if (pid < 0 || pid >= (int)mcParts.size()) { + ++tIdx; + continue; + } + + const auto& mc = mcParts[pid]; + + // reco params (already in fitter convention) + TrackingValidationHelpers::HelixParams reco; + reco.D0 = float(stReco.D0); + reco.Z0 = float(stReco.Z0); + reco.phi = float(stReco.phi); + reco.omega = float(stReco.omega); + reco.tanLambda = float(stReco.tanLambda); + reco.pT = TrackingValidationHelpers::ptFromState(stReco, m_Bz.value()); + reco.p = TrackingValidationHelpers::momentumFromState(stReco, m_Bz.value()); + + // ref from MC using the SAME convention as fitter (PCA + phi0 + ZPCA + omega=a*B/pT) + const TrackingValidationHelpers::HelixParams refMC = TrackingValidationHelpers::truthFromMC_GenfitConvention(mc, m_Bz.value(), m_refX.value(), m_refY.value(), m_refZ.value()); + + + // --- vs MC --- + m_fit_vs_mc.track_index.push_back(tIdx); + m_fit_vs_mc.track_location.push_back(int(stReco.location)); + m_fit_vs_mc.resD0.push_back(reco.D0 - refMC.D0); + m_fit_vs_mc.resZ0.push_back(reco.Z0 - refMC.Z0); + m_fit_vs_mc.resPhi.push_back(TrackingValidationHelpers::wrapDeltaPhi(reco.phi, refMC.phi)); + m_fit_vs_mc.resOmega.push_back(reco.omega - refMC.omega); + m_fit_vs_mc.resTanL.push_back(reco.tanLambda - refMC.tanLambda); + m_fit_vs_mc.p_reco.push_back(reco.p); + m_fit_vs_mc.p_ref.push_back(refMC.p); + m_fit_vs_mc.pT_reco.push_back(reco.pT); + m_fit_vs_mc.pT_ref.push_back(refMC.pT); + + // --- vs perfect-fitted --- + if (doPerfect) { + auto it = perfectAtIPByPid.find(pid); + if (it != perfectAtIPByPid.end()) { + const auto& stPerf = it->second.st; + + TrackingValidationHelpers::HelixParams refP; + refP.D0 = float(stPerf.D0); + refP.Z0 = float(stPerf.Z0); + refP.phi = float(stPerf.phi); + refP.omega = float(stPerf.omega); + refP.tanLambda = float(stPerf.tanLambda); + refP.pT = TrackingValidationHelpers::ptFromState(stPerf, m_Bz.value()); + refP.p = TrackingValidationHelpers::momentumFromState(stPerf, m_Bz.value()); + + m_fit_vs_perfect.track_index.push_back(tIdx); + m_fit_vs_perfect.track_location.push_back(int(stReco.location)); + m_fit_vs_perfect.resD0.push_back(reco.D0 - refP.D0); + m_fit_vs_perfect.resZ0.push_back(reco.Z0 - refP.Z0); + m_fit_vs_perfect.resPhi.push_back(TrackingValidationHelpers::wrapDeltaPhi(reco.phi, refP.phi)); + m_fit_vs_perfect.resOmega.push_back(reco.omega - refP.omega); + m_fit_vs_perfect.resTanL.push_back(reco.tanLambda - refP.tanLambda); + m_fit_vs_perfect.p_reco.push_back(reco.p); + m_fit_vs_perfect.p_ref.push_back(refP.p); + m_fit_vs_perfect.pT_reco.push_back(reco.pT); + m_fit_vs_perfect.pT_ref.push_back(refP.pT); + } else { + m_fit_vs_perfect.track_index.push_back(tIdx); + m_fit_vs_perfect.track_location.push_back(int(stReco.location)); + m_fit_vs_perfect.resD0.push_back(NaN); + m_fit_vs_perfect.resZ0.push_back(NaN); + m_fit_vs_perfect.resPhi.push_back(NaN); + m_fit_vs_perfect.resOmega.push_back(NaN); + m_fit_vs_perfect.resTanL.push_back(NaN); + m_fit_vs_perfect.p_reco.push_back(reco.p); + m_fit_vs_perfect.p_ref.push_back(NaN); + m_fit_vs_perfect.pT_reco.push_back(reco.pT); + m_fit_vs_perfect.pT_ref.push_back(NaN); + } + } + + ++tIdx; + } + + if (m_fit_vs_mc.tree) m_fit_vs_mc.tree->Fill(); + if (m_fit_vs_perfect.tree) m_fit_vs_perfect.tree->Fill(); + } + +private: + mutable int m_evt = 0; + mutable bool m_warnedMissingPerfectInput = false; + + std::unique_ptr m_outFile; + + mutable AssocTree m_finder_p2t; + mutable AssocTree m_finder_t2p; + mutable AssocTree m_perf_p2t; + mutable AssocTree m_perf_t2p; + + mutable FitterTree m_fit_vs_mc; + mutable FitterTree m_fit_vs_perfect; +}; + +DECLARE_COMPONENT(TrackingValidation) \ No newline at end of file diff --git a/TrackingPerformance/test/inputFiles/SimpleGatrIDEAv3o1.onnx.md5 b/TrackingPerformance/test/inputFiles/SimpleGatrIDEAv3o1.onnx.md5 new file mode 100644 index 0000000..3df28d7 --- /dev/null +++ b/TrackingPerformance/test/inputFiles/SimpleGatrIDEAv3o1.onnx.md5 @@ -0,0 +1 @@ +303e5169a22383d74627a3c412707835 \ No newline at end of file diff --git a/TrackingPerformance/test/runTrackingValidation.py b/TrackingPerformance/test/runTrackingValidation.py new file mode 100644 index 0000000..3ac555a --- /dev/null +++ b/TrackingPerformance/test/runTrackingValidation.py @@ -0,0 +1,336 @@ +import os +import math + +from Gaudi.Configuration import INFO +from Configurables import EventDataSvc, GeoSvc, UniqueIDGenSvc, RndmGenSvc +from k4FWCore import IOSvc, ApplicationMgr +from k4FWCore.parseArgs import parser + +# -------------------- +# Arguments +# -------------------- +parser.add_argument("--inputFile", required=True, + help="Input EDM4hep ROOT file") +parser.add_argument("--modelPath", default="", + help="Path to the GGTF ONNX model (required only if --runFinder 1)") +parser.add_argument("--outputFile", default="out_reco.root", + help="Output EDM4hep ROOT file with reconstructed collections") +parser.add_argument("--validationFile", default="validation.root", + help="Output ROOT file written by TrackingValidation") + +parser.add_argument("--geom", + default=os.path.join(os.environ["K4GEO"], "FCCee/IDEA/compact/IDEA_o1_v03/IDEA_o1_v03.xml"), + help="Detector geometry XML file") + +# Pipeline control +parser.add_argument("--runDigi", type=int, default=1, + help="0=skip digitization, 1=run digitization") +parser.add_argument("--runFinder", type=int, default=1, + help="0=skip track finder, 1=run track finder") +parser.add_argument("--runFitter", type=int, default=1, + help="0=skip reco fitter, 1=run reco fitter") +parser.add_argument("--runPerfectTracking", type=int, default=1, + help="0=skip perfect tracking/perfect fitter, 1=run them") +parser.add_argument("--runValidation", type=int, default=1, + help="0=skip validation, 1=run TrackingValidation") +parser.add_argument("--useDCH", type=int, default=1, + help="0=disable DCH collections, 1=use DCH collections") + +# Validation control +parser.add_argument("--mode", type=int, default=0, + help="Validation mode: 0=Full, 1=FinderOnly, 2=FitterOnly") +parser.add_argument("--doPerfectFit", type=int, default=1, + help="0=do not fill fitter_vs_perfect, 1=fill fitter_vs_perfect") +parser.add_argument("--finderEfficiencyDefinition", type=int, default=1, + help="1=purity-based definition, 2=purity+efficiency >= 0.5 definition") +parser.add_argument("--finderPurityThreshold", type=float, default=0.75, + help="Purity threshold used when FinderEfficiencyDefinition = 1") + +args = parser.parse_args() + +if args.runFinder == 1 and not args.modelPath: + parser.error("--modelPath is required when --runFinder 1") + +if args.runValidation == 0 and args.doPerfectFit == 1: + print("WARNING: --doPerfectFit is ignored when --runValidation 0") + +if args.runDigi == 0 and args.runFinder == 1: + print("WARNING: --runFinder 1 with --runDigi 0 assumes digi collections are already present in the input file") + +if args.runFinder == 0 and args.runFitter == 1: + print("WARNING: --runFitter 1 with --runFinder 0 assumes finder-track collections are already present in the input file") + +if args.runPerfectTracking == 0 and args.doPerfectFit == 1: + print("WARNING: --doPerfectFit 1 with --runPerfectTracking 0 assumes PerfectFittedTracks is already present in the input file") + +if all(flag == 0 for flag in [args.runDigi, args.runFinder, args.runFitter, args.runPerfectTracking, args.runValidation]): + parser.error("Nothing to do: all run flags are set to 0") + +# -------------------- +# Fixed collection names +# -------------------- +MC_COLLECTION = "MCParticles" + +PLANAR_LINK_COLLECTIONS = [ + "SiWrBSimDigiLinks", + "SiWrDSimDigiLinks", + "VTXBSimDigiLinks", + "VTXDSimDigiLinks", +] + +DCH_LINK_COLLECTIONS = ["DCH_DigiSimAssociationCollection"] if args.useDCH == 1 else [] + +PLANAR_DIGI_COLLECTIONS = [ + "VTXBDigis", + "VTXDDigis", + "SiWrBDigis", + "SiWrDDigis", +] + +DCH_DIGI_COLLECTIONS = ["DCH_DigiCollection"] if args.useDCH == 1 else [] + +FINDER_TRACK_COLLECTION = "GGTFTracks" +FITTED_TRACK_COLLECTION = "FittedTracks" +PERFECT_TRACK_COLLECTION = "PerfectTracks" +PERFECT_FITTED_TRACK_COLLECTION = "PerfectFittedTracks" + +# -------------------- +# IO +# -------------------- +io = IOSvc("IOSvc") +io.Input = args.inputFile +io.Output = args.outputFile + +# -------------------- +# Geometry +# -------------------- +geoservice = GeoSvc("GeoSvc") +geoservice.detectors = [args.geom] +geoservice.EnableGeant4Geo = False +geoservice.OutputLevel = INFO + +# -------------------- +# Algorithm sequence +# -------------------- +TopAlg = [] + +# -------------------- +# Digitizers +# -------------------- +if args.runDigi == 1: + from Configurables import DDPlanarDigi, DCHdigi_v02 + + innerVertexResolution_x = 0.003 + innerVertexResolution_y = 0.003 + innerVertexResolution_t = 1000 + + outerVertexResolution_x = 0.050 / math.sqrt(12) + outerVertexResolution_y = 0.150 / math.sqrt(12) + outerVertexResolution_t = 1000 + + siWrapperResolution_x = 0.050 / math.sqrt(12) + siWrapperResolution_y = 1.0 / math.sqrt(12) + siWrapperResolution_t = 0.040 + + vtxb_digitizer = DDPlanarDigi("VTXBdigitizer") + vtxb_digitizer.SubDetectorName = "Vertex" + vtxb_digitizer.IsStrip = False + vtxb_digitizer.ResolutionU = [ + innerVertexResolution_x, + innerVertexResolution_x, + innerVertexResolution_x, + outerVertexResolution_x, + outerVertexResolution_x, + ] + vtxb_digitizer.ResolutionV = [ + innerVertexResolution_y, + innerVertexResolution_y, + innerVertexResolution_y, + outerVertexResolution_y, + outerVertexResolution_y, + ] + vtxb_digitizer.ResolutionT = [ + innerVertexResolution_t, + innerVertexResolution_t, + innerVertexResolution_t, + outerVertexResolution_t, + outerVertexResolution_t, + ] + vtxb_digitizer.SimTrackHitCollectionName = ["VertexBarrelCollection"] + vtxb_digitizer.SimTrkHitRelCollection = ["VTXBSimDigiLinks"] + vtxb_digitizer.TrackerHitCollectionName = ["VTXBDigis"] + vtxb_digitizer.ForceHitsOntoSurface = True + + vtxd_digitizer = DDPlanarDigi("VTXDdigitizer") + vtxd_digitizer.SubDetectorName = "Vertex" + vtxd_digitizer.IsStrip = False + vtxd_digitizer.ResolutionU = [ + outerVertexResolution_x, + outerVertexResolution_x, + outerVertexResolution_x, + ] + vtxd_digitizer.ResolutionV = [ + outerVertexResolution_y, + outerVertexResolution_y, + outerVertexResolution_y, + ] + vtxd_digitizer.ResolutionT = [ + outerVertexResolution_t, + outerVertexResolution_t, + outerVertexResolution_t, + ] + vtxd_digitizer.SimTrackHitCollectionName = ["VertexEndcapCollection"] + vtxd_digitizer.SimTrkHitRelCollection = ["VTXDSimDigiLinks"] + vtxd_digitizer.TrackerHitCollectionName = ["VTXDDigis"] + vtxd_digitizer.ForceHitsOntoSurface = True + + siwrb_digitizer = DDPlanarDigi("SiWrBdigitizer") + siwrb_digitizer.SubDetectorName = "SiWrB" + siwrb_digitizer.IsStrip = False + siwrb_digitizer.ResolutionU = [siWrapperResolution_x, siWrapperResolution_x] + siwrb_digitizer.ResolutionV = [siWrapperResolution_y, siWrapperResolution_y] + siwrb_digitizer.ResolutionT = [siWrapperResolution_t, siWrapperResolution_t] + siwrb_digitizer.SimTrackHitCollectionName = ["SiWrBCollection"] + siwrb_digitizer.SimTrkHitRelCollection = ["SiWrBSimDigiLinks"] + siwrb_digitizer.TrackerHitCollectionName = ["SiWrBDigis"] + siwrb_digitizer.ForceHitsOntoSurface = True + + siwrd_digitizer = DDPlanarDigi("SiWrDdigitizer") + siwrd_digitizer.SubDetectorName = "SiWrD" + siwrd_digitizer.IsStrip = False + siwrd_digitizer.ResolutionU = [siWrapperResolution_x, siWrapperResolution_x] + siwrd_digitizer.ResolutionV = [siWrapperResolution_y, siWrapperResolution_y] + siwrd_digitizer.ResolutionT = [siWrapperResolution_t, siWrapperResolution_t] + siwrd_digitizer.SimTrackHitCollectionName = ["SiWrDCollection"] + siwrd_digitizer.SimTrkHitRelCollection = ["SiWrDSimDigiLinks"] + siwrd_digitizer.TrackerHitCollectionName = ["SiWrDDigis"] + siwrd_digitizer.ForceHitsOntoSurface = True + + TopAlg += [vtxb_digitizer, vtxd_digitizer, siwrb_digitizer, siwrd_digitizer] + + if args.useDCH == 1: + dch_digitizer = DCHdigi_v02( + "DCHdigi2", + InputSimHitCollection=["DCHCollection"], + OutputDigihitCollection=["DCH_DigiCollection"], + OutputLinkCollection=["DCH_DigiSimAssociationCollection"], + DCH_name="DCH_v2", + zResolution_mm=30.0, + xyResolution_mm=0.1, + Deadtime_ns=400.0, + GasType=0, + ReadoutWindowStartTime_ns=1.0, + ReadoutWindowDuration_ns=450.0, + DriftVelocity_um_per_ns=-1.0, + SignalVelocity_mm_per_ns=200.0, + OutputLevel=INFO, + ) + TopAlg += [dch_digitizer] + +# -------------------- +# Track finder +# -------------------- +if args.runFinder == 1: + from Configurables import GGTFTrackFinder + + ggtf = GGTFTrackFinder( + "GGTFTrackFinder", + InputPlanarHitCollections=PLANAR_DIGI_COLLECTIONS, + InputWireHitCollections=DCH_DIGI_COLLECTIONS, + OutputTracksGGTF=[FINDER_TRACK_COLLECTION], + ModelPath=args.modelPath, + Tbeta=0.6, + Td=0.3, + OutputLevel=INFO, + ) + TopAlg += [ggtf] + +# -------------------- +# Reco fitter +# -------------------- +if args.runFitter == 1: + from Configurables import GenfitTrackFitter + + reco_fitter = GenfitTrackFitter("RecoTrackFitter") + reco_fitter.InputTracks = [FINDER_TRACK_COLLECTION] + reco_fitter.OutputFittedTracks = [FITTED_TRACK_COLLECTION] + reco_fitter.RunSingleEvaluation = True + reco_fitter.UseBrems = False + reco_fitter.BetaInit = 100.0 + reco_fitter.BetaFinal = 0.1 + reco_fitter.BetaSteps = 10 + reco_fitter.InitializationType = 1 + reco_fitter.SkipTrackOrdering = False + reco_fitter.SkipUnmatchedTracks = False + reco_fitter.OutputLevel = INFO + + TopAlg += [reco_fitter] + +# -------------------- +# Perfect tracking + perfect fitter +# -------------------- +if args.runPerfectTracking == 1: + from Configurables import PerfectTrackFinder, GenfitTrackFitter + + perfect = PerfectTrackFinder("PerfectTrackFinder") + perfect.InputMCParticles = [MC_COLLECTION] + perfect.InputPlanarHitCollections = PLANAR_LINK_COLLECTIONS + perfect.InputWireHitCollections = DCH_LINK_COLLECTIONS + perfect.OutputPerfectTracks = [PERFECT_TRACK_COLLECTION] + perfect.OutputLevel = INFO + + perfect_fitter = GenfitTrackFitter("PerfectTrackFitter") + perfect_fitter.InputTracks = [PERFECT_TRACK_COLLECTION] + perfect_fitter.OutputFittedTracks = [PERFECT_FITTED_TRACK_COLLECTION] + perfect_fitter.RunSingleEvaluation = True + perfect_fitter.UseBrems = False + perfect_fitter.BetaInit = 100.0 + perfect_fitter.BetaFinal = 0.1 + perfect_fitter.BetaSteps = 10 + perfect_fitter.InitializationType = 1 + perfect_fitter.SkipTrackOrdering = False + perfect_fitter.SkipUnmatchedTracks = False + perfect_fitter.OutputLevel = INFO + + TopAlg += [perfect, perfect_fitter] + +# -------------------- +# Validation consumer +# -------------------- +if args.runValidation == 1: + from Configurables import TrackingValidation + + val = TrackingValidation("TrackingValidation") + val.OutputFile = args.validationFile + val.Mode = args.mode + + # Fixed fitter/truth convention settings + val.Bz = 2.0 + val.RefPointX = 0.0 + val.RefPointY = 0.0 + val.RefPointZ = 0.0 + + val.DoPerfectFit = bool(args.doPerfectFit) + val.FinderEfficiencyDefinition = args.finderEfficiencyDefinition + val.FinderPurityThreshold = args.finderPurityThreshold + + val.MCParticles = [MC_COLLECTION] + val.PlanarLinks = PLANAR_LINK_COLLECTIONS + val.DCHLinks = DCH_LINK_COLLECTIONS + val.FinderTracks = [FINDER_TRACK_COLLECTION] + val.FittedTracks = [FITTED_TRACK_COLLECTION] + val.PerfectFittedTracks = [PERFECT_FITTED_TRACK_COLLECTION] if args.doPerfectFit == 1 else [] + val.OutputLevel = INFO + + TopAlg += [val] + +# -------------------- +# AppMgr +# -------------------- +ApplicationMgr( + TopAlg=TopAlg, + EvtSel="NONE", + EvtMax=-1, + ExtSvc=[geoservice, EventDataSvc("EventDataSvc"), UniqueIDGenSvc("uidSvc"), RndmGenSvc()], + OutputLevel=INFO, +) \ No newline at end of file diff --git a/TrackingPerformance/test/testTrackingValidation.sh b/TrackingPerformance/test/testTrackingValidation.sh new file mode 100755 index 0000000..379b282 --- /dev/null +++ b/TrackingPerformance/test/testTrackingValidation.sh @@ -0,0 +1,188 @@ +#!/bin/bash +set -euo pipefail + +SCRIPT_DIR="$(cd "$(dirname "$0")" && pwd)" +MODEL_FILE="${1:-}" + +# Optional shell-level control of the simulation step +TRACKINGPERF_RUN_SIM="${TRACKINGPERF_RUN_SIM:-1}" +TRACKINGPERF_INPUT_FILE_OVERRIDE="${TRACKINGPERF_INPUT_FILE_OVERRIDE:-}" + +if [ -z "${MODEL_FILE}" ]; then + echo "ERROR: missing ONNX model path argument" + echo "Usage: $0 /full/path/to/SimpleGatrIDEAv3o1.onnx" + exit 1 +fi + +if [ ! -f "${MODEL_FILE}" ]; then + echo "ERROR: ONNX model file not found: ${MODEL_FILE}" + exit 1 +fi + +if [ -z "${K4GEO:-}" ]; then + echo "ERROR: K4GEO is not set" + exit 1 +fi + +if [ ! -d "${K4GEO}" ]; then + echo "ERROR: K4GEO does not point to a valid directory: ${K4GEO}" + exit 1 +fi + +if ! command -v k4run >/dev/null 2>&1; then + echo "ERROR: k4run not found in PATH" + exit 1 +fi + +if [ "${TRACKINGPERF_RUN_SIM}" -eq 1 ]; then + if ! command -v ddsim >/dev/null 2>&1; then + echo "ERROR: ddsim not found in PATH" + exit 1 + fi + + if ! command -v curl >/dev/null 2>&1; then + echo "ERROR: curl not found in PATH" + exit 1 + fi +fi + +XML_FILE="${K4GEO}/FCCee/IDEA/compact/IDEA_o1_v03/IDEA_o1_v03.xml" +RUN_FILE="${SCRIPT_DIR}/runTrackingValidation.py" +VAL_FILE="${SCRIPT_DIR}/validation_output_test.root" + +if [ ! -f "${XML_FILE}" ]; then + echo "ERROR: geometry XML file not found: ${XML_FILE}" + exit 1 +fi + +if [ ! -f "${RUN_FILE}" ]; then + echo "ERROR: tracking validation run file not found: ${RUN_FILE}" + exit 1 +fi + +TMPDIR="$(mktemp -d "${TMPDIR:-/tmp}/tracking_validation.XXXXXX")" +trap 'rm -rf "${TMPDIR}"' EXIT + +STEERING_FILE="${TMPDIR}/SteeringFile_IDEA_o1_v03.py" +SIM_FILE="${TMPDIR}/out_sim_edm4hep.root" +RECO_FILE="${TMPDIR}/out_reco.root" + +rm -f "${VAL_FILE}" + +N_EVENTS=5 +SEED=42 + +# Decide which input file to pass to k4run +if [ "${TRACKINGPERF_RUN_SIM}" -eq 1 ]; then + INPUT_FILE="${SIM_FILE}" +else + if [ -z "${TRACKINGPERF_INPUT_FILE_OVERRIDE}" ]; then + echo "ERROR: TRACKINGPERF_RUN_SIM=0 but TRACKINGPERF_INPUT_FILE_OVERRIDE is empty" + echo "Please provide an existing EDM4hep file, e.g." + echo " TRACKINGPERF_RUN_SIM=0 TRACKINGPERF_INPUT_FILE_OVERRIDE=/path/to/input.root $0 /path/to/model.onnx" + exit 1 + fi + + if [ ! -f "${TRACKINGPERF_INPUT_FILE_OVERRIDE}" ]; then + echo "ERROR: TRACKINGPERF_INPUT_FILE_OVERRIDE does not exist: ${TRACKINGPERF_INPUT_FILE_OVERRIDE}" + exit 1 + fi + + INPUT_FILE="${TRACKINGPERF_INPUT_FILE_OVERRIDE}" +fi + +echo "=== Test configuration ===" +echo "Script dir: ${SCRIPT_DIR}" +echo "Temporary dir: ${TMPDIR}" +echo "Geometry XML: ${XML_FILE}" +echo "Run script: ${RUN_FILE}" +echo "ONNX model: ${MODEL_FILE}" +echo "TRACKINGPERF_RUN_SIM: ${TRACKINGPERF_RUN_SIM}" +echo "Input file: ${INPUT_FILE}" +echo "Reco file: ${RECO_FILE}" +echo "Validation file: ${VAL_FILE}" +echo "Events: ${N_EVENTS}" +echo "Seed: ${SEED}" +echo "runDigi: 1" +echo "runFinder: 1" +echo "runFitter: 1" +echo "runPerfectTracking: 1" +echo "runValidation: 1" +echo "useDCH: 1" +echo "mode: 0" +echo "doPerfectFit: 1" +echo "finderEfficiencyDefinition: 1" +echo "finderPurityThreshold: 0.75" + +if [ "${TRACKINGPERF_RUN_SIM}" -eq 1 ]; then + echo "=== Downloading DDSim steering file ===" + curl -L \ + -o "${STEERING_FILE}" \ + https://raw.githubusercontent.com/key4hep/k4geo/master/example/SteeringFile_IDEA_o1_v03.py + + if [ ! -f "${STEERING_FILE}" ]; then + echo "ERROR: failed to download DDSim steering file" + exit 1 + fi + + echo "=== Step 1: DDSim ===" + ddsim \ + --steeringFile "${STEERING_FILE}" \ + --compactFile "${XML_FILE}" \ + -G \ + --gun.particle mu- \ + --gun.energy "5*GeV" \ + --gun.distribution uniform \ + --gun.thetaMin "89*deg" \ + --gun.thetaMax "89*deg" \ + --gun.phiMin "0*deg" \ + --gun.phiMax "0*deg" \ + --random.seed "${SEED}" \ + --numberOfEvents "${N_EVENTS}" \ + --outputFile "${SIM_FILE}" + + if [ ! -f "${SIM_FILE}" ]; then + echo "ERROR: simulation output was not created: ${SIM_FILE}" + exit 1 + fi +else + echo "=== Step 1: DDSim skipped (TRACKINGPERF_RUN_SIM=0) ===" +fi + +echo "=== Step 2: digitization + tracking + validation ===" +k4run "${RUN_FILE}" \ + --inputFile "${INPUT_FILE}" \ + --modelPath "${MODEL_FILE}" \ + --outputFile "${RECO_FILE}" \ + --validationFile "${VAL_FILE}" \ + --geom "${XML_FILE}" \ + --runDigi 1 \ + --runFinder 1 \ + --runFitter 1 \ + --runPerfectTracking 1 \ + --runValidation 1 \ + --useDCH 1 \ + --mode 0 \ + --doPerfectFit 1 \ + --finderEfficiencyDefinition 1 \ + --finderPurityThreshold 0.75 + +if [ ! -f "${RECO_FILE}" ]; then + echo "ERROR: reconstruction output was not created: ${RECO_FILE}" + exit 1 +fi + +if [ ! -f "${VAL_FILE}" ]; then + echo "ERROR: validation output was not created: ${VAL_FILE}" + exit 1 +fi + +echo "=== Step 3: check outputs ===" +if [ "${TRACKINGPERF_RUN_SIM}" -eq 1 ]; then + test -f "${SIM_FILE}" +fi +test -f "${RECO_FILE}" +test -f "${VAL_FILE}" + +echo "Test completed successfully." +echo "Validation file: ${VAL_FILE}" \ No newline at end of file