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()