diff --git a/README.md b/README.md index 08b1c5bf42..7fceba2698 100644 --- a/README.md +++ b/README.md @@ -8,11 +8,11 @@ Demonstrator tracking chain for accelerators. | ------------------ | ---------------------- | --- | ---- | ---- | ------- | | **Clusterization** | CCL | ✅ | ✅ | ✅ | ✅ | | | Measurement creation | ✅ | ✅ | ✅ | ✅ | -| | Spacepoint formation | ✅ | ✅ | ✅ | ⚪ | -| **Track finding** | Spacepoint binning | ✅ | ✅ | ✅ | ⚪ | +| **Seeding** | Spacepoint formation | ✅ | ✅ | ✅ | ⚪ | +| | Spacepoint binning | ✅ | ✅ | ✅ | ⚪ | | | Seed finding | ✅ | ✅ | ✅ | ⚪ | | | Track param estimation | ✅ | ✅ | ✅ | ⚪ | -| | Combinatorial KF | ✅ | ✅ | 🟡 | ⚪ | +| **Track finding** | Combinatorial KF | ✅ | ✅ | 🟡 | ⚪ | | **Track fitting** | KF | ✅ | ✅ | ✅ | ⚪ | ✅: exists, 🟡: work started, ⚪: work not started yet @@ -255,7 +255,7 @@ cmake --build /bin/traccc_seeding_example_cuda --input_directory=detray_simulation/toy_detector/n_particles_2000/ --check_performance=true --detector_file=/toy_detector_geometry.json --material_file=/toy_detector_homogeneous_material.json --grid_file=/toy_detector_surface_grids.json --event=1 --track_candidates_range=3:10 --constraint-step-size-mm=1000 --run_cpu=1 -/bin/traccc_truth_finding_example_cuda --input_directory=detray_simulation/toy_detector/n_particles_2000/ --check_performance=true --detector_file=/toy_detector_geometry.json --material_file=/toy_detector_homogeneous_material.json --grid_file=/toy_detector_surface_grids.json --event=1 --track_candidates_range=3:10 --constraint-step-size-mm=1 --run_cpu=1 +/bin/traccc_truth_finding_example_cuda --input_directory=detray_simulation/toy_detector/n_particles_2000/ --check_performance=true --detector_file=/toy_detector_geometry.json --material_file=/toy_detector_homogeneous_material.json --grid_file=/toy_detector_surface_grids.json --event=1 --track_candidates_range=3:30 --constraint-step-size-mm=1 --run_cpu=1 ``` ```sh @@ -280,6 +280,9 @@ cmake --build Users can generate muon-like particle simulation data by running following example commands: ```sh +# Generate telescope geometry data +/bin/traccc_simulate_telescope --gen-vertex-xyz-mm=0:0:0 --gen-vertex-xyz-std-mm=0:0:0 --gen-mom-gev=100:100 --gen-phi-degree=0:0 --events=10 --gen-nparticles=2000 --output_directory=detray_simulation/telescope_detector/n_particles_2000/ --gen-eta=1:3 + # Generate toy geometry data /bin/traccc_simulate_toy_detector --gen-vertex-xyz-mm=0:0:0 --gen-vertex-xyz-std-mm=0:0:0 --gen-mom-gev=100:100 --gen-phi-degree=0:360 --events=10 --gen-nparticles=2000 --output_directory=detray_simulation/toy_detector/n_particles_2000/ --gen-eta=-3:3 --constraint-step-size-mm=1 @@ -287,11 +290,17 @@ Users can generate muon-like particle simulation data by running following examp /bin/traccc_simulate_wire_chamber --gen-vertex-xyz-mm=0:0:0 --gen-vertex-xyz-std-mm=0:0:0 --gen-mom-gev=2:2 --gen-phi-degree=0:360 --events=10 --gen-nparticles=100 --output_directory=detray_simulation/wire_chamber/n_particles_100/ --gen-eta=-1:1 --constraint-step-size-mm=1 ``` -The simulation will also generate the detector json files (geometry, material and surface_grid) in the current directory. It is user's responsibility to move them to an appropriate place (say, ) and match it to the input file arguments of the reconstruction examples. +The simulation will also generate the detector json files (geometry, material and surface_grid) in the current directory. It is user's responsibility to move them to an appropriate place (e.g. ``) and match them to the input file arguments of reconstruction chains. + +There are three types of partial reconstruction chain users can operate: `seeding_example`, `truth_finding_example`, and `truth_fitting_example` where their algorithm coverages are shown in the table below. Each of them starts from truth measurements, truth seeds, and truth tracks, respectively. -Currently, there are two types of partial reconstruction chain users can operate: seeding_example and truth_finding_example. seeding_example takes the truth measurement input w/o clusterization and it goes through seeding, finding and fitting, which generate performance root files, respectively. On the other hand, truth_finding_examples starts from the truth initial parameter of particles (no duplicate seeds; i.e. the number of seeds for CKF = the number of truth particles) and run track finding and track fitting. +| Category | Clusterization | Seeding | Track finding | Track fitting | +| ----------------------- | -------------- | ------- | ------------- | ------------- | +| `seeding_example` | | ✅ | ✅ | ✅ | +| `truth_finding_example` | | | ✅ | ✅ | +| `truth_fitting_example` | | | | ✅ | -The dirft chamber will not produce meaningful results with seeding_examples as the current seeding algorithm is only designed for 2D measurement objects. Truth finding works OK in general but the combinatoric explosion can occur for a few unlucky events, leading to poor pull value distributions. +The dirft chamber will not produce meaningful results with `seeding_example` as the current seeding algorithm is only designed for 2D measurement objects. Truth finding works OK in general but the combinatoric explosion can occur for a few unlucky events, leading to poor pull value distributions. The followings are example commands: ```sh # Run cuda seeding example for toy geometry @@ -308,7 +317,12 @@ The dirft chamber will not produce meaningful results with seeding_examples as t /bin/traccc_truth_finding_example_cuda --input_directory=detray_simulation/wire_chamber/n_particles_100/ --check_performance=true --detector_file=/wire_chamber_geometry.json --material_file=/wire_chamber_homogeneous_material.json --grid_file=/wire_chamber_surface_grids.json --event=10 --track_candidates_range=6:30 --constraint-step-size-mm=1 --run_cpu=1 ``` -Users can open the performance root files (with --check_performance=true) and draw the histograms. +```sh +# Run cpu truth fitting example for drift chamber +/bin/traccc_truth_fitting_example --input_directory=detray_simulation/wire_chamber/n_particles_2000_100GeV/ --check_performance=true --detector_file=/wire_chamber_geometry.json --material_file=/wire_chamber_homogeneous_material.json --grid_file=/wire_chamber_surface_grids.json --event=10 --constraint-step-size-mm=1 +``` + +Users can open the performance root files (with `--check_performance=true`) and draw the histograms. ```sh $ root -l performance_track_finding.root diff --git a/core/include/traccc/fitting/kalman_filter/kalman_fitter.hpp b/core/include/traccc/fitting/kalman_filter/kalman_fitter.hpp index 2244a25428..fac8e42429 100644 --- a/core/include/traccc/fitting/kalman_filter/kalman_fitter.hpp +++ b/core/include/traccc/fitting/kalman_filter/kalman_fitter.hpp @@ -211,6 +211,7 @@ class kalman_fitter { auto& last = track_states.back(); last.smoothed().set_vector(last.filtered().vector()); last.smoothed().set_covariance(last.filtered().covariance()); + last.smoothed_chi2() = last.filtered_chi2(); for (typename vector_type< track_state>::reverse_iterator it = diff --git a/examples/options/include/traccc/options/telescope_detector_options.hpp b/examples/options/include/traccc/options/telescope_detector_options.hpp new file mode 100644 index 0000000000..f6ece378b3 --- /dev/null +++ b/examples/options/include/traccc/options/telescope_detector_options.hpp @@ -0,0 +1,62 @@ +/** TRACCC library, part of the ACTS project (R&D line) + * + * (c) 2023 CERN for the benefit of the ACTS project + * + * Mozilla Public License Version 2.0 + */ + +// Project include(s). +#include "traccc/options/options.hpp" + +// Detray include(s). +#include "detray/definitions/units.hpp" + +// Boost +#include + +namespace traccc { + +namespace po = boost::program_options; + +template +struct telescope_detector_options { + bool empty_material; + unsigned int n_planes; + scalar_t thickness; + scalar_t spacing; + scalar_t smearing; + scalar_t half_length; + + telescope_detector_options(po::options_description& desc) { + desc.add_options()("empty_material", + po::value()->default_value(false), + "Build detector without materials"); + desc.add_options()("n_planes", + po::value()->default_value(9), + "Number of planes"); + desc.add_options()("thickness", + po::value()->default_value(0.5f), + "Slab thickness in [mm]"); + desc.add_options()("spacing", + po::value()->default_value(20.f), + "Space between planes in [mm]"); + desc.add_options()("smearing", + po::value()->default_value(50.f), + "Measurement smearing in [um]"); + desc.add_options()("half_length", + po::value()->default_value(1000000.f), + "Half length of plane [mm]"); + } + + void read(const po::variables_map& vm) { + empty_material = vm["empty_material"].as(); + n_planes = vm["n_planes"].as(); + thickness = vm["thickness"].as() * detray::unit::mm; + spacing = vm["spacing"].as() * detray::unit::mm; + smearing = vm["smearing"].as() * detray::unit::um; + half_length = + vm["half_length"].as() * detray::unit::mm; + } +}; + +} // namespace traccc \ No newline at end of file diff --git a/examples/run/cpu/CMakeLists.txt b/examples/run/cpu/CMakeLists.txt index 11e0331a1b..89e0c6e989 100644 --- a/examples/run/cpu/CMakeLists.txt +++ b/examples/run/cpu/CMakeLists.txt @@ -16,6 +16,10 @@ traccc_add_executable( truth_finding_example "truth_finding_example.cpp" LINK_LIBRARIES vecmem::core detray::utils traccc::core traccc::io traccc::performance traccc::options) +traccc_add_executable( truth_fitting_example "truth_fitting_example.cpp" + LINK_LIBRARIES vecmem::core detray::io detray::utils traccc::core + traccc::io traccc::performance traccc::options) + traccc_add_executable( ccl_example "ccl_example.cpp" LINK_LIBRARIES vecmem::core traccc::core traccc::io) diff --git a/examples/run/cpu/truth_fitting_example.cpp b/examples/run/cpu/truth_fitting_example.cpp new file mode 100644 index 0000000000..0bfe028fdc --- /dev/null +++ b/examples/run/cpu/truth_fitting_example.cpp @@ -0,0 +1,165 @@ +/** TRACCC library, part of the ACTS project (R&D line) + * + * (c) 2023 CERN for the benefit of the ACTS project + * + * Mozilla Public License Version 2.0 + */ + +// Project include(s). +#include "traccc/definitions/common.hpp" +#include "traccc/definitions/primitives.hpp" +#include "traccc/fitting/fitting_algorithm.hpp" +#include "traccc/fitting/kalman_filter/kalman_fitter.hpp" +#include "traccc/io/read_geometry.hpp" +#include "traccc/io/read_measurements.hpp" +#include "traccc/io/utils.hpp" +#include "traccc/options/common_options.hpp" +#include "traccc/options/handle_argument_errors.hpp" +#include "traccc/options/propagation_options.hpp" +#include "traccc/resolution/fitting_performance_writer.hpp" +#include "traccc/utils/seed_generator.hpp" + +// Detray include(s). +#include "detray/core/detector.hpp" +#include "detray/core/detector_metadata.hpp" +#include "detray/detectors/bfield.hpp" +#include "detray/io/common/detector_reader.hpp" +#include "detray/propagator/navigator.hpp" +#include "detray/propagator/propagator.hpp" +#include "detray/propagator/rk_stepper.hpp" + +// VecMem include(s). +#include + +// System include(s). +#include +#include +#include + +using namespace traccc; +namespace po = boost::program_options; + +// The main routine +// +int main(int argc, char* argv[]) { + // Set up the program options + po::options_description desc("Allowed options"); + + // Add options + desc.add_options()("help,h", "Give some help with the program's options"); + traccc::common_options common_opts(desc); + traccc::propagation_options propagation_opts(desc); + + po::variables_map vm; + po::store(po::parse_command_line(argc, argv, desc), vm); + + // Check errors + traccc::handle_argument_errors(vm, desc); + + // Read options + common_opts.read(vm); + propagation_opts.read(vm); + + std::cout << "Running " << argv[0] << " " << common_opts.input_directory + << " " << common_opts.events << std::endl; + + /// Type declarations + using host_detector_type = detray::detector; + + using b_field_t = covfie::field; + using rk_stepper_type = + detray::rk_stepper>; + + using host_navigator_type = detray::navigator; + using host_fitter_type = + traccc::kalman_fitter; + + // Memory resources used by the application. + vecmem::host_memory_resource host_mr; + + // Performance writer + traccc::fitting_performance_writer fit_performance_writer( + traccc::fitting_performance_writer::config{}); + + /***************************** + * Build a geometry + *****************************/ + + // B field value and its type + // @TODO: Set B field as argument + const traccc::vector3 B{0, 0, 2 * detray::unit::T}; + auto field = detray::bfield::create_const_field(B); + + // Read the detector + detray::io::detector_reader_config reader_cfg{}; + reader_cfg + .add_file(traccc::io::data_directory() + common_opts.detector_file) + .add_file(traccc::io::data_directory() + common_opts.material_file) + .add_file(traccc::io::data_directory() + common_opts.grid_file); + + const auto [host_det, names] = + detray::io::read_detector(host_mr, reader_cfg); + + /***************************** + * Do the reconstruction + *****************************/ + + /// Standard deviations for seed track parameters + static constexpr std::array stddevs = { + 0.03 * detray::unit::mm, + 0.03 * detray::unit::mm, + 0.017, + 0.017, + 0.01 / detray::unit::GeV, + 1 * detray::unit::ns}; + + // Fitting algorithm object + typename traccc::fitting_algorithm::config_type fit_cfg; + fit_cfg.step_constraint = propagation_opts.step_constraint; + traccc::fitting_algorithm host_fitting(fit_cfg); + + // Seed generator + traccc::seed_generator sg(host_det, stddevs); + + // Iterate over events + for (unsigned int event = common_opts.skip; + event < common_opts.events + common_opts.skip; ++event) { + + // Truth Track Candidates + traccc::event_map2 evt_map2(event, common_opts.input_directory, + common_opts.input_directory, + common_opts.input_directory); + + traccc::track_candidate_container_types::host truth_track_candidates = + evt_map2.generate_truth_candidates(sg, host_mr); + + // Run fitting + auto track_states = + host_fitting(host_det, field, truth_track_candidates); + + std::cout << "Number of fitted tracks: " << track_states.size() + << std::endl; + + const unsigned int n_fitted_tracks = track_states.size(); + + if (common_opts.check_performance) { + + for (unsigned int i = 0; i < n_fitted_tracks; i++) { + const auto& trk_states_per_track = track_states.at(i).items; + + const auto& fit_info = track_states[i].header; + + fit_performance_writer.write(trk_states_per_track, fit_info, + host_det, evt_map2); + } + } + } + + if (common_opts.check_performance) { + fit_performance_writer.finalize(); + } + + return 1; +} diff --git a/examples/run/cuda/CMakeLists.txt b/examples/run/cuda/CMakeLists.txt index 444167c314..7e5560c537 100644 --- a/examples/run/cuda/CMakeLists.txt +++ b/examples/run/cuda/CMakeLists.txt @@ -22,6 +22,10 @@ traccc_add_executable( truth_finding_example_cuda "truth_finding_example_cuda.cp LINK_LIBRARIES vecmem::core vecmem::cuda traccc::io traccc::performance traccc::core traccc::device_common traccc::cuda traccc::options ) +traccc_add_executable( truth_fitting_example_cuda "truth_fitting_example_cuda.cpp" + LINK_LIBRARIES vecmem::core vecmem::cuda traccc::io traccc::performance + traccc::core traccc::device_common traccc::cuda + traccc::options ) # # Set up the "throughput applications". # diff --git a/examples/run/cuda/truth_finding_example_cuda.cpp b/examples/run/cuda/truth_finding_example_cuda.cpp index 9dac4390a1..06b09b460b 100644 --- a/examples/run/cuda/truth_finding_example_cuda.cpp +++ b/examples/run/cuda/truth_finding_example_cuda.cpp @@ -297,6 +297,13 @@ int seq_run(const traccc::finding_input_config& i_cfg, std::cout << "Track candidate matching Rate: " << float(n_matches) / track_candidates.size() << std::endl; + + // Compare the track parameters made on the host and on the device. + traccc::collection_comparator> + compare_fitter_infos{"fitted tracks"}; + compare_fitter_infos( + vecmem::get_data(track_states.get_headers()), + vecmem::get_data(track_states_cuda.get_headers())); } /// Statistics diff --git a/examples/run/cuda/truth_fitting_example_cuda.cpp b/examples/run/cuda/truth_fitting_example_cuda.cpp new file mode 100644 index 0000000000..51127409a2 --- /dev/null +++ b/examples/run/cuda/truth_fitting_example_cuda.cpp @@ -0,0 +1,266 @@ +/** TRACCC library, part of the ACTS project (R&D line) + * + * (c) 2023 CERN for the benefit of the ACTS project + * + * Mozilla Public License Version 2.0 + */ + +// Project include(s). +#include "traccc/cuda/fitting/fitting_algorithm.hpp" +#include "traccc/cuda/utils/stream.hpp" +#include "traccc/definitions/common.hpp" +#include "traccc/definitions/primitives.hpp" +#include "traccc/device/container_d2h_copy_alg.hpp" +#include "traccc/device/container_h2d_copy_alg.hpp" +#include "traccc/fitting/fitting_algorithm.hpp" +#include "traccc/fitting/kalman_filter/kalman_fitter.hpp" +#include "traccc/io/read_geometry.hpp" +#include "traccc/io/read_measurements.hpp" +#include "traccc/io/utils.hpp" +#include "traccc/options/common_options.hpp" +#include "traccc/options/handle_argument_errors.hpp" +#include "traccc/options/propagation_options.hpp" +#include "traccc/performance/collection_comparator.hpp" +#include "traccc/performance/container_comparator.hpp" +#include "traccc/performance/timer.hpp" +#include "traccc/resolution/fitting_performance_writer.hpp" +#include "traccc/utils/seed_generator.hpp" + +// detray include(s). +#include "detray/core/detector.hpp" +#include "detray/core/detector_metadata.hpp" +#include "detray/detectors/bfield.hpp" +#include "detray/io/common/detector_reader.hpp" +#include "detray/propagator/navigator.hpp" +#include "detray/propagator/propagator.hpp" +#include "detray/propagator/rk_stepper.hpp" + +// VecMem include(s). +#include +#include +#include +#include +#include + +// System include(s). +#include +#include +#include + +using namespace traccc; +namespace po = boost::program_options; + +// The main routine +// +int main(int argc, char* argv[]) { + // Set up the program options + po::options_description desc("Allowed options"); + + // Add options + desc.add_options()("help,h", "Give some help with the program's options"); + traccc::common_options common_opts(desc); + traccc::propagation_options propagation_opts(desc); + desc.add_options()("run_cpu", po::value()->default_value(false), + "run cpu tracking as well"); + + po::variables_map vm; + po::store(po::parse_command_line(argc, argv, desc), vm); + + // Check errors + traccc::handle_argument_errors(vm, desc); + + // Read options + common_opts.read(vm); + propagation_opts.read(vm); + auto run_cpu = vm["run_cpu"].as(); + + std::cout << "Running " << argv[0] << " " << common_opts.input_directory + << " " << common_opts.events << std::endl; + + /// Type declarations + using host_detector_type = detray::detector; + using device_detector_type = + detray::detector; + + using b_field_t = covfie::field; + using rk_stepper_type = + detray::rk_stepper>; + using host_navigator_type = detray::navigator; + using host_fitter_type = + traccc::kalman_fitter; + using device_navigator_type = detray::navigator; + using device_fitter_type = + traccc::kalman_fitter; + + // Memory resources used by the application. + vecmem::host_memory_resource host_mr; + vecmem::cuda::host_memory_resource cuda_host_mr; + vecmem::cuda::managed_memory_resource mng_mr; + vecmem::cuda::device_memory_resource device_mr; + traccc::memory_resource mr{device_mr, &cuda_host_mr}; + + // Performance writer + traccc::fitting_performance_writer fit_performance_writer( + traccc::fitting_performance_writer::config{}); + + // Output Stats + std::size_t n_fitted_tracks = 0; + std::size_t n_fitted_tracks_cuda = 0; + + /***************************** + * Build a geometry + *****************************/ + + // B field value and its type + // @TODO: Set B field as argument + const traccc::vector3 B{0, 0, 2 * detray::unit::T}; + auto field = detray::bfield::create_const_field(B); + + // Read the detector + detray::io::detector_reader_config reader_cfg{}; + reader_cfg + .add_file(traccc::io::data_directory() + common_opts.detector_file) + .add_file(traccc::io::data_directory() + common_opts.material_file) + .add_file(traccc::io::data_directory() + common_opts.grid_file); + + auto [host_det, names] = + detray::io::read_detector(mng_mr, reader_cfg); + + // Detector view object + auto det_view = detray::get_data(host_det); + + /***************************** + * Do the reconstruction + *****************************/ + + // Stream object + traccc::cuda::stream stream; + + // Copy object + vecmem::cuda::async_copy async_copy{stream.cudaStream()}; + + traccc::device::container_h2d_copy_alg< + traccc::track_candidate_container_types> + track_candidate_h2d{mr, async_copy}; + + traccc::device::container_d2h_copy_alg + track_state_d2h{mr, async_copy}; + + /// Standard deviations for seed track parameters + static constexpr std::array stddevs = { + 0.03 * detray::unit::mm, + 0.03 * detray::unit::mm, + 0.017, + 0.017, + 0.01 / detray::unit::GeV, + 1 * detray::unit::ns}; + + // Fitting algorithm object + typename traccc::fitting_algorithm::config_type fit_cfg; + fit_cfg.step_constraint = propagation_opts.step_constraint; + traccc::fitting_algorithm host_fitting(fit_cfg); + traccc::cuda::fitting_algorithm device_fitting( + fit_cfg, mr, async_copy, stream); + + // Seed generator + traccc::seed_generator sg(host_det, stddevs); + + traccc::performance::timing_info elapsedTimes; + + // Iterate over events + for (unsigned int event = common_opts.skip; + event < common_opts.events + common_opts.skip; ++event) { + + // Truth Track Candidates + traccc::event_map2 evt_map2(event, common_opts.input_directory, + common_opts.input_directory, + common_opts.input_directory); + + traccc::track_candidate_container_types::host truth_track_candidates = + evt_map2.generate_truth_candidates(sg, host_mr); + + // track candidates buffer + const traccc::track_candidate_container_types::buffer + truth_track_candidates_cuda_buffer = + track_candidate_h2d(traccc::get_data(truth_track_candidates)); + + // Navigation buffer + auto navigation_buffer = detray::create_candidates_buffer( + host_det, truth_track_candidates.size(), mr.main, mr.host); + + // Instantiate cuda containers/collections + traccc::track_state_container_types::buffer track_states_cuda_buffer{ + {{}, *(mr.host)}, {{}, *(mr.host), mr.host}}; + + { + traccc::performance::timer t("Track fitting (cuda)", elapsedTimes); + + // Run fitting + track_states_cuda_buffer = + device_fitting(det_view, field, navigation_buffer, + truth_track_candidates_cuda_buffer); + } + + traccc::track_state_container_types::host track_states_cuda = + track_state_d2h(track_states_cuda_buffer); + + // CPU container(s) + traccc::fitting_algorithm::output_type track_states; + + if (run_cpu) { + + { + traccc::performance::timer t("Track fitting (cpu)", + elapsedTimes); + + // Run fitting + track_states = + host_fitting(host_det, field, truth_track_candidates); + } + } + + if (run_cpu) { + // Show which event we are currently presenting the results for. + std::cout << "===>>> Event " << event << " <<<===" << std::endl; + + // Compare the track parameters made on the host and on the device. + traccc::collection_comparator> + compare_fitter_infos{"fitted tracks"}; + compare_fitter_infos( + vecmem::get_data(track_states.get_headers()), + vecmem::get_data(track_states_cuda.get_headers())); + } + + // Statistics + n_fitted_tracks += track_states.size(); + n_fitted_tracks_cuda += track_states_cuda.size(); + + if (common_opts.check_performance) { + for (unsigned int i = 0; i < track_states_cuda.size(); i++) { + const auto& trk_states_per_track = + track_states_cuda.at(i).items; + + const auto& fit_info = track_states_cuda[i].header; + + fit_performance_writer.write(trk_states_per_track, fit_info, + host_det, evt_map2); + } + } + } + + if (common_opts.check_performance) { + fit_performance_writer.finalize(); + } + + std::cout << "==> Statistics ... " << std::endl; + std::cout << "- created (cuda) " << n_fitted_tracks_cuda << " fitted tracks" + << std::endl; + std::cout << "- created (cpu) " << n_fitted_tracks << " fitted tracks" + << std::endl; + std::cout << "==>Elapsed times...\n" << elapsedTimes << std::endl; + + return 1; +} diff --git a/examples/simulation/simulate_telescope.cpp b/examples/simulation/simulate_telescope.cpp index 719f98e205..614d36051c 100644 --- a/examples/simulation/simulate_telescope.cpp +++ b/examples/simulation/simulate_telescope.cpp @@ -13,6 +13,8 @@ #include "traccc/options/handle_argument_errors.hpp" #include "traccc/options/options.hpp" #include "traccc/options/particle_gen_options.hpp" +#include "traccc/options/propagation_options.hpp" +#include "traccc/options/telescope_detector_options.hpp" #include "traccc/simulation/measurement_smearer.hpp" #include "traccc/simulation/simulator.hpp" #include "traccc/simulation/smearing_writer.hpp" @@ -34,7 +36,9 @@ using namespace traccc; namespace po = boost::program_options; int simulate(std::string output_directory, unsigned int events, - const traccc::particle_gen_options& pg_opts) { + const traccc::particle_gen_options& pg_opts, + const traccc::propagation_options& propagation_opts, + const traccc::telescope_detector_options& telescope_opts) { // Use deterministic random number generator for testing using uniform_gen_t = @@ -48,25 +52,33 @@ int simulate(std::string output_directory, unsigned int events, * Build a telescope geometry *****************************/ - // Plane alignment direction (aligned to x-axis) - detray::detail::ray traj{{0, 0, 0}, 0, {1, 0, 0}, -1}; + // Plane alignment direction (aligned to z-axis) + detray::detail::ray traj{{0, 0, 0}, 0, {0, 0, 1}, -1}; // Position of planes (in mm unit) - std::vector plane_positions = {20., 40., 60., 80., 100., - 120., 140, 160, 180.}; + std::vector plane_positions; + + for (unsigned int i = 0; i < telescope_opts.n_planes; i++) { + plane_positions.push_back((i + 1) * telescope_opts.spacing); + } // B field value and its type using b_field_t = covfie::field; - const vector3 B{2 * detray::unit::T, 0, 0}; + const vector3 B{0, 0, 2 * detray::unit::T}; auto field = detray::bfield::create_const_field(B); - // Create the detector - const auto mat = detray::silicon_tml(); - const scalar thickness = 0.5 * detray::unit::mm; + // Set material and thickness + detray::material mat; + if (!telescope_opts.empty_material) { + mat = detray::silicon(); + } else { + mat = detray::vacuum(); + } + + const scalar thickness = telescope_opts.thickness; // Use rectangle surfaces detray::mask> rectangle{ - 0u, 10000.f * detray::unit::mm, - 10000.f * detray::unit::mm}; + 0u, telescope_opts.half_length, telescope_opts.half_length}; detray::tel_det_config<> tel_cfg{rectangle}; tel_cfg.positions(plane_positions); @@ -95,7 +107,7 @@ int simulate(std::string output_directory, unsigned int events, // Smearing value for measurements traccc::measurement_smearer meas_smearer( - 50 * detray::unit::um, 50 * detray::unit::um); + telescope_opts.smearing, telescope_opts.smearing); // Type declarations using detector_type = decltype(det); @@ -115,6 +127,9 @@ int simulate(std::string output_directory, unsigned int events, writer_type>( events, det, field, std::move(generator), std::move(smearer_writer_cfg), full_path); + sim.get_config().step_constraint = propagation_opts.step_constraint; + sim.get_config().overstep_tolerance = propagation_opts.overstep_tolerance; + sim.run(); // Create detector file @@ -139,6 +154,8 @@ int main(int argc, char* argv[]) { desc.add_options()("events", po::value()->required(), "number of events"); traccc::particle_gen_options pg_opts(desc); + traccc::propagation_options propagation_opts(desc); + traccc::telescope_detector_options telescope_opts(desc); po::variables_map vm; po::store(po::parse_command_line(argc, argv, desc), vm); @@ -150,9 +167,12 @@ int main(int argc, char* argv[]) { auto output_directory = vm["output_directory"].as(); auto events = vm["events"].as(); pg_opts.read(vm); + propagation_opts.read(vm); + telescope_opts.read(vm); std::cout << "Running " << argv[0] << " " << output_directory << " " << events << std::endl; - return simulate(output_directory, events, pg_opts); + return simulate(output_directory, events, pg_opts, propagation_opts, + telescope_opts); } diff --git a/extern/detray/CMakeLists.txt b/extern/detray/CMakeLists.txt index 942401997d..af094a3fbb 100644 --- a/extern/detray/CMakeLists.txt +++ b/extern/detray/CMakeLists.txt @@ -18,7 +18,7 @@ message( STATUS "Building Detray as part of the TRACCC project" ) # Declare where to get Detray from. set( TRACCC_DETRAY_SOURCE -"URL;https://github.com/acts-project/detray/archive/refs/tags/v0.46.0.tar.gz;URL_MD5;db85b907f5906972a0b8025ef135a830" +"URL;https://github.com/acts-project/detray/archive/refs/tags/v0.47.0.tar.gz;URL_MD5;200f7fb00ac1d16b22e7330da62d03fb" CACHE STRING "Source for Detray, when built as part of this project" ) mark_as_advanced( TRACCC_DETRAY_SOURCE ) diff --git a/performance/include/traccc/performance/details/is_same_object.hpp b/performance/include/traccc/performance/details/is_same_object.hpp index b3cd4a2dc4..dd81cf9609 100644 --- a/performance/include/traccc/performance/details/is_same_object.hpp +++ b/performance/include/traccc/performance/details/is_same_object.hpp @@ -45,6 +45,7 @@ class is_same_object { #include "traccc/performance/impl/is_same_object.ipp" // Include specialisations for the core library types. +#include "traccc/performance/impl/is_same_fitter_info.ipp" #include "traccc/performance/impl/is_same_measurement.ipp" #include "traccc/performance/impl/is_same_seed.ipp" #include "traccc/performance/impl/is_same_spacepoint.ipp" diff --git a/performance/include/traccc/performance/impl/is_same_fitter_info.ipp b/performance/include/traccc/performance/impl/is_same_fitter_info.ipp new file mode 100644 index 0000000000..f26e7f2d97 --- /dev/null +++ b/performance/include/traccc/performance/impl/is_same_fitter_info.ipp @@ -0,0 +1,35 @@ +/** TRACCC library, part of the ACTS project (R&D line) + * + * (c) 2023 CERN for the benefit of the ACTS project + * + * Mozilla Public License Version 2.0 + */ + +#pragma once + +// Project include(s). +#include "traccc/edm/track_state.hpp" + +namespace traccc::details { + +/// @c traccc::is_same_object specialisation for @c traccc::fitter_info +template <> +class is_same_object> { + + public: + /// Constructor with a reference object, and an allowed uncertainty + is_same_object(const fitter_info& ref, + scalar unc = float_epsilon); + + /// Specialised implementation for @c traccc::measurement + bool operator()(const fitter_info& obj) const; + + private: + /// The reference object + std::reference_wrapper> m_ref; + /// The uncertainty + scalar m_unc; + +}; // class is_same_object> + +} // namespace traccc::details diff --git a/performance/include/traccc/resolution/fitting_performance_writer.hpp b/performance/include/traccc/resolution/fitting_performance_writer.hpp index 0fe901459c..8f59acde90 100644 --- a/performance/include/traccc/resolution/fitting_performance_writer.hpp +++ b/performance/include/traccc/resolution/fitting_performance_writer.hpp @@ -12,6 +12,7 @@ #include "traccc/resolution/stat_plot_tool_config.hpp" // Project include(s). +#include "traccc/edm/particle.hpp" #include "traccc/edm/track_parameters.hpp" #include "traccc/edm/track_state.hpp" #include "traccc/io/event_map2.hpp" @@ -94,8 +95,10 @@ class fitting_performance_writer { ptc.charge / getter::norm(global_mom); // For the moment, only fill with the first measurements - write_res(truth_param, trk_state.smoothed()); - write_stat(fit_info); + if (fit_info.ndf > 0 && !trk_state.is_hole) { + write_res(truth_param, trk_state.smoothed(), ptc); + } + write_stat(fit_info, track_states_per_track); } /// Writing caches into the file @@ -104,10 +107,12 @@ class fitting_performance_writer { private: /// Non-templated part of the @c write(...) function void write_res(const bound_track_parameters& truth_param, - const bound_track_parameters& fit_param); + const bound_track_parameters& fit_param, + const particle& ptc); /// Non-templated part of the @c write(...) function - void write_stat(const fitter_info& fit_info); + void write_stat(const fitter_info& fit_info, + const track_state_collection_types::host& track_states); /// Configuration for the tool config m_cfg; diff --git a/performance/include/traccc/resolution/res_plot_tool_config.hpp b/performance/include/traccc/resolution/res_plot_tool_config.hpp index e8dc8cbcf8..1d64cb65d0 100644 --- a/performance/include/traccc/resolution/res_plot_tool_config.hpp +++ b/performance/include/traccc/resolution/res_plot_tool_config.hpp @@ -21,8 +21,8 @@ namespace traccc { struct res_plot_tool_config { /// parameter sets to do plots - std::vector param_names = {"d0", "z0", "phi", - "theta", "qop", "t"}; + std::vector param_names = {"d0", "z0", "phi", "theta", + "qop", "t", "qopT"}; /// Binning setups std::map var_binning = { @@ -33,9 +33,27 @@ struct res_plot_tool_config { plot_helpers::binning("r_{#phi} [rad]", 100, -0.01, 0.01)}, {"residual_theta", plot_helpers::binning("r_{#theta} [rad]", 100, -0.01, 0.01)}, - {"residual_qop", - plot_helpers::binning("r_{q/p} [c/GeV]", 100, -0.1, 0.1)}, - {"residual_t", plot_helpers::binning("r_{t} [s]", 100, -1000, 1000)}}; + {"residual_qop", plot_helpers::binning("r_{q/p} [c/GeV]", 100, -1, 1)}, + {"residual_t", plot_helpers::binning("r_{t} [s]", 100, -1000, 1000)}, + {"residual_qopT", + plot_helpers::binning("r_{q/p_{T}} [c/GeV]", 100, -1, 1)}, + {"resolution_d0", + plot_helpers::binning("#sigma_{d0} [mm]", 100, 0, 0.5)}, + {"resolution_z0", + plot_helpers::binning("#sigma_{z0} [mm]", 100, 0, 0.5)}, + {"resolution_phi", + plot_helpers::binning("#sigma_{#phi} [rad]", 100, 0, 0.01)}, + {"resolution_theta", + plot_helpers::binning("#sigma_{#theta} [rad]", 100, 0, 0.01)}, + {"resolution_qop", + plot_helpers::binning("#sigma_{q/p} [c/GeV]", 100, 0, 0.1)}, + {"resolution_t", + plot_helpers::binning("#sigma_{t} [s]", 100, -1000, 1000)}, + {"resolution_qopT", + plot_helpers::binning("#sigma_{q/p_{T}} [c/GeV]", 100, 0, 0.1)}, + {"Eta", plot_helpers::binning("#eta", 40, -4, 4)}, + {"Pt", plot_helpers::binning("p_{T} [GeV/c]", 40, 0, 100)}, + }; }; } // namespace traccc diff --git a/performance/include/traccc/resolution/stat_plot_tool_config.hpp b/performance/include/traccc/resolution/stat_plot_tool_config.hpp index b724274e02..a1a22bd4d2 100644 --- a/performance/include/traccc/resolution/stat_plot_tool_config.hpp +++ b/performance/include/traccc/resolution/stat_plot_tool_config.hpp @@ -24,7 +24,8 @@ struct stat_plot_tool_config { {"ndf", plot_helpers::binning("ndf", 35, -5., 30.)}, {"chi2", plot_helpers::binning("chi2", 100, 0., 50.)}, {"reduced_chi2", plot_helpers::binning("chi2/ndf", 100, 0., 10.)}, - {"pval", plot_helpers::binning("pval", 100, 0., 1.)}}; + {"pval", plot_helpers::binning("pval", 100, 0., 1.)}, + {"chi2_local", plot_helpers::binning("chi2", 100, 0., 10.)}}; }; } // namespace traccc diff --git a/performance/src/performance/details/is_same_object.cpp b/performance/src/performance/details/is_same_object.cpp index e45c04afce..dc81cf581a 100644 --- a/performance/src/performance/details/is_same_object.cpp +++ b/performance/src/performance/details/is_same_object.cpp @@ -128,4 +128,23 @@ bool is_same_object::operator()( /// @} +/// @name Implementation for +/// @c traccc::details::is_same_object +/// @{ + +is_same_object>::is_same_object( + const fitter_info& ref, scalar unc) + : m_ref(ref), m_unc(unc) {} + +bool is_same_object>::operator()( + const fitter_info& obj) const { + + return (is_same_object(m_ref.get().fit_params, + m_unc)(obj.fit_params) && + is_same_angle(obj.ndf, m_ref.get().ndf, m_unc) && + is_same_angle(obj.chi2, m_ref.get().chi2, m_unc)); +} + +/// @} + } // namespace traccc::details diff --git a/performance/src/resolution/fitting_performance_writer.cpp b/performance/src/resolution/fitting_performance_writer.cpp index 5fff2d9d4b..a2ce9eec56 100644 --- a/performance/src/resolution/fitting_performance_writer.cpp +++ b/performance/src/resolution/fitting_performance_writer.cpp @@ -74,16 +74,21 @@ void fitting_performance_writer::finalize() { void fitting_performance_writer::write_res( const bound_track_parameters& truth_param, - const bound_track_parameters& fit_param) { + const bound_track_parameters& fit_param, const particle& ptc) { m_data->m_res_plot_tool.fill(m_data->m_res_plot_cache, truth_param, - fit_param); + fit_param, ptc); } void fitting_performance_writer::write_stat( - const fitter_info& fit_info) { + const fitter_info& fit_info, + const track_state_collection_types::host& track_states) { m_data->m_stat_plot_tool.fill(m_data->m_stat_plot_cache, fit_info); + + for (const auto& trk_state : track_states) { + m_data->m_stat_plot_tool.fill(m_data->m_stat_plot_cache, trk_state); + } } } // namespace traccc diff --git a/performance/src/resolution/res_plot_tool.cpp b/performance/src/resolution/res_plot_tool.cpp index 10ec4cc4ec..edb4241c5a 100644 --- a/performance/src/resolution/res_plot_tool.cpp +++ b/performance/src/resolution/res_plot_tool.cpp @@ -8,6 +8,15 @@ // Library include(s). #include "res_plot_tool.hpp" +// Detray include(s). +#include "detray/utils/statistics.hpp" + +// ROOT include(s). +#ifdef TRACCC_HAVE_ROOT +#include +#include +#endif // TRACCC_HAVE_ROOT + namespace traccc { res_plot_tool::res_plot_tool(const res_plot_tool_config& cfg) : m_cfg(cfg) {} @@ -15,12 +24,16 @@ res_plot_tool::res_plot_tool(const res_plot_tool_config& cfg) : m_cfg(cfg) {} void res_plot_tool::book(res_plot_cache& cache) const { plot_helpers::binning b_pull = m_cfg.var_binning.at("pull"); + plot_helpers::binning b_eta = m_cfg.var_binning.at("Eta"); + plot_helpers::binning b_pT = m_cfg.var_binning.at("Pt"); // Avoid unused variable warnings when building the code without ROOT. (void)cache; (void)b_pull; + (void)b_eta; + (void)b_pT; - for (std::size_t idx = 0; idx < e_bound_size; idx++) { + for (std::size_t idx = 0; idx < m_cfg.param_names.size(); idx++) { std::string par_name = m_cfg.param_names.at(idx); // Binning for residual is parameter dependent @@ -40,38 +53,101 @@ void res_plot_tool::book(res_plot_cache& cache) const { cache.pulls[par_name] = plot_helpers::book_histo( Form("pull_%s", par_name.c_str()), Form("Pull of %s", par_name.c_str()), b_pull); + + // resolution vs. eta plots + cache.resolutions_eta[par_name] = plot_helpers::book_histo( + Form("resolution_%s_vs_eta", par_name.c_str()), + Form("Resolution of %s vs. eta", par_name.c_str()), b_eta); + + // resolution vs. pT plots + cache.resolutions_pT[par_name] = plot_helpers::book_histo( + Form("resolution_%s_vs_pT", par_name.c_str()), + Form("Resolution of %s vs. pT", par_name.c_str()), b_pT); + + // residuals vs. eta plots + cache.residuals_eta[par_name] = plot_helpers::book_histo( + Form("residual_%s_vs_eta", par_name.c_str()), + Form("Residual of %s vs. eta", par_name.c_str()), b_eta, + b_residual); + + // residuals vs. pT plots + cache.residuals_pT[par_name] = plot_helpers::book_histo( + Form("residual_%s_vs_pT", par_name.c_str()), + Form("Residual of %s vs. pT", par_name.c_str()), b_pT, b_residual); + + for (int i = 0; i < b_eta.n_bins; i++) { + cache.residuals_per_eta[par_name][i] = plot_helpers::book_histo( + Form("residual_%s_at_%d_th_eta", par_name.c_str(), i), + Form("Residual of %s at %d th eta", par_name.c_str(), i), + b_residual); + } + + for (int i = 0; i < b_pT.n_bins; i++) { + cache.residuals_per_pT[par_name][i] = plot_helpers::book_histo( + Form("residual_%s_at_%d_th_pT", par_name.c_str(), i), + Form("Residual of %s at %d th pT", par_name.c_str(), i), + b_residual); + } #endif // TRACCC_HAVE_ROOT } } void res_plot_tool::fill(res_plot_cache& cache, const bound_track_parameters& truth_param, - const bound_track_parameters& fit_param) const { + const bound_track_parameters& fit_param, + const particle& ptc) const { + + // Find index of eta and pT for resolution histogram + const scalar eta = getter::eta(ptc.mom); + const scalar pT = std::hypot(ptc.mom[0], ptc.mom[1]); // Avoid unused variable warnings when building the code without ROOT. (void)cache; - for (std::size_t idx = 0; idx < e_bound_size; idx++) { + for (std::size_t idx = 0; idx < m_cfg.param_names.size(); idx++) { std::string par_name = m_cfg.param_names.at(idx); - const auto residual = getter::element(fit_param.vector(), idx, 0) - - getter::element(truth_param.vector(), idx, 0); + scalar residual = 0.f; + scalar pull = 0.f; + + if (idx < e_bound_size) { + residual = getter::element(fit_param.vector(), idx, 0) - + getter::element(truth_param.vector(), idx, 0); + pull = residual / + std::sqrt(getter::element(fit_param.covariance(), idx, idx)); + } else if (par_name == "qopT") { + residual = fit_param.qopT() - truth_param.qopT(); - const auto pull = residual / std::sqrt(getter::element( - fit_param.covariance(), idx, idx)); + const scalar sigma = + std::sqrt(getter::element(fit_param.covariance(), + e_bound_qoverp, e_bound_qoverp)) / + std::sin(fit_param.theta()); + pull = residual / sigma; + } // Avoid unused variable warnings when building the code without ROOT. (void)residual; (void)pull; + const auto eta_idx = + std::min(cache.resolutions_eta[par_name]->FindBin(eta) - 1, + cache.resolutions_eta[par_name]->GetNbinsX() - 1); + const auto pT_idx = + std::min(cache.resolutions_pT[par_name]->FindBin(pT) - 1, + cache.resolutions_pT[par_name]->GetNbinsX() - 1); + #ifdef TRACCC_HAVE_ROOT cache.residuals.at(par_name)->Fill(residual); cache.pulls.at(par_name)->Fill(pull); + cache.residuals_eta.at(par_name)->Fill(eta, residual); + cache.residuals_pT.at(par_name)->Fill(pT, residual); + cache.residuals_per_eta.at(par_name).at(eta_idx)->Fill(residual); + cache.residuals_per_pT.at(par_name).at(pT_idx)->Fill(residual); #endif // TRACCC_HAVE_ROOT } } -void res_plot_tool::write(const res_plot_cache& cache) const { +void res_plot_tool::write(res_plot_cache& cache) const { // Avoid unused variable warnings when building the code without ROOT. (void)cache; @@ -83,6 +159,83 @@ void res_plot_tool::write(const res_plot_cache& cache) const { for (const auto& pull : cache.pulls) { pull.second->Write(); } + for (const auto& residual : cache.residuals_eta) { + residual.second->Write(); + } + for (const auto& residual : cache.residuals_pT) { + residual.second->Write(); + } + + auto create_resolution_graph = [](std::shared_ptr H, auto residuals, + const char* x_axis_title, + const char* y_axis_title) { + const auto n_bins = H->GetNbinsX(); + + std::vector sigmas; + + for (int i = 0; i < n_bins; i++) { + auto& data = residuals[i]; + + // When there is no entry + if (data->GetEntries() == 0) { + H->SetBinContent(i, 0); + sigmas.push_back(0); + continue; + } + + // Function used for the fit. + TF1 gaus{"gaus", "gaus", -1., 1.}; + double fit_par[3]; + auto res = data->Fit("gaus", "Q0S"); + gaus.GetParameters(&fit_par[0]); + H->SetBinContent(i + 1, fit_par[2]); + sigmas.push_back(gaus.GetParError(2)); + } + + std::unique_ptr G = + std::make_unique(H.get()); + + for (int i = 0; i < n_bins; i++) { + G->SetPointError(i, H->GetBinWidth(i) / 2.f, sigmas[i]); + } + G->SetName(H->GetName()); + G->SetMinimum(1e-5); + G->SetMaximum(1e-1); + G->GetHistogram()->GetYaxis()->SetMaxDigits(2); + G->GetHistogram()->GetXaxis()->SetTitle(x_axis_title); + G->GetHistogram()->GetYaxis()->SetTitle(y_axis_title); + + return G; + }; + + for (const auto& resolution : cache.resolutions_eta) { + const auto& par_name = resolution.first; + auto& H = resolution.second; + + plot_helpers::binning binning = + m_cfg.var_binning.at("resolution_" + par_name); + + auto G = + create_resolution_graph(H, cache.residuals_per_eta.at(par_name), + "#eta", binning.title.c_str()); + + G->Write(); + } + + for (const auto& resolution : cache.resolutions_pT) { + const auto& par_name = resolution.first; + auto& H = resolution.second; + + plot_helpers::binning binning = + m_cfg.var_binning.at("resolution_" + par_name); + + auto G = + create_resolution_graph(H, cache.residuals_per_pT.at(par_name), + "p_{T} [GeV/c]", binning.title.c_str()); + + G->Write(); + } + #endif // TRACCC_HAVE_ROOT } diff --git a/performance/src/resolution/res_plot_tool.hpp b/performance/src/resolution/res_plot_tool.hpp index 3d1ab13fff..1564dc8829 100644 --- a/performance/src/resolution/res_plot_tool.hpp +++ b/performance/src/resolution/res_plot_tool.hpp @@ -12,6 +12,7 @@ #include "traccc/resolution/res_plot_tool_config.hpp" // Project include(s). +#include "traccc/edm/particle.hpp" #include "traccc/edm/track_parameters.hpp" // System include(s). @@ -28,8 +29,16 @@ class res_plot_tool { struct res_plot_cache { #ifdef TRACCC_HAVE_ROOT // Residuals and pulls for parameters - std::map > residuals; - std::map > pulls; + std::map> residuals; + std::map> pulls; + std::map> resolutions_eta; + std::map> resolutions_pT; + std::map> residuals_eta; + std::map> residuals_pT; + std::map>> + residuals_per_eta; + std::map>> + residuals_per_pT; #endif // TRACCC_HAVE_ROOT }; @@ -49,12 +58,13 @@ class res_plot_tool { /// @param truth_param truth track parameter /// @param fit_param fitted track parameter void fill(res_plot_cache& cache, const bound_track_parameters& truth_param, - const bound_track_parameters& fit_param) const; + const bound_track_parameters& fit_param, + const particle& ptc) const; /// @brief write the resolution plots into ROOT /// /// @param cache the cache for resolution plots - void write(const res_plot_cache& cache) const; + void write(res_plot_cache& cache) const; private: res_plot_tool_config m_cfg; ///< The Config class diff --git a/performance/src/resolution/stat_plot_tool.cpp b/performance/src/resolution/stat_plot_tool.cpp index 7c2cb33cbe..c22b6b2bb3 100644 --- a/performance/src/resolution/stat_plot_tool.cpp +++ b/performance/src/resolution/stat_plot_tool.cpp @@ -27,11 +27,26 @@ void stat_plot_tool::book(stat_plot_cache& cache) const { plot_helpers::binning b_chi2 = m_cfg.var_binning.at("chi2"); plot_helpers::binning b_reduced_chi2 = m_cfg.var_binning.at("reduced_chi2"); plot_helpers::binning b_pval = m_cfg.var_binning.at("pval"); + plot_helpers::binning b_chi2_local = m_cfg.var_binning.at("chi2_local"); cache.ndf_hist = plot_helpers::book_histo("ndf", "NDF", b_ndf); cache.chi2_hist = plot_helpers::book_histo("chi2", "Chi2", b_chi2); cache.reduced_chi2_hist = plot_helpers::book_histo("reduced_chi2", "Chi2/NDF", b_reduced_chi2); cache.pval_hist = plot_helpers::book_histo("pval", "p value", b_pval); + for (unsigned int D = 1u; D <= 2u; D++) { + cache.chi2_filtered_hist[D] = plot_helpers::book_histo( + Form("chi2_%dD_filtered", D), + Form("chi2 of %dD filtered parameters", D), b_chi2_local); + cache.chi2_smoothed_hist[D] = plot_helpers::book_histo( + Form("chi2_%dD_smoothed", D), + Form("chi2 of %dD smoothed parameters", D), b_chi2_local); + cache.pval_filtered_hist[D] = plot_helpers::book_histo( + Form("pval_%dD_filtered", D), + Form("p value of %dD filtered parameters", D), b_pval); + cache.pval_smoothed_hist[D] = plot_helpers::book_histo( + Form("pval_%dD_smoothed", D), + Form("p value of %dD smoothed parameters", D), b_pval); + } #endif // TRACCC_HAVE_ROOT } @@ -52,6 +67,26 @@ void stat_plot_tool::fill(stat_plot_cache& cache, #endif // TRACCC_HAVE_ROOT } +void stat_plot_tool::fill(stat_plot_cache& cache, + const track_state& trk_state) const { + + // Avoid unused variable warnings when building the code without ROOT. + (void)cache; + (void)trk_state; + +#ifdef TRACCC_HAVE_ROOT + const unsigned int D = trk_state.get_measurement().meas_dim; + const auto filtered_chi2 = trk_state.filtered_chi2(); + const auto smoothed_chi2 = trk_state.smoothed_chi2(); + cache.chi2_filtered_hist[D]->Fill(filtered_chi2); + cache.chi2_smoothed_hist[D]->Fill(smoothed_chi2); + cache.pval_filtered_hist[D]->Fill( + ROOT::Math::chisquared_cdf_c(filtered_chi2, D)); + cache.pval_smoothed_hist[D]->Fill( + ROOT::Math::chisquared_cdf_c(smoothed_chi2, D)); +#endif // TRACCC_HAVE_ROOT +} + void stat_plot_tool::write(const stat_plot_cache& cache) const { // Avoid unused variable warnings when building the code without ROOT. @@ -62,6 +97,24 @@ void stat_plot_tool::write(const stat_plot_cache& cache) const { cache.chi2_hist->Write(); cache.reduced_chi2_hist->Write(); cache.pval_hist->Write(); + cache.ndf_hist->Write(); + cache.chi2_hist->Write(); + cache.reduced_chi2_hist->Write(); + cache.pval_hist->Write(); + + for (const auto& flt : cache.chi2_filtered_hist) { + flt.second->Write(); + } + for (const auto& smt : cache.chi2_smoothed_hist) { + smt.second->Write(); + } + for (const auto& flt : cache.pval_filtered_hist) { + flt.second->Write(); + } + for (const auto& smt : cache.pval_smoothed_hist) { + smt.second->Write(); + } + #endif // TRACCC_HAVE_ROOT } diff --git a/performance/src/resolution/stat_plot_tool.hpp b/performance/src/resolution/stat_plot_tool.hpp index e380d0c25b..ad53e8f439 100644 --- a/performance/src/resolution/stat_plot_tool.hpp +++ b/performance/src/resolution/stat_plot_tool.hpp @@ -30,6 +30,14 @@ class stat_plot_tool { std::unique_ptr reduced_chi2_hist; // Histogram for the pvalue std::unique_ptr pval_hist; + // Histogram for chi2 of filtered states + std::map> chi2_filtered_hist; + // Histogram for chi2 of smoothed states + std::map> chi2_smoothed_hist; + // Histogram for p-value of filtered states + std::map> pval_filtered_hist; + // Histogram for p-value of smoothed states + std::map> pval_smoothed_hist; #endif // TRACCC_HAVE_ROOT }; @@ -50,6 +58,13 @@ class stat_plot_tool { void fill(stat_plot_cache& cache, const fitter_info& fit_info) const; + /// @brief fill the cache + /// + /// @param cache the cache for statistics plots + /// @param trk_state track state at local measurements + void fill(stat_plot_cache& cache, + const track_state& trk_state) const; + /// @brief write the statistics plots into ROOT /// /// @param cache the cache for statistics plots diff --git a/performance/src/utils/helpers.cpp b/performance/src/utils/helpers.cpp index c7247cf789..6389e6ab3d 100644 --- a/performance/src/utils/helpers.cpp +++ b/performance/src/utils/helpers.cpp @@ -24,6 +24,22 @@ std::unique_ptr book_histo(std::string_view hist_name, return result; } +std::unique_ptr book_histo(std::string_view hist_name, + std::string_view hist_title, + const binning& var_x_binning, + const binning& var_y_binning) { + + auto result = std::make_unique( + hist_name.data(), hist_title.data(), var_x_binning.n_bins, + var_x_binning.min, var_x_binning.max, var_y_binning.n_bins, + var_y_binning.min, var_y_binning.max); + + result->GetXaxis()->SetTitle(var_x_binning.title.c_str()); + result->GetYaxis()->SetTitle(var_y_binning.title.c_str()); + result->Sumw2(); + return result; +} + std::unique_ptr book_eff(std::string_view eff_name, std::string_view eff_title, const binning& var_binning) { diff --git a/performance/src/utils/helpers.hpp b/performance/src/utils/helpers.hpp index ba5c69cd1b..3d2955dfd8 100644 --- a/performance/src/utils/helpers.hpp +++ b/performance/src/utils/helpers.hpp @@ -14,6 +14,7 @@ #ifdef TRACCC_HAVE_ROOT #include #include +#include #include #endif // TRACCC_HAVE_ROOT @@ -35,6 +36,18 @@ std::unique_ptr book_histo(std::string_view hist_name, std::string_view hist_title, const binning& var_binning); +/// @brief book a 1D histogram +/// @param hist_name the name of histogram +/// @param hist_title the title of histogram +/// @param var_x_binning the binning info of x variable +/// @param var_y_binning the binning info of y variable +/// @return histogram pointer +/// +std::unique_ptr book_histo(std::string_view hist_name, + std::string_view hist_title, + const binning& var_x_binning, + const binning& var_y_binning); + /// @brief book a 1D efficiency plot /// @param effName the name of plot /// @param effTitle the title of plot diff --git a/tests/cpu/test_simulation.cpp b/tests/cpu/test_simulation.cpp index a6f14fa0a9..368d05284c 100644 --- a/tests/cpu/test_simulation.cpp +++ b/tests/cpu/test_simulation.cpp @@ -121,7 +121,7 @@ GTEST_TEST(detray_simulation, toy_geometry_simulation) { std::vector particles; auto particle_reader = traccc::io::csv::make_particle_reader( - detail::get_event_filename(i_event, "-particles.csv")); + traccc::io::get_event_filename(i_event, "-particles.csv")); traccc::io::csv::particle io_particle; while (particle_reader.read(io_particle)) { particles.push_back(io_particle); @@ -129,7 +129,7 @@ GTEST_TEST(detray_simulation, toy_geometry_simulation) { std::vector hits; auto hit_reader = traccc::io::csv::make_hit_reader( - detail::get_event_filename(i_event, "-hits.csv")); + traccc::io::get_event_filename(i_event, "-hits.csv")); traccc::io::csv::hit io_hit; while (hit_reader.read(io_hit)) { hits.push_back(io_hit); @@ -137,7 +137,7 @@ GTEST_TEST(detray_simulation, toy_geometry_simulation) { std::vector measurements; auto measurement_reader = traccc::io::csv::make_measurement_reader( - detail::get_event_filename(i_event, "-measurements.csv")); + traccc::io::get_event_filename(i_event, "-measurements.csv")); traccc::io::csv::measurement io_measurement; while (measurement_reader.read(io_measurement)) { measurements.push_back(io_measurement); @@ -146,8 +146,8 @@ GTEST_TEST(detray_simulation, toy_geometry_simulation) { std::vector meas_hit_ids; auto measurement_hit_id_reader = traccc::io::csv::make_measurement_hit_id_reader( - detail::get_event_filename(i_event, - "-measurement-simhit-map.csv")); + traccc::io::get_event_filename(i_event, + "-measurement-simhit-map.csv")); traccc::io::csv::measurement_hit_id io_meas_hit_id; while (measurement_hit_id_reader.read(io_meas_hit_id)) { meas_hit_ids.push_back(io_meas_hit_id); @@ -267,7 +267,7 @@ TEST_P(TelescopeDetectorSimulation, telescope_detector_simulation) { std::vector measurements; auto measurement_reader = traccc::io::csv::make_measurement_reader( directory + - detail::get_event_filename(i_event, "-measurements.csv")); + traccc::io::get_event_filename(i_event, "-measurements.csv")); traccc::io::csv::measurement io_measurement; while (measurement_reader.read(io_measurement)) { measurements.push_back(io_measurement);