diff --git a/Rec/Allen/src/GaudiAllenTrackViewsToV3Tracks.cpp b/Rec/Allen/src/GaudiAllenTrackViewsToV3Tracks.cpp index 9100d887e9108f266c773b8ef3a7dac5382634e2..fcf38b7e112c334b7184093a4b51b912299a5b9e 100644 --- a/Rec/Allen/src/GaudiAllenTrackViewsToV3Tracks.cpp +++ b/Rec/Allen/src/GaudiAllenTrackViewsToV3Tracks.cpp @@ -16,9 +16,13 @@ // LHCb #include "Event/Track.h" #include "Event/Track_v3.h" +#include "Event/PrVeloTracks.h" #include "Event/TrackEnums.h" #include "Event/UniqueIDGenerator.h" #include "Event/StateParameters.h" +#include "LHCbDet/InteractionRegion.h" +#include "DetDesc/GenericConditionAccessorHolder.h" +#include "LHCbAlgs/Transformer.h" // Allen #include "Logger.h" @@ -810,6 +814,195 @@ namespace GaudiAllen::Converters::v3 { } }; + using AllenOutType = std::tuple< + std::vector, + std::vector, + std::vector>; + + /** + * + * Converter from PrVeloTracks to Allen::Views::Velo::Consolidated::Tracks and beamline states. + * The converter is a MultiTransformer that takes two PrVeloTracks as input and returns + * Allen::Views::Velo::Consolidated::Tracks and beamline states. + * + */ + class PrVeloToGaudiAllenVeloTracks final + : public Gaudi::Functional::MultiTransformer< + AllenOutType( + LHCb::Pr::Velo::Tracks const&, + LHCb::Pr::Velo::Tracks const&, + const LHCb::Conditions::InteractionRegion&), + LHCb::Algorithm::Traits::usesConditions> { + public: + //============================================================================= + // Standard constructor, initializes variables + //============================================================================= + PrVeloToGaudiAllenVeloTracks(const std::string& name, ISvcLocator* pSvcLocator) : + MultiTransformer( + name, + pSvcLocator, + {KeyValue {"PrVeloTracksForward", ""}, + KeyValue {"PrVeloTracksBackward", ""}, + {KeyValue {"InteractionRegionCache", "AlgorithmSpecific-" + name + "-InteractionRegion"}}}, + {KeyValue {"AllenVeloTracksView", ""}, + KeyValue {"AllenBeamlineStatesView", ""}, + KeyValue {"NumberOfReconstructedVeloTracks", ""}}) + {} + + //============================================================================= + // initialize() + //============================================================================= + StatusCode initialize() override + { + auto sc = MultiTransformer::initialize().andThen([&]() { + LHCb::Conditions::InteractionRegion::addConditionDerivation( + this, inputLocation()); + }); + return sc; + } + + //============================================================================= + // execute() + //============================================================================= + AllenOutType operator()( + LHCb::Pr::Velo::Tracks const& pr_tracks_forward, + LHCb::Pr::Velo::Tracks const& pr_tracks_backward, + LHCb::Conditions::InteractionRegion const& region) const override + { + const unsigned n_events = 1; // only one event is processed at a time + const auto simd = SIMDWrapper::InstructionSet::Scalar; // scalar instruction set + const unsigned n_tracks = (pr_tracks_forward.size() + pr_tracks_backward.size()) / + SIMDWrapper::type_map_t::size; // total number of tracks + + auto offset_tracks = std::make_unique(n_events + 1); // memory for the track offsets + offset_tracks[n_events - 1] = 0; + offset_tracks[n_events] = n_tracks; + + auto ptr_tracks = new Allen::Views::Velo::Consolidated::Track[n_tracks]; // memory for the tracks + unsigned nb_elements_state = 6; // number of parameters in the state + unsigned nb_elements_cov = 8; // number of elements in the covariance matrix + auto ptr_states = new float[n_tracks * (nb_elements_state + nb_elements_cov)]; // memory for the states + std::vector ptr_hits; // pointer to hits + // it is necessary to instantiate the Consolidated::Hits object + auto hits = new Allen::Views::Velo::Consolidated::Hits[n_tracks]; // memory for the hits + auto offset_hits = std::make_unique(n_tracks + 1); // memory for the hit offsets + + unsigned i_track = 0; // track index + unsigned tot_hits = 0; // total number of hits + // Loop over the tracks + for (const auto& tracks : {&pr_tracks_forward, &pr_tracks_backward}) { + // Get scalar view of the tracks, to cast to the correct type + for (auto const& track : tracks->scalar()) { + // Get the state parameters + auto loc = LHCb::Event::Enum::State::Location::ClosestToBeam; + auto pos = track.StatePos(loc); + auto dir = track.StateDir(loc); + auto covX = track.StateCovX(loc); + auto covY = track.StateCovY(loc); + + // Get the beamline position and direction + const auto beamspot = region.avgPosition; + const auto beamlineX = beamspot.x(); + const auto beamlineY = beamspot.y(); + const auto beamlineTx = region.tX(); + const auto beamlineTy = region.tY(); + + /// Extrapolate the state at the beamline + auto const dx = beamlineX - pos.x().cast(); + auto const dy = beamlineY - pos.y().cast(); + auto const tx = dir.x().cast() - beamlineTx; + auto const ty = dir.y().cast() - beamlineTy; + auto const dz = tx * ty != 0 ? (tx * dx + ty * dy) / (tx * tx + ty * ty) : 0.f; + auto const dz2 = dz * dz; + auto const xbeam = pos.x().cast() + tx * dz; + auto const ybeam = pos.y().cast() + ty * dz; + auto const zbeam = pos.z().cast() + dz; + auto const Vx = (covX.x().cast() + 2.f * dz * covX.y() + dz2 * covX.z().cast()).cast(); + auto const Vy = (covY.x().cast() + 2.f * dz * covY.y() + dz2 * covY.z().cast()).cast(); + auto const Vxy = covX.y().cast() + dz * covX.z().cast(); + auto const Vyx = covY.y().cast() + dz * covY.z().cast(); + auto const Wx = Vx != 0 ? 1. / Vx : 1.f; + auto const Wy = Vy != 0 ? 1. / Vy : 1.f; + auto const blchi2 = dx * Wx * dx + dy * Wy * dy; + + // Fill the state parameters + ptr_states[i_track * nb_elements_state] = xbeam; + ptr_states[i_track * nb_elements_state + 1] = ybeam; + ptr_states[i_track * nb_elements_state + 2] = zbeam; + ptr_states[i_track * nb_elements_state + 3] = tx; + ptr_states[i_track * nb_elements_state + 4] = ty; + // no qop in the PrVeloTracks, it should be extracted from the first hit + // we can use the function get_qop to get the value, but it is not used anywhere + // in the Allen PV finding algorithm, so set it to 0 + ptr_states[i_track * nb_elements_state + 5] = 0; + ptr_states[n_tracks * nb_elements_state + i_track * nb_elements_cov] = Vx; + ptr_states[n_tracks * nb_elements_state + i_track * nb_elements_cov + 1] = Vxy; + ptr_states[n_tracks * nb_elements_state + i_track * nb_elements_cov + 2] = covX.z().cast(); + ptr_states[n_tracks * nb_elements_state + i_track * nb_elements_cov + 3] = Vy; + ptr_states[n_tracks * nb_elements_state + i_track * nb_elements_cov + 4] = Vyx; + ptr_states[n_tracks * nb_elements_state + i_track * nb_elements_cov + 5] = covY.z().cast(); + ptr_states[n_tracks * nb_elements_state + i_track * nb_elements_cov + 6] = blchi2; + ptr_states[n_tracks * nb_elements_state + i_track * nb_elements_cov + 7] = 5; // default ndf (x,y,z,tx,ty) + + // Get the number of hits + unsigned n_hits = track.nHits().cast(); + + offset_hits[i_track] = tot_hits; + offset_hits[i_track + 1] = tot_hits + n_hits; + + // Loop over the hits to get the LHCbID + for (uint i_hit = 0; i_hit < n_hits; i_hit++) { + ptr_hits.push_back(track.vp_lhcbID(i_hit).LHCbID().lhcbID()); + } + + tot_hits += n_hits; + i_track++; + } + } + + for (UInt_t i_track = 0; i_track < n_tracks; i_track++) { + hits[i_track] = Allen::Views::Velo::Consolidated::Hits( + reinterpret_cast(ptr_hits.data()), + offset_tracks.get(), + offset_hits.get(), + n_events - 1, + n_events); // Create a new Allen hit container + + ptr_tracks[i_track] = Allen::Views::Velo::Consolidated::Track( + hits, offset_tracks.get(), offset_hits.get(), i_track, n_events - 1); // Create a new Allen track + } + + // Create the output container + std::vector allen_states_view; + allen_states_view.emplace_back(reinterpret_cast(ptr_states), offset_tracks.get(), n_events - 1, n_events); + std::vector allen_tracks_view; + allen_tracks_view.emplace_back(ptr_tracks, offset_tracks.get(), n_events - 1); + + std::vector number_of_tracks; + number_of_tracks.push_back(static_cast(allen_tracks_view[0].size())); + + assert(number_of_tracks[0] == n_tracks); + + return std::tuple {allen_tracks_view, allen_states_view, number_of_tracks}; + } + + private: + template + float get_qop(const PrVeloTrack& track) const + { + auto dir = track.closestToBeamStateDir(); + // no qop in the PrVeloTracks, needs to be extracted from the PrVeloTrack + const float tx = dir.x().cast(); + const float ty = dir.y().cast(); + const float slope2 = std::max(tx * tx + ty * ty, 1.e-20f); + const int firstRow = track.vp_lhcbID(0).LHCbID().channelID(); + const float charge = (firstRow % 2 == 0 ? -1.f : 1.f); + const float pT = 400 * Gaudi::Units::MeV; // default pT for Velo tracks + const float qop = charge / (pT * std::sqrt(1.f + 1.f / slope2)); + return qop; + } + }; + using GaudiAllenMEBasicParticlesToV3Tracks = GaudiAllenTrackViewsToV3Tracks; DECLARE_COMPONENT_WITH_ID(GaudiAllenMEBasicParticlesToV3Tracks, "GaudiAllenMEBasicParticlesToV3Tracks") @@ -823,4 +1016,5 @@ namespace GaudiAllen::Converters::v3 { beamline_states, endvelo_states>; DECLARE_COMPONENT_WITH_ID(GaudiAllenUTToV3Tracks, "GaudiAllenUTToV3Tracks") + DECLARE_COMPONENT_WITH_ID(PrVeloToGaudiAllenVeloTracks, "PrVeloToGaudiAllenVeloTracks") } // namespace GaudiAllen::Converters::v3 diff --git a/configuration/python/AllenConf/primary_vertex_reconstruction.py b/configuration/python/AllenConf/primary_vertex_reconstruction.py index e6f21861d2b402ce5675ed1ffc51f99e58f70d78..e430bfeadbda5ced8c954a40b227117f8d80336c 100644 --- a/configuration/python/AllenConf/primary_vertex_reconstruction.py +++ b/configuration/python/AllenConf/primary_vertex_reconstruction.py @@ -25,7 +25,8 @@ def make_pvs(velo_tracks, zmin=-541., zmax=307., SMOG2_pp_separation=-334., - Nbins=3392): + Nbins=3392, + use_converted_tracks=False): dz = 0.25 pp_maxTrackZ0Err = 1.5 @@ -42,16 +43,18 @@ def make_pvs(velo_tracks, number_of_events = initialize_number_of_events() host_number_of_events = number_of_events["host_number_of_events"] - dev_number_of_events = number_of_events["dev_number_of_events"] host_number_of_reconstructed_velo_tracks = velo_tracks[ "host_number_of_reconstructed_velo_tracks"] - dev_offsets_all_velo_tracks = velo_tracks["dev_offsets_all_velo_tracks"] - dev_offsets_velo_track_hit_number = velo_tracks[ - "dev_offsets_velo_track_hit_number"] - dev_velo_track_hits = velo_tracks["dev_velo_track_hits"] - velo_states = run_velo_kalman_filter(velo_tracks, pv_name) + if use_converted_tracks: + velo_states_view = velo_tracks["dev_velo_kalman_beamline_states_view"] + else: + velo_states = run_velo_kalman_filter(velo_tracks, pv_name) + velo_states_view = velo_states["dev_velo_kalman_beamline_states_view"] + + if use_converted_tracks: + pv_name += "_{hash}" pv_beamline_extrapolate = make_algorithm( pv_beamline_extrapolate_t, @@ -59,8 +62,7 @@ def make_pvs(velo_tracks, host_number_of_reconstructed_velo_tracks_t= host_number_of_reconstructed_velo_tracks, dev_velo_tracks_view_t=velo_tracks["dev_velo_tracks_view"], - dev_velo_states_view_t=velo_states[ - "dev_velo_kalman_beamline_states_view"]) + dev_velo_states_view_t=velo_states_view) pv_beamline_histo = make_algorithm( pv_beamline_histo_t,