From 77fffe17e87e4237d9aaa106d9c343fc7f805947 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Mon, 16 Sep 2024 15:48:50 -0600 Subject: [PATCH] Update tracking to use Acts v36.0 (#1454) * move acts to v36 and required ldmx::tracking code changes, ldmx-sw builds but does not link yet * remove now non-existent acts component from linking The ActsPluginIdentification does not exist within ACTS v36 and so asking for it to be built just gets ignored since its not an option ACTS looks for; however, asking to link to it errors since this component is not a valid library that can be linked to (it is never built). We can remove reference to this plugin and point out that the default-build of acts appears to satisfy our needs. * v36 is building, linking, running and give resonable results. Still need to update GSF track fitting and restore a few features. * forgot a couple of changes to MeasurementCalibrator * fix up CKF so that it makes everything for ldmx::track; GSF builds but not sure if working yet * Once in a while CKFProcessor crashes when it cant extrapolate to target state. Put in a catch and some info statements * Add back a block I accidently deleted * adapted GSF for v36 and added tagger/recoil flexibility * move acts to v36 and required ldmx::tracking code changes, ldmx-sw builds but does not link yet * remove now non-existent acts component from linking The ActsPluginIdentification does not exist within ACTS v36 and so asking for it to be built just gets ignored since its not an option ACTS looks for; however, asking to link to it errors since this component is not a valid library that can be linked to (it is never built). We can remove reference to this plugin and point out that the default-build of acts appears to satisfy our needs. * v36 is building, linking, running and give resonable results. Still need to update GSF track fitting and restore a few features. * forgot a couple of changes to MeasurementCalibrator * fix up CKF so that it makes everything for ldmx::track; GSF builds but not sure if working yet * Once in a while CKFProcessor crashes when it cant extrapolate to target state. Put in a catch and some info statements * Add back a block I accidently deleted * adapted GSF for v36 and added tagger/recoil flexibility * add new plots and ranges to dqm * update the example reco so that it actually works right * clang_format and some other cleanup * clang-format checker didnt like my comment for some reason * damn formatting * Update Tracking/include/Tracking/Reco/DigitizationProcessor.h Co-authored-by: Tamas Vami * Update Tracking/include/Tracking/Reco/GSFProcessor.h Co-authored-by: Tamas Vami * Update Tracking/include/Tracking/Sim/MeasurementCalibrator.h Co-authored-by: Tamas Vami * Update Tracking/python/tracking.py Co-authored-by: Tamas Vami * Update Tracking/include/Tracking/Reco/TruthSeedProcessor.h Co-authored-by: Tamas Vami * Update Tracking/src/Tracking/Reco/CKFProcessor.cxx Co-authored-by: Tamas Vami * Update Tracking/include/Tracking/Sim/MeasurementCalibrator.h Co-authored-by: Tamas Vami * Update Tracking/include/Tracking/Sim/MeasurementCalibrator.h Co-authored-by: Tamas Vami * Update Tracking/include/Tracking/Sim/MeasurementCalibrator.h Co-authored-by: Tamas Vami * Update Tracking/include/Tracking/Sim/MeasurementCalibrator.h Co-authored-by: Tamas Vami * Update Tracking/src/Tracking/Reco/CKFProcessor.cxx Co-authored-by: Tamas Vami * Update Tracking/src/Tracking/Reco/SeedFinderProcessor.cxx Co-authored-by: Tamas Vami * Update Tracking/include/Tracking/Reco/TrackExtrapolatorTool.h Co-authored-by: Tamas Vami * Update Tracking/include/Tracking/Reco/TrackExtrapolatorTool.h Co-authored-by: Tamas Vami * Update Tracking/include/Tracking/Reco/TrackExtrapolatorTool.h Co-authored-by: Tamas Vami * Update Tracking/include/Tracking/Reco/TrackExtrapolatorTool.h Co-authored-by: Tamas Vami * Update Tracking/include/Tracking/Reco/TrackExtrapolatorTool.h Co-authored-by: Tamas Vami * Update Tracking/src/Tracking/Reco/Vertexer.cxx Co-authored-by: Tamas Vami * Update Tracking/src/Tracking/Reco/VertexProcessor.cxx Co-authored-by: Tamas Vami * Update Tracking/src/Tracking/Reco/Vertexer.cxx Co-authored-by: Tamas Vami * Update Tracking/src/Tracking/Reco/CKFProcessor.cxx Co-authored-by: Tamas Vami * Update Tracking/src/Tracking/Reco/CKFProcessor.cxx Co-authored-by: Tamas Vami * Update Tracking/src/Tracking/Reco/CKFProcessor.cxx Co-authored-by: Tamas Vami * add typename to make it clearer * clean up some comments from tvami. Probably some more clang-format issues though * a few more changes for clang-format * add specific measurement_collection to pn validation config --------- Co-authored-by: tomeichlersmith Co-authored-by: EBerzin Co-authored-by: Tamas Vami --- .github/validation_samples/ecal_pn/config.py | 1 + Tracking/CMakeLists.txt | 13 +- Tracking/acts | 2 +- Tracking/exampleConfigs/reco.py | 162 +++--- Tracking/include/Tracking/Reco/CKFProcessor.h | 24 +- .../Tracking/Reco/CustomStatePropagator.h | 98 ---- .../Tracking/Reco/DigitizationProcessor.h | 3 - Tracking/include/Tracking/Reco/GSFProcessor.h | 9 +- .../Tracking/Reco/SeedFinderProcessor.h | 2 - .../Tracking/Reco/TrackExtrapolatorTool.h | 64 +-- .../Tracking/Reco/TruthSeedProcessor.h | 1 - .../include/Tracking/Sim/BFieldXYZUtils.h | 30 +- .../include/Tracking/Sim/LdmxSpacePoint.h | 6 +- .../Tracking/Sim/MeasurementCalibrator.h | 65 ++- .../Tracking/Sim/PropagatorStepWriter.h | 2 +- .../Tracking/Sim/SeedToTrackParamMaker.ipp | 4 +- Tracking/include/Tracking/Sim/TrackingUtils.h | 39 +- .../include/Tracking/dqm/TrackingRecoDQM.h | 4 + Tracking/python/dqm.py | 43 +- Tracking/python/geo.py | 3 +- Tracking/python/tracking.py | 3 + .../Tracking/Reco/AlignmentTestProcessor.cxx | 3 +- Tracking/src/Tracking/Reco/CKFProcessor.cxx | 484 +++++++++--------- .../Tracking/Reco/CustomStatePropagator.cxx | 271 ---------- .../Tracking/Reco/DigitizationProcessor.cxx | 40 +- Tracking/src/Tracking/Reco/GSFProcessor.cxx | 136 +++-- .../src/Tracking/Reco/SeedFinderProcessor.cxx | 53 +- .../src/Tracking/Reco/TruthSeedProcessor.cxx | 19 +- .../src/Tracking/Reco/VertexProcessor.cxx | 30 +- Tracking/src/Tracking/Reco/Vertexer.cxx | 22 +- .../src/Tracking/Sim/PropagatorStepWriter.cxx | 26 +- Tracking/src/Tracking/dqm/TrackingRecoDQM.cxx | 51 +- .../src/Tracking/geo/EcalTrackingGeometry.cxx | 3 +- .../Tracking/geo/TrackersTrackingGeometry.cxx | 4 +- .../src/Tracking/geo/TrackingGeometry.cxx | 3 +- 35 files changed, 756 insertions(+), 967 deletions(-) delete mode 100644 Tracking/include/Tracking/Reco/CustomStatePropagator.h delete mode 100644 Tracking/src/Tracking/Reco/CustomStatePropagator.cxx diff --git a/.github/validation_samples/ecal_pn/config.py b/.github/validation_samples/ecal_pn/config.py index 0c8d1c988..cc45376b1 100644 --- a/.github/validation_samples/ecal_pn/config.py +++ b/.github/validation_samples/ecal_pn/config.py @@ -176,6 +176,7 @@ seed_recoil_dqm.title = "" recoil_dqm = tkdqm.TrackingRecoDQM("RecoilTrackerDQM") +recoil_dqm.measurement_collection=digi_recoil.out_collection recoil_dqm.buildHistograms() recoil_dqm.track_collection = tracking_recoil.out_trk_collection recoil_dqm.truth_collection = "RecoilTruthTracks" diff --git a/Tracking/CMakeLists.txt b/Tracking/CMakeLists.txt index 6a9f427dd..11a45e730 100644 --- a/Tracking/CMakeLists.txt +++ b/Tracking/CMakeLists.txt @@ -51,15 +51,11 @@ if(BUILD_EVENT_ONLY) endif() # since ACTS changes frequently enough, we build it as a submodule of Tracking -# here we 'set' some CMake variables that we would have passed to 'cmake' on +# here we could 'set' some CMake variables that we would have passed to 'cmake' on # the command line with '-D' if we were building ACTS directly. -#set(ACTS_BUILD_PLUGIN_DD4HEP ON) # representing geometry with DD4hep right now -#set(ACTS_BUILD_PLUGIN_TGEO ON) # representing geometry with DD4hep right now -set(ACTS_BUILD_PLUGIN_IDENTIFICATION ON) #This turns on the identiers -set(ACTS_BUILD_EXAMPLES OFF) # don't waste time building examples -set(CMAKE_CXX_STANDARD 17) # set C++ standard we use within ldmx-sw -add_subdirectory(acts) # now build acts -#add_subdirectory(acts-dd4hep) +# for example, disable the building of examples (which is already disabled by default) +#set(ACTS_BUILD_EXAMPLES OFF) # don't waste time building examples +add_subdirectory(acts) # now build acts # luckily for us, adding ACTS as a subdirectory produces the same CMake targets # as find_package, so we don't need to do anything else from here on out @@ -79,7 +75,6 @@ setup_library(module Tracking dependencies Framework::Configure Framework::Framework ActsCore - ActsPluginIdentification Geant4::Interface ROOT::Physics Tracking::Event diff --git a/Tracking/acts b/Tracking/acts index ff2f97619..6eca77c45 160000 --- a/Tracking/acts +++ b/Tracking/acts @@ -1 +1 @@ -Subproject commit ff2f976198d08aa9a3d70133051f2097f62e587e +Subproject commit 6eca77c45b136861272694edbb61bb77200948a5 diff --git a/Tracking/exampleConfigs/reco.py b/Tracking/exampleConfigs/reco.py index 3a67b288b..8664b00bc 100644 --- a/Tracking/exampleConfigs/reco.py +++ b/Tracking/exampleConfigs/reco.py @@ -4,6 +4,8 @@ import os,math from LDMX.Framework import ldmxcfg +from LDMX.SimCore import generators +from LDMX.SimCore import simulator p = ldmxcfg.Process("TrackerReco") @@ -14,21 +16,33 @@ # From the conditions from LDMX.Tracking import geo +n_evts=1000 #number of events to gen/reco + +# set up a simple particle gun for this example # +# just 8gev electrons started upstream of tagger and first ts # +partGunString='single_8gev_e_upstream_tagger' +detector = 'ldmx-det-v14-8gev-no-cals' +#### set up beam simulation +sim = simulator.simulator('inclusive_single_8gev') +sim.setDetector(detector,include_scoring_planes = True) +sim.description = 'A single 8gev electron shot from upstream of the 8gev tagger.' +sim.beamSpotSmear = [20., 80., 0] +particle_gun = generators.single_8gev_e_upstream_tagger() +sim.generators.append(particle_gun) +#### end beam simulation # Truth seeder # Runs truth tracking producing tracks from target scoring plane hits for Recoil # and generated electros for Tagger. # Truth tracks can be used for assessing tracking performance or using as seeds -truth_tracking = tracking.TruthSeedProcessor() -truth_tracking.debug = True -truth_tracking.trk_coll_name = "RecoilTruthSeeds" -truth_tracking.pdgIDs = [11] -truth_tracking.scoring_hits = "TargetScoringPlaneHits" -truth_tracking.z_min = 0. -truth_tracking.track_id = -1 -truth_tracking.p_cut = 0.05 # In MeV -truth_tracking.pz_cut = 0.03 -truth_tracking.p_cutEcal = 0. # In MeV +truth_tracking = tracking.TruthSeedProcessor("TruthSeeds") +truth_tracking.debug = False +#truth_tracking.trk_coll_name = "RecoilTruthSeeds" +#truth_tracking.pdgIDs = [11] +#truth_tracking.scoring_hits = "TargetScoringPlaneHits" +#truth_tracking.z_min = 0. +#truth_tracking.track_id = -1 +#truth_tracking.p_cut = 0.05 # In MeV # These smearing quantities are default. We expect around 6um hit resolution in bending plane @@ -41,20 +55,20 @@ # Smearing Processor - Tagger # Runs G4 hit smearing producing measurements in the Tagger tracker. # Hits that belong to the same sensor with the same trackID are merged together to reduce combinatorics -digiTagger = tracking.DigitizationProcessor("DigitizationProcessor") -digiTagger.hit_collection = "TaggerSimHits" -digiTagger.out_collection = "DigiTaggerSimHits" -digiTagger.merge_hits = True -digiTagger.sigma_u = uSmearing -digiTagger.sigma_v = vSmearing +digi_tagger = tracking.DigitizationProcessor("DigitizationProcessor") +digi_tagger.hit_collection = "TaggerSimHits" +digi_tagger.out_collection = "DigiTaggerSimHits" +digi_tagger.merge_hits = True +digi_tagger.sigma_u = uSmearing +digi_tagger.sigma_v = vSmearing # Smearing Processor - Recoil -digiRecoil = tracking.DigitizationProcessor("DigitizationProcessorRecoil") -digiRecoil.hit_collection = "RecoilSimHits" -digiRecoil.out_collection = "DigiRecoilSimHits" -digiRecoil.merge_hits = True -digiRecoil.sigma_u = uSmearing -digiRecoil.sigma_v = vSmearing +digi_recoil = tracking.DigitizationProcessor("DigitizationProcessorRecoil") +digi_recoil.hit_collection = "RecoilSimHits" +digi_recoil.out_collection = "DigiRecoilSimHits" +digi_recoil.merge_hits = True +digi_recoil.sigma_u = uSmearing +digi_recoil.sigma_v = vSmearing # Seed Finder Tagger @@ -63,25 +77,26 @@ # parameters and the impact parameters at the target or generation point. For the tagger one should look # for compatibility with the beam orbit / beam spot -seederTagger = tracking.SeedFinderProcessor() -seederTagger.input_hits_collection = digiTagger.out_collection -seederTagger.out_seed_collection = "TaggerRecoSeeds" -seederTagger.pmin = 2. -seederTagger.pmax = 8. -seederTagger.d0min = -60. -seederTagger.d0max = 0. +seeder_tagger = tracking.SeedFinderProcessor("SeedTagger") +seeder_tagger.input_hits_collection = digi_tagger.out_collection +seeder_tagger.out_seed_collection = "TaggerRecoSeeds" +seeder_tagger.pmin = 0.1 +seeder_tagger.pmax = 10.0 +seeder_tagger.d0min = -45. +seeder_tagger.d0max = 45. +seeder_tagger.z0max = 60. #Seed finder processor - Recoil -seederRecoil = tracking.SeedFinderProcessor("SeedRecoil") -seederRecoil.perigee_location = [0.,0.,0.] -seederRecoil.input_hits_collection = digiRecoil.out_collection -seederRecoil.out_seed_collection = "RecoilRecoSeeds" -seederRecoil.bfield = 1.5 -seederRecoil.pmin = 0.1 -seederRecoil.pmax = 4. -seederRecoil.d0min = -0.5 -seederRecoil.d0max = 0.5 -seederRecoil.z0max = 10. +seeder_recoil = tracking.SeedFinderProcessor("SeedRecoil") +seeder_recoil.perigee_location = [0.,0.,0.] +seeder_recoil.input_hits_collection = digi_recoil.out_collection +seeder_recoil.out_seed_collection = "RecoilRecoSeeds" +seeder_recoil.bfield = 1.5 +seeder_recoil.pmin = 0.1 +seeder_recoil.pmax = 10.0 +seeder_recoil.d0min = -40.0 +seeder_recoil.d0max = 40.0 +seeder_recoil.z0max = 50. # Producer for running the CKF track finding starting from the found seeds. @@ -93,14 +108,14 @@ tracking_tagger.const_b_field = False #Target location for the CKF extrapolation -tracking_tagger.seed_coll_name = seederTagger.out_seed_collection +tracking_tagger.seed_coll_name = seeder_tagger.out_seed_collection tracking_tagger.out_trk_collection = "TaggerTracks" #smear the hits used for finding/fitting tracking_tagger.trackID = -1 #1 tracking_tagger.pdgID = -9999 #11 -tracking_tagger.measurement_collection = digiTagger.out_collection -tracking_tagger.min_hits = 5 +tracking_tagger.measurement_collection = digi_tagger.out_collection +tracking_tagger.min_hits = 6 #CKF Options @@ -110,32 +125,51 @@ tracking_recoil.propagator_step_size = 1000. #mm tracking_recoil.bfield = -1.5 #in T #From looking at the BField map tracking_recoil.const_b_field = False - +tracking_recoil.taggerTracking=False #Target location for the CKF extrapolation -#tracking_recoil.seed_coll_name = seederRecoil.out_seed_collection +#tracking_recoil.seed_coll_name = seeder_recoil.out_seed_collection tracking_recoil.seed_coll_name = "RecoilTruthSeeds" tracking_recoil.out_trk_collection = "RecoilTracks" #smear the hits used for finding/fitting tracking_recoil.trackID = -1 #1 tracking_recoil.pdgID = -9999 #11 -tracking_recoil.measurement_collection = digiRecoil.out_collection -tracking_recoil.min_hits = 5 +tracking_recoil.measurement_collection = digi_recoil.out_collection +tracking_recoil.min_hits = 6 from LDMX.Tracking import dqm -digi_dqm = dqm.TrackerDigiDQM() -tracking_dqm = dqm.TrackingRecoDQM() + +seed_tagger_dqm = dqm.TrackingRecoDQM("SeedTaggerDQM") +seed_tagger_dqm.track_collection = seeder_tagger.out_seed_collection +seed_tagger_dqm.truth_collection = "TaggerTruthTracks" +seed_tagger_dqm.title = "" +seed_tagger_dqm.buildHistograms() + +tagger_dqm = dqm.TrackingRecoDQM("TaggerDQM") +tagger_dqm.track_collection = tracking_tagger.out_trk_collection +tagger_dqm.truth_collection = "TaggerTruthTracks" +tagger_dqm.trackStates = ["target"] +tagger_dqm.title = "" +tagger_dqm.measurement_collection=digi_tagger.out_collection +tagger_dqm.truth_hit_collection="TaggerSimHits" +tagger_dqm.buildHistograms() + seed_recoil_dqm = dqm.TrackingRecoDQM("SeedRecoilDQM") -seed_recoil_dqm.track_collection = seederRecoil.out_seed_collection +seed_recoil_dqm.track_collection = seeder_recoil.out_seed_collection seed_recoil_dqm.truth_collection = "RecoilTruthTracks" seed_recoil_dqm.title = "" +seed_recoil_dqm.buildHistograms() recoil_dqm = dqm.TrackingRecoDQM("RecoilDQM") recoil_dqm.track_collection = tracking_recoil.out_trk_collection recoil_dqm.truth_collection = "RecoilTruthTracks" +recoil_dqm.trackStates = ["ecal","target"] recoil_dqm.title = "" +recoil_dqm.measurement_collection=digi_recoil.out_collection +recoil_dqm.truth_hit_collection="RecoilSimHits" +recoil_dqm.buildHistograms() # This sequence runs the digitization in the tagger and recoil # Then the truth tracking to have TruthTracks in the final state @@ -143,40 +177,22 @@ # Track finding is then raun in the tagger and in the recoil # Finally two dqm examples are run in the recoil tracks and using the seed tracks -p.sequence = [digiTagger, digiRecoil, - truth_tracking, - seederTagger, seederRecoil, +p.sequence = [sim,truth_tracking,digi_tagger, digi_recoil, + seeder_tagger, seeder_recoil, tracking_tagger, tracking_recoil, - recoil_dqm, seed_recoil_dqm] - -# The input file to be added. -# for now, we just take the first argument on the command line -import sys -p.inputFiles = [sys.argv[1]] - -# Example about how to drop some of the collections in the output file. -p.keep = [ - # "drop .*SimHits.*", #drop all sim hits - "drop .*Ecal.*", #drop all ecal (Digis are not removed) - # "drop .*Magnet*", - "drop .*Hcal.*", - "drop .*Scoring.*", - "drop .*SimParticles.*", - "drop .*TriggerPad.*", - "drop .*trig.*" -] + recoil_dqm, seed_recoil_dqm, + tagger_dqm, seed_tagger_dqm] # Output name # just append '_withTracking' to the name of the input file from pathlib import Path -input_filepath = Path(p.inputFiles[0]) -p.outputFiles = [str(input_filepath.with_suffix(''))+'_withTracking.root'] +p.outputFiles = ['test_8gev_electrons_withTracking.root'] # lower log level so 'info' and above messages can be printed p.termLogLevel=1 # Number of events -p.maxEvents = 10000 +p.maxEvents = n_evts # Where to store DQM plots p.histogramFile = "test_dqmMonitoringFile.root" diff --git a/Tracking/include/Tracking/Reco/CKFProcessor.h b/Tracking/include/Tracking/Reco/CKFProcessor.h index 6ac96299e..41ba929fc 100644 --- a/Tracking/include/Tracking/Reco/CKFProcessor.h +++ b/Tracking/include/Tracking/Reco/CKFProcessor.h @@ -19,7 +19,6 @@ #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/Definitions/Units.hpp" #include "Acts/EventData/TrackParameters.hpp" -#include "Acts/EventData/detail/TransformationFreeToBound.hpp" #include "Acts/Utilities/Logger.hpp" // geometry @@ -76,22 +75,17 @@ //--- Interpolated magnetic field ---// #include "Tracking/Sim/BFieldXYZUtils.h" +// mg Aug 2024 not sure if these are needed... +using Updater = Acts::GainMatrixUpdater; +using Smoother = Acts::GainMatrixSmoother; using ActionList = Acts::ActionList; using AbortList = Acts::AbortList; using CkfPropagator = Acts::Propagator, Acts::Navigator>; -using GsfPropagator = Acts::Propagator< - Acts::MultiEigenStepperLoop< - Acts::StepperExtensionList< - Acts::detail::GenericDefaultExtension>, - Acts::WeightedComponentReducerLoop, Acts::detail::VoidAuctioneer>, - Acts::Navigator>; - -//?! -// using PropagatorOptions = -// Acts::DenseStepperPropagatorOptions; +using TrackContainer = Acts::TrackContainer; namespace tracking { namespace reco { @@ -161,8 +155,10 @@ class CKFProcessor final : public TrackingGeometryUser { // time profiling std::map profiling_map_; - bool debug_{false}; + bool debug_acts_{false}; + std::shared_ptr target_surface; + Acts::RotationMatrix3 surf_rotation; // Constant BField double bfield_{0}; // Use constant bfield @@ -209,8 +205,8 @@ class CKFProcessor final : public TrackingGeometryUser { std::unique_ptr propagator_; // The CKF - std::unique_ptr> + std::unique_ptr< + const Acts::CombinatorialKalmanFilter> ckf_; // Track Extrapolator Tool diff --git a/Tracking/include/Tracking/Reco/CustomStatePropagator.h b/Tracking/include/Tracking/Reco/CustomStatePropagator.h deleted file mode 100644 index fa11d3bb9..000000000 --- a/Tracking/include/Tracking/Reco/CustomStatePropagator.h +++ /dev/null @@ -1,98 +0,0 @@ -#pragma once - -//--- Framework ---// -#include "Framework/Configure/Parameters.h" -#include "Framework/EventProcessor.h" - -//--- ACTS ---// -#include "Acts/Definitions/Units.hpp" -#include "Acts/Geometry/GeometryContext.hpp" -#include "Acts/MagneticField/MagneticFieldContext.hpp" -#include "Acts/MagneticField/MagneticFieldProvider.hpp" -#include "Acts/Propagator/AbortList.hpp" -#include "Acts/Propagator/EigenStepper.hpp" -#include "Acts/Propagator/Propagator.hpp" -#include "Acts/Propagator/StandardAborters.hpp" -#include "Acts/Surfaces/PerigeeSurface.hpp" -#include "Acts/Surfaces/PlaneSurface.hpp" - -//--- Tracking ---// -#include "Tracking/Sim/BFieldXYZUtils.h" -#include "Tracking/Sim/TrackingUtils.h" - -//--- C++ ---// -#include - -//--- ROOT ---// -#include - -#include "TFile.h" -#include "TH1F.h" -#include "TTree.h" - -using AbortList = Acts::AbortList; - -namespace tracking::reco { - -class CustomStatePropagator : public framework::Producer { - public: - CustomStatePropagator(const std::string& name, framework::Process& process); - ~CustomStatePropagator(); - - void onProcessStart() override; - void onProcessEnd() override; - - void configure(framework::config::Parameters& parameters) override; - - void produce(framework::Event& event) override{}; - - void fillTree(int state, int q, const Acts::Vector3 gen_pos, - const Acts::Vector3 gen_mom, - const Acts::BoundTrackParameters& endParams); - - Acts::GeometryContext gctx_; - Acts::MagneticFieldContext bctx_; - - // The interpolated bfield - std::string field_map_{""}; - double surf_location_{0.}; - int nstates_{0}; - std::vector bs_size_; - std::vector prange_; - std::vector thetarange_; - std::vector phirange_; - - // Output ntuple - TFile* outFile_; - TTree* outTree_; - std::shared_ptr histo_end_px; - std::shared_ptr histo_end_py; - std::shared_ptr histo_end_pz; - std::shared_ptr histo_gen_px; - std::shared_ptr histo_gen_py; - std::shared_ptr histo_gen_pz; - std::shared_ptr histo_gen_p; - std::shared_ptr histo_end_p; - std::shared_ptr histo_loc01; - - double state_nr{0.}; - int charge_{0}; - double gen_x{0.}; - double gen_y{0.}; - double gen_z{0.}; - double gen_px{0.}; - double gen_py{0.}; - double gen_pz{0.}; - - double end_x{0.}; - double end_y{0.}; - double end_z{0.}; - double end_loc0{0.}; - double end_loc1{0.}; - - double end_px{0.}; - double end_py{0.}; - double end_pz{0.}; -}; - -} // namespace tracking::reco diff --git a/Tracking/include/Tracking/Reco/DigitizationProcessor.h b/Tracking/include/Tracking/Reco/DigitizationProcessor.h index 548143178..2378b26b9 100644 --- a/Tracking/include/Tracking/Reco/DigitizationProcessor.h +++ b/Tracking/include/Tracking/Reco/DigitizationProcessor.h @@ -5,9 +5,6 @@ //--- ACTS ---// #include "Acts/Definitions/Units.hpp" -#include "Acts/Digitization/CartesianSegmentation.hpp" -#include "Acts/Digitization/DigitizationModule.hpp" -#include "Acts/Digitization/PlanarModuleStepper.hpp" #include "Acts/Surfaces/RectangleBounds.hpp" #include "Acts/Surfaces/Surface.hpp" diff --git a/Tracking/include/Tracking/Reco/GSFProcessor.h b/Tracking/include/Tracking/Reco/GSFProcessor.h index dcf9f1da1..773afb541 100644 --- a/Tracking/include/Tracking/Reco/GSFProcessor.h +++ b/Tracking/include/Tracking/Reco/GSFProcessor.h @@ -19,7 +19,6 @@ #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/Definitions/Units.hpp" #include "Acts/EventData/TrackParameters.hpp" -#include "Acts/EventData/detail/TransformationFreeToBound.hpp" #include "Acts/Utilities/Logger.hpp" // geometry @@ -66,6 +65,7 @@ #include "Acts/Propagator/MultiEigenStepperLoop.hpp" #include "Acts/TrackFitting/BetheHeitlerApprox.hpp" #include "Acts/TrackFitting/GaussianSumFitter.hpp" +#include "Acts/TrackFitting/GsfMixtureReduction.hpp" //--- Tracking ---// #include "Tracking/Event/Measurement.h" @@ -93,7 +93,7 @@ using AbortList = Acts::AbortList; using MultiStepper = Acts::MultiEigenStepperLoop<>; using Propagator = Acts::Propagator, Acts::Navigator>; using GsfPropagator = Acts::Propagator; -using BetheHeitlerApprox = Acts::Experimental::AtlasBetheHeitlerApprox<6, 5>; +using BetheHeitlerApprox = Acts::AtlasBetheHeitlerApprox<6, 5>; namespace tracking { namespace reco { @@ -216,7 +216,7 @@ class GSFProcessor final : public TrackingGeometryUser { std::string seed_coll_name_{"seedTracks"}; // The GSF Fitter - std::unique_ptr> gsf_; @@ -237,6 +237,9 @@ class GSFProcessor final : public TrackingGeometryUser { bool usePerigee_{false}; // bool usePlaneSurface_{false}; + // Keep track on which system this processor is running on + bool taggerTracking_{true}; + // The Propagators std::unique_ptr propagator_; diff --git a/Tracking/include/Tracking/Reco/SeedFinderProcessor.h b/Tracking/include/Tracking/Reco/SeedFinderProcessor.h index 0503f431c..d41759ec1 100644 --- a/Tracking/include/Tracking/Reco/SeedFinderProcessor.h +++ b/Tracking/include/Tracking/Reco/SeedFinderProcessor.h @@ -20,8 +20,6 @@ //---< ACTS >---// #include "Acts/Definitions/Algebra.hpp" #include "Acts/MagneticField/MagneticFieldContext.hpp" -#include "Acts/Seeding/BinFinder.hpp" -#include "Acts/Seeding/BinnedSPGroup.hpp" #include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp" #include "Acts/Seeding/Seed.hpp" #include "Acts/Seeding/SeedFilter.hpp" diff --git a/Tracking/include/Tracking/Reco/TrackExtrapolatorTool.h b/Tracking/include/Tracking/Reco/TrackExtrapolatorTool.h index a00783256..d06d06f71 100644 --- a/Tracking/include/Tracking/Reco/TrackExtrapolatorTool.h +++ b/Tracking/include/Tracking/Reco/TrackExtrapolatorTool.h @@ -51,16 +51,18 @@ class TrackExtrapolatorTool { @return optional with BoundTrackParameters */ + using PropagatorOptions = + typename propagator_t::template Options; + std::optional extrapolate( const Acts::BoundTrackParameters pars, const std::shared_ptr& target_surface) { - // Just to make it explicit - bool boundaryCheck = false; - auto intersection = target_surface->intersect( - gctx_, pars.position(gctx_), pars.unitDirection(), boundaryCheck); + auto intersection = target_surface->intersect(gctx_, pars.position(gctx_), + pars.direction()); + + PropagatorOptions pOptions(gctx_, mctx_); - Acts::PropagatorOptions pOptions(gctx_, mctx_); - pOptions.direction = intersection.intersection.pathLength >= 0 + pOptions.direction = intersection.intersections()[0].pathLength() >= 0 ? Acts::Direction::Forward : Acts::Direction::Backward; @@ -98,15 +100,15 @@ class TrackExtrapolatorTool { template std::optional extrapolate( track_t track, const std::shared_ptr& target_surface) { - auto outermost = *(track.trackStates().begin()); - auto begin = track.trackStates().begin(); + // get first and last track state on surface + auto outermost = *(track.trackStatesReversed().begin()); + auto begin = track.trackStatesReversed().begin(); std::advance(begin, track.nTrackStates() - 1); auto innermost = *begin; // I'm checking which track state is closer to the origin of the target // surface to decide from where to start the extrapolation to the surface. I // use the coordinate along the beam axis. - double first_dis = std::abs( innermost.referenceSurface().transform(gctx_).translation()(0) - target_surface->transform(gctx_).translation()(0)); @@ -119,26 +121,26 @@ class TrackExtrapolatorTool { const auto& ts = first_dis < last_dis ? innermost : outermost; - // std::cout<<"Selected track state for extrapolation"< 0 - ? 1 * Acts::UnitConstants::e - : -1 * Acts::UnitConstants::e; - - Acts::BoundTrackParameters sp(surface.getSharedPtr(), smoothed, q, cov); - + if (debug_) { + std::cout << "Surface::" << surface.transform(gctx_).translation() + << std::endl; + if (hasSmoothed) + std::cout << "Smoothed::" << smoothed.transpose() << std::endl; + std::cout << "HasSmoothed::" << hasSmoothed << std::endl; + std::cout << "Filtered::" << filtered.transpose() << std::endl; + } + // mg Aug 2024 ... v36 takes the particle...assume electron + auto partHypo{Acts::SinglyChargedParticleHypothesis::electron()}; + Acts::BoundTrackParameters sp(surface.getSharedPtr(), smoothed, cov, + partHypo); return extrapolate(sp, target_surface); } @@ -149,6 +151,7 @@ class TrackExtrapolatorTool { // Now.. I'm taking whatever it is. I'm not checking here if it is a // measurement. + auto& tsc = track.container().trackStateContainer(); auto begin = track.trackStates().begin(); auto ts_last = *begin; const auto& surface = (ts_last).referenceSurface(); @@ -156,19 +159,14 @@ class TrackExtrapolatorTool { const auto& cov = (ts_last).smoothedCovariance(); // Get the BoundTrackStateParameters - - Acts::ActsScalar q = smoothed[Acts::eBoundQOverP] > 0 - ? 1 * Acts::UnitConstants::e - : -1 * Acts::UnitConstants::e; + // assume electron for now + auto partHypo{Acts::SinglyChargedParticleHypothesis::electron()}; Acts::BoundTrackParameters state_parameters(surface.getSharedPtr(), - smoothed, q, cov); + smoothed, cov, partHypo); // One can also use directly the extrapolate method - - Acts::PropagatorOptions pOptions(gctx_, mctx_); - pOptions.direction = Acts::Direction::Forward; - + PropagatorOptions pOptions(gctx_, mctx_); auto result = propagator_.propagate(state_parameters, *target_surface, pOptions); @@ -178,6 +176,8 @@ class TrackExtrapolatorTool { return std::nullopt; } + /** + /** Create an ldmx::TrackState to the extrapolated position @param Acts::Track @param extrapolation surface diff --git a/Tracking/include/Tracking/Reco/TruthSeedProcessor.h b/Tracking/include/Tracking/Reco/TruthSeedProcessor.h index 2c570fa78..8cb1ab561 100644 --- a/Tracking/include/Tracking/Reco/TruthSeedProcessor.h +++ b/Tracking/include/Tracking/Reco/TruthSeedProcessor.h @@ -21,7 +21,6 @@ #include "Acts/Definitions/Algebra.hpp" #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/EventData/TrackParameters.hpp" -#include "Acts/EventData/detail/TransformationFreeToBound.hpp" #include "Acts/Propagator/Navigator.hpp" #include "Acts/Surfaces/PerigeeSurface.hpp" diff --git a/Tracking/include/Tracking/Sim/BFieldXYZUtils.h b/Tracking/include/Tracking/Sim/BFieldXYZUtils.h index 83145344a..da3df3022 100644 --- a/Tracking/include/Tracking/Sim/BFieldXYZUtils.h +++ b/Tracking/include/Tracking/Sim/BFieldXYZUtils.h @@ -8,17 +8,17 @@ #include "Acts/MagneticField/BFieldMapUtils.hpp" #include "Acts/MagneticField/InterpolatedBFieldMap.hpp" #include "Acts/MagneticField/MagneticFieldContext.hpp" +#include "Acts/Utilities/AxisFwd.hpp" +#include "Acts/Utilities/Grid.hpp" #include "Acts/Utilities/Interpolation.hpp" #include "Acts/Utilities/Result.hpp" -#include "Acts/Utilities/detail/AxisFwd.hpp" -#include "Acts/Utilities/detail/Grid.hpp" static const double DIPOLE_OFFSET = 400.; // 400 mm -using InterpolatedMagneticField3 = - Acts::InterpolatedBFieldMap>; +using InterpolatedMagneticField3 = Acts::InterpolatedBFieldMap< + Acts::Grid, + Acts::Axis, + Acts::Axis>>; using GenericTransformPos = std::function; using GenericTransformBField = @@ -113,17 +113,17 @@ inline InterpolatedMagneticField3 rotateFieldMapXYZ( nBinsY = 2 * nBinsY - 1; nBinsZ = 2 * nBinsZ - 1; } - Acts::detail::EquidistantAxis xAxis(xMin * lengthUnit, xMax * lengthUnit, - nBinsX); - Acts::detail::EquidistantAxis yAxis(yMin * lengthUnit, yMax * lengthUnit, - nBinsY); - Acts::detail::EquidistantAxis zAxis(zMin * lengthUnit, zMax * lengthUnit, - nBinsZ); + Acts::Axis xAxis(xMin * lengthUnit, + xMax * lengthUnit, nBinsX); + Acts::Axis yAxis(yMin * lengthUnit, + yMax * lengthUnit, nBinsY); + Acts::Axis zAxis(zMin * lengthUnit, + zMax * lengthUnit, nBinsZ); // Create the grid using Grid_t = - Acts::detail::Grid; + Acts::Grid, + Acts::Axis, + Acts::Axis>; Grid_t grid( std::make_tuple(std::move(xAxis), std::move(yAxis), std::move(zAxis))); diff --git a/Tracking/include/Tracking/Sim/LdmxSpacePoint.h b/Tracking/include/Tracking/Sim/LdmxSpacePoint.h index 7bfffb702..f74d343d2 100644 --- a/Tracking/include/Tracking/Sim/LdmxSpacePoint.h +++ b/Tracking/include/Tracking/Sim/LdmxSpacePoint.h @@ -8,7 +8,7 @@ // sort of sense // TODO:: The projector should support time as well - +// MG Aug 2024 Change all SymMatrix2 to SquareMatrix2 to support ACTs v36 #include //--- Acts ---// @@ -102,14 +102,14 @@ class LdmxSpacePoint { projector_(1, 1) = 1; } - const Acts::SymMatrix2 getLocalCovariance() const { return local_cov_; }; + const Acts::SquareMatrix2 getLocalCovariance() const { return local_cov_; }; const Acts::Vector3 getGlobalPosition() const { return global_pos_; }; const Acts::Vector2 getLocalPosition() const { return local_pos_; }; Acts::Vector3 global_pos_; Acts::Vector2 local_pos_; // TODO:: not sure about this - Acts::SymMatrix2 local_cov_; + Acts::SquareMatrix2 local_cov_; // Projection matrix from the full space to the (u,v) space. // This can be expanded to (u,v,t) space in the case time needs to be added. diff --git a/Tracking/include/Tracking/Sim/MeasurementCalibrator.h b/Tracking/include/Tracking/Sim/MeasurementCalibrator.h index f66bb2454..9e08ef0f8 100644 --- a/Tracking/include/Tracking/Sim/MeasurementCalibrator.h +++ b/Tracking/include/Tracking/Sim/MeasurementCalibrator.h @@ -7,6 +7,7 @@ #include "Acts/EventData/MultiTrajectory.hpp" #include "Acts/EventData/SourceLink.hpp" #include "Acts/EventData/VectorMultiTrajectory.hpp" +#include "Acts/Utilities/CalibrationContext.hpp" #include "Tracking/Event/Measurement.h" #include "Tracking/Sim/IndexSourceLink.h" #include "Tracking/Sim/LdmxSpacePoint.h" @@ -53,14 +54,13 @@ class LdmxMeasurementCalibrator { /// @tparam parameters_t Track parameters type /// @param gctx The geometry context (unused) /// @param trackState The track state to calibrate - void calibrate( - const Acts::GeometryContext& /*gctx*/, - Acts::MultiTrajectory::TrackStateProxy - trackState) const { - ActsExamples::IndexSourceLink sourceLink = - trackState.getUncalibratedSourceLink() - .get(); - + template + void calibrate(const Acts::GeometryContext& /*gctx*/, + const Acts::CalibrationContext& /*cctx*/, + const Acts::SourceLink& genericSourceLink /*sourceLink*/, + typename traj_t::TrackStateProxy trackState) const { + ActsExamples::IndexSourceLink sourceLink{ + genericSourceLink.get()}; assert(m_measurements and "Undefined measurement container in LdmxMeasurementCalibrator"); assert((sourceLink.index() < m_measurements->size()) and @@ -68,20 +68,19 @@ class LdmxMeasurementCalibrator { "LdmxMeasurementCalibrator"); auto meas = m_measurements->at(sourceLink.index()); - - trackState.calibrated<2>().setZero(); - Acts::Vector2 local_pos{meas.getLocalPosition()[0], meas.getLocalPosition()[1]}; - trackState.calibrated<2>().head<2>() = local_pos; - // trackState.data().measdim = 2; - trackState.calibratedCovariance<2>().setZero(); - - Acts::SymMatrix2 local_cov; + auto tsCal{trackState.template calibrated<2>()}; + auto tsCalCov{trackState.template calibratedCovariance<2>()}; + tsCal.setZero(); + tsCal.template head<2>() = local_pos; + Acts::SquareMatrix2 local_cov; local_cov.setZero(); local_cov(0, 0) = meas.getLocalCovariance()[0]; local_cov(1, 1) = meas.getLocalCovariance()[1]; - trackState.calibratedCovariance<2>().block<2, 2>(0, 0) = local_cov; + tsCalCov.setZero(); + // make tsCalCov 2x2 block the local_cov we just set + tsCalCov.block(0, 0, 2, 2) = local_cov; Acts::ActsMatrix<2, 6> projector; projector.setZero(); @@ -97,14 +96,13 @@ class LdmxMeasurementCalibrator { /// @tparam parameters_t Track parameters type /// @param gctx The geometry context (unused) /// @param trackState The track state to calibrate - void calibrate_1d( - const Acts::GeometryContext& /*gctx*/, - Acts::MultiTrajectory::TrackStateProxy - trackState) const { - // Use by value - life management is not working properly - ActsExamples::IndexSourceLink sourceLink = - trackState.getUncalibratedSourceLink() - .get(); + template + void calibrate_1d(const Acts::GeometryContext& /*gctx*/, + const Acts::CalibrationContext& /*cctx*/, + const Acts::SourceLink& genericSourceLink /*sourceLink*/, + typename traj_t::TrackStateProxy trackState) const { + ActsExamples::IndexSourceLink sourceLink{ + genericSourceLink.get()}; assert(m_measurements and "Undefined measurement container in LdmxMeasurementCalibrator"); @@ -112,22 +110,23 @@ class LdmxMeasurementCalibrator { "Source link index is outside the container bounds in " "LdmxMeasurementCalibrator"); - // std::cout<<"calibrate_1d ==> NMeasurements available=" << - // m_measurements->size()<at(sourceLink.index()); - // You need to explicitly allocate measurements here trackState.allocateCalibrated(1); - trackState.calibrated<1>().setZero(); - trackState.calibrated<1>()(0) = (meas.getLocalPosition())[0]; - trackState.calibratedCovariance<1>().setZero(); - trackState.calibratedCovariance<1>()(0, 0) = (meas.getLocalCovariance())[0]; + auto tsCal{trackState.template calibrated<1>()}; + auto tsCalCov{trackState.template calibratedCovariance<1>()}; + + tsCal.setZero(); + tsCal(0) = (meas.getLocalPosition())[0]; + tsCalCov.setZero(); + tsCalCov(0, 0) = (meas.getLocalCovariance())[0]; Acts::ActsMatrix<2, 6> projector; projector.setZero(); projector(0, 0) = 1.; projector(1, 1) = 1.; trackState.setProjector(projector.row(0)); + trackState.setUncalibratedSourceLink(genericSourceLink); } // Function to test the measurement calibrator @@ -136,8 +135,6 @@ class LdmxMeasurementCalibrator { void test(const Acts::GeometryContext& /*gctx*/, const ActsExamples::IndexSourceLink& sourceLink) const { auto meas = m_measurements->at(sourceLink.index()); - // get the measurement - std::cout << "Measurement layer::\n" << meas.getLayer() << std::endl; Acts::Vector3 global_pos{meas.getGlobalPosition()[0], meas.getGlobalPosition()[1], diff --git a/Tracking/include/Tracking/Sim/PropagatorStepWriter.h b/Tracking/include/Tracking/Sim/PropagatorStepWriter.h index ad9e0c4a8..794b90d36 100644 --- a/Tracking/include/Tracking/Sim/PropagatorStepWriter.h +++ b/Tracking/include/Tracking/Sim/PropagatorStepWriter.h @@ -58,7 +58,7 @@ class PropagatorStepWriter { TTree* m_outputTree; ///< the output tree int m_eventNr; ///< the event number of - std::vector m_volumeID; ///< volume identifier + // std::vector m_volumeID; ///< volume identifier std::vector m_boundaryID; ///< boundary identifier std::vector m_layerID; ///< layer identifier if std::vector m_approachID; ///< surface identifier diff --git a/Tracking/include/Tracking/Sim/SeedToTrackParamMaker.ipp b/Tracking/include/Tracking/Sim/SeedToTrackParamMaker.ipp index 8fb3281ff..6c17a85e8 100644 --- a/Tracking/include/Tracking/Sim/SeedToTrackParamMaker.ipp +++ b/Tracking/include/Tracking/Sim/SeedToTrackParamMaker.ipp @@ -129,7 +129,7 @@ bool SeedToTrackParamMaker::KarimakiFit( double Vphid_m1 = u * (rho * Sgamma - u * Sbeta); double Vdd_m1 = rho * (rho * Saa - 2 * u * Salpha) + u * u * Sw; - Acts::SymMatrix3 kariCov_inv; + Acts::SquareMatrix3 kariCov_inv; // TODO:: This one should be fixed kariCov_inv(0, 0) = Vrr_m1; @@ -145,7 +145,7 @@ bool SeedToTrackParamMaker::KarimakiFit( kariCov_inv(0, 2) = Vrhod_m1; kariCov_inv(2, 2) = Vdd_m1; - Acts::SymMatrix3 kariCov = kariCov_inv.inverse(); + Acts::SquareMatrix3 kariCov = kariCov_inv.inverse(); // Compute the corrections to the estimated parameters double sigma = -rho * Sdelta + 2 * u * Saa - d * (1 + u) * Salpha; diff --git a/Tracking/include/Tracking/Sim/TrackingUtils.h b/Tracking/include/Tracking/Sim/TrackingUtils.h index ae0b65ea9..9ebab69b9 100644 --- a/Tracking/include/Tracking/Sim/TrackingUtils.h +++ b/Tracking/include/Tracking/Sim/TrackingUtils.h @@ -33,10 +33,12 @@ // --- < ACTS > --- // #include "Acts/Definitions/Algebra.hpp" +#include "Acts/Definitions/PdgParticle.hpp" #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/Definitions/Units.hpp" #include "Acts/EventData/TrackParameters.hpp" #include "Acts/Surfaces/PerigeeSurface.hpp" +#include "Acts/Surfaces/PlaneSurface.hpp" #include "Tracking/Event/Measurement.h" namespace tracking { @@ -138,15 +140,18 @@ inline ldmx::LdmxSpacePoint* convertSimHitToLdmxSpacePoint( sigma_v * sigma_v, hit.getID()); } -inline void flatCov(Acts::BoundSymMatrix cov, std::vector& v_cov) { +// BoundSymMatrix doesn't exist in v36 .. use BoundSquareMatrix +// have to change this everywhere .. I think using BoundSysMatrix was defined +// exactly the same as BoundSquareMatrix is now in ACTs +inline void flatCov(Acts::BoundSquareMatrix cov, std::vector& v_cov) { v_cov.clear(); v_cov.reserve(cov.rows() * (cov.rows() + 1) / 2); for (int i = 0; i < cov.rows(); i++) for (int j = i; j < cov.cols(); j++) v_cov.push_back(cov(i, j)); } -inline Acts::BoundSymMatrix unpackCov(const std::vector& v_cov) { - Acts::BoundSymMatrix cov; +inline Acts::BoundSquareMatrix unpackCov(const std::vector& v_cov) { + Acts::BoundSquareMatrix cov; int e{0}; for (int i = 0; i < cov.rows(); i++) for (int j = i; j < cov.cols(); j++) { @@ -167,7 +172,7 @@ inline Acts::BoundSymMatrix unpackCov(const std::vector& v_cov) { inline Acts::Vector3 Ldmx2Acts(Acts::Vector3 ldmx_v) { // TODO::Move it to a static member - Acts::SymMatrix3 acts_rot; + Acts::SquareMatrix3 acts_rot; acts_rot << 0., 0., 1., 1., 0., 0., 0., 1, 0.; return acts_rot * ldmx_v; @@ -219,21 +224,31 @@ inline Acts::BoundVector boundState(const ldmx::Track::TrackState& ts) { inline Acts::BoundTrackParameters boundTrackParameters( const ldmx::Track& trk, std::shared_ptr perigee) { Acts::BoundVector paramVec = boundState(trk); - Acts::BoundSymMatrix covMat = unpackCov(trk.getPerigeeCov()); - return Acts::BoundTrackParameters(perigee, paramVec, std::move(covMat)); + Acts::BoundSquareMatrix covMat = unpackCov(trk.getPerigeeCov()); + auto partHypo{Acts::SinglyChargedParticleHypothesis::electron()}; + // auto + // part{Acts::GenericParticleHypothesis(Acts::ParticleHypothesis(Acts::PdgParticle(trk.getPdgID())))}; + // return Acts::BoundTrackParameters(perigee, paramVec, std::move(covMat)); + // need to add particle hypothesis + return Acts::BoundTrackParameters(perigee, paramVec, std::move(covMat), + partHypo); } inline Acts::BoundTrackParameters btp(const ldmx::Track::TrackState& ts, - std::shared_ptr surf) { + std::shared_ptr surf, + int pdgid) { Acts::BoundVector paramVec = boundState(ts); - Acts::BoundSymMatrix covMat = unpackCov(ts.cov); - return Acts::BoundTrackParameters(surf, paramVec, std::move(covMat)); + Acts::BoundSquareMatrix covMat = unpackCov(ts.cov); + auto partHypo{Acts::SinglyChargedParticleHypothesis::electron()}; + // auto + // part{Acts::GenericParticleHypothesis(Acts::ParticleHypothesis(Acts::PdgParticle(pdgid)))}; + return Acts::BoundTrackParameters(surf, paramVec, std::move(covMat), + partHypo); } // Return an unbound surface -inline const std::shared_ptr unboundSurface(double xloc, - double yloc = 0., - double zloc = 0.) { +inline const std::shared_ptr unboundSurface( + double xloc, double yloc = 0., double zloc = 0.) { // Define the target surface - be careful: // x - downstream // y - left (when looking along x) diff --git a/Tracking/include/Tracking/dqm/TrackingRecoDQM.h b/Tracking/include/Tracking/dqm/TrackingRecoDQM.h index 8912b29e2..4373b8368 100644 --- a/Tracking/include/Tracking/dqm/TrackingRecoDQM.h +++ b/Tracking/include/Tracking/dqm/TrackingRecoDQM.h @@ -4,6 +4,7 @@ #include "Framework/Event.h" #include "Framework/EventProcessor.h" #include "SimCore/Event/SimTrackerHit.h" +#include "Tracking/Event/Measurement.h" #include "Tracking/Event/Track.h" #include "Tracking/Event/TruthTrack.h" @@ -31,10 +32,12 @@ class TrackingRecoDQM : public framework::Analyzer { void analyze(const framework::Event& event) override; void TrackMonitoring(const std::vector& tracks, + const std::vector& measurements, const std::string title, const bool& doDetail, const bool& doTruth); void EfficiencyPlots(const std::vector& tracks, + const std::vector& measurements, const std::string& title); /** Monitoring plots for tracks extrapolated to the ECAL Scoring plane. @@ -65,6 +68,7 @@ class TrackingRecoDQM : public framework::Analyzer { private: std::string trackCollection_{"TruthTracks"}; std::string truthCollection_{"TaggerTruthTracks"}; + std::string measurementCollection_{"DigiTaggerSimHits"}; std::string title_{"tagger_trk_"}; double trackProb_cut_{0.5}; std::string subdetector_{"Tagger"}; diff --git a/Tracking/python/dqm.py b/Tracking/python/dqm.py index 957d3d531..1f399a00c 100644 --- a/Tracking/python/dqm.py +++ b/Tracking/python/dqm.py @@ -35,13 +35,13 @@ def __init__(self, name="TrackingRecoDQM"): self.d0max = 50# 15. self.z0min = -200#-45. self.z0max = 200# 45. - self.phimin = -2#-0.2 - self.phimax = 2#0.2 - self.thetamin = -2#1.4 - self.thetamax = 4#1.7 - self.qopmin = -50#-10 - self.qopmax = 50#10 - self.pmax = 8. + self.phimin = -0.1#-0.2 + self.phimax = 0.1#0.2 + self.thetamin = 1.4#1.4 + self.thetamax = 1.7#1.7 + self.qopmin = -10#-10 + self.qopmax = 10#10 + self.pmax = 10. self.trackStates = ["ecal","target"] self.doTruth= True @@ -85,12 +85,14 @@ def buildHistograms(self) : "qOverP [GeV^{-1}]",nbins,qopmin,qopmax) self.build1DHistogram("pt_bending", - "pT bending plane [GeV]",nbins,-pmax,pmax) + "pT bending plane [GeV]",nbins,0.,pmax) self.build1DHistogram("pt_beam", "pT beam axis [GeV]",nbins,0,0.5) self.build1DHistogram("nHits", "nHits",15,0,15) + self.build1DHistogram("LayersHit", + "LayersHit",15,0,15) self.build1DHistogram("Chi2", "Chi2",nbins,0,100) self.build1DHistogram("ndf", @@ -102,11 +104,11 @@ def buildHistograms(self) : self.build1DHistogram("nHoles", "nHoles",5,0,5) self.build1DHistogram("px", - "pX [GeV]",nbins,-pmax,pmax) + "pX [GeV]",nbins,0.0,pmax) self.build1DHistogram("py", - "pY [GeV]",nbins,-pmax,pmax) + "pY [GeV]",nbins,-pmax/20.0,pmax/20.0) self.build1DHistogram("pz", - "pZ [GeV]",nbins,-pmax,pmax) + "pZ [GeV]",nbins,-pmax/20.0,pmax/20.0) self.build1DHistogram("d0_err", "#sigma_{d0} [mm]",nbins,0,0.2) self.build1DHistogram("z0_err", @@ -209,6 +211,12 @@ def buildHistograms(self) : if self.doTruth: + self.build1DHistogram("truth_N_tracks", + "truth_N tracks",10,0,10) + self.build1DHistogram("truth_nHits", + "truth nHits", 15, 0,15) + self.build1DHistogram("truth_LayersHit", + "truth_LayersHit",15,0,15) self.build1DHistogram("truth_d0", "truth d0 [mm]", nbins, d0min,d0max) self.build1DHistogram("truth_z0", @@ -281,6 +289,8 @@ def buildHistograms(self) : #Efficiency plots + self.build1DHistogram("match_prob", + "reco truth match probability", nbins, 0.0,1.1) self.build1DHistogram("match_d0", "reco match d0 [mm]", nbins, d0min,d0max) self.build1DHistogram("match_z0", @@ -297,6 +307,11 @@ def buildHistograms(self) : "angle wrt beam axis", 20, 0, 2) self.build1DHistogram("match_PID", "Particles",8,-4,4) + self.build1DHistogram("match_nHits", + "match nHits",15,0,15) + self.build1DHistogram("match_LayersHit", + "match_LayersHit",15,0,15) + self.build1DHistogram("match_kminus_p", @@ -335,6 +350,8 @@ def buildHistograms(self) : "fake qOverP [GeV^{-1}]",nbins,-40,40) self.build1DHistogram("fake_nHits", "fake nHits",15,0,15) + self.build1DHistogram("fake_LayersHit", + "fake_LayersHit",15,0,15) self.build1DHistogram("fake_Chi2", "fake Chi2",100,0,chi2Fake_max) self.build1DHistogram("fake_Chi2_per_ndf", @@ -358,9 +375,11 @@ def buildHistograms(self) : self.build1DHistogram("dup_p", "dup p [GeV]",100,0,pmax) self.build1DHistogram("dup_qop", - "dup qOverP [GeV^{-1}]",nbins,-20,20) + "dup qOverP [GeV^{-1}]",nbins,qopmin,qopmax) self.build1DHistogram("dup_nHits", "dup nHits",15,0,15) + self.build1DHistogram("dup_LayersHit", + "dup_LayersHit",15,0,15) self.build1DHistogram("dup_Chi2", "dup Chi2",100,0,100) self.build1DHistogram("dup_Chi2_per_ndf", diff --git a/Tracking/python/geo.py b/Tracking/python/geo.py index 3c3e2dd2e..e3545c85b 100644 --- a/Tracking/python/geo.py +++ b/Tracking/python/geo.py @@ -38,6 +38,7 @@ def setDetector(self, det_name): """ from LDMX.Detectors import makePath as mP + print("Setting detector for tracking to "+det_name) self.detector = mP.makeDetectorPath( det_name ) def __init__(self): @@ -46,7 +47,7 @@ def __init__(self): else: super().__init__('TrackersTrackingGeometry', 'tracking::geo::TrackersTrackingGeometryProvider', 'Tracking') self.debug = False - self.setDetector('ldmx-det-v14') + self.setDetector('ldmx-det-v14-8gev-no-cals') TrackersTrackingGeometryProvider.__instance = self TrackersTrackingGeometryProvider.get_instance() diff --git a/Tracking/python/tracking.py b/Tracking/python/tracking.py index 26ed2f275..f5ddf711c 100644 --- a/Tracking/python/tracking.py +++ b/Tracking/python/tracking.py @@ -231,10 +231,13 @@ def __init__(self, instance_name='GSFProcessor'): self.abortOnError = False self.disableAllMaterialHandling = False self.weightCutoff = 1.0e-4 + self.debug = False self.propagator_step_size = 200. self.propagator_maxSteps = 1000 self.field_map = makeFieldMapPath() + self.taggerTracking = True + self.out_trk_collection = "GSFTracks" diff --git a/Tracking/src/Tracking/Reco/AlignmentTestProcessor.cxx b/Tracking/src/Tracking/Reco/AlignmentTestProcessor.cxx index 4aff256d4..8475e361c 100644 --- a/Tracking/src/Tracking/Reco/AlignmentTestProcessor.cxx +++ b/Tracking/src/Tracking/Reco/AlignmentTestProcessor.cxx @@ -52,7 +52,8 @@ void AlignmentTestProcessor::produce(framework::Event& event) { std::cout << entry.first << std::endl; std::cout << "Dumping surfaces information" << std::endl; - (entry.second)->toStream(align_gctx, std::cout); + // (entry.second)->toStream(align_gctx, std::cout); + (entry.second)->toStream(align_gctx); } } diff --git a/Tracking/src/Tracking/Reco/CKFProcessor.cxx b/Tracking/src/Tracking/Reco/CKFProcessor.cxx index 5941f1c35..d570943d5 100644 --- a/Tracking/src/Tracking/Reco/CKFProcessor.cxx +++ b/Tracking/src/Tracking/Reco/CKFProcessor.cxx @@ -1,6 +1,7 @@ #include "Tracking/Reco/CKFProcessor.h" -#include "Acts/EventData/TrackHelpers.hpp" +#include "Acts/EventData/TrackContainer.hpp" +#include "Acts/Utilities/TrackHelpers.hpp" #include "SimCore/Event/SimParticle.h" #include "Tracking/Reco/TruthMatchingTool.h" #include "Tracking/Sim/GeometryContainers.h" @@ -8,7 +9,7 @@ //--- C++ StdLib ---// #include //std::vector reverse #include - +#include // eN files #include @@ -34,6 +35,28 @@ void CKFProcessor::onNewRun(const ldmx::RunHeader& rh) { // Setup a constant magnetic field const auto constBField = std::make_shared(b_field); + // Define the target surface - be careful: + // x - downstream + // y - left (when looking along x) + // z - up + // Passing identity here means that your target surface is oriented in the + // same way + surf_rotation = Acts::RotationMatrix3::Zero(); + // u direction along +Y + surf_rotation(1, 0) = 1; + // v direction along +Z + surf_rotation(2, 1) = 1; + // w direction along +X + surf_rotation(0, 2) = 1; + + Acts::Vector3 target_pos(0., 0., 0.); + Acts::Translation3 target_translation(target_pos); + Acts::Transform3 target_transform(target_translation * surf_rotation); + + // Unbounded surface + target_surface = + Acts::Surface::makeShared(target_transform); + // Custom transformation of the interpolated bfield map bool debugTransform = false; auto transformPos = [this, debugTransform](const Acts::Vector3& pos) { @@ -99,7 +122,7 @@ void CKFProcessor::onNewRun(const ldmx::RunHeader& rh) { transformPos, transformBField)); auto acts_loggingLevel = Acts::Logging::FATAL; - if (debug_) acts_loggingLevel = Acts::Logging::VERBOSE; + if (debug_acts_) acts_loggingLevel = Acts::Logging::VERBOSE; // Setup the steppers const auto stepper = Acts::EigenStepper<>{map}; @@ -111,7 +134,6 @@ void CKFProcessor::onNewRun(const ldmx::RunHeader& rh) { navCfg.resolveMaterial = true; navCfg.resolvePassive = true; navCfg.resolveSensitive = true; - navCfg.boundaryCheckLayerResolving = false; const Acts::Navigator navigator(navCfg); // Setup the propagators @@ -122,8 +144,6 @@ void CKFProcessor::onNewRun(const ldmx::RunHeader& rh) { stepper, navigator, Acts::getDefaultLogger("ACTS_PROP", acts_loggingLevel)); - // auto gsf_propagator = GsfPropagator(multi_stepper, navigator); - // Setup the finder / fitters ckf_ = std::make_unique>( *propagator_, Acts::getDefaultLogger("CKF", acts_loggingLevel)); @@ -150,11 +170,11 @@ void CKFProcessor::produce(framework::Event& event) { Acts::getDefaultLogger("LDMX Tracking Geometry Maker", loggingLevel)); // Move this at the start of the producer - Acts::PropagatorOptions propagator_options( - geometry_context(), magnetic_field_context()); + Acts::PropagatorOptions + propagator_options(geometry_context(), magnetic_field_context()); propagator_options.pathLimit = std::numeric_limits::max(); - // Activate loop protection at some pt value propagator_options.loopProtection = false; //(startParameters.transverseMomentum() < cfg.ptLoopers); @@ -171,13 +191,10 @@ void CKFProcessor::produce(framework::Event& event) { propagator_options.actionList.get(); sLogger.sterile = true; // Set a maximum step size - propagator_options.maxStepSize = + propagator_options.stepping.maxStepSize = propagator_step_size_ * Acts::UnitConstants::mm; propagator_options.maxSteps = propagator_maxSteps_; - // Electron hypothesis - propagator_options.mass = 0.511 * Acts::UnitConstants::MeV; - // #######################// // Kalman Filter algorithm// // #######################// @@ -221,6 +238,8 @@ void CKFProcessor::produce(framework::Event& event) { const std::vector seed_tracks = event.getCollection(seed_coll_name_); + ldmx_log(debug) << "Number of seeds::" << seed_tracks.size(); + // Run the CKF on each seed and produce a track candidate std::vector startParameters; @@ -229,20 +248,15 @@ void CKFProcessor::produce(framework::Event& event) { for (auto& seed : seed_tracks) { // Transform the seed track to bound parameters - // MG ... make these into pararamers at target surface std::shared_ptr perigeeSurface = Acts::Surface::makeShared(Acts::Vector3( seed.getPerigeeX(), seed.getPerigeeY(), seed.getPerigeeZ())); - // MG ... make these into pararamers at target surface - // MG ... do I need to extrapolate from perigee to target z? Acts::BoundVector paramVec; paramVec << seed.getD0(), seed.getZ0(), seed.getPhi(), seed.getTheta(), seed.getQoP(), seed.getT(); - // MG ... make these into bound pararamers at target surface - // MG ... do I need to extrapolate from perigee to target z? - Acts::BoundSymMatrix covMat = + Acts::BoundSquareMatrix covMat = tracking::sim::utils::unpackCov(seed.getPerigeeCov()); ldmx_log(debug) << "perigee" << std::endl @@ -253,11 +267,10 @@ void CKFProcessor::produce(framework::Event& event) { ldmx_log(debug) << "cov matrix" << std::endl << covMat << std::endl; - Acts::ActsScalar q = seed.getQoP() < 0 ? -1 * Acts::UnitConstants::e - : Acts::UnitConstants::e; - // MG ... if the above is at target surface, this should just work + // need to set particle hypothesis...set to electron for now... + auto partHypo{Acts::SinglyChargedParticleHypothesis::electron()}; startParameters.push_back( - Acts::BoundTrackParameters(perigeeSurface, paramVec, q, covMat)); + Acts::BoundTrackParameters(perigeeSurface, paramVec, covMat, partHypo)); seedPDGID.push_back(seed.getPdgID()); @@ -275,7 +288,6 @@ void CKFProcessor::produce(framework::Event& event) { std::chrono::duration(seeds - hits).count(); Acts::GainMatrixUpdater kfUpdater; - Acts::GainMatrixSmoother kfSmoother; // configuration for the measurement selector. Empty geometry identifier means // applicable to all the detector elements @@ -289,25 +301,21 @@ void CKFProcessor::produce(framework::Event& event) { tracking::sim::LdmxMeasurementCalibrator calibrator{measurements}; - Acts::CombinatorialKalmanFilterExtensions - ckf_extensions; + Acts::CombinatorialKalmanFilterExtensions ckf_extensions; if (use1Dmeasurements_) ckf_extensions.calibrator - .connect<&tracking::sim::LdmxMeasurementCalibrator::calibrate_1d>( - &calibrator); + .connect<&tracking::sim::LdmxMeasurementCalibrator::calibrate_1d< + Acts::VectorMultiTrajectory>>(&calibrator); else ckf_extensions.calibrator - .connect<&tracking::sim::LdmxMeasurementCalibrator::calibrate>( - &calibrator); + .connect<&tracking::sim::LdmxMeasurementCalibrator::calibrate< + Acts::VectorMultiTrajectory>>(&calibrator); ckf_extensions.updater.connect< &Acts::GainMatrixUpdater::operator()>( &kfUpdater); - ckf_extensions.smoother.connect< - &Acts::GainMatrixSmoother::operator()>( - &kfSmoother); ckf_extensions.measurementSelector .connect<&Acts::MeasurementSelector::select>( @@ -356,30 +364,9 @@ void CKFProcessor::produce(framework::Event& event) { ldmx_log(debug) << "Surfaces..." << std::endl; - Acts::RotationMatrix3 target_surf_rotation = Acts::RotationMatrix3::Zero(); - // u direction along +Y - target_surf_rotation(1, 0) = 1; - // v direction along +Z - target_surf_rotation(2, 1) = 1; - // w direction along +X - target_surf_rotation(0, 2) = 1; - // - Acts::Vector3 target_center_unbound(0., 0., 0.); - Acts::Translation3 target_center_translation(target_center_unbound); - Acts::Transform3 target_center_transform(target_center_translation * - target_surf_rotation); - - // Unbounded surface...different from above tgt_surface is that this is a - // plane and not perigee surface (cylinder) - const std::shared_ptr tgt_surface_unbound = - Acts::Surface::makeShared(target_center_transform); - - auto extr_surface = &(*tgt_surface_unbound); - - // std::shared_ptr origin_surface = - // Acts::Surface::makeShared( - // Acts::Vector3(0., 0., 0.)); - // auto extr_surface = &(*origin_surface); + std::shared_ptr origin_surface = + Acts::Surface::makeShared( + Acts::Vector3(0., 0., 0.)); ldmx_log(debug) << "About to run CKF..." << std::endl; @@ -399,212 +386,223 @@ void CKFProcessor::produce(framework::Event& event) { for (size_t trackId = 0u; trackId < startParameters.size(); ++trackId) { // The seed has a track PdgID associated if (seedPDGID.at(trackId) != 0) { - int pdgID = seedPDGID.at(trackId); - - if (pdgID == 2212 || pdgID == -2212) - propagator_options.mass = 938 * Acts::UnitConstants::MeV; + // int pdgID = seedPDGID.at(trackId); } // Define the CKF options here: const Acts::CombinatorialKalmanFilterOptions - ckfOptions(geometry_context(), magnetic_field_context(), - calibration_context(), sourceLinkAccessorDelegate, - ckf_extensions, propagator_options, &(*extr_surface)); + TrackContainer> + ckfOptions(TrackingGeometryUser::geometry_context(), + TrackingGeometryUser::magnetic_field_context(), + TrackingGeometryUser::calibration_context(), + sourceLinkAccessorDelegate, ckf_extensions, + propagator_options, true /* multiple scattering */, + false /* energy loss */); ldmx_log(debug) << "Running CKF on seed params " << startParameters.at(trackId).parameters().transpose() << std::endl; - + ldmx_log(debug) << "Checking options: multiple scattering = " + << ckfOptions.multipleScattering + << " energy loss = " << ckfOptions.energyLoss; auto results = ckf_->findTracks(startParameters.at(trackId), ckfOptions, tc); - + ldmx_log(debug) << "findTracks returned ... checking if ok"; if (not results.ok()) { ldmx_log(debug) << "CKF Fit failed" << std::endl; continue; } // No track found - if (tc.size() < trackId + 1) continue; - - ldmx_log(debug) << "Filling track info" << std::endl; - - // The track tips are the last measurement index - // Acts::MultiTrajectory mj = - // tc.getTrack(trackId); - // //.container() - // //.trackStateContainer(); - - auto track = tc.getTrack(trackId); - calculateTrackQuantities(track); - // MG ... if I converted above to target surface, these should be parameters - // at target (should change names) - const Acts::BoundVector& perigee_pars = track.parameters(); - const Acts::BoundMatrix& trk_cov = track.covariance(); - const Acts::Surface& perigee_surface = track.referenceSurface(); - - ldmx_log(debug) - << "Found track: nMeas " << track.nMeasurements() << std::endl - << "Track states " << track.nTrackStates() << std::endl - << perigee_pars[Acts::eBoundLoc0] << " " - << perigee_pars[Acts::eBoundLoc1] << " " - << perigee_pars[Acts::eBoundPhi] << " " - << perigee_pars[Acts::eBoundTheta] << " " - << perigee_pars[Acts::eBoundQOverP] << std::endl - << "Reference Surface" << std::endl - << " " << perigee_surface.transform(geometry_context()).translation()(0) - << " " << perigee_surface.transform(geometry_context()).translation()(1) - << " " << perigee_surface.transform(geometry_context()).translation()(2) - << std::endl - << "nHoles " << track.nHoles(); - - ldmx::Track trk = ldmx::Track(); - trk.setPerigeeLocation( - perigee_surface.transform(geometry_context()).translation()(0), - perigee_surface.transform(geometry_context()).translation()(1), - perigee_surface.transform(geometry_context()).translation()(2)); - - trk.setChi2(track.chi2()); - trk.setNhits(track.nMeasurements()); - // trk.setNdf(track.nDoF()); - // TODO Switch back to nDoF when Acts is fixed. - trk.setNdf(track.nMeasurements() - 5); - trk.setNsharedHits(track.nSharedHits()); - - trk.setPerigeeParameters( - tracking::sim::utils::convertActsToLdmxPars(perigee_pars)); - std::vector v_trk_cov; - tracking::sim::utils::flatCov(trk_cov, v_trk_cov); - trk.setPerigeeCov(v_trk_cov); - - Acts::Vector3 trk_momentum = track.momentum(); - trk.setMomentum(trk_momentum(0), trk_momentum(1), trk_momentum(2)); - - // Add measurements on track - for (auto ts : track.trackStates()) { - // Check TrackStates Quality - ldmx_log(debug) << "Checking Track State at location " - << ts.referenceSurface() - .transform(geometry_context()) - .translation() - .transpose() - << std::endl; - - ldmx_log(debug) << "Smoothed? \n" << ts.hasSmoothed() << std::endl; - if (ts.hasSmoothed()) { - ldmx_log(debug) << "Parameters \n" - << ts.smoothed().transpose() << std::endl; - ldmx_log(debug) << "Covariance \n" - << ts.smoothedCovariance() << std::endl; + // if (tc.size() < trackId + 1) continue; + + auto& tracksFromSeed = results.value(); + + ldmx_log(debug) << "number of entries in results " << tracksFromSeed.size(); + for (auto& track : tracksFromSeed) { + // do the track smoothing...this is not done in the CKF code anymore + Acts::smoothTrack(geometry_context(), track); // from TrackHelpers + // make the empty ldmx::Track() and track state at target + ldmx::Track trk = ldmx::Track(); + ldmx::Track::TrackState tsAtTarget; + ldmx_log(debug) << "Found track: nMeas " << track.nMeasurements(); + ldmx_log(debug) << "Track states " << track.nTrackStates(); + ldmx_log(debug) << "chi2 " << track.chi2(); + + for (const auto ts : track.trackStatesReversed()) { + // Check TrackStates Quality + ldmx_log(debug) << "Checking Track State at location " + << ts.referenceSurface() + .transform(geometry_context()) + .translation() + .transpose() + << std::endl; + + ldmx_log(debug) << "Smoothed? " << ts.hasSmoothed() << std::endl; + if (ts.hasSmoothed()) { + ldmx_log(debug) << "Parameters \n" + << ts.smoothed().transpose() << std::endl; + ldmx_log(debug) << "Covariance \n" + << ts.smoothedCovariance() << std::endl; + } + + // Check if the track state is a measurement + auto typeFlags = ts.typeFlags(); + + if (typeFlags.test(Acts::TrackStateFlag::MeasurementFlag) && + ts.hasUncalibratedSourceLink()) { + ldmx_log(debug) << " getting source link for this measurement"; + + const ActsExamples::IndexSourceLink sl = + ts.getUncalibratedSourceLink() + .template get(); + + ldmx_log(debug) << " looking up this index in measurements list"; + ldmx::Measurement ldmx_meas = measurements.at(sl.index()); + ldmx_log(debug) << "SourceLink Index::" << sl.index(); + ldmx_log(debug) << "Measurement:\n" << ldmx_meas << "\n"; + ldmx_log(debug) << " adding measurement to ldmx::track"; + trk.addMeasurementIndex(sl.index()); + } } - - // Check if the track state is a measurement - auto typeFlags = ts.typeFlags(); - if (typeFlags.test(Acts::TrackStateFlag::MeasurementFlag)) { - ActsExamples::IndexSourceLink sl = - ts.getUncalibratedSourceLink().get(); - ldmx::Measurement ldmx_meas = measurements.at(sl.index()); - ldmx_log(debug) << "SourceLink Index::" << sl.index(); - ldmx_log(debug) << "Measurement:\n" << ldmx_meas << "\n"; - trk.addMeasurementIndex(sl.index()); + bool success = trk_extrap_->TrackStateAtSurface( + track, target_surface, tsAtTarget, ldmx::TrackStateType::AtTarget); + ldmx_log(debug) << "target extrapolation success??? " << success; + if (success) { + ldmx_log(debug) << "Successfully obtained TS at target"; + ldmx_log(debug) << "Parameters At Target: \n" + << tsAtTarget.params[0] << " " << tsAtTarget.params[1] + << " " << tsAtTarget.params[2] << " " + << tsAtTarget.params[3] << " " << tsAtTarget.params[4]; + + trk.addTrackState(tsAtTarget); + } else { + ldmx_log(info) + << "Could not extrapolate to target? Printing track states: "; + ldmx_log(info) << " nhits = " << track.nMeasurements(); + for (const auto ts : track.trackStatesReversed()) { + ldmx_log(info) << "Smoothed? " << ts.hasSmoothed() << std::endl; + if (ts.hasSmoothed()) { + ldmx_log(info) << "momentum for track state = " + << 1 / ts.smoothed()[Acts::eBoundQOverP]; + ldmx_log(info) << "Parameters \n" + << ts.smoothed().transpose() << std::endl; + } else { + ldmx_log(info) << "Track state not smoothed?"; + } + } + ldmx_log(info) << "...skipping this track..."; + continue; } - } - // Extrapolations - - // Define the target surface - be careful: - // x - downstream - // y - left (when looking along x) - // z - up - // Passing identity here means that your target surface is oriented in the - // same way - Acts::RotationMatrix3 surf_rotation = Acts::RotationMatrix3::Zero(); - // u direction along +Y - surf_rotation(1, 0) = 1; - // v direction along +Z - surf_rotation(2, 1) = 1; - // w direction along +X - surf_rotation(0, 2) = 1; - - const double ECAL_SCORING_PLANE = 240.5; - Acts::Vector3 pos(ECAL_SCORING_PLANE, 0., 0.); - Acts::Translation3 surf_translation(pos); - Acts::Transform3 surf_transform(surf_translation * surf_rotation); - - // Unbounded surface - const std::shared_ptr ecal_surface = - Acts::Surface::makeShared(surf_transform); - - Acts::Vector3 target_pos(0., 0., 0.); - Acts::Translation3 target_translation(target_pos); - Acts::Transform3 target_transform(target_translation * surf_rotation); - - // Unbounded surface - const std::shared_ptr target_surface = - Acts::Surface::makeShared(target_transform); - - // Beam Origin unbounded surface - const std::shared_ptr beamOrigin_surface = - tracking::sim::utils::unboundSurface(-700); - - ldmx_log(debug) << "Target extrapolation ... this should not change " - "anything since track is stored at target plane"; - ldmx::Track::TrackState tsAtTarget; - bool successAtTarget = trk_extrap_->TrackStateAtSurface( - track, target_surface, tsAtTarget, ldmx::TrackStateType::AtTarget); - if (successAtTarget) { - ldmx_log(debug) << "Successfully obtained TS at target"; - ldmx_log(debug) << "Parameters At Target: \n" - << tsAtTarget.params[0] << " " << tsAtTarget.params[1] - << " " << tsAtTarget.params[2] << " " - << tsAtTarget.params[3] << " " << tsAtTarget.params[4]; - - trk.addTrackState(tsAtTarget); - } - if (taggerTracking_) { - ldmx_log(debug) << "Beam Origin Extrapolation"; - ldmx::Track::TrackState tsAtBeamOrigin; - bool successAtBeam = trk_extrap_->TrackStateAtSurface( - track, beamOrigin_surface, tsAtBeamOrigin, - ldmx::TrackStateType::AtBeamOrigin); - - if (successAtBeam) { - trk.addTrackState(tsAtBeamOrigin); - ldmx_log(debug) << "Successfully obtained TS at beam origin"; + // get the BoundTrackParameters at the target + // ...use to fill in the Acts::TrackProxy object + // This isn't really necessary, since we can take + // most everything for making the ldmx::track + // from tsAtTarget...maybe useful for something? + // -->one thing this does is allow Acts to + // calculate the momentum 3-vector for you + Acts::BoundTrackParameters boundStateAtTarget = + tracking::sim::utils::btp(tsAtTarget, target_surface, 11); + track.setReferenceSurface(target_surface); + track.parameters() = boundStateAtTarget.parameters(); + + ldmx_log(debug) << typeid(track).name(); + // These are the parameters at the target surface + const Acts::BoundVector& track_pars = track.parameters(); + // const Acts::BoundMatrix& trk_cov = track.covariance(); + const Acts::Surface& track_surface = track.referenceSurface(); + ldmx_log(debug) << "Got the parameters, covariance, and perigee surface"; + + ldmx_log(debug) << track_pars[Acts::eBoundLoc0]; + ldmx_log(debug) << track_pars[Acts::eBoundLoc1]; + ldmx_log(debug) << track_pars[Acts::eBoundTheta]; + ldmx_log(debug) << track_pars[Acts::eBoundPhi]; + ldmx_log(debug) + << "Reference Surface" << std::endl + << " " << track_surface.transform(geometry_context()).translation()(0) + << " " << track_surface.transform(geometry_context()).translation()(1) + << " " + << track_surface.transform(geometry_context()).translation()(2); + + trk.setPerigeeLocation( + 0, 0, 0); // the target...it's not really perigee anymore. + trk.setPerigeeParameters(tsAtTarget.params); + trk.setPerigeeCov(tsAtTarget.cov); + + ldmx_log(debug) << "setting chi2 and nHits: " << track.chi2() << " " + << track.nMeasurements(); + trk.setChi2(track.chi2()); + trk.setNhits(track.nMeasurements()); + // trk.setNdf(track.nDoF()); + // TODO Switch back to nDoF when Acts is fixed. + trk.setNdf(track.nMeasurements() - 5); + trk.setNsharedHits(track.nSharedHits()); + + ldmx_log(debug) << "setting track momentum: " << track.momentum(); + trk.setMomentum(track.momentum()[0], track.momentum()[1], + track.momentum()[2]); + + ldmx_log(debug) << "starting extrapolations"; + // Extrapolations + + const double ECAL_SCORING_PLANE = 240.5; + Acts::Vector3 pos(ECAL_SCORING_PLANE, 0., 0.); + Acts::Translation3 surf_translation(pos); + Acts::Transform3 surf_transform(surf_translation * surf_rotation); + + // Unbounded surface + const std::shared_ptr ecal_surface = + Acts::Surface::makeShared(surf_transform); + + // Beam Origin unbounded surface + const std::shared_ptr beamOrigin_surface = + tracking::sim::utils::unboundSurface(-700); + + if (taggerTracking_) { + ldmx_log(debug) << "Beam Origin Extrapolation"; + ldmx::Track::TrackState tsAtBeamOrigin; + success = trk_extrap_->TrackStateAtSurface( + track, beamOrigin_surface, tsAtBeamOrigin, + ldmx::TrackStateType::AtBeamOrigin); + + if (success) { + trk.addTrackState(tsAtBeamOrigin); + ldmx_log(debug) << "Successfully obtained TS at beam origin"; + } } - } - // Recoil Extrapolation to ECAL only - if (!taggerTracking_) { - ldmx_log(debug) << "Ecal Extrapolation"; - ldmx::Track::TrackState tsAtEcal; - bool successAtEcal = trk_extrap_->TrackStateAtSurface( - track, ecal_surface, tsAtEcal, ldmx::TrackStateType::AtECAL); - - if (successAtEcal) { - trk.addTrackState(tsAtEcal); - ldmx_log(debug) << "Successfully obtained TS at Ecal"; - ldmx_log(debug) << "Parameters At Ecal: \n" - << tsAtEcal.params[0] << " " << tsAtEcal.params[1] - << " " << tsAtEcal.params[2] << " " - << tsAtEcal.params[3] << " " << tsAtEcal.params[4]; + // Recoil Extrapolation to ECAL only + if (!taggerTracking_) { + ldmx_log(debug) << "Ecal Extrapolation"; + ldmx::Track::TrackState tsAtEcal; + success = trk_extrap_->TrackStateAtSurface( + track, ecal_surface, tsAtEcal, ldmx::TrackStateType::AtECAL); + + if (success) { + trk.addTrackState(tsAtEcal); + ldmx_log(debug) << "Successfully obtained TS at Ecal"; + ldmx_log(debug) << "Parameters At Ecal: \n" + << tsAtEcal.params[0] << " " << tsAtEcal.params[1] + << " " << tsAtEcal.params[2] << " " + << tsAtEcal.params[3] << " " << tsAtEcal.params[4]; + } } - } - // Truth matching - if (truthMatchingTool) { - auto truthInfo = truthMatchingTool->TruthMatch(trk); - trk.setTrackID(truthInfo.trackID); - trk.setPdgID(truthInfo.pdgID); - trk.setTruthProb(truthInfo.truthProb); - } + // Truth matching + if (truthMatchingTool) { + auto truthInfo = truthMatchingTool->TruthMatch(trk); + trk.setTrackID(truthInfo.trackID); + trk.setPdgID(truthInfo.pdgID); + trk.setTruthProb(truthInfo.truthProb); + } - // At least 8 hits and p > 50 MeV - if (trk.getNhits() > min_hits_ && abs(1. / trk.getQoP()) > 0.05) { - tracks.push_back(trk); - ntracks_++; + // At least min_hits_ hits and p > 50 MeV + if (trk.getNhits() > min_hits_ && abs(1. / trk.getQoP()) > 0.05) { + tracks.push_back(trk); + ntracks_++; + } } - } // loop seed track parameters auto result_loop = std::chrono::high_resolution_clock::now(); @@ -670,6 +668,8 @@ void CKFProcessor::configure(framework::config::Parameters& parameters) { "measurement_collection", "TaggerMeasurements"); outlier_pval_ = parameters.getParameter("outlier_pval_", 3.84); + debug_acts_ = parameters.getParameter("debug_acts", false); + remove_stereo_ = parameters.getParameter("remove_stereo", false); use1Dmeasurements_ = parameters.getParameter("use1Dmeasurements", true); min_hits_ = parameters.getParameter("min_hits", 7); @@ -710,7 +710,7 @@ auto CKFProcessor::makeGeoIdSourceLinkMap( ActsExamples::IndexSourceLink> geoId_sl_map; - ldmx_log(debug) << "makeGeoIdSourceLinkMap::Available measurements" + ldmx_log(debug) << "makeGeoIdSourceLinkMap::Available measurements = " << measurements.size(); // Check the hits associated to the surfaces @@ -725,7 +725,9 @@ auto CKFProcessor::makeGeoIdSourceLinkMap( // information ActsExamples::IndexSourceLink idx_sl(hit_surface->geometryId(), i_meas); - + // mg aug 2024 ... these don't print statements + // don't compile using v36 in Acts...figure out later + /* ldmx_log(debug) << "Insert measurement on surface located at::" << hit_surface->transform(geometry_context()).translation(); @@ -734,7 +736,7 @@ auto CKFProcessor::makeGeoIdSourceLinkMap( ldmx_log(debug) << "Surface info::" << std::tie(*hit_surface, geometry_context()); - + */ geoId_sl_map.insert(std::make_pair(hit_surface->geometryId(), idx_sl)); } else diff --git a/Tracking/src/Tracking/Reco/CustomStatePropagator.cxx b/Tracking/src/Tracking/Reco/CustomStatePropagator.cxx deleted file mode 100644 index f1bab53bf..000000000 --- a/Tracking/src/Tracking/Reco/CustomStatePropagator.cxx +++ /dev/null @@ -1,271 +0,0 @@ -#include "Tracking/Reco/CustomStatePropagator.h" - -#include -#include - -#include "Acts/Utilities/Logger.hpp" - -namespace tracking { -namespace reco { - -CustomStatePropagator::CustomStatePropagator(const std::string& name, - framework::Process& process) - : framework::Producer(name, process) {} - -CustomStatePropagator::~CustomStatePropagator() {} - -void CustomStatePropagator::onProcessStart() { - outFile_ = new TFile("prop_states.root", "RECREATE"); - outTree_ = new TTree("prop_states", "prop_states"); - histo_gen_p = - std::make_shared("histo_gen_p", "histo_gen_p", 100, 0, 10); - histo_end_p = std::make_shared("histo_end_p", "histo_end_p", 100, 0, - 10); // smart pointer - histo_end_px = - std::make_shared("histo_end_px", "histo_end_px", 100, -50, 50); - histo_end_py = - std::make_shared("histo_end_py", "histo_end_py", 100, -50, 50); - histo_end_pz = - std::make_shared("histo_end_pz", "histo_end_pz", 100, -50, 50); - histo_gen_px = - std::make_shared("histo_gen_px", "histo_gen_px", 100, -50, 50); - histo_gen_py = - std::make_shared("histo_gen_py", "histo_gen_py", 100, -50, 50); - histo_gen_pz = - std::make_shared("histo_gen_pz", "histo_gen_pz", 100, -50, 50); - histo_loc01 = std::make_shared("histo_loc01", "end_loc0-vs-end_loc1", - 100, -50, 50, 100, -50, 50); - - outTree_->Branch("state_nr", &state_nr); - outTree_->Branch("charge", &charge_); - outTree_->Branch("gen_x", &gen_x); - outTree_->Branch("gen_y", &gen_y); - outTree_->Branch("gen_z", &gen_z); - outTree_->Branch("gen_px", &gen_px); - outTree_->Branch("gen_py", &gen_py); - outTree_->Branch("gen_pz", &gen_pz); - - outTree_->Branch("end_x", &end_x); - outTree_->Branch("end_y", &end_y); - outTree_->Branch("end_z", &end_z); - outTree_->Branch("end_loc0", &end_loc0); - outTree_->Branch("end_loc1", &end_loc1); - - outTree_->Branch("end_px", &end_px); - outTree_->Branch("end_py", &end_py); - outTree_->Branch("end_pz", &end_pz); - - // Setup a interpolated bfield map - const auto map = - std::make_shared(loadDefaultBField( - field_map_, default_transformPos, default_transformBField)); - - const auto stepper = Acts::EigenStepper<>{map}; - - // Set up propagator with void navigator (No material) - auto propagator = - std::make_shared>>(stepper); - - Acts::Vector3 gen_pos{0., 0., 0.}; - Acts::Vector3 gen_mom{0., 0., 0.}; - - std::default_random_engine generator; - generator.seed(1); - - // Beamspot - std::uniform_real_distribution bY(-bs_size_[0], bs_size_[0]); - std::uniform_real_distribution bX(-bs_size_[1], bs_size_[1]); - - // 3-momentum in cartesian coordinates - std::uniform_real_distribution PX(-4, 4); - std::uniform_real_distribution PY(-4, 4); - std::uniform_real_distribution PZ(0.0, 4); - - // 3-momentum in polar coordinates - std::uniform_real_distribution P(prange_[0], prange_[1]); - std::uniform_real_distribution THETA(thetarange_[0], thetarange_[1]); - std::uniform_real_distribution PHI(phirange_[0], phirange_[1]); - - std::uniform_real_distribution CHARGE(-1, 1); - - for (int i_state = 0; i_state < nstates_; i_state++) { - double p = P(generator); - double theta = THETA(generator); - double phi = PHI(generator); - int charge = CHARGE(generator) > 0 ? 1 : -1; - - double px = p * cos(theta); - double py = p * sin(theta) * cos(phi); - double pz = p * sin(theta) * sin(phi); - - double by = bY(generator); - double bx = bX(generator); - - gen_pos(0) = 0.; - gen_pos(1) = bx; - gen_pos(2) = by; - - // Transform to MeV because that's what TrackUtils assumes - gen_mom(0) = px / Acts::UnitConstants::MeV; - gen_mom(1) = py / Acts::UnitConstants::MeV; - gen_mom(2) = pz / Acts::UnitConstants::MeV; - - // std::cout<<"Generated position"< gen_surface = - Acts::Surface::makeShared(Acts::Vector3( - part_free[Acts::eFreePos0], part_free[Acts::eFreePos1], - part_free[Acts::eFreePos2])); - - // Transform the free parameters to the bound parameters - auto bound_params = Acts::detail::transformFreeToBoundParameters( - part_free, *gen_surface, gctx_) - .value(); - - Acts::BoundTrackParameters startParams(gen_surface, std::move(bound_params), - std::move(std::nullopt)); - - // const auto pLogger = Acts::getDefaultLogger("Propagator", - // Acts::Logging::INFO); - Acts::PropagatorOptions<> propagator_options( - gctx_, bctx_); //, Acts::LoggerWrapper{*pLogger}); - - // propagator_options.direction = Acts::Direction::Forward; // should be the - // default - propagator_options.pathLimit = std::numeric_limits::max(); - propagator_options.loopProtection = true; - propagator_options.maxStepSize = 1 * Acts::UnitConstants::mm; - propagator_options.maxSteps = 2000; - - // Define the target surface - be careful: - // x - downstream - // y - left (when looking along x) - // z - up - // Passing identity here means that your target surface is oriented in the - // same way - Acts::RotationMatrix3 surf_rotation = Acts::RotationMatrix3::Zero(); - // u direction along +Y - surf_rotation(1, 0) = 1; - // v direction along +Z - surf_rotation(2, 1) = 1; - // w direction along +X - surf_rotation(0, 2) = 1; - - Acts::Vector3 pos(surf_location_, 0., 0.); - Acts::Translation3 surf_translation(pos); - Acts::Transform3 surf_transform(surf_translation * surf_rotation); - - // Unbounded surface - const std::shared_ptr target_surface = - Acts::Surface::makeShared(surf_transform); - - // Do the propagation to the surface - - auto result = - propagator->propagate(startParams, *target_surface, propagator_options); - - if (not result.ok()) { - return; - } - - const auto& endParams = *result->endParameters; - - fillTree(i_state, q, gen_pos, gen_mom, endParams); - - // loc0 // loc1 will give you the u-v location of the hit on the ecal face - - } // state propagation - -} // on Process Start - -void CustomStatePropagator::configure( - framework::config::Parameters& parameters) { - const double PIo2 = 1.57079632679; - - field_map_ = parameters.getParameter("field_map", ""); - surf_location_ = parameters.getParameter("surf_location", 350.); - nstates_ = parameters.getParameter("nstates", 10); - bs_size_ = - parameters.getParameter>("bs_size", {40., 10.}); - prange_ = parameters.getParameter>("prange", {0.05, 4.}); - thetarange_ = - parameters.getParameter>("thetarange", {0, PIo2}); - phirange_ = - parameters.getParameter>("phirange", {0, PIo2 * 4.}); -} - -void CustomStatePropagator::fillTree( - int state, int q, const Acts::Vector3 gen_pos, const Acts::Vector3 gen_mom, - const Acts::BoundTrackParameters& endParams) { - state_nr = state; - charge_ = q; - gen_x = gen_pos(0); - gen_y = gen_pos(1); - gen_z = gen_pos(2); - - gen_px = gen_mom(0); - gen_py = gen_mom(1); - gen_pz = gen_mom(2); - - Acts::Vector3 end_pos = endParams.position(gctx_); - end_x = end_pos(0); - end_y = end_pos(1); - end_z = end_pos(2); - - Acts::BoundVector bound_parameters = endParams.parameters(); - end_loc0 = bound_parameters[Acts::eBoundLoc0]; - end_loc1 = bound_parameters[Acts::eBoundLoc1]; - - Acts::Vector3 end_mom = endParams.momentum(); - end_px = end_mom(0); - end_py = end_mom(1); - end_pz = end_mom(2); - - // Calculate magnatitude of generating & end momentum - const double gen_mag_p = - sqrt(pow(gen_px, 2) + pow(gen_py, 2) + pow(gen_pz, 2)); - const double end_mag_p = - sqrt(pow(end_px, 2) + pow(end_py, 2) + pow(end_pz, 2)); - outTree_->Fill(); - - histo_end_px->Fill(end_px); - histo_end_py->Fill(end_py); - histo_end_pz->Fill(end_pz); - histo_gen_px->Fill(gen_px / 1000); - histo_gen_py->Fill(gen_py / 1000); - histo_gen_pz->Fill(gen_pz / 1000); - histo_gen_p->Fill(gen_mag_p / 1000); - histo_end_p->Fill(end_mag_p); - histo_loc01->Fill(end_loc0, end_loc1); -} - -void CustomStatePropagator::onProcessEnd() { - outFile_->cd(); - outTree_->Write(); - histo_end_px->Write(); - histo_end_py->Write(); - histo_end_pz->Write(); - histo_gen_px->Write(); - histo_gen_py->Write(); - histo_gen_pz->Write(); - histo_gen_p->Write(); - histo_end_p->Write(); - histo_loc01->Write(); - outFile_->Close(); - delete outFile_; -} - -} // namespace reco -} // namespace tracking - -DECLARE_PRODUCER_NS(tracking::reco, CustomStatePropagator) diff --git a/Tracking/src/Tracking/Reco/DigitizationProcessor.cxx b/Tracking/src/Tracking/Reco/DigitizationProcessor.cxx index 6dc1eac36..276f5c254 100644 --- a/Tracking/src/Tracking/Reco/DigitizationProcessor.cxx +++ b/Tracking/src/Tracking/Reco/DigitizationProcessor.cxx @@ -5,6 +5,10 @@ #include "Tracking/Event/Measurement.h" #include "Tracking/Sim/TrackingUtils.h" +// Use this instead of CartesianSegmeter I think +// BinUtility(std::size_t bins, float min, float max, BinningOption opt = open, +// BinningValue value = BinningValue::binX, +// const Transform3& tForm = Transform3::Identity()) using namespace framework; namespace tracking::reco { @@ -15,42 +19,6 @@ DigitizationProcessor::DigitizationProcessor(const std::string& name, void DigitizationProcessor::onProcessStart() { normal_ = std::make_shared>(0., 1.); - - ldmx_log(info) << "Loading the tracking geometry"; - - // Module Bounds => Take them from the tracking geometry TODO - auto moduleBounds = std::make_shared( - 20.17 * Acts::UnitConstants::mm, 50 * Acts::UnitConstants::mm); - - // I assume 5 APVs - int nbinsx = 128 * 5; - - // Strips - int nbinsy = 1; - - // Thickness = 0.320 mm - double thickness = 0.320 * Acts::UnitConstants::mm; - - // Lorentz angle - double lAngle = 0.01; - - // Energy threshold - double eThresh = 0.; - - // Analogue readout - bool isAnalog = true; - - // Cartesian segmentation - auto cSegmentation = std::make_shared( - moduleBounds, nbinsx, nbinsy); - - // Negative side readout => TODO Make sure this is correct! - // - Ask Paul what does this mean: depending on how local w is oriented - // TODO: load proper lorentz angle - - Acts::DigitizationModule ndModule(cSegmentation, thickness * 0.5, -1, lAngle, - eThresh, isAnalog); - ldmx_log(info) << "Initialization done" << std::endl; } diff --git a/Tracking/src/Tracking/Reco/GSFProcessor.cxx b/Tracking/src/Tracking/Reco/GSFProcessor.cxx index 3197af489..85775ba8b 100644 --- a/Tracking/src/Tracking/Reco/GSFProcessor.cxx +++ b/Tracking/src/Tracking/Reco/GSFProcessor.cxx @@ -81,12 +81,13 @@ void GSFProcessor::onNewRun(const ldmx::RunHeader& rh) { // Acts::MixtureReductionMethod finalReductionMethod; // const auto multi_stepper = Acts::MultiEigenStepperLoop{map}; - Acts::MixtureReductionMethod reductionMethod = - Acts::MixtureReductionMethod::eMaxWeight; - Acts::MultiEigenStepperLoop multi_stepper( - map, reductionMethod, - Acts::getDefaultLogger("GSF_STEP", acts_loggingLevel)); + Acts::ComponentMergeMethod reductionMethod = + Acts::ComponentMergeMethod::eMaxWeight; + // Acts::MultiEigenStepperLoop multi_stepper( + // map, reductionMethod, + // Acts::getDefaultLogger("GSF_STEP", acts_loggingLevel)); + Acts::MultiEigenStepperLoop multi_stepper(map); // Detailed Stepper // Acts::MultiEigenStepperLoop multi_stepper(map, finalReductionMethod); @@ -96,15 +97,13 @@ void GSFProcessor::onNewRun(const ldmx::RunHeader& rh) { navCfg.resolveMaterial = true; navCfg.resolvePassive = false; navCfg.resolveSensitive = true; - navCfg.boundaryCheckLayerResolving = false; const Acts::Navigator navigator(navCfg); auto gsf_propagator = GsfPropagator(std::move(multi_stepper), std::move(navigator), Acts::getDefaultLogger("GSF_PROP", acts_loggingLevel)); - BetheHeitlerApprox betheHeitler = - Acts::Experimental::makeDefaultBetheHeitlerApprox(); + BetheHeitlerApprox betheHeitler = Acts::makeDefaultBetheHeitlerApprox(); const auto gsfLogger = Acts::getDefaultLogger("GSF", acts_loggingLevel); @@ -124,6 +123,11 @@ void GSFProcessor::configure(framework::config::Parameters& parameters) { out_trk_collection_ = parameters.getParameter("out_trk_collection", "GSFTracks"); + trackCollection_ = + parameters.getParameter("trackCollection", "TaggerTracks"); + measCollection_ = parameters.getParameter("measCollection", + "DigiTaggerSimHits"); + maxComponents_ = parameters.getParameter("maxComponents", 4); abortOnError_ = parameters.getParameter("abortOnError", false); disableAllMaterialHandling_ = @@ -138,6 +142,7 @@ void GSFProcessor::configure(framework::config::Parameters& parameters) { usePerigee_ = parameters.getParameter("usePerigee", false); debug_ = parameters.getParameter("debug", false); + taggerTracking_ = parameters.getParameter("taggerTracking", true); // finalReductionMethod_ = // parameters.getParameter("finalReductionMethod",); @@ -159,21 +164,38 @@ void GSFProcessor::produce(framework::Event& event) { tracking::sim::LdmxMeasurementCalibrator calibrator{measurements}; // GSF Setup - Acts::GainMatrixUpdater updater; - Acts::Experimental::GsfExtensions gsf_extensions; + Acts::GsfExtensions gsf_extensions; gsf_extensions.updater.connect< &Acts::GainMatrixUpdater::operator()>( &updater); gsf_extensions.calibrator - .connect<&tracking::sim::LdmxMeasurementCalibrator::calibrate_1d>( - &calibrator); + .connect<&tracking::sim::LdmxMeasurementCalibrator::calibrate_1d< + Acts::VectorMultiTrajectory>>(&calibrator); + + // Surface Accessor + struct SurfaceAccessor { + const Acts::TrackingGeometry* trackingGeometry; + + const Acts::Surface* operator()(const Acts::SourceLink& sourceLink) const { + const auto& indexSourceLink = + sourceLink.get(); + return trackingGeometry->findSurface(indexSourceLink.geometryId()); + } + }; + + SurfaceAccessor m_slSurfaceAccessor{tg.getTG().get()}; + // m_slSurfaceAccessor.trackingGeometry = tg.getTG(); + gsf_extensions.surfaceAccessor.connect<&SurfaceAccessor::operator()>( + &m_slSurfaceAccessor); + gsf_extensions.mixtureReducer.connect<&Acts::reduceMixtureLargestWeights>(); // Propagator Options // Move this at the start of the producer - Acts::PropagatorOptions propagator_options( - geometry_context(), magnetic_field_context()); + Acts::PropagatorOptions + propagator_options(geometry_context(), magnetic_field_context()); propagator_options.pathLimit = std::numeric_limits::max(); @@ -193,12 +215,12 @@ void GSFProcessor::produce(framework::Event& event) { propagator_options.actionList.get(); sLogger.sterile = true; // Set a maximum step size - propagator_options.maxStepSize = + propagator_options.stepping.maxStepSize = propagator_step_size_ * Acts::UnitConstants::mm; propagator_options.maxSteps = propagator_maxSteps_; // Electron hypothesis - propagator_options.mass = 0.511 * Acts::UnitConstants::MeV; + // propagator_options.mass = 0.511 * Acts::UnitConstants::MeV; // GSF Options, to be moved out of produce loop @@ -206,12 +228,34 @@ void GSFProcessor::produce(framework::Event& event) { Acts::Surface::makeShared( Acts::Vector3(0., 0., 0.)); - Acts::Experimental::GsfOptions gsfOptions{ + std::shared_ptr tagger_layer_surface = + Acts::Surface::makeShared( + Acts::Vector3(-700., 0., 0.)); + + std::shared_ptr reference_surface = + origin_surface; + if (taggerTracking_) { + reference_surface = tagger_layer_surface; + } + + /* + Acts::GsfOptions gsfOptions{ geometry_context(), magnetic_field_context(), calibration_context(), gsf_extensions, propagator_options, &(*origin_surface), maxComponents_, weightCutoff_, abortOnError_, disableAllMaterialHandling_}; + */ + Acts::GsfOptions gsfOptions{ + geometry_context(), magnetic_field_context(), calibration_context()}; + + gsfOptions.extensions = gsf_extensions; + gsfOptions.propagatorPlainOptions = propagator_options; + gsfOptions.referenceSurface = &(*reference_surface); + gsfOptions.maxComponents = maxComponents_; + gsfOptions.weightCutoff = weightCutoff_; + gsfOptions.abortOnError = abortOnError_; + gsfOptions.disableAllMaterialHandling = disableAllMaterialHandling_; // Output track container std::vector out_tracks; @@ -266,16 +310,36 @@ void GSFProcessor::produce(framework::Event& event) { std::shared_ptr beamOrigin_surface = tracking::sim::utils::unboundSurface(-700); - if (!track.getTrackState(ldmx::TrackStateType::AtBeamOrigin).has_value()) { - ldmx_log(warn) - << "Failed retreiving AtBeamOrigin TrackState for track. Skipping.."; - continue; - } + const std::shared_ptr target_surface = + tracking::sim::utils::unboundSurface(0.); + + const std::shared_ptr ecal_surface = + tracking::sim::utils::unboundSurface(240.5); - auto ts = track.getTrackState(ldmx::TrackStateType::AtBeamOrigin).value(); Acts::BoundTrackParameters trk_btp_bO = - tracking::sim::utils::btp(ts, beamOrigin_surface); + tracking::sim::utils::boundTrackParameters(track, perigee); + if (taggerTracking_) { + if (!track.getTrackState(ldmx::TrackStateType::AtBeamOrigin) + .has_value()) { + ldmx_log(warn) << "Failed retreiving AtBeamOrigin TrackState for " + "track. Skipping.."; + continue; + } + + auto ts = track.getTrackState(ldmx::TrackStateType::AtBeamOrigin).value(); + trk_btp_bO = tracking::sim::utils::btp( + ts, beamOrigin_surface, + 11); // 11 == electron PDGid...hardcode for now + } else { + if (!track.getTrackState(ldmx::TrackStateType::AtTarget).has_value()) { + ldmx_log(warn) + << "Failed retreiving AtTarget TrackState for track. Skipping.."; + continue; + } + auto ts = track.getTrackState(ldmx::TrackStateType::AtTarget).value(); + trk_btp_bO = tracking::sim::utils::btp(ts, target_surface, 11); + } const Acts::BoundVector& trkpars = trk_btp.parameters(); ldmx_log(debug) << "CKF Track parameters" << std::endl << trkpars[0] << " " << trkpars[1] << " " << trkpars[2] @@ -334,17 +398,23 @@ void GSFProcessor::produce(framework::Event& event) { ldmx::Track trk = ldmx::Track(); - // Unbounded surface - const std::shared_ptr target_surface = - tracking::sim::utils::unboundSurface(0.); + bool success = false; + if (taggerTracking_) { + ldmx_log(debug) << "Target extrapolation"; + ldmx::Track::TrackState tsAtTarget; - ldmx_log(debug) << "Target extrapolation"; - ldmx::Track::TrackState tsAtTarget; + success = trk_extrap_->TrackStateAtSurface( + gsftrk, target_surface, tsAtTarget, ldmx::TrackStateType::AtTarget); - bool success = trk_extrap_->TrackStateAtSurface( - gsftrk, target_surface, tsAtTarget, ldmx::TrackStateType::AtTarget); + if (success) trk.addTrackState(tsAtTarget); + } else { + ldmx_log(debug) << "Ecal Extrapolation"; + ldmx::Track::TrackState tsAtEcal; + success = trk_extrap_->TrackStateAtSurface(gsftrk, ecal_surface, tsAtEcal, + ldmx::TrackStateType::AtECAL); - if (success) trk.addTrackState(tsAtTarget); + if (success) trk.addTrackState(tsAtEcal); + } trk.setPerigeeLocation( perigee_surface.transform(geometry_context()).translation()(0), @@ -352,6 +422,8 @@ void GSFProcessor::produce(framework::Event& event) { perigee_surface.transform(geometry_context()).translation()(2)); trk.setChi2(gsftrk.chi2()); + trk.setNhits(gsftrk.nMeasurements()); + trk.setNdf(gsftrk.nMeasurements() - 5); trk.setPerigeeParameters( tracking::sim::utils::convertActsToLdmxPars(perigee_pars)); std::vector v_trk_cov; diff --git a/Tracking/src/Tracking/Reco/SeedFinderProcessor.cxx b/Tracking/src/Tracking/Reco/SeedFinderProcessor.cxx index dbd42c166..d2e62fdaf 100644 --- a/Tracking/src/Tracking/Reco/SeedFinderProcessor.cxx +++ b/Tracking/src/Tracking/Reco/SeedFinderProcessor.cxx @@ -121,7 +121,7 @@ void SeedFinderProcessor::produce(framework::Event& event) { if (ts.has_value()) { auto trackState = ts.value(); - Acts::BoundSymMatrix cov = + Acts::BoundSquareMatrix cov = tracking::sim::utils::unpackCov(trackState.cov); double locu = trackState.params[0]; double locv = trackState.params[1]; @@ -154,18 +154,24 @@ void SeedFinderProcessor::produce(framework::Event& event) { ldmx_log(debug) << "Preparing the strategies"; groups_map.clear(); + // set the seeding strategy + // strategy is a list of layers from which to make the seed + // this must include 5 layers; layer numbering starts at 0. + // std::vector strategy = {9,10,11,12,13}; std::vector strategy = {0, 1, 2, 3, 4}; bool success = GroupStrips(measurements, strategy); if (success) FindSeedsFromMap(seed_tracks, target_pseudo_meas); + // currently, we only use a single strategy but eventually + // we will use more. Below is an example of how to add them /* groups_map.clear(); - strategy = {3,4,5,6,7}; + strategy = {9,10,11,12,13}; success = GroupStrips(measurements,strategy); if (success) - FindSeedsFromMap(seed_tracks); - + FindSeedsFromMap(seed_tracks, target_pseudo_meas); */ + groups_map.clear(); // outputTree_->Fill(); ntracks_ += seed_tracks.size(); @@ -209,6 +215,9 @@ void SeedFinderProcessor::produce(framework::Event& event) { // yOrigin is the location along the beam about which we fit the seed helix // perigee_location is where the track parameters will be extracted +// while this takes in a target measurement (from tagger, this is pmeas_tgt) +// this code doesn't do anything with it yet. + ldmx::Track SeedFinderProcessor::SeedTracker( const ldmx::Measurements& vmeas, double xOrigin, const Acts::Vector3& perigee_location, @@ -322,13 +331,21 @@ ldmx::Track SeedFinderProcessor::SeedTracker( // This is mainly necessary for the perigee surface, where // the mean might not fulfill the perigee condition. + // mg Aug 2024 .. interect has changed, but just remove boundary check + // and change intersection to intersections + // auto intersection = + // (*seed_perigee).intersect(geometry_context(), seed_pos, dir, false); + + // Acts::FreeVector seed_free = tracking::sim::utils::toFreeParameters( + // intersection.intersection.position, seed_mom, q); + auto intersection = - (*seed_perigee).intersect(geometry_context(), seed_pos, dir, false); + (*seed_perigee).intersect(geometry_context(), seed_pos, dir); Acts::FreeVector seed_free = tracking::sim::utils::toFreeParameters( - intersection.intersection.position, seed_mom, q); + intersection.intersections()[0].position(), seed_mom, q); - auto bound_params = Acts::detail::transformFreeToBoundParameters( + auto bound_params = Acts::transformFreeToBoundParameters( seed_free, *seed_perigee, geometry_context()) .value(); @@ -336,6 +353,7 @@ ldmx::Track SeedFinderProcessor::SeedTracker( << bound_params; Acts::BoundVector stddev; + // sigma set to 75% of momentum double sigma_p = 0.75 * p * Acts::UnitConstants::GeV; stddev[Acts::eBoundLoc0] = inflate_factors_[Acts::eBoundLoc0] * 2 * Acts::UnitConstants::mm; @@ -350,7 +368,11 @@ ldmx::Track SeedFinderProcessor::SeedTracker( stddev[Acts::eBoundTime] = inflate_factors_[Acts::eBoundTime] * 1000 * Acts::UnitConstants::ns; - Acts::BoundSymMatrix bound_cov = stddev.cwiseProduct(stddev).asDiagonal(); + ldmx_log(debug) + << "Making covariance matrix as diagonal matrix with inflated terms"; + Acts::BoundSquareMatrix bound_cov = stddev.cwiseProduct(stddev).asDiagonal(); + + ldmx_log(debug) << "...now putting together the seed track ..."; ldmx::Track trk = ldmx::Track(); trk.setPerigeeLocation(perigee_location(0), perigee_location(1), @@ -375,9 +397,15 @@ ldmx::Track SeedFinderProcessor::SeedTracker( trk.setMomentum(seed_free[Acts::eFreeDir0], seed_free[Acts::eFreeDir1], seed_free[Acts::eFreeDir2]); - Acts::BoundTrackParameters seedParameters(seed_perigee, - std::move(bound_params), bound_cov); + ldmx_log(debug) + << "...making the ParticleHypothesis ...assume electron for now"; + auto partHypo{Acts::SinglyChargedParticleHypothesis::electron()}; + + ldmx_log(debug) << "Making BoundTrackParameters seedParameters"; + Acts::BoundTrackParameters seedParameters( + seed_perigee, std::move(bound_params), bound_cov, partHypo); + ldmx_log(debug) << "Returning seed track"; return trk; } @@ -414,10 +442,11 @@ bool SeedFinderProcessor::GroupStrips( // std::cout< 0) pdgid = -11; // positron + auto part{Acts::GenericParticleHypothesis( + Acts::ParticleHypothesis(Acts::PdgParticle(pdgid)))}; Acts::BoundTrackParameters boundTrkPars(gen_surface, bound_params, - std::nullopt); + std::nullopt, part); // CAUTION:: The target surface should be close to the gen surface // Linear propagation to the target surface. I assume 1mm of tolerance @@ -389,7 +392,8 @@ ldmx::Track TruthSeedProcessor::seedFromTruth(const ldmx::Track& tt, (bound_params).data(), bound_params.data() + bound_params.rows() * bound_params.cols()); - Acts::BoundSymMatrix bound_cov = stddev.cwiseProduct(stddev).asDiagonal(); + Acts::BoundSquareMatrix bound_cov = + stddev.cwiseProduct(stddev).asDiagonal(); std::vector v_seed_cov; tracking::sim::utils::flatCov(bound_cov, v_seed_cov); seed.setPerigeeParameters(v_seed_params); @@ -414,7 +418,8 @@ ldmx::Track TruthSeedProcessor::seedFromTruth(const ldmx::Track& tt, stddev[Acts::eBoundTheta] = 5 * Acts::UnitConstants::degree; stddev[Acts::eBoundQOverP] = (1. / p) * (1. / p) * sigma_p; - Acts::BoundSymMatrix bound_cov = stddev.cwiseProduct(stddev).asDiagonal(); + Acts::BoundSquareMatrix bound_cov = + stddev.cwiseProduct(stddev).asDiagonal(); std::vector v_seed_cov; tracking::sim::utils::flatCov(bound_cov, v_seed_cov); seed.setPerigeeParameters(v_seed_params); diff --git a/Tracking/src/Tracking/Reco/VertexProcessor.cxx b/Tracking/src/Tracking/Reco/VertexProcessor.cxx index feb794402..54f57f1fa 100644 --- a/Tracking/src/Tracking/Reco/VertexProcessor.cxx +++ b/Tracking/src/Tracking/Reco/VertexProcessor.cxx @@ -60,26 +60,27 @@ void VertexProcessor::produce(framework::Event &event) { propagator_ = std::make_shared(stepper); // Track linearizer in the proximity of the vertex location - using Linearizer = Acts::HelicalTrackLinearizer; - Linearizer::Config linearizerConfig(sp_interpolated_bField_, propagator_); + using Linearizer = Acts::HelicalTrackLinearizer; + Linearizer::Config linearizerConfig; + linearizerConfig.bField = sp_interpolated_bField_; + linearizerConfig.propagator = propagator_; Linearizer linearizer(linearizerConfig); // Set up Billoir Vertex Fitter - using VertexFitter = - Acts::FullBilloirVertexFitter; + using VertexFitter = Acts::FullBilloirVertexFitter; VertexFitter::Config vertexFitterCfg; VertexFitter billoirFitter(vertexFitterCfg); - VertexFitter::State state(sp_interpolated_bField_->makeCache(bctx_)); + // VertexFitter::State state(sp_interpolated_bField_->makeCache(bctx_)); // Unconstrained fit // See // https://github.com/acts-project/acts/blob/main/Tests/UnitTests/Core/Vertexing/FullBilloirVertexFitterTests.cpp#L149 // For constraint implementation - Acts::VertexingOptions vfOptions(gctx_, bctx_); + Acts::VertexingOptions vfOptions(gctx_, bctx_); // Retrieve the track collection const std::vector tracks = @@ -109,11 +110,12 @@ void VertexProcessor::produce(framework::Event &event) { tracks.at(iTrack).getPhi(), tracks.at(iTrack).getTheta(), tracks.at(iTrack).getQoP(), tracks.at(iTrack).getT(); - Acts::BoundSymMatrix covMat = + Acts::BoundSquareMatrix covMat = tracking::sim::utils::unpackCov(tracks.at(iTrack).getPerigeeCov()); - + auto part{Acts::GenericParticleHypothesis(Acts::ParticleHypothesis( + Acts::PdgParticle(tracks.at(iTrack).getPdgID())))}; billoir_tracks.push_back(Acts::BoundTrackParameters( - perigeeSurface, paramVec, std::move(covMat))); + perigeeSurface, paramVec, std::move(covMat), part)); } // Select exactly 2 tracks @@ -149,11 +151,15 @@ void VertexProcessor::produce(framework::Event &event) { seeds.at(iSeed).getPhi(), seeds.at(iSeed).getTheta(), seeds.at(iSeed).getQoP(), seeds.at(iSeed).getT(); - Acts::BoundSymMatrix covMat = + Acts::BoundSquareMatrix covMat = tracking::sim::utils::unpackCov(seeds.at(iSeed).getPerigeeCov()); - + int pionPdgId = 211; // pi+ + if (seeds.at(iSeed).q() < 0) pionPdgId = -211; + // BoundTrackParameters needs the particle hypothesis + auto part{Acts::GenericParticleHypothesis( + Acts::ParticleHypothesis(Acts::PdgParticle(pionPdgId)))}; auto boundSeedParams = Acts::BoundTrackParameters( - perigeeSurface2, paramVec, std::move(covMat)); + perigeeSurface, paramVec, std::move(covMat), part); TLorentzVector pion4v; pion4v.SetXYZM(boundSeedParams.momentum()(0), diff --git a/Tracking/src/Tracking/Reco/Vertexer.cxx b/Tracking/src/Tracking/Reco/Vertexer.cxx index df909d59e..fbaced288 100644 --- a/Tracking/src/Tracking/Reco/Vertexer.cxx +++ b/Tracking/src/Tracking/Reco/Vertexer.cxx @@ -92,14 +92,14 @@ void Vertexer::produce(framework::Event& event) { // auto start = std::chrono::high_resolution_clock::now(); // Track linearizer in the proximity of the vertex location - using Linearizer = Acts::HelicalTrackLinearizer; - // Linearizer::Config linearizerConfig(sp_interpolated_bField_,propagator_); - Linearizer::Config linearizerConfig(bField_, propagator_); + using Linearizer = Acts::HelicalTrackLinearizer; + Linearizer::Config linearizerConfig; + linearizerConfig.bField = bField_; + linearizerConfig.propagator = propagator_; Linearizer linearizer(linearizerConfig); // Set up Billoir Vertex Fitter - using VertexFitter = - Acts::FullBilloirVertexFitter; + using VertexFitter = Acts::FullBilloirVertexFitter; // Alternatively one can use // using VertexFitter = @@ -107,15 +107,18 @@ void Vertexer::produce(framework::Event& event) { VertexFitter::Config vertexFitterCfg; VertexFitter billoirFitter(vertexFitterCfg); - - VertexFitter::State state(sp_interpolated_bField_->makeCache(bctx_)); + // mg Aug 2024 .. State doesn't exist in v36 and isn't used here anyway + // VertexFitter::State state(sp_interpolated_bField_->makeCache(bctx_)); // Unconstrained fit // See // https://github.com/acts-project/acts/blob/main/Tests/UnitTests/Core/Vertexing/FullBilloirVertexFitterTests.cpp#L149 // For constraint implementation - Acts::VertexingOptions vfOptions(gctx_, bctx_); + // Acts::VertexingOptions vfOptions(gctx_, + // bctx_); + // mg Aug 2024 ... VertexingOptions template change in v36 + Acts::VertexingOptions vfOptions(gctx_, bctx_); // Retrive the two track collections @@ -156,7 +159,8 @@ void Vertexer::produce(framework::Event& event) { tracking::sim::utils::boundTrackParameters(trk, perigeeSurface)); } - std::vector > fit_vertices; + // std::vector > fit_vertices; + std::vector fit_vertices; for (auto& b_trk_1 : billoir_tracks_1) { std::vector fit_tracks_ptr; diff --git a/Tracking/src/Tracking/Sim/PropagatorStepWriter.cxx b/Tracking/src/Tracking/Sim/PropagatorStepWriter.cxx index e1e89a6a9..cb0066886 100644 --- a/Tracking/src/Tracking/Sim/PropagatorStepWriter.cxx +++ b/Tracking/src/Tracking/Sim/PropagatorStepWriter.cxx @@ -7,8 +7,9 @@ #include #include "Acts/Geometry/DetectorElementBase.hpp" -#include "Acts/Plugins/Identification/IdentifiedDetectorElement.hpp" -#include "Acts/Plugins/Identification/Identifier.hpp" +// mg ... I don't think these are used, and they are not defined in acts v36 +//#include "Acts/Plugins/Identification/IdentifiedDetectorElement.hpp" +//#include "Acts/Plugins/Identification/Identifier.hpp" namespace tracking { namespace sim { @@ -35,7 +36,7 @@ PropagatorStepWriter::PropagatorStepWriter( // Set the branches m_outputTree->Branch("event_nr", &m_eventNr); - m_outputTree->Branch("volume_id", &m_volumeID); + // m_outputTree->Branch("volume_id", &m_volumeID); m_outputTree->Branch("boundary_id", &m_boundaryID); m_outputTree->Branch("layer_id", &m_layerID); m_outputTree->Branch("approach_id", &m_approachID); @@ -105,7 +106,7 @@ bool PropagatorStepWriter::WriteSteps( // loop over the step vector of each test propagation in this for (auto& steps : stepCollection) { // clear the vectors for each collection - m_volumeID.clear(); + // m_volumeID.clear(); m_boundaryID.clear(); m_layerID.clear(); m_approachID.clear(); @@ -125,7 +126,7 @@ bool PropagatorStepWriter::WriteSteps( // loop over single steps for (auto& step : steps) { // the identification of the step - Acts::GeometryIdentifier::Value volumeID = 0; + // Acts::GeometryIdentifier::Value volumeID = 0; Acts::GeometryIdentifier::Value boundaryID = 0; Acts::GeometryIdentifier::Value layerID = 0; Acts::GeometryIdentifier::Value approachID = 0; @@ -133,22 +134,23 @@ bool PropagatorStepWriter::WriteSteps( // get the identification from the surface first if (step.surface) { auto geoID = step.surface->geometryId(); - volumeID = geoID.volume(); + // volumeID = geoID.volume(); boundaryID = geoID.boundary(); layerID = geoID.layer(); approachID = geoID.approach(); sensitiveID = geoID.sensitive(); } // a current volume overwrites the surface tagged one - if (step.volume) { - volumeID = step.volume->geometryId().volume(); - } + // mg ... v36 Step does not include volume + // if (step.volume) { + // volumeID = step.volume->geometryId().volume(); + // } // now fill m_sensitiveID.push_back(sensitiveID); m_approachID.push_back(approachID); m_layerID.push_back(layerID); m_boundaryID.push_back(boundaryID); - m_volumeID.push_back(volumeID); + // m_volumeID.push_back(volumeID); // kinematic information m_x.push_back(step.position.x()); @@ -159,7 +161,9 @@ bool PropagatorStepWriter::WriteSteps( m_dy.push_back(direction.y()); m_dz.push_back(direction.z()); - double accuracy = step.stepSize.value(Acts::ConstrainedStep::accuracy); + // double accuracy = + // step.stepSize.value(Acts::ConstrainedStep::accuracy()); + double accuracy = step.stepSize.accuracy(); double actor = step.stepSize.value(Acts::ConstrainedStep::actor); double aborter = step.stepSize.value(Acts::ConstrainedStep::aborter); double user = step.stepSize.value(Acts::ConstrainedStep::user); diff --git a/Tracking/src/Tracking/dqm/TrackingRecoDQM.cxx b/Tracking/src/Tracking/dqm/TrackingRecoDQM.cxx index 2f2b5d033..8cc52873c 100644 --- a/Tracking/src/Tracking/dqm/TrackingRecoDQM.cxx +++ b/Tracking/src/Tracking/dqm/TrackingRecoDQM.cxx @@ -17,6 +17,8 @@ void TrackingRecoDQM::configure(framework::config::Parameters& parameters) { subdetector_ = parameters.getParameter("subdetector", "Tagger"); trackStates_ = parameters.getParameter>("trackStates", {}); + measurementCollection_ = parameters.getParameter( + "measurement_collection", "DigiTaggerSimHits"); ldmx_log(info) << "Track Collection " << trackCollection_ << std::endl; ldmx_log(info) << "Truth Collection " << truthCollection_ << std::endl; @@ -40,7 +42,8 @@ void TrackingRecoDQM::analyze(const framework::Event& event) { return; } auto tracks{event.getCollection(trackCollection_)}; - + auto measurements{ + event.getCollection(measurementCollection_)}; // The truth track collection if (event.exists(truthCollection_)) { truthTrackCollection_ = std::make_shared( @@ -73,12 +76,14 @@ void TrackingRecoDQM::analyze(const framework::Event& event) { histograms_.fill(title_ + "N_tracks", tracks.size()); ldmx_log(debug) << "Track Monitoring on Unique Tracks" << std::endl; - TrackMonitoring(uniqueTracks_, title_, true, true); + + TrackMonitoring(uniqueTracks_, measurements, title_, true, true); ldmx_log(debug) << "Track Monitoring on duplicates and fakes" << std::endl; // Fakes and duplicates - TrackMonitoring(duplicateTracks_, title_ + "dup_", false, false); - TrackMonitoring(fakeTracks_, title_ + "fake_", false, false); + TrackMonitoring(duplicateTracks_, measurements, title_ + "dup_", false, + false); + TrackMonitoring(fakeTracks_, measurements, title_ + "fake_", false, false); // Track Extrapolation to Ecal Monitoring @@ -99,7 +104,7 @@ void TrackingRecoDQM::analyze(const framework::Event& event) { } // Technical Efficiency plots - EfficiencyPlots(tracks, title_); + EfficiencyPlots(tracks, measurements, title_); // Tagger Recoil Matching @@ -113,10 +118,13 @@ void TrackingRecoDQM::onProcessEnd() { // Produce the efficiency plots. (TODO::Switch to TEfficiency instead) } -void TrackingRecoDQM::EfficiencyPlots(const std::vector& tracks, - const std::string& title) { +void TrackingRecoDQM::EfficiencyPlots( + const std::vector& tracks, + const std::vector& measurements, + const std::string& title) { // Do all truth track plots - denominator + histograms_.fill(title + "truth_N_tracks", truthTrackCollection_->size()); for (auto& truth_trk : *(truthTrackCollection_)) { double truth_phi = truth_trk.getPhi(); double truth_d0 = truth_trk.getD0(); @@ -124,6 +132,8 @@ void TrackingRecoDQM::EfficiencyPlots(const std::vector& tracks, double truth_theta = truth_trk.getTheta(); double truth_qop = truth_trk.getQoP(); double truth_p = 1. / abs(truth_trk.getQoP()); + int truth_nHits = truth_trk.getNhits(); + std::vector truth_mom = truth_trk.getMomentum(); // Polar angle @@ -133,6 +143,7 @@ void TrackingRecoDQM::EfficiencyPlots(const std::vector& tracks, double truth_beam_angle = std::atan2(truth_pt_beam, truth_mom[0]); + histograms_.fill(title + "truth_nHits", truth_nHits); histograms_.fill(title + "truth_d0", truth_d0); histograms_.fill(title + "truth_z0", truth_z0); histograms_.fill(title + "truth_phi", truth_phi); @@ -211,7 +222,7 @@ void TrackingRecoDQM::EfficiencyPlots(const std::vector& tracks, double truth_beam_angle = std::atan2(truth_pt_beam, truth_mom[0]); // Fill reco plots for efficiencies - numerator. The quantities are truth - + histograms_.fill(title + "match_prob", trackTruthProb); histograms_.fill(title + "match_d0", truth_d0); histograms_.fill(title + "match_z0", truth_z0); histograms_.fill(title + "match_phi", truth_phi); @@ -219,6 +230,11 @@ void TrackingRecoDQM::EfficiencyPlots(const std::vector& tracks, histograms_.fill(title + "match_p", truth_p); histograms_.fill(title + "match_qop", truth_qop); histograms_.fill(title + "match_beam_angle", truth_beam_angle); + histograms_.fill(title + "match_nHits", measurements.size()); + for (auto trackHit : track.getMeasurementsIdxs()) { + histograms_.fill(title + "match_LayersHit", + measurements.at(trackHit).getLayer()); + } // For some particles @@ -259,19 +275,22 @@ void TrackingRecoDQM::EfficiencyPlots(const std::vector& tracks, } // Efficiency plots -void TrackingRecoDQM::TrackMonitoring(const std::vector& tracks, - const std::string title, - const bool& doDetail, - const bool& doTruth) { +void TrackingRecoDQM::TrackMonitoring( + const std::vector& tracks, + const std::vector& measurements, const std::string title, + const bool& doDetail, const bool& doTruth) { for (auto& track : tracks) { // Perigee track parameters - double trk_d0 = track.getD0(); double trk_z0 = track.getZ0(); double trk_qop = track.getQoP(); double trk_theta = track.getTheta(); double trk_phi = track.getPhi(); double trk_p = 1. / abs(trk_qop); + for (auto trackHit : track.getMeasurementsIdxs()) { + histograms_.fill(title + "LayersHit", + measurements.at(trackHit).getLayer()); + } std::vector trk_mom = track.getMomentum(); @@ -284,7 +303,7 @@ void TrackingRecoDQM::TrackMonitoring(const std::vector& tracks, std::sqrt(trk_mom[1] * trk_mom[1] + trk_mom[2] * trk_mom[2]); // Covariance matrix - Acts::BoundSymMatrix cov = + Acts::BoundSquareMatrix cov = tracking::sim::utils::unpackCov(track.getPerigeeCov()); double sigmad0 = sqrt( @@ -367,7 +386,6 @@ void TrackingRecoDQM::TrackMonitoring(const std::vector& tracks, double truth_theta = truth_trk->getTheta(); double truth_qop = truth_trk->getQoP(); double truth_p = 1. / abs(truth_trk->getQoP()); - std::vector truth_mom = truth_trk->getMomentum(); // Polar angle // The momentum in the plane transverse wrt the beam axis @@ -477,7 +495,8 @@ void TrackingRecoDQM::TrackStateMonitoring(const ldmx::Tracks& tracks, ldmx::Track::TrackState& TargetState = trk_ts.value(); ldmx_log(debug) << "Unpacking covariance matrix" << std::endl; - Acts::BoundSymMatrix cov = tracking::sim::utils::unpackCov(TargetState.cov); + Acts::BoundSquareMatrix cov = + tracking::sim::utils::unpackCov(TargetState.cov); [[maybe_unused]] double sigmaloc0 = sqrt( cov(Acts::BoundIndices::eBoundLoc0, Acts::BoundIndices::eBoundLoc0)); diff --git a/Tracking/src/Tracking/geo/EcalTrackingGeometry.cxx b/Tracking/src/Tracking/geo/EcalTrackingGeometry.cxx index 65b1aed16..fa50f8dce 100644 --- a/Tracking/src/Tracking/geo/EcalTrackingGeometry.cxx +++ b/Tracking/src/Tracking/geo/EcalTrackingGeometry.cxx @@ -353,7 +353,8 @@ EcalTrackingGeometry::convertHexToActsSurface(const G4VPhysicalVolume& phex) { surface->assignSurfaceMaterial( std::make_shared(silicon_slab)); - surface->toStream(*gctx_, std::cout); + // surface->toStreamImpl(*gctx_, std::cout); + surface->toStream(*gctx_); return surface; } diff --git a/Tracking/src/Tracking/geo/TrackersTrackingGeometry.cxx b/Tracking/src/Tracking/geo/TrackersTrackingGeometry.cxx index 5d976bd48..5651a2353 100644 --- a/Tracking/src/Tracking/geo/TrackersTrackingGeometry.cxx +++ b/Tracking/src/Tracking/geo/TrackersTrackingGeometry.cxx @@ -190,7 +190,9 @@ TrackersTrackingGeometry::buildTrackerVolume() { if (debug_) { std::cout << layer.first << " : surfaces==>" << layer.second.size() << std::endl; - for (auto& surface : layer.second) surface->toStream(gctx_, std::cout); + // for (auto& surface : layer.second) surface->toStreamImpl(gctx_, + // std::cout); + for (auto& surface : layer.second) surface->toStream(gctx_); } Acts::CuboidVolumeBuilder::LayerConfig lcfg; diff --git a/Tracking/src/Tracking/geo/TrackingGeometry.cxx b/Tracking/src/Tracking/geo/TrackingGeometry.cxx index 2b4670349..3b9b1f491 100644 --- a/Tracking/src/Tracking/geo/TrackingGeometry.cxx +++ b/Tracking/src/Tracking/geo/TrackingGeometry.cxx @@ -154,7 +154,8 @@ void TrackingGeometry::dumpGeometry(const std::string& outputDir, for (auto const& surfaceId : layer_surface_map_) { std::cout << " " << surfaceId.first << std::endl; std::cout << " Check the surface" << std::endl; - surfaceId.second->toStream(gctx, std::cout); + // surfaceId.second->toStream(gctx, std::cout); + surfaceId.second->toStream(gctx); std::cout << " GeometryID::" << surfaceId.second->geometryId() << std::endl; std::cout << " GeometryID::" << surfaceId.second->geometryId().value()