diff --git a/Trigger/EFTracking/FPGATrackSim/FPGATrackSimHough/CMakeLists.txt b/Trigger/EFTracking/FPGATrackSim/FPGATrackSimHough/CMakeLists.txt index 6631195c2487144c6f63e80fc329b86730dfa49e..aead05d87b75bd973ecf919e2422769a7c8854fd 100644 --- a/Trigger/EFTracking/FPGATrackSim/FPGATrackSimHough/CMakeLists.txt +++ b/Trigger/EFTracking/FPGATrackSim/FPGATrackSimHough/CMakeLists.txt @@ -12,7 +12,7 @@ atlas_add_library( FPGATrackSimHoughLib src/*.cxx PUBLIC_HEADERS FPGATrackSimHough INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} - LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps GaudiKernel FPGATrackSimMapsLib FPGATrackSimObjectsLib + LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps GaudiKernel FourMomUtils FPGATrackSimMapsLib FPGATrackSimObjectsLib PRIVATE_INCLUDE_DIRS ${Boost_INCLUDE_DIRS} PRIVATE_LINK_LIBRARIES ${Boost_LIBRARIES} FPGATrackSimBanksLib FPGATrackSimConfToolsLib ) diff --git a/Trigger/EFTracking/FPGATrackSim/FPGATrackSimHough/src/FPGATrackSimGenScanBinning.cxx b/Trigger/EFTracking/FPGATrackSim/FPGATrackSimHough/src/FPGATrackSimGenScanBinning.cxx new file mode 100644 index 0000000000000000000000000000000000000000..8b3b907acdef39e499b5ec2ba05304a0c6a0e0e6 --- /dev/null +++ b/Trigger/EFTracking/FPGATrackSim/FPGATrackSimHough/src/FPGATrackSimGenScanBinning.cxx @@ -0,0 +1,396 @@ +// Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration + +/** + * @file FPGATrackSimGenScanBinning.cxx + * @author Elliot Lipeles + * @date Sept 10th, 2024 + * @brief See header file. + */ + +#include "FPGATrackSimGenScanBinning.h" + +//-------------------------------------------------------------------------------------------------- +// +// FPGATrackSimGenScanBinningBase base class implementation +// +//-------------------------------------------------------------------------------------------------- + +bool FPGATrackSimGenScanBinningBase::inRange(const ParSet &pars) const +{ + for (unsigned i = 0; i < NPars; i++) + { + if (!inRange(i, pars[i])) + return false; + } + return true; +} + +FPGATrackSimGenScanBinningBase::IdxSet FPGATrackSimGenScanBinningBase::binIdx(const ParSet &pars) const +{ + IdxSet retv(NPars); + for (unsigned i = 0; i < NPars; i++) + { + retv[i] = binIdx(i, pars[i]); + } + return retv; +} +const FPGATrackSimGenScanBinningBase::IdxSet FPGATrackSimGenScanBinningBase::parsToBin(FPGATrackSimTrackPars &pars) const +{ + ParSet parset = trackParsToParSet(pars); + if (inRange(parset)) + { + return binIdx(parset); + } + return m_invalidBin; +} + +FPGATrackSimGenScanBinningBase::ParSet FPGATrackSimGenScanBinningBase::center() const +{ + ParSet parset(NPars); + for (unsigned i = 0; i < NPars; i++) + { + parset[i] = (m_parMin[i] + m_parMin[i]) / 2.0; + } + return parset; +} + +FPGATrackSimGenScanBinningBase::ParSet FPGATrackSimGenScanBinningBase::binLowEdge(const IdxSet &idx) const +{ + ParSet parset(NPars); + for (unsigned i = 0; i < NPars; i++) + { + parset[i] = binLowEdge(i, idx[i]); + } + return parset; +} + +FPGATrackSimGenScanBinningBase::ParSet FPGATrackSimGenScanBinningBase::binCenter(const IdxSet &idx) const +{ + ParSet parset(NPars); + for (unsigned i = 0; i < NPars; i++) + { + parset[i] = binCenter(i, idx[i]); + } + return parset; +} + +std::vector<unsigned> FPGATrackSimGenScanBinningBase::subVec(const std::vector<unsigned>& elems, const std::vector<unsigned>& invec) const +{ + std::vector<unsigned> retv; + for (auto elem : elems) + { + retv.push_back(invec[elem]); + } + return retv; +} + +StatusCode FPGATrackSimGenScanBinningBase::setIdxSubVec(IdxSet &idx, const std::vector<unsigned>& subvecelems, const std::vector<unsigned>& subvecidx) const +{ + if (subvecelems.size()!=subvecidx.size()) { + return StatusCode::FAILURE; + } + for (unsigned i = 0; i < subvecelems.size(); i++) + { + if (subvecelems[i] >= idx.size()) { + return StatusCode::FAILURE; + } + idx[subvecelems[i]] = subvecidx[i]; + } + return StatusCode::SUCCESS; +} + + +// This gives a list tracks parameters for the corners of bin of dimensions scanpars.size() +std::vector<FPGATrackSimGenScanBinningBase::ParSet> FPGATrackSimGenScanBinningBase::makeVariationSet(const std::vector<unsigned> &scanpars, const IdxSet &idx) const +{ + std::vector<ParSet> retv; + for (unsigned corners = 0; corners < unsigned((1 << scanpars.size())); corners++) + { + IdxSet newidx = idx; + int scanDimCnt = 0; + for (auto &par : scanpars) + { + newidx[par] = idx[par] + ((corners >> scanDimCnt) & 1); + scanDimCnt++; + } + retv.push_back(binLowEdge(newidx)); + } + return retv; +} + +// find valid range for sliceVar +bool FPGATrackSimGenScanBinningBase::hitInSlice(const FPGATrackSimGenScanBinningBase::IdxSet &idx, FPGATrackSimHit const *hit) const +{ + // Get float values + double slicevarmin = std::numeric_limits<double>::max(); + double slicevarmax = std::numeric_limits<double>::min(); + + auto variations = makeVariationSet(slicePars(), idx); + + for (auto &pars : variations) + { + double slicevarhit = sliceVarExpected(pars, hit); + slicevarmin = std::min(slicevarmin, slicevarhit); + slicevarmax = std::max(slicevarmax, slicevarhit); + } + double slicevar = sliceVar(hit); + return ( slicevar > slicevarmin) && (slicevar < slicevarmax); +} + + + +std::pair<unsigned, unsigned> FPGATrackSimGenScanBinningBase::idxsetToRowParBinRange(const FPGATrackSimGenScanBinningBase::IdxSet &idx, FPGATrackSimHit const *hit) const +{ + // Get float values + double rowparmin = std::numeric_limits<double>::max(); + double rowparmax = std::numeric_limits<double>::min(); + + auto variations = makeVariationSet(scanPars(), idx); + + for (auto &pars : variations) + { + double rowparhit = rowPar(pars, hit); + rowparmin = std::min(rowparmin, rowparhit); + rowparmax = std::max(rowparmax, rowparhit); + } + + // only consider parameter bins in range + unsigned lowbin = std::max(rowParBinIdx(rowparmin), unsigned(0)); + unsigned highbin = std::min(rowParBinIdx(rowparmax) + 1, m_parBins[rowParIdx()] ); + + return {lowbin, highbin}; +} + +// Default implementations +double FPGATrackSimGenScanBinningBase::sliceVar([[maybe_unused]] FPGATrackSimHit const *hit) const { return 0.0; }; +double FPGATrackSimGenScanBinningBase::sliceVarExpected([[maybe_unused]] const ParSet &pars, [[maybe_unused]] FPGATrackSimHit const *hit) const {return 0.0; }; +double FPGATrackSimGenScanBinningBase::rowPar([[maybe_unused]] const ParSet &pars, [[maybe_unused]] FPGATrackSimHit const *hit) const {return 0.0; }; + + +//-------------------------------------------------------------------------------------------------- +// +// Geometry Helpers +// +//-------------------------------------------------------------------------------------------------- + +double FPGATrackSimGenScanGeomHelpers::ThetaFromEta(double eta) +{ + return 2.0 * atan(exp(-1.0 * eta)); +} + +double FPGATrackSimGenScanGeomHelpers::EtaFromTheta(double theta) +{ + return -log(tan(theta / 2.0)); +} + +double FPGATrackSimGenScanGeomHelpers::zFromPars(double r, const FPGATrackSimTrackPars &pars) +{ + double theta = ThetaFromEta(pars.eta); + double zhit = pars.z0 + r / tan(theta); + if (std::abs(pars.qOverPt) > 0) + { + zhit = pars.z0 + 1.0 / tan(theta) * asin(r * CurvatureConstant * pars.qOverPt) / (CurvatureConstant * pars.qOverPt); + } + return zhit; +} + +double FPGATrackSimGenScanGeomHelpers::phiFromPars(double r, const FPGATrackSimTrackPars &pars) +{ + double phi_hit = xAOD::P4Helpers::deltaPhi(pars.phi,asin(r * CurvatureConstant * pars.qOverPt - pars.d0 / r)); + + return phi_hit; +} + +double FPGATrackSimGenScanGeomHelpers::parsToTrkPhi(const FPGATrackSimTrackPars &pars, FPGATrackSimHit const *hit) +{ + double r = hit->getR(); // mm + double phi_hit = hit->getGPhi(); // radians + double phi_trk = xAOD::P4Helpers::deltaPhi(phi_hit,asin(r * CurvatureConstant * pars.qOverPt - pars.d0 / r)); + return phi_trk; +} + +// +// +// Tool for doing Key Layer Math +// +// + +std::pair<double, double> FPGATrackSimGenScanKeyLyrHelper::getXY(double r, double phi) const +{ + return std::pair<double, double>(r * cos(phi), r * sin(phi)); +} + +std::pair<double, double> FPGATrackSimGenScanKeyLyrHelper::getRotationAngles(const std::pair<double, double>& xy1, const std::pair<double, double>& xy2) const +{ + double dr12x = xy2.first - xy1.first; + double dr12y = xy2.second - xy1.second; + double dr12mag = std::hypot(dr12x, dr12y); + double cos_rot = dr12y / dr12mag; + double sin_rot = -dr12x / dr12mag; + return std::pair<double, double>(cos_rot, sin_rot); +} + +std::pair<double, double> FPGATrackSimGenScanKeyLyrHelper::rotateXY(const std::pair<double, double>& xy, const std::pair<double, double>& ang) const +{ + return std::pair<double, double>(xy.first * ang.first + xy.second * ang.second, + -xy.first * ang.second + xy.second * ang.first); +} + +FPGATrackSimGenScanKeyLyrHelper::rotatedConfig FPGATrackSimGenScanKeyLyrHelper::getRotatedConfig(const KeyLyrPars &keypars) const +{ + rotatedConfig result; + + auto xy1 = getXY(m_R1, keypars.phi1); + auto xy2 = getXY(m_R2, keypars.phi2); + + result.rotang = getRotationAngles(xy1, xy2); + result.xy1p = rotateXY(xy1, result.rotang); + result.xy2p = rotateXY(xy2, result.rotang); + result.y = result.xy2p.second - result.xy1p.second; + + return result; +} + +std::pair<double, double> FPGATrackSimGenScanKeyLyrHelper::getRotatedHit(const std::pair<double, double>& rotang, const FPGATrackSimHit *hit) const +{ + auto xyh = getXY(hit->getR(), hit->getGPhi()); + return rotateXY(xyh, rotang); +} + +FPGATrackSimGenScanKeyLyrHelper::KeyLyrPars FPGATrackSimGenScanKeyLyrHelper::trackParsToKeyPars(const FPGATrackSimTrackPars &pars) const +{ + KeyLyrPars retv; + retv.z1 = FPGATrackSimGenScanGeomHelpers::zFromPars(m_R1, pars); + retv.z2 = FPGATrackSimGenScanGeomHelpers::zFromPars(m_R2, pars); + retv.phi1 = FPGATrackSimGenScanGeomHelpers::phiFromPars(m_R1, pars); + retv.phi2 = FPGATrackSimGenScanGeomHelpers::phiFromPars(m_R2, pars); + + double Rinv = pars[FPGATrackSimTrackPars::IHIP] * (2.0 * FPGATrackSimGenScanGeomHelpers::CurvatureConstant); + if (Rinv != 0) + { + double R = 1 / Rinv; + auto xy1 = getXY(m_R1, retv.phi1); + auto xy2 = getXY(m_R2, retv.phi2); + double ysqr = (xy2.first - xy1.first) * (xy2.first - xy1.first) + (xy2.second - xy1.second) * (xy2.second - xy1.second); + double sign = (R > 0) - (R < 0); + retv.xm = -1 * sign * (std::abs(R) - sqrt(R * R - ysqr / 4.0)); + } + else + { + retv.xm = 0; + } + + return retv; +} + +FPGATrackSimTrackPars FPGATrackSimGenScanKeyLyrHelper::keyParsToTrackPars(const KeyLyrPars &keypars) const +{ + + FPGATrackSimTrackPars pars; + pars[FPGATrackSimTrackPars::IZ0] = (keypars.z1 * m_R2 - keypars.z2 * m_R1) / (m_R2 - m_R1); + pars[FPGATrackSimTrackPars::IETA] = FPGATrackSimGenScanGeomHelpers::EtaFromTheta(atan2(m_R2 - m_R1, keypars.z2 - keypars.z1)); + + // This is the exact math which is a bit messy + // if you want the math contact lipeles@sas.upenn.edu + + double xm = keypars.xm; + + auto rotated_coords = getRotatedConfig(keypars); + auto xy1p = rotated_coords.xy1p; + auto rotang = rotated_coords.rotang; + auto y = rotated_coords.y; + + // reverse rotation + auto revang = rotang; + revang.second = -rotang.second; + + if (xm != 0) + { + double Rinv = 2 * xm / (xm * xm + (y / 2) * (y / 2)); + double d = (y * y / 4.0 - xm * xm) / (2.0 * xm); + double sign = (xm > 0) - (xm < 0); + + pars[FPGATrackSimTrackPars::IHIP] = -1 * Rinv / (2.0 * FPGATrackSimGenScanGeomHelpers::CurvatureConstant); + + std::pair<double, double> xycp(-d + xy1p.first, y / 2.0 + xy1p.second); + + pars[FPGATrackSimTrackPars::ID0] = -1 * sign * (std::abs(1 / Rinv) - std::hypot(xycp.first, xycp.second)); + + auto xyc = rotateXY(xycp, revang); + + pars[FPGATrackSimTrackPars::IPHI] = atan2(sign * -xyc.first, sign * xyc.second); + } + else + { + pars[FPGATrackSimTrackPars::IHIP] = 0.0; + pars[FPGATrackSimTrackPars::ID0] = -1 * xy1p.first; + pars[FPGATrackSimTrackPars::IPHI] = atan2(rotang.first, -rotang.second); + } + + return pars; +} + +double FPGATrackSimGenScanKeyLyrHelper::zExpected(KeyLyrPars keypars, double r) const +{ + return (keypars.z2 - keypars.z1) / (m_R2 - m_R1) * (r - m_R1) + keypars.z1; +} + +// takes hit position and calculated what xm should be for that hit +double FPGATrackSimGenScanKeyLyrHelper::xmForHit(KeyLyrPars keypars, const FPGATrackSimHit *hit) const +{ + auto rotated_coords = getRotatedConfig(keypars); + auto xy1p = rotated_coords.xy1p; + auto rotang = rotated_coords.rotang; + auto y = rotated_coords.y; + + auto xyhp = getRotatedHit(rotang, hit); + double xh = xyhp.first - xy1p.first; + double yh = xyhp.second - xy1p.second; + + // use taylor expanded xm calculation + double sign = ((yh > 0) && (yh < y)) ? 1 : -1; + if ((std::abs(yh) < std::abs(xh)) || (std::abs(y - yh) < std::abs(xh))) + { + return ((xh > 0) ? 1 : -1) * 100000; + } + + double xm_taylor = sign * y * y / (y * y - 4 * xh * xh - 4 * (yh - y / 2.0) * (yh - y / 2.0)) * xh; + + return xm_taylor; +} + +// Find shift from nominal track to hit in the "x" direction +double FPGATrackSimGenScanKeyLyrHelper::deltaX(KeyLyrPars keypars, const FPGATrackSimHit *hit) const +{ + auto rotated_coords = getRotatedConfig(keypars); + + auto rotang = rotated_coords.rotang; + + auto xyhp = getRotatedHit(rotang, hit); + + double xh = xyhp.first-rotated_coords.xy1p.first; + + return xh - xExpected(keypars,hit); +} + +// Find shift from nominal track to hit in the "x" direction +double FPGATrackSimGenScanKeyLyrHelper::xExpected(KeyLyrPars keypars, const FPGATrackSimHit *hit) const +{ + auto rotated_coords = getRotatedConfig(keypars); + + auto rotang = rotated_coords.rotang; + auto y = rotated_coords.y; + + auto xyhp = getRotatedHit(rotang, hit); + + double yh = xyhp.second-rotated_coords.xy1p.second; + + return keypars.xm * (keypars.xm * keypars.xm + yh * (y - yh)) / (keypars.xm * keypars.xm + (y / 2) * (y / 2)); +} + +double FPGATrackSimGenScanPhiSlicedKeyLyrBinning::phiResidual(const ParSet &parset, FPGATrackSimHit const *hit, [[maybe_unused]] bool debug) const +{ + return m_keylyrtool.deltaX(parSetToKeyPars(parset), hit); +} + + diff --git a/Trigger/EFTracking/FPGATrackSim/FPGATrackSimHough/src/FPGATrackSimGenScanBinning.h b/Trigger/EFTracking/FPGATrackSim/FPGATrackSimHough/src/FPGATrackSimGenScanBinning.h new file mode 100644 index 0000000000000000000000000000000000000000..93041b77fe4679ec751a20f6740384475a7ca0d8 --- /dev/null +++ b/Trigger/EFTracking/FPGATrackSim/FPGATrackSimHough/src/FPGATrackSimGenScanBinning.h @@ -0,0 +1,530 @@ +// Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration + +#ifndef FPGATrackSimGenScanBinning_H +#define FPGATrackSimGenScanBinning_H + +/** + * @file FPGATrackSimGenScanBinning.h + * @author Elliot Lipeles + * @date Sept 10th, 2024 + * @brief Binning Classes for GenScanTool + * + * Declarations in this file (there are a series of small classes): + * class FPGATrackSimGenScanBinningBase + * base class that defines interface to FPGATrackSimGenScanTool + * + * class FPGATrackSimGenScanGeomHelpers + * very small class to implement common track parameter calculations + * that are shared by clases + * + * class GenScanStdBinning : FPGATrackSimGenScanBinningBase + * implements standard pT, eta, phi, d0, z0 binning + * + * class GenScanKeyLyrTool + * implements key layer math for binning using where in phi and z a track + * crosses two "keylayers" plus the distance the track has bent from a + * straight line midway between the layers (called xm in the code, but is + * also known as sagitta). This common helper tool avoid code duplication. + * + * class GenScanKeyLyrBinning : FPGATrackSimGenScanBinningBase + * implements keylyr binning with slice=z values at two radii, scan=phi + * values at two radii, and row=xm (sagitta) + * + * class GenScanKeyLyrPhiSlicedBinning : FPGATrackSimGenScanBinningBase + * implements keylyr binning with slice=phi values at two radii and xm, + * scan =z values at two radii, and row=xm. The track rates will be the + * same as GenScanKeyLyrBinning, but the internal throughput rates are + * different and are needed to understand firmware requirements + * + * + * Overview of stucture: + * -- This is to be used by the FPGATrackSimGenScanTool to define the binning + * (read that header first) + * -- A binning is defined by inheriting from FPGATrackSimGenScanBinningBase + * -- There are 3 implementations currently in this file, see descriptions in + * "Declarations" section above + * -- In addition, there are two helper tools for the math for standard track + * parameters and keylayer track parameters + * + + * + * References: + * + */ + +#include "FPGATrackSimObjects/FPGATrackSimTrackPars.h" +#include "FPGATrackSimObjects/FPGATrackSimHit.h" +#include "FPGATrackSimObjects/FPGATrackSimConstants.h" + +#include "FourMomUtils/xAODP4Helpers.h" +#include "GaudiKernel/StatusCode.h" + +#include <string> +#include <vector> + +//------------------------------------------------------------------------------------------------------- +// Binning base class +// Nomenclature: +// slice = first larger region binning (z0-eta binning in typical Hough) +// scan, row are the independent and dependent variables in a standard Hough ( e.g. pT and phi_track) +// slicePar = variables used to slice into overlapping regions (e.g. z0,eta) +// sliceVar = variable that slicePars specify a valid range of (e.g. hit_z) +// scanPar = scan variables of Hough-like scan (e.g. pT,d0) +// rowPar = variable that scanPars specify a valid range of for a given hit (e.g. phi_track) +//------------------------------------------------------------------------------------------------------- +class FPGATrackSimGenScanBinningBase +{ +public: + typedef std::vector<double> ParSet; + typedef std::vector<unsigned int> IdxSet; + + //-------------------------------------------------------------------------------------------------- + // + // Virtual methods that are overloaded to define the binning + // + //-------------------------------------------------------------------------------------------------- + + // Specification of parameters + virtual const std::string &parNames(unsigned i) const = 0; + + // overloaded to define which parameters are slice, scan, and row + // they return vectors which are the numbers of the parameters in the slice or scan + // and just an unsigned for row since there can be only one row parmeter + virtual std::vector<unsigned> slicePars() const = 0; + virtual std::vector<unsigned> scanPars() const = 0; + virtual unsigned rowParIdx() const = 0; + + // convert back and forth from pT, eta, phi, d0, z0 and internal paramater set + virtual const ParSet trackParsToParSet(const FPGATrackSimTrackPars &pars) const = 0; + virtual const FPGATrackSimTrackPars parSetToTrackPars(const ParSet &parset) const = 0; + + // calculate the distance in phi or eta from a track defined by parset to a hit + // these can be implemented as any variable in the r-phi or r-eta plane (not necessarily eta and phi). + virtual double phiResidual(const ParSet &parset, FPGATrackSimHit const *hit, [[maybe_unused]] bool debug=false) const = 0; + virtual double etaResidual(const ParSet &parset, FPGATrackSimHit const *hit, [[maybe_unused]] bool debug=false) const = 0; + + // Used to scale the histograms binning appropriately (very approximate, but at least the default ranges come up reasonable) + virtual double phiHistScale() const = 0; + virtual double etaHistScale() const = 0; + + // To determine if a hit is in a slice or row there are two options + // 1) Directly implement hitInSlice and idxsetToRowParBinRange. + // 2) Use the default implementaions of those in which case sliceVar, sliceVarExpected, and rowPar should be implemented + // Tells if a hit is in a slice + virtual bool hitInSlice(const IdxSet &idx, FPGATrackSimHit const *hit) const; + + // Tells the range of row parameters a hit is conistent with assuming the other parameters (slice and scan) + // are set by the given bin idx + virtual std::pair<unsigned, unsigned> idxsetToRowParBinRange(const IdxSet &idx, FPGATrackSimHit const *hit) const; + + private: + // Hit variable used in slicing (e.g. for eta-z0 slicing it would be the z of the hit) + virtual double sliceVar([[maybe_unused]] FPGATrackSimHit const *hit) const; + + // Find expected sliceVar for a track with the given parameters (e.g. or eta-z0 slicing this is the expect z for a given hit radius) + virtual double sliceVarExpected([[maybe_unused]] const ParSet &pars, [[maybe_unused]] FPGATrackSimHit const *hit) const; + + // Find the parameter for a given hit and parameters (e.g. given pT and d0 what is track phi) + virtual double rowPar([[maybe_unused]] const ParSet &pars, [[maybe_unused]] FPGATrackSimHit const *hit) const; + + + //-------------------------------------------------------------------------------------------------- + // + // Functional Implementation + // + //-------------------------------------------------------------------------------------------------- + + public: + // merge of both slice and scan par lists + std::vector<unsigned> sliceAndScanPars() const { + std::vector<unsigned> retv = slicePars(); + std::vector<unsigned> scan = scanPars(); + retv.insert(retv.end(), scan.begin(), scan.end()); + return retv; + }; + + // Find distance from bin center to hit in either slice direction or the row direction + double phiShift(const IdxSet &idx, FPGATrackSimHit const *hit, bool debug=false) const { return phiResidual(binCenter(idx),hit,debug);} + double etaShift(const IdxSet &idx, FPGATrackSimHit const *hit, bool debug=false) const { return etaResidual(binCenter(idx),hit,debug);} + + // accessors to get the numbers of bins in just the slice, scan, sliceAndScan, or row (subset of the 5-d space) + virtual std::vector<unsigned> sliceBins() const { return subVec(slicePars(), m_parBins); } + virtual std::vector<unsigned> scanBins() const { return subVec(scanPars(), m_parBins); } + virtual std::vector<unsigned> sliceAndScanBins() const { return subVec(sliceAndScanPars(), m_parBins); } + virtual unsigned rowBin() const { return m_parBins[rowParIdx()]; } + + // extract just the slice, scan, sliceAndScan, or row part of a 5-d index + virtual std::vector<unsigned> sliceIdx(const IdxSet &idx) const { return subVec(slicePars(), idx); } + virtual std::vector<unsigned> scanIdx(const IdxSet &idx) const { return subVec(scanPars(), idx); } + virtual std::vector<unsigned> sliceAndScanIdx(const IdxSet &idx) const { return subVec(sliceAndScanPars(), idx); } + virtual unsigned rowIdx(const IdxSet &idx) const { return idx[rowParIdx()]; } + + // Bin boundary utilities + double binCenter(unsigned par, unsigned bin) const { return m_parMin[par] + m_parStep[par] * (double(bin) + 0.5); } + double binLowEdge(unsigned par, unsigned bin) const { return m_parMin[par] + m_parStep[par] * (double(bin)); } + double binHighEdge(unsigned par, unsigned bin) const { return m_parMin[par] + m_parStep[par] * (double(bin) + 1.0); } + ParSet binLowEdge(const IdxSet &idx) const; + ParSet binCenter(const IdxSet &idx) const; + + // center of whole region + ParSet center() const; + + // get bin value for a specific parameter value + unsigned binIdx(unsigned par, double val) const { return (val > m_parMin[par]) ? unsigned(floor((val - m_parMin[par]) / m_parStep[par])) : 0; } + unsigned rowParBinIdx(double val) const { return binIdx(rowParIdx(), val); } + + // check if 1-d or 5-d parameter is within the range of the binning + bool inRange(unsigned par, double val) const { return ((val < m_parMax[par]) && (val > m_parMin[par])); } + bool inRange(const ParSet &pars) const; + + // convert parset (the binning parameters) to a 5-d bin + IdxSet binIdx(const ParSet &pars) const; + + // convert FPGATrackSimTrackPars to 5-d bin + const IdxSet parsToBin(FPGATrackSimTrackPars &pars) const; + + // Generic Utility for splitting in vector (e.g. idx or #bins 5-d vectors) + // into subvectors (e.g. idx for just the scan parameters). Technically, for + // a list of parameter indices (elems) gives the subvector of the invec with just those indices + std::vector<unsigned> subVec(const std::vector<unsigned>& elems, const std::vector<unsigned>& invec) const; + + // Opposite of above subVec, this sets the subvector + StatusCode setIdxSubVec(IdxSet &idx, const std::vector<unsigned>& subvecelems, const std::vector<unsigned>& subvecidx) const; + + // Makes are set of parameters corresponding to the corners specified by scanpars of the bin specified by idx + // e.g. if scan pars is (pT,d0) then the set is (low pT,low d0), (low pT, high d0), (high pT,low d0), (high pT, high d0) + std::vector<ParSet> makeVariationSet(const std::vector<unsigned> &scanpars, const IdxSet &idx) const; + + // + // Internal data + // + static constexpr unsigned NPars = 5; + std::vector<double> m_parMin; + std::vector<double> m_parMax; + std::vector<double> m_parStep; + std::vector<unsigned> m_parBins; + + // invalid bin value, there is no way a true bin could be there + const std::vector<unsigned> m_invalidBin = {std::numeric_limits<unsigned>::max(),std::numeric_limits<unsigned>::max(), + std::numeric_limits<unsigned>::max(),std::numeric_limits<unsigned>::max(),std::numeric_limits<unsigned>::max()}; +}; + +//------------------------------------------------------------------------------------------------------- +// +// Geometry Helpers -- does basic helix calculations +// +//------------------------------------------------------------------------------------------------------- +class FPGATrackSimGenScanGeomHelpers +{ +public: + // This is the constant needed to relate hit phi to track phi due to curvature + static constexpr double CurvatureConstant = fpgatracksim::A; + + // standard eta to theta calculation + static double ThetaFromEta(double eta); + + // standard theta to eta calculation + static double EtaFromTheta(double theta); + + // find the expected z position from the radius and track parameters + static double zFromPars(double r, const FPGATrackSimTrackPars &pars); + + // find the expected z position from the radius and track parameters + static double phiFromPars(double r, const FPGATrackSimTrackPars &pars); + + // find the track phi that would be consistent with the other track parameters and the hit (r,phi) + static double parsToTrkPhi(const FPGATrackSimTrackPars &pars, FPGATrackSimHit const *hit); +}; + +//------------------------------------------------------------------------------------------------------- +// +// Standard pT, d0, phi, eta, z0 implmentation of FPGATrackSimGenScanStdTrkBinning +// -- this is then effectively a "standard" Hough transform except that it slices +// in (z0,eta) and then the accumulator scans is in both pT and d0 with phi_track +// as the row of the accumulator +// +//------------------------------------------------------------------------------------------------------- +class FPGATrackSimGenScanStdTrkBinning : public FPGATrackSimGenScanBinningBase +{ +public: + FPGATrackSimGenScanStdTrkBinning() : m_parNames({"z0", "eta", "qOverPt", "d0", "phi"}) {} + virtual const std::string &parNames(unsigned i) const override { return m_parNames[i]; } + virtual unsigned rowParIdx() const override { return 4;} + + virtual std::vector<unsigned> slicePars() const override { return std::vector<unsigned>({0,1}); } + virtual std::vector<unsigned> scanPars() const override { return std::vector<unsigned>({2,3}); } + + virtual double etaHistScale() const {return 200.0;} + virtual double phiHistScale() const {return 0.1;} + + virtual const ParSet trackParsToParSet(const FPGATrackSimTrackPars &pars) const override + { + return std::vector<double>({pars[FPGATrackSimTrackPars::IZ0], pars[FPGATrackSimTrackPars::IETA], + pars[FPGATrackSimTrackPars::IHIP], pars[FPGATrackSimTrackPars::ID0], pars[FPGATrackSimTrackPars::IPHI]}); + } + virtual const FPGATrackSimTrackPars parSetToTrackPars(const ParSet &parset) const override + { + FPGATrackSimTrackPars pars; + pars[FPGATrackSimTrackPars::IZ0] = parset[0]; + pars[FPGATrackSimTrackPars::IETA] = parset[1]; + pars[FPGATrackSimTrackPars::IHIP] = parset[2]; + pars[FPGATrackSimTrackPars::ID0] = parset[3]; + pars[FPGATrackSimTrackPars::IPHI] = parset[4]; + return pars; + } + + virtual double sliceVarExpected(const ParSet &pars, FPGATrackSimHit const *hit) const override + { + return FPGATrackSimGenScanGeomHelpers::zFromPars(hit->getR(),parSetToTrackPars(pars)); + } + virtual double sliceVar(FPGATrackSimHit const *hit) const override + { + return hit->getZ(); + } + virtual double rowPar(const ParSet &pars, FPGATrackSimHit const *hit) const override + { + return FPGATrackSimGenScanGeomHelpers::parsToTrkPhi(parSetToTrackPars(pars),hit); + } + + virtual double etaResidual(const ParSet &parset ,FPGATrackSimHit const * hit, [[maybe_unused]] bool debug=false) const { + // this uses a shift in "eta*radius" instead of "z" + double theta = FPGATrackSimGenScanGeomHelpers::ThetaFromEta(parset[1]); // 1 = eta + return (hit->getZ() - parset[0]) * sin(theta) - hit->getR() * cos(theta); // 0=z0 + } + virtual double phiResidual(const ParSet &parset,FPGATrackSimHit const *hit, [[maybe_unused]] bool debug=false) const { + return xAOD::P4Helpers::deltaPhi(hit->getGPhi(), FPGATrackSimGenScanGeomHelpers::phiFromPars(hit->getR(),parSetToTrackPars(parset))); + } + +private: + const std::vector<std::string> m_parNames; +}; + +//------------------------------------------------------------------------------------------------------- +// +// Tool for doing Key Layer Math +// -- that is converting be the track parameters set by the standard (pT,eta,phi,d0,z0) and +// the "KeyLayer" parametes set by the phi and z at to radii (R1,R2) and a deviation from +// straight line between the two points (xm) +// -- All the math for the r-phi plane is based finding the rotated coordinate system where +// the points (R1,phiR1) and (R2,phiR2) both lie on the x=0 axis. Then the x of a track halfway +// between the points (called xm) is the sagitta. +// +//------------------------------------------------------------------------------------------------------- +class FPGATrackSimGenScanKeyLyrHelper +{ +public: + FPGATrackSimGenScanKeyLyrHelper(double r1, double r2) : m_R1(r1), m_R2(r2) {} + + struct KeyLyrPars { + double z1; + double z2; + double phi1; + double phi2; + double xm; + }; + + // convert (r,phi) to (x,y) + std::pair<double, double> getXY(double r, double phi) const; + + // Get rotation angles need to rotate the two given points to be only seperated in the y direction (both have same x) + // results is the sin and cos of rotation + std::pair<double, double> getRotationAngles(const std::pair<double, double>& xy1, const std::pair<double, double>& xy2) const; + + // Apply a rotation angles (sin and cos) to a point xy + std::pair<double, double> rotateXY(const std::pair<double, double>& xy, const std::pair<double, double>& ang) const; + + // Simple struct and function to get the rotation angle and full rotated information in one call + struct rotatedConfig + { + std::pair<double, double> xy1p; // point 1 after rotation + std::pair<double, double> xy2p; // point 2 after rotation + std::pair<double, double> rotang; // rotation angle + double y; // y seperation after rotation + }; + rotatedConfig getRotatedConfig(const KeyLyrPars &keypars) const; + + // get the x,y coordinates of a hit in the rotated coordinate system specified by rotang + std::pair<double, double> getRotatedHit(const std::pair<double, double>& rotang, const FPGATrackSimHit *hit) const; + + // Convert back and forth to standard track parameters + KeyLyrPars trackParsToKeyPars(const FPGATrackSimTrackPars &pars) const; + FPGATrackSimTrackPars keyParsToTrackPars(const KeyLyrPars &keypars) const; + + // find expected z hit position given a radius + double zExpected(KeyLyrPars keypars, double r) const; + + // find expected x position of a hit at given a radius in the rotated coordinate system + double xExpected(KeyLyrPars keypars, const FPGATrackSimHit *hit) const; + + // takes hit position and calculated what xm would be for track going through that hit + double xmForHit(KeyLyrPars keypars, const FPGATrackSimHit *hit) const; + + // Find shift from nominal track to hit in the "x" direction + double deltaX(KeyLyrPars keypars, const FPGATrackSimHit *hit) const; + + // accessors + double R1() const {return m_R1;} + double R2() const {return m_R2;} + + private: + double m_R1; + double m_R2; +}; + +//------------------------------------------------------------------------------------------------------- +// +// Key Layer Binning implementation (z slice then phi scan) +// -- Using the FPGATrackSimGenScanKeyLyrHelper to implement an instance of +// FPGATrackSimGenScanBinningBase which slices in z and then scans in (phiR1,phiR2), +// and then uses xm as the final row variable +// +//------------------------------------------------------------------------------------------------------- +class FPGATrackSimGenScanKeyLyrBinning : public FPGATrackSimGenScanBinningBase +{ +public: + + FPGATrackSimGenScanKeyLyrBinning(double r_in, double r_out) : + m_keylyrtool(r_in,r_out), m_parNames({"zR1", "zR2", "phiR1", "phiR2", "xm"}) {} + virtual const std::string &parNames(unsigned i) const override { return m_parNames[i]; } + virtual unsigned rowParIdx() const override { return 4;} + virtual std::vector<unsigned> slicePars() const override { return std::vector<unsigned>({0,1}); } + virtual std::vector<unsigned> scanPars() const override { return std::vector<unsigned>({2,3}); } + + virtual double etaHistScale() const {return 300.0;} + virtual double phiHistScale() const {return 60.0;} + + ParSet keyparsToParSet(const FPGATrackSimGenScanKeyLyrHelper::KeyLyrPars& keypars) const { + return ParSet({keypars.z1,keypars.z2,keypars.phi1,keypars.phi2,keypars.xm}); + } + + FPGATrackSimGenScanKeyLyrHelper::KeyLyrPars parSetToKeyPars(const ParSet &parset) const { + FPGATrackSimGenScanKeyLyrHelper::KeyLyrPars keypars; + keypars.z1 = parset[0]; + keypars.z2 = parset[1]; + keypars.phi1 = parset[2]; + keypars.phi2 = parset[3]; + keypars.xm = parset[4]; + return keypars; + } + + virtual const ParSet trackParsToParSet(const FPGATrackSimTrackPars &pars) const override { + return keyparsToParSet(m_keylyrtool.trackParsToKeyPars(pars)); + } + virtual const FPGATrackSimTrackPars parSetToTrackPars(const ParSet &parset) const override { + return m_keylyrtool.keyParsToTrackPars(parSetToKeyPars(parset)); + } + virtual double sliceVarExpected(const ParSet &pars, FPGATrackSimHit const *hit) const override + { + return m_keylyrtool.zExpected(parSetToKeyPars(pars),hit->getR()); + } + virtual double sliceVar(FPGATrackSimHit const *hit) const override { + return hit->getZ(); + } + virtual double rowPar(const ParSet &parset, FPGATrackSimHit const *hit) const override + { + return m_keylyrtool.xmForHit(parSetToKeyPars(parset),hit); + } + virtual double etaResidual(const ParSet &parset,FPGATrackSimHit const * hit, [[maybe_unused]] bool debug=false) const override { + return hit->getZ()- m_keylyrtool.zExpected(parSetToKeyPars(parset),hit->getR()); + } + + virtual double phiResidual(const ParSet &parset, FPGATrackSimHit const *hit, [[maybe_unused]] bool debug=false) const override { + return m_keylyrtool.deltaX(parSetToKeyPars(parset),hit); + } + +private: + FPGATrackSimGenScanKeyLyrHelper m_keylyrtool; + const std::vector<std::string> m_parNames; + + }; + +//------------------------------------------------------------------------------------------------------- +// +// Key Layer Binning Phi Sliced (then z scanned) +// -- This scans slices in all phi parameters ("phiR1", "phiR2", "xm") +// and then scans in in ("zR1", "zR2") +// -- This was used to explore the throughput of inverting these opperations +// -- Also the hitInSlice and idxsetToRowParBinRange functions were switched +// to be simple windows so that its more like real firmware +// +//------------------------------------------------------------------------------------------------------- +class FPGATrackSimGenScanPhiSlicedKeyLyrBinning : public FPGATrackSimGenScanBinningBase +{ +public: + + FPGATrackSimGenScanPhiSlicedKeyLyrBinning(double r_in, double r_out) : + m_keylyrtool(r_in,r_out), m_parNames({"zR1", "zR2", "phiR1", "phiR2", "xm"}) {} + virtual const std::string &parNames(unsigned i) const override { return m_parNames[i]; } + virtual unsigned rowParIdx() const override { return 4;} + virtual std::vector<unsigned> slicePars() const override { return std::vector<unsigned>({2,3,4}); } + virtual std::vector<unsigned> scanPars() const override { return std::vector<unsigned>({0,1}); } + + virtual double etaHistScale() const {return 60.0;} + virtual double phiHistScale() const {return 30.0;} + + ParSet keyparsToParSet(const FPGATrackSimGenScanKeyLyrHelper::KeyLyrPars& keypars) const { + return ParSet({keypars.z1,keypars.z2,keypars.phi1,keypars.phi2,keypars.xm}); + } + + FPGATrackSimGenScanKeyLyrHelper::KeyLyrPars parSetToKeyPars(const ParSet &parset) const { + FPGATrackSimGenScanKeyLyrHelper::KeyLyrPars keypars; + keypars.z1 = parset[0]; + keypars.z2 = parset[1]; + keypars.phi1 = parset[2]; + keypars.phi2 = parset[3]; + keypars.xm = parset[4]; + return keypars; + } + + virtual const ParSet trackParsToParSet(const FPGATrackSimTrackPars &pars) const override { + return keyparsToParSet(m_keylyrtool.trackParsToKeyPars(pars)); + } + virtual const FPGATrackSimTrackPars parSetToTrackPars(const ParSet &parset) const override { + return m_keylyrtool.keyParsToTrackPars(parSetToKeyPars(parset)); + } + + virtual double phiResidual(const ParSet &parset, FPGATrackSimHit const *hit, [[maybe_unused]] bool debug) const override; + + virtual double etaResidual(const ParSet &parset, FPGATrackSimHit const *hit, [[maybe_unused]] bool debug) const override + { + return hit->getZ()- m_keylyrtool.zExpected(parSetToKeyPars(parset),hit->getR()); + } + + virtual bool hitInSlice(const IdxSet &idx, FPGATrackSimHit const *hit) const override { + double r1 = m_keylyrtool.R1(); + double r2 = m_keylyrtool.R2(); + auto keypars = parSetToKeyPars(binCenter(idx)); + auto tmppars = keypars; + tmppars.xm = m_parStep[4] / 2.0; + double xrange = m_keylyrtool.xExpected(tmppars, hit) + (r1*m_parStep[2] + (r2*m_parStep[3] - r1*m_parStep[2]) / (r2 - r1) * (hit->getR() - r1))/2.0; + return (std::abs(phiShift(idx, hit)) < xrange); + } + + virtual std::pair<unsigned, unsigned> idxsetToRowParBinRange(const IdxSet &idx, [[maybe_unused]] FPGATrackSimHit const *hit) const override + { + double r1 = m_keylyrtool.R1(); + double r2 = m_keylyrtool.R2(); + double lowz_in = binLowEdge(0,idx[0]); + double highz_in = binHighEdge(0,idx[0]); + double lowz_out = binLowEdge(1,idx[1]); + double highz_out = binHighEdge(1,idx[1]); + + + // simple box cut + if (hit->getZ() < lowz_in + (lowz_out-lowz_in) * (hit->getR()-r1)/(r2-r1)) + return std::pair<unsigned, unsigned>(0, 0); // empty range + + if (hit->getZ() > highz_in + (highz_out-highz_in) * (hit->getR()-r1)/(r2-r1)) + return std::pair<unsigned, unsigned>(0, 0); // empty range + + return std::pair<unsigned, unsigned>(rowIdx(idx), rowIdx(idx) + 1); // range covers just 1 row bin + } + +private: + FPGATrackSimGenScanKeyLyrHelper m_keylyrtool; + const std::vector<std::string> m_parNames; + + }; + + +#endif // FPGATrackSimGenScanBinning_H