GaussianTrackDensity.h

Go to the documentation of this file.

/*
  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
*/
 
#ifndef TRKVERTEXSEEDFINDERUTILIS_GAUSSIANTRACKDENSITY_H
#define TRKVERTEXSEEDFINDERUTILIS_GAUSSIANTRACKDENSITY_H
 
#include "AthenaBaseComps/AthAlgTool.h"
#include "TrkVertexFitterInterfaces/IVertexTrackDensityEstimator.h"
 
#include <map>
 
#include "AthAllocators/ArenaPoolSTLAllocator.h"
namespace Trk
{
 
  class Track;
  class GaussianTrackDensity;

 
  class GaussianTrackDensity final: public extends<AthAlgTool, IVertexTrackDensityEstimator>
  {
  public:
    using base_class::base_class;

    virtual double
    globalMaximum (const std::vector<const Track*>& vectorTrk) const override final;

    virtual double
    globalMaximum (const std::vector<const Track*>& vectorTrk,
                   std::unique_ptr<ITrackDensity>& density) const override final;

    virtual double
    globalMaximum (const std::vector<const TrackParameters*>& perigeeList) const override final;

    virtual double
    globalMaximum (const std::vector<const TrackParameters*>& perigeeList,
                   std::unique_ptr<ITrackDensity>& density) const override final;
 
    virtual std::pair<double, double> globalMaximumWithWidth(
        const std::vector<const TrackParameters*>& perigeeList)
        const override final;
 
   private:
    struct TrackEntry {
      TrackEntry() = default;
      TrackEntry(const TrackEntry&) = default;
      TrackEntry(TrackEntry&&) = default;
      TrackEntry& operator=(const TrackEntry&) = default;
      TrackEntry& operator=(TrackEntry&&) = default;
      ~TrackEntry() = default;
 
      TrackEntry(double c0, double c1, double c2, double zMin, double zMax)
          : c_0(c0), c_1(c1), c_2(c2), lowerBound(zMin), upperBound(zMax) {}
      explicit TrackEntry(double z)
          : c_0(0), c_1(0), c_2(0), lowerBound(z), upperBound(z) {}
      // Cached information for a single track
      double c_0 = 0;  // z-independent term in exponent
      double c_1 = 1;  // linear coefficient in exponent
      double c_2 = 0;  // quadratic coefficient in exponent
      double lowerBound = 0;
      double upperBound = 0;
    };
 
    // helper to handle the evaluation of the parametrised track density 
    struct TrackDensityEval{
      public:
        // initialise with the z coordinate at which the density is to be evaluated 
        TrackDensityEval(double z_coordinate): m_z(z_coordinate){} 
        // add the contribution for one track to the density. 
        // will internally check if the z coordinate is within 
        // the bounds of the track 
        void addTrack (const TrackEntry & entry);
        // retrieve the density and its derivatives
        inline double density() const {return m_density;}
        inline double firstDerivative() const {return m_firstDerivative;}
        inline double secondDerivative() const {return m_secondDerivative;}
      private: 
        double m_z; 
        double m_density{0};
        double m_firstDerivative{0};
        double m_secondDerivative{0}; 
    };
 
    class TrackDensity final : public ITrackDensity {
       public:
        explicit TrackDensity(bool gaussStep) : m_gaussStep(gaussStep) {}
        virtual ~TrackDensity() = default;

        virtual double trackDensity(double z) const override final;

        virtual void trackDensity(
            double z, double& density, double& firstDerivative,
            double& secondDerivative) const override final;
        double globalMaximum() const;
        std::pair<double, double> globalMaximumWithWidth() const;
        void addTrack(const Perigee& itrk, const double d0SignificanceCut,
                      const double z0SignificanceCut);
 
       private:
        inline void updateMaximum(double trialZ, double trialValue,
                                  double secondDerivative, double& maxZ,
                                  double& maxValue,
                                  double& maxSecondDerivative) const {
        if (trialValue > maxValue) {
          maxZ = trialZ;
          maxValue = trialValue;
          maxSecondDerivative = secondDerivative;
        }
        }
 
        inline double stepSize(double y, double dy, double ddy) const {
        return (m_gaussStep ? (y * dy) / (dy * dy - y * ddy) : -dy / ddy);
        }
        bool m_gaussStep;
        double m_maxRange = 0;
        //  Cache for track information
        // functor to compare two Perigee values
        struct pred_perigee {
        inline bool operator()(const Perigee& left,
                               const Perigee& right) const {
          return left.parameters()[Trk::z0] < right.parameters()[Trk::z0];
        }
        };
        // functor to compare two TrackEntry values based on their upper limits (low
        // to high)
        struct pred_entry_by_max {
        inline bool operator()(const TrackEntry& left,
                               const TrackEntry& right) const {
          return left.upperBound < right.upperBound;
        }
        };
 
        using trackMap =
            std::map<Perigee, GaussianTrackDensity::TrackEntry, pred_perigee,
                     SG::ArenaPoolSTLAllocator<std::pair<
                         const Perigee, GaussianTrackDensity::TrackEntry>>>;
 
        using lowerMap = std::map<
            GaussianTrackDensity::TrackEntry, Perigee, pred_entry_by_max,
            SG::ArenaPoolSTLAllocator<
                std::pair<const GaussianTrackDensity::TrackEntry, Perigee>>>;
 
        trackMap m_trackMap;
        lowerMap m_lowerMap;
    };

    double
    globalMaximumImpl (const std::vector<const TrackParameters*>& perigeeList,
                       TrackDensity& density) const;

    std::pair<double,double>
    globalMaximumWithWidthImpl (const std::vector<const TrackParameters*>& perigeeList,
                       TrackDensity& density) const;
 

    void addTracks(const std::vector<const TrackParameters*>& perigeeList,
                   TrackDensity& density) const;
 
 
 
    //  Cuts set by configurable properties
    
    //  Maximum allowed d0 significance to use (in sigma)
    Gaudi::Property<double> m_d0MaxSignificance{
      this,
      "MaxD0Significance",
      3.5,
      "Maximum radial impact parameter significance to use track"
    };
 
    // Tracks within this many sigma(z) are added to weight; increasing cut
    // trades CPU efficiency for improved smoothness in tails
    Gaudi::Property<double> m_z0MaxSignificance{
      this,
      "MaxZ0Significance",
      12.0,
      "Maximum longitudinal impact parameter significance to include track in "
      "weight"
    };
 
    // Assumed shape of density function near peak; Gaussian (true) or parabolic
    // (false)
    Gaudi::Property<bool> m_gaussStep{
      this,
      "GaussianStep",
      true,
      "Peak search: True means assume Gaussian behavior, False means "
      "Newton/parabolic"
    };
  };
}
#endif