HLTriggerOffline/Scouting/plugins/ScoutingTrackMonitor.cc - master - githubmirror/cmssw

HLTriggerOffline/Scouting/plugins/ScoutingTrackMonitor.cc
1 057 строк43 KB
Marco Musich
miscellaneous improvements to ScoutingTrackMonitor
31 май 2026, 17:01
31 май 2026, 17:01907459d
// system includes
#include <cmath>
#include <vector>
#include <numbers>
#include <fmt/format.h>
#include <boost/range/adaptor/indexed.hpp>

// ROOT includes
#include "TMath.h"

// user includes
#include "DQMServices/Core/interface/DQMEDAnalyzer.h"
#include "DQMServices/Core/interface/MonitorElement.h"
#include "DataFormats/Scouting/interface/Run3ScoutingTrack.h"
#include "DataFormats/Scouting/interface/Run3ScoutingVertex.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
#include "FWCore/Utilities/interface/InputTag.h"

namespace sctTrackMonitor {
  // same logic used for the MTV:
  // cf https://github.com/cms-sw/cmssw/blob/master/Validation/RecoTrack/src/MTVHistoProducerAlgoForTracker.cc
  typedef dqm::reco::DQMStore DQMStore;

  inline void setBinLog(TAxis* axis) {
    int bins = axis->GetNbins();
    float from = axis->GetXmin();
    float to = axis->GetXmax();
    float width = (to - from) / bins;
    std::vector<float> new_bins(bins + 1, 0);
    for (int i = 0; i <= bins; i++) {
      new_bins[i] = TMath::Power(10, from + i * width);
    }
    axis->Set(bins, new_bins.data());
  }

  inline void setBinLogX(TH1* h) {
    TAxis* axis = h->GetXaxis();
    setBinLog(axis);
  }
  inline void setBinLogY(TH1* h) {
    TAxis* axis = h->GetYaxis();
    setBinLog(axis);
  }

  template <typename... Args>
  dqm::reco::MonitorElement* makeProfileIfLog(DQMStore::IBooker& ibook, bool logx, bool logy, Args&&... args) {
    auto prof = std::make_unique<TProfile>(std::forward<Args>(args)...);
    if (logx)
      setBinLogX(prof.get());
    if (logy)
      setBinLogY(prof.get());
    const auto& name = prof->GetName();
    return ibook.bookProfile(name, prof.release());
  }

  template <typename... Args>
  dqm::reco::MonitorElement* makeTH1IfLog(DQMStore::IBooker& ibook, bool logx, bool logy, Args&&... args) {
    auto h1 = std::make_unique<TH1F>(std::forward<Args>(args)...);
    if (logx)
      setBinLogX(h1.get());
    if (logy)
      setBinLogY(h1.get());
    const auto& name = h1->GetName();
    return ibook.book1D(name, h1.release());
  }
}  // namespace sctTrackMonitor

class ScoutingTrackMonitor : public DQMEDAnalyzer {
public:
  explicit ScoutingTrackMonitor(const edm::ParameterSet&);
  ~ScoutingTrackMonitor() override = default;
  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

  struct IPMonitoring {
    std::string varname_;
    float pTcut_;
    dqm::reco::MonitorElement *IP_, *IPErr_, *IPPull_;
    dqm::reco::MonitorElement *IPVsPhi_, *IPVsEta_, *IPVsPt_;
    dqm::reco::MonitorElement *IPErrVsPhi_, *IPErrVsEta_, *IPErrVsPt_;
    dqm::reco::MonitorElement *IPVsEtaVsPhi_, *IPErrVsEtaVsPhi_;

    void bookIPMonitor(DQMStore::IBooker&, const edm::ParameterSet&);

  private:
    int PhiBin_, EtaBin_, PtBin_;
    double PhiMin_, PhiMax_, EtaMin_, EtaMax_, PtMin_, PtMax_;
  };

  struct ProfileConfig {
    std::string name;   // Base name
    std::string title;  // Human-readable title
    double ymin;        // Minimum Y value
    double ymax;        // Maximum Y value

    // MonitorElements
    MonitorElement* p2_eta_phi = nullptr;
    MonitorElement* p_eta = nullptr;
    MonitorElement* p_phi = nullptr;
  };

protected:
  void analyze(const edm::Event&, const edm::EventSetup&) override;
  void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;

private:
  static inline std::pair<float, float> trk_vtx_offSet(const Run3ScoutingTrack& tk, const Run3ScoutingVertex& vtx) {
    const auto pt = tk.tk_pt();
    const auto phi = tk.tk_phi();
    const auto eta = tk.tk_eta();

    const auto px = pt * std::cos(phi);
    const auto py = pt * std::sin(phi);
    const auto pz = pt * std::sinh(eta);
    const auto pt2 = pt * pt;

    const auto dx = tk.tk_vx() - vtx.x();
    const auto dy = tk.tk_vy() - vtx.y();
    const auto dz = tk.tk_vz() - vtx.z();

    const auto tk_dxyPV = (-dx * py + dy * px) / pt;
    const auto tk_dzPV = dz - (dx * px + dy * py) * pz / pt2;

    return {tk_dxyPV, tk_dzPV};
  }

  // configuration
  const edm::ParameterSet conf_;

  // tokens
  const edm::EDGetTokenT<std::vector<Run3ScoutingTrack>> tracksToken_;
  const edm::EDGetTokenT<std::vector<Run3ScoutingVertex>> verticesToken_;
  const edm::EDGetTokenT<reco::BeamSpot> beamSpotToken_;

  const std::string topFolderName_;  // top folder name where to book histograms

  // histograms
  MonitorElement* h_dxy;
  MonitorElement* h_dz;
  MonitorElement* h_vtx_idx;
  MonitorElement* h2_eta_phi;

  MonitorElement* p_dxy_eta;
  MonitorElement* p_dxy_phi;
  MonitorElement* p2_dxy_eta_phi;

  MonitorElement* p_dz_eta;
  MonitorElement* p_dz_phi;
  MonitorElement* p2_dz_eta_phi;

  MonitorElement* h_vr;  // production radius for all tracks
  MonitorElement* h_vz;  // productio  z for all tracks

  // beamspot histograms
  MonitorElement *h_bsX, *h_bsY, *h_bsZ, *h_bsSigmaZ, *h_bsDxdz, *h_bsDydz, *h_bsBeamWidthX, *h_bsBeamWidthY, *h_bsType;

  // vertex histograms
  MonitorElement *h_vtx_chi2ndf, *h_vtx_prob;
  MonitorElement *h_vtx_sumPt2, *h_vtx_nTracks;

  // hit profiles configuration
  std::vector<ProfileConfig> hitProfiles_;
  std::vector<ProfileConfig> impactParameterProfiles_;

  static constexpr int cmToUm = 10000;

  // IP monitoring structs
  IPMonitoring dxy_pt1;
  IPMonitoring dxy_pt10;

  IPMonitoring dz_pt1;
  IPMonitoring dz_pt10;

  // profiles
  std::vector<MonitorElement*> vTrackProfiles_;

  // helpers
  reco::Track makeRecoTrack(const Run3ScoutingTrack& sTrack) const;
  reco::Vertex makeRecoVertex(const Run3ScoutingVertex& sVertex) const;
  std::pair<unsigned int, const Run3ScoutingVertex*> findClosestScoutingVertex(
      const reco::Track* track, const std::vector<Run3ScoutingVertex>& vertices);

  template <typename T>
  bool getValidHandle(const edm::Event& iEvent,
                      const edm::EDGetTokenT<T>& token,
                      edm::Handle<T>& handle,
                      const std::string& label);

  template <class OBJECT_TYPE>
  int index(const std::vector<OBJECT_TYPE*>& vec, const TString& name) {
    for (const auto& iter : vec | boost::adaptors::indexed(0)) {
      if (iter.value() && iter.value()->getName() == name) {
        return iter.index();
      }
    }
    edm::LogError("ScoutingTrackMonitor") << "@SUB=ScoutingTrackMonitor::index"
                                          << " could not find " << name;
    return -1;
  }
};

// constructor
ScoutingTrackMonitor::ScoutingTrackMonitor(const edm::ParameterSet& iConfig)
    : conf_(iConfig),
      tracksToken_{consumes<std::vector<Run3ScoutingTrack>>(iConfig.getParameter<edm::InputTag>("tracks"))},
      verticesToken_{consumes<std::vector<Run3ScoutingVertex>>(iConfig.getParameter<edm::InputTag>("vertices"))},
      beamSpotToken_{consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpotLabel"))},
      topFolderName_{iConfig.getParameter<std::string>("topFolderName")},
      h_bsX(nullptr),
      h_bsY(nullptr),
      h_bsZ(nullptr),
      h_bsSigmaZ(nullptr),
      h_bsDxdz(nullptr),
      h_bsDydz(nullptr),
      h_bsBeamWidthX(nullptr),
      h_bsBeamWidthY(nullptr),
      h_bsType(nullptr) {
  hitProfiles_ = {{"nValidPixelHits", "nValidPixelHits", 0., 10.},
                  {"nTrackerLayersWithMeasurement", "nTrackerLayersWithMeasurement", 0., 20.},
                  {"nValidStripHits", "nValidStripHits", 0., 30.}};
  impactParameterProfiles_ = {{"dxy", "d_{xy}", -0.15 * cmToUm, 0.15 * cmToUm},
                              {"dz", "d_{z}", -0.35 * cmToUm, 0.35 * cmToUm}};
}

// histogram booking
void ScoutingTrackMonitor::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const&, edm::EventSetup const&) {
  ibooker.setCurrentFolder(topFolderName_);

  h_dxy = ibooker.book1DD("dxy", "d_{xy};d_{xy} [#mum];Tracks", 100, -0.15 * cmToUm, 0.15 * cmToUm);
  h_dz = ibooker.book1DD("dz", "d_{z};d_{z} [#mum];Tracks", 100, -0.35 * cmToUm, 0.35 * cmToUm);

  h_vtx_idx = ibooker.book1DD("vertexIndex", "tracks Vertex Index;Vertex index;Tracks", 17, -1.5, 15.5);

  h_vtx_sumPt2 =
      ibooker.book1DD("vtxSumPt2", "Vertex #Sigma p_{T}^{2};#Sigma p_{T}^{2} [GeV^{2}];Vertices", 100, 0., 10000.);

  h_vtx_nTracks = ibooker.book1DD("vtxNTracks", "Tracks per vertex;N_{tracks};Vertices", 50, -0.5, 49.5);

  // 2D eta-phi occupancy histograms
  h2_eta_phi = ibooker.book2I(
      "eta_vs_phi", "Track occupancy;#eta;#phi [rad]", 50, -3.0, 3.0, 50, -std::numbers::pi, std::numbers::pi);
  h2_eta_phi->setOption("colz");

  h_vtx_chi2ndf = ibooker.book1DD("vtxChi2ndf", "PV #chi^{2}/ndof", 100, 0., 20.);
  h_vtx_prob = ibooker.book1DD("vtxChi2prob", "PV #chi^{2} probability", 100, 0., 1.);

  h_vr = ibooker.book1DD("vr", "radius at DCA;r_{DCA} (cm);Tracks", 100, 0., 1.);
  h_vz = ibooker.book1DD("vz", "z at DCA;z_{DCA} (cm);Tracks", 100, -30., 30.);

  // Profiles
  constexpr int nEtaBins = 50;
  constexpr double etaMin = -3.0;
  constexpr double etaMax = 3.0;

  constexpr int nPhiBins = 50;
  constexpr double phiMin = -std::numbers::pi;
  constexpr double phiMax = std::numbers::pi;

  // n. hits profiles
  for (auto& cfg : hitProfiles_) {
    const std::string& base = cfg.name;
    const std::string& title = cfg.title;

    // 2D Profile: vs eta-phi
    cfg.p2_eta_phi = ibooker.bookProfile2D(base + "_vs_eta_phi_prof",
                                           title + " vs #eta-#phi;#eta;#phi [rad];#LT" + title + "#GT",
                                           nEtaBins,
                                           etaMin,
                                           etaMax,
                                           nPhiBins,
                                           phiMin,
                                           phiMax,
                                           cfg.ymin,
                                           cfg.ymax,
                                           "");
    cfg.p2_eta_phi->setOption("colz");

    // 1D Profile: vs eta
    cfg.p_eta = ibooker.bookProfile(base + "_vs_eta_prof",
                                    title + " vs #eta;#eta;#LT" + title + "#GT",
                                    nEtaBins,
                                    etaMin,
                                    etaMax,
                                    cfg.ymin,
                                    cfg.ymax,
                                    "");

    // 1D Profile: vs phi
    cfg.p_phi = ibooker.bookProfile(base + "_vs_phi_prof",
                                    title + " vs #phi;#phi [rad];#LT" + title + "#GT",
                                    nPhiBins,
                                    phiMin,
                                    phiMax,
                                    cfg.ymin,
                                    cfg.ymax,
                                    "");
  }

  // IP profiles
  for (auto& cfg : impactParameterProfiles_) {
    // Profile vs eta
    cfg.p_eta = ibooker.bookProfile(cfg.name + "_vs_eta",
                                    cfg.title + " vs #eta;#eta;#LT" + cfg.title + "#GT [#mum]",
                                    nEtaBins,
                                    etaMin,
                                    etaMax,
                                    cfg.ymin,
                                    cfg.ymax,
                                    "");

    // Profile vs phi
    cfg.p_phi = ibooker.bookProfile(cfg.name + "_vs_phi",
                                    cfg.title + " vs #phi;#phi [rad];#LT" + cfg.title + "#GT [#mum]",
                                    nPhiBins,
                                    phiMin,
                                    phiMax,
                                    cfg.ymin,
                                    cfg.ymax,
                                    "");

    // 2D profile vs eta-phi
    cfg.p2_eta_phi = ibooker.bookProfile2D(cfg.name + "_vs_eta_phi",
                                           cfg.title + " vs #eta-#phi;#eta;#phi [rad];#LT" + cfg.title + "#GT [#mum]",
                                           nEtaBins,
                                           etaMin,
                                           etaMax,
                                           nPhiBins,
                                           phiMin,
                                           phiMax,
                                           cfg.ymin,
                                           cfg.ymax,
                                           "");
    cfg.p2_eta_phi->setOption("colz");
  }

  // intialize the profiles
  double xBins[19] = {0., 0.15, 0.5, 1., 1.5, 2., 2.5, 3., 3.5, 4., 4.5, 5., 7., 10., 15., 25., 40., 100., 200.};
  vTrackProfiles_.push_back(ibooker.bookProfile("p_d0_vs_phi",
                                                "Transverse Impact Parameter vs. #phi;#phi_{Track};#LT d_{0} #GT [cm]",
                                                100,
                                                -std::numbers::pi,
                                                std::numbers::pi,
                                                -0.15 * cmToUm,
                                                0.15 * cmToUm,
                                                ""));
  vTrackProfiles_.push_back(
      ibooker.bookProfile("p_dz_vs_phi",
                          "Longitudinal Impact Parameter vs. #phi;#phi_{Track};#LT d_{z} #GT [cm]",
                          100,
                          -std::numbers::pi,
                          std::numbers::pi,
                          -0.35 * cmToUm,
                          0.35 * cmToUm,
                          ""));
  vTrackProfiles_.push_back(ibooker.bookProfile("p_d0_vs_eta",
                                                "Transverse Impact Parameter vs. #eta;#eta_{Track};#LT d_{0} #GT [cm]",
                                                100,
                                                -3.,
                                                3.,
                                                -0.15 * cmToUm,
                                                0.15 * cmToUm,
                                                ""));
  vTrackProfiles_.push_back(
      ibooker.bookProfile("p_dz_vs_eta",
                          "Longitudinal Impact Parameter vs. #eta;#eta_{Track};#LT d_{z} #GT [cm]",
                          100,
                          -3.,
                          3.,
                          -0.35 * cmToUm,
                          0.35 * cmToUm,
                          ""));
  vTrackProfiles_.push_back(ibooker.bookProfile("p_chi2_vs_phi",
                                                "#chi^{2} vs. #phi;#phi_{Track};#LT #chi^{2} #GT",
                                                100,
                                                -std::numbers::pi,
                                                std::numbers::pi,
                                                0,
                                                100,
                                                ""));
  vTrackProfiles_.push_back(
      ibooker.bookProfile("p_chi2Prob_vs_phi",
                          "#chi^{2} probablility vs. #phi;#phi_{Track};#LT #chi^{2} probability#GT",
                          100,
                          -std::numbers::pi,
                          std::numbers::pi,
                          0.,
                          1.,
                          ""));
  vTrackProfiles_.push_back(
      ibooker.bookProfile("p_chi2Prob_vs_d0",
                          "#chi^{2} probablility vs. |d_{0}|;|d_{0}|[cm];#LT #chi^{2} probability#GT",
                          100,
                          0.,
                          80.,
                          0.,
                          1.,
                          ""));
  vTrackProfiles_.push_back(ibooker.bookProfile("p_chi2Prob_vs_dz",
                                                "#chi^{2} probablility vs. dz;d_{z} [cm];#LT #chi^{2} probability#GT",
                                                100,
                                                -30,
                                                30,
                                                0.,
                                                1.,
                                                ""));
  vTrackProfiles_.push_back(ibooker.bookProfile("p_normchi2_vs_phi",
                                                "#chi^{2}/ndof vs. #phi;#phi_{Track};#LT #chi^{2}/ndof #GT",
                                                100,
                                                -std::numbers::pi,
                                                std::numbers::pi,
                                                0.,
                                                5.,
                                                ""));
  vTrackProfiles_.push_back(ibooker.bookProfile(
      "p_chi2_vs_eta", "#chi^{2} vs. #eta;#eta_{Track};#LT #chi^{2} #GT", 100, -3., 3., 0., 100., ""));
  vTrackProfiles_.push_back(ibooker.bookProfile("p_normchi2_vs_pt",
                                                "norm #chi^{2} vs. p_{T}_{Track}; p_{T}_{Track};#LT #chi^{2}/ndof #GT",
                                                18,
                                                xBins,
                                                0.,
                                                5.,
                                                ""));
  vTrackProfiles_.push_back(ibooker.bookProfile(
      "p_normchi2_vs_p", "#chi^{2}/ndof vs. p_{Track};p_{Track};#LT #chi^{2}/ndof #GT", 18, xBins, 0., 5., ""));
  vTrackProfiles_.push_back(
      ibooker.bookProfile("p_chi2Prob_vs_eta",
                          "#chi^{2} probability vs. #eta;#eta_{Track};#LT #chi^{2} probability #GT",
                          100,
                          -3.,
                          3.,
                          0.,
                          1.,
                          ""));
  vTrackProfiles_.push_back(ibooker.bookProfile(
      "p_normchi2_vs_eta", "#chi^{2}/ndof vs. #eta;#eta_{Track};#LT #chi^{2}/ndof #GT", 100, -3., 3., 0., 5, ""));
  vTrackProfiles_.push_back(ibooker.bookProfile(
      "p_kappa_vs_phi", "#kappa vs. #phi;#phi_{Track};#kappa", 100, -std::numbers::pi, std::numbers::pi, -5., 5., ""));
  vTrackProfiles_.push_back(
      ibooker.bookProfile("p_kappa_vs_eta", "#kappa vs. #eta;#eta_{Track};#kappa", 100, -3., 3., -5., 5., ""));
  vTrackProfiles_.push_back(
      ibooker.bookProfile("p_ptResolution_vs_phi",
                          "#delta_{p_{T}}/p_{T}^{track};#phi^{track};#delta_{p_{T}}/p_{T}^{track}",
                          100,
                          -std::numbers::pi,
                          std::numbers::pi,
                          0.,
                          1.,
                          ""));
  vTrackProfiles_.push_back(
      ibooker.bookProfile("p_ptResolution_vs_eta",
                          "#delta_{p_{T}}/p_{T}^{track};#eta^{track};#delta_{p_{T}}/p_{T}^{track}",
                          100,
                          -3.,
                          3.,
                          0.,
                          1.,
                          ""));
  vTrackProfiles_.push_back(
      ibooker.bookProfile("p_ptResolution_vs_pt",
                          "#delta_{p_{T}}/p_{T}^{track};p_{T}^{track};#delta_{p_{T}}/p_{T}^{track}",
                          100,
                          0.,
                          100.,
                          0.,
                          1.,
                          ""));

  // initialize and book the monitors;
  dxy_pt1.varname_ = "xy";
  dxy_pt1.pTcut_ = 1.f;
  dxy_pt1.bookIPMonitor(ibooker, conf_);

  dxy_pt10.varname_ = "xy";
  dxy_pt10.pTcut_ = 10.f;
  dxy_pt10.bookIPMonitor(ibooker, conf_);

  dz_pt1.varname_ = "z";
  dz_pt1.pTcut_ = 1.f;
  dz_pt1.bookIPMonitor(ibooker, conf_);

  dz_pt10.varname_ = "z";
  dz_pt10.pTcut_ = 10.f;
  dz_pt10.bookIPMonitor(ibooker, conf_);

  // BeamSpot
  auto vposx = conf_.getParameter<double>("Xpos");
  auto vposy = conf_.getParameter<double>("Ypos");

  h_bsX = ibooker.book1D("bsX", "BeamSpot x0", 100, vposx - 0.1, vposx + 0.1);
  h_bsY = ibooker.book1D("bsY", "BeamSpot y0", 100, vposy - 0.1, vposy + 0.1);
  h_bsZ = ibooker.book1D("bsZ", "BeamSpot z0", 100, -2., 2.);
  h_bsSigmaZ = ibooker.book1D("bsSigmaZ", "BeamSpot sigmaZ", 100, 0., 10.);
  h_bsDxdz = ibooker.book1D("bsDxdz", "BeamSpot dxdz", 100, -0.0003, 0.0003);
  h_bsDydz = ibooker.book1D("bsDydz", "BeamSpot dydz", 100, -0.0003, 0.0003);
  h_bsBeamWidthX = ibooker.book1D("bsBeamWidthX", "BeamSpot BeamWidthX", 500, 0., 15.);
  h_bsBeamWidthY = ibooker.book1D("bsBeamWidthY", "BeamSpot BeamWidthY", 500, 0., 15.);
  h_bsType = ibooker.book1D("bsType", "BeamSpot type", 4, -1.5, 2.5);
  h_bsType->setBinLabel(1, "Unknown");
  h_bsType->setBinLabel(2, "Fake");
  h_bsType->setBinLabel(3, "LHC");
  h_bsType->setBinLabel(4, "Tracker");
}

void ScoutingTrackMonitor::IPMonitoring::bookIPMonitor(DQMStore::IBooker& iBooker, const edm::ParameterSet& config) {
  int VarBin = config.getParameter<int>(fmt::format("D{}Bin", varname_));
  double VarMin = config.getParameter<double>(fmt::format("D{}Min", varname_));
  double VarMax = config.getParameter<double>(fmt::format("D{}Max", varname_));

  PhiBin_ = config.getParameter<int>("PhiBin");
  PhiMin_ = config.getParameter<double>("PhiMin");
  PhiMax_ = config.getParameter<double>("PhiMax");
  int PhiBin2D = config.getParameter<int>("PhiBin2D");

  EtaBin_ = config.getParameter<int>("EtaBin");
  EtaMin_ = config.getParameter<double>("EtaMin");
  EtaMax_ = config.getParameter<double>("EtaMax");
  int EtaBin2D = config.getParameter<int>("EtaBin2D");

  PtBin_ = config.getParameter<int>("PtBin");
  PtMin_ = config.getParameter<double>("PtMin") * pTcut_;
  PtMax_ = config.getParameter<double>("PtMax") * pTcut_;

  // 1D variables
  IP_ = iBooker.book1DD(fmt::format("d{}_pt{}", varname_, pTcut_),
                        fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} (#mum)", pTcut_, varname_),
                        VarBin,
                        VarMin,
                        VarMax);

  IPErr_ = iBooker.book1DD(fmt::format("d{}Err_pt{}", varname_, pTcut_),
                           fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} error (#mum)", pTcut_, varname_),
                           100,
                           0.,
                           (varname_.find("xy") != std::string::npos) ? 2000. : 10000.);

  IPPull_ = iBooker.book1DD(
      fmt::format("d{}Pull_pt{}", varname_, pTcut_),
      fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}}/#sigma_{{d_{{{}}}}}", pTcut_, varname_, varname_),
      100,
      -5.,
      5.);

  // IP profiles
  IPVsPhi_ = iBooker.bookProfile(fmt::format("d{}VsPhi_pt{}", varname_, pTcut_),
                                 fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} VS track #phi", pTcut_, varname_),
                                 PhiBin_,
                                 PhiMin_,
                                 PhiMax_,
                                 VarBin,
                                 VarMin,
                                 VarMax,
                                 "");
  IPVsPhi_->setAxisTitle("PV track (p_{T} > 1 GeV) #phi", 1);
  IPVsPhi_->setAxisTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} (#mum)", pTcut_, varname_), 2);

  IPVsEta_ = iBooker.bookProfile(fmt::format("d{}VsEta_pt{}", varname_, pTcut_),
                                 fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} VS track #eta", pTcut_, varname_),
                                 EtaBin_,
                                 EtaMin_,
                                 EtaMax_,
                                 VarBin,
                                 VarMin,
                                 VarMax,
                                 "");
  IPVsEta_->setAxisTitle("PV track (p_{T} > 1 GeV) #eta", 1);
  IPVsEta_->setAxisTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} (#mum)", pTcut_, varname_), 2);

  IPVsPt_ = sctTrackMonitor::makeProfileIfLog(
      iBooker,
      true,  /* x-axis */
      false, /* y-axis */
      fmt::format("d{}VsPt_pt{}", varname_, pTcut_).c_str(),
      fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} VS track p_{{T}}", pTcut_, varname_).c_str(),
      PtBin_,
      log10(PtMin_),
      log10(PtMax_),
      VarMin,
      VarMax,
      "");
  IPVsPt_->setAxisTitle("PV track (p_{T} > 1 GeV) p_{T} [GeV]", 1);
  IPVsPt_->setAxisTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} (#mum)", pTcut_, varname_), 2);

  // IP error profiles
  IPErrVsPhi_ =
      iBooker.bookProfile(fmt::format("d{}ErrVsPhi_pt{}", varname_, pTcut_),
                          fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} error VS track #phi", pTcut_, varname_),
                          PhiBin_,
                          PhiMin_,
                          PhiMax_,
                          VarBin,
                          0.,
                          (varname_.find("xy") != std::string::npos) ? 100. : 200.,
                          "");
  IPErrVsPhi_->setAxisTitle("PV track (p_{T} > 1 GeV) #phi", 1);
  IPErrVsPhi_->setAxisTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} error (#mum)", pTcut_, varname_), 2);

  IPErrVsEta_ =
      iBooker.bookProfile(fmt::format("d{}ErrVsEta_pt{}", varname_, pTcut_),
                          fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} error VS track #eta", pTcut_, varname_),
                          EtaBin_,
                          EtaMin_,
                          EtaMax_,
                          VarBin,
                          0.,
                          (varname_.find("xy") != std::string::npos) ? 100. : 200.,
                          "");
  IPErrVsEta_->setAxisTitle("PV track (p_{T} > 1 GeV) #eta", 1);
  IPErrVsEta_->setAxisTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} error (#mum)", pTcut_, varname_), 2);

  IPErrVsPt_ = sctTrackMonitor::makeProfileIfLog(
      iBooker,
      true,  /* x-axis */
      false, /* y-axis */
      fmt::format("d{}ErrVsPt_pt{}", varname_, pTcut_).c_str(),
      fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} error VS track p_{{T}}", pTcut_, varname_).c_str(),
      PtBin_,
      log10(PtMin_),
      log10(PtMax_),
      VarMin,
      VarMax,
      "");
  IPErrVsPt_->setAxisTitle("PV track (p_{T} > 1 GeV) p_{T} [GeV]", 1);
  IPErrVsPt_->setAxisTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} error (#mum)", pTcut_, varname_), 2);

  // 2D profiles
  IPVsEtaVsPhi_ = iBooker.bookProfile2D(
      fmt::format("d{}VsEtaVsPhi_pt{}", varname_, pTcut_),
      fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} VS track #eta VS track #phi", pTcut_, varname_),
      EtaBin2D,
      EtaMin_,
      EtaMax_,
      PhiBin2D,
      PhiMin_,
      PhiMax_,
      VarBin,
      VarMin,
      VarMax,
      "");
  IPVsEtaVsPhi_->setAxisTitle("PV track (p_{T} > 1 GeV) #eta", 1);
  IPVsEtaVsPhi_->setAxisTitle("PV track (p_{T} > 1 GeV) #phi", 2);
  IPVsEtaVsPhi_->setAxisTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} (#mum)", pTcut_, varname_), 3);
  IPVsEtaVsPhi_->setOption("colz");

  IPErrVsEtaVsPhi_ = iBooker.bookProfile2D(
      fmt::format("d{}ErrVsEtaVsPhi_pt{}", varname_, pTcut_),
      fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} error VS track #eta VS track #phi", pTcut_, varname_),
      EtaBin2D,
      EtaMin_,
      EtaMax_,
      PhiBin2D,
      PhiMin_,
      PhiMax_,
      VarBin,
      0.,
      (varname_.find("xy") != std::string::npos) ? 100. : 200.,
      "");
  IPErrVsEtaVsPhi_->setAxisTitle("PV track (p_{T} > 1 GeV) #eta", 1);
  IPErrVsEtaVsPhi_->setAxisTitle("PV track (p_{T} > 1 GeV) #phi", 2);
  IPErrVsEtaVsPhi_->setAxisTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} error (#mum)", pTcut_, varname_),
                                 3);
  IPErrVsEtaVsPhi_->setOption("colz");
}

template <typename T>
bool ScoutingTrackMonitor::getValidHandle(const edm::Event& iEvent,
                                          const edm::EDGetTokenT<T>& token,
                                          edm::Handle<T>& handle,
                                          const std::string& label) {
  iEvent.getByToken(token, handle);
  if (!handle.isValid()) {
    edm::LogWarning("ScoutingTrackMonitor") << "Invalid handle for " << label;
    return false;
  }
  return true;
}

// main event loop
void ScoutingTrackMonitor::analyze(const edm::Event& iEvent, const edm::EventSetup&) {
  edm::Handle<std::vector<Run3ScoutingVertex>> primaryVerticesH;
  edm::Handle<std::vector<Run3ScoutingTrack>> tracksH;

  if (!getValidHandle(iEvent, verticesToken_, primaryVerticesH, "primary vertices") ||
      !getValidHandle(iEvent, tracksToken_, tracksH, "tracks")) {
    return;
  }

  // derefernce handles when it's safe to do so.
  auto const& tracks = *tracksH;
  auto const& vertices = *primaryVerticesH;

  if (vertices.empty())
    return;

  for (const auto& vtx : vertices) {
    h_vtx_chi2ndf->Fill(vtx.chi2() / vtx.ndof());
    h_vtx_prob->Fill(TMath::Prob(vtx.chi2(), vtx.ndof()));
    const int scoutingTracksSize = static_cast<int>(vtx.tracksSize());
    h_vtx_nTracks->Fill(scoutingTracksSize);
  }

  // --- Per-vertex accumulators (indexed by vertex index) ---
  const unsigned int nVtx = vertices.size();
  std::vector<double> vtxSumPt2(nVtx, 0.);
  std::vector<unsigned int> vtxNTracks(nVtx, 0);

  for (const auto& trk : tracks) {
    // --- build reco track ---
    reco::Track recoTrk = makeRecoTrack(trk);

    //auto [vtxIndex, closestVtx] = findClosestScoutingVertex(&recoTrk, vertices);
    //if (!closestVtx)
    //  continue;

    // initialize the impact parameters to large values
    // std::pair<float, float> best_offset{9999.f, 99999.f};

    // // loop on all the vertices and find the closest one
    // unsigned int vtxIndex = 9999;
    // unsigned int idx = 0;
    // for (const auto& vtx : vertices) {
    //   const auto offset = trk_vtx_offSet(trk, vtx);
    //   if (std::abs(offset.second) < std::abs(best_offset.second)) {
    //     best_offset = offset;
    //     vtxIndex = idx;  // save the index of the best vertex
    //   }
    //   idx++;
    // }

    const unsigned int vtxIndex = static_cast<unsigned int>(trk.tk_vtxInd());
    if (vtxIndex >= vertices.size())
      continue;

    const Run3ScoutingVertex* closestVtx = &vertices[vtxIndex];
    h_vtx_idx->Fill(vtxIndex);

    // accumulate per-vertex quantities before filling per-track histos
    if (vtxIndex < nVtx) {
      const float pt = trk.tk_pt();
      vtxSumPt2[vtxIndex] += pt * pt;
      vtxNTracks[vtxIndex]++;
    }

    const float eta = trk.tk_eta();
    const float phi = trk.tk_phi();
    const float pt = trk.tk_pt();

    const int nValidPixelHits = trk.tk_nValidPixelHits();
    const int nTrackerLayers = trk.tk_nTrackerLayersWithMeasurement();
    const int nValidStripHits = trk.tk_nValidStripHits();

    // --- fill 2D eta-phi occupancy histograms ---
    h2_eta_phi->Fill(eta, phi);
    const std::array<double, 3> values = {static_cast<double>(nValidPixelHits),
                                          static_cast<double>(nTrackerLayers),
                                          static_cast<double>(nValidStripHits)};

    for (size_t i = 0; i < hitProfiles_.size(); ++i) {
      const auto& cfg = hitProfiles_[i];
      const double value = values[i];
      cfg.p2_eta_phi->Fill(eta, phi, value);
      cfg.p_eta->Fill(eta, value);
      cfg.p_phi->Fill(phi, value);
    }

    // --- build reco vertex ---
    reco::Vertex recoVtx = makeRecoVertex(*closestVtx);

    // --- impact parameters (standard CMSSW definitions) ---
    float dxy = recoTrk.dxy(recoVtx.position()) * cmToUm;
    float dz = recoTrk.dz(recoVtx.position()) * cmToUm;

    float dzErr = std::sqrt(recoTrk.dzError() * recoTrk.dzError() + recoVtx.zError() * recoVtx.zError()) * cmToUm;
    float dxyErr = std::sqrt(recoTrk.dxyError() * recoTrk.dxyError() + recoVtx.xError() * recoVtx.xError() +
                             recoVtx.yError() * recoVtx.yError()) *
                   cmToUm;

    // production radius from stored reference point
    float vr = std::sqrt(trk.tk_vx() * trk.tk_vx() + trk.tk_vy() * trk.tk_vy());
    h_vr->Fill(vr);
    h_vz->Fill(trk.tk_vz());

    // ---- impact parameters using offsets
    //float dxy = best_offset.first;
    //float dz = best_offset.second;

    // --- impact parameters (directly from scouting) ---
    //float dxy = trk.tk_dxy() * cmToUm;
    //float dz = trk.tk_dz() * cmToUm;
    //float dxyErr = trk.tk_dxy_Error() * cmToUm;
    //float dzErr = trk.tk_dz_Error() * cmToUm;

    //float dxyErr = recoTrk.dxyError() * cmToUm;
    //float dzErr = recoTrk.dzError() * cmToUm;

    // --- fill histograms ---
    h_dxy->Fill(dxy);
    h_dz->Fill(dz);

    const std::array<double, 2> ipvalues = {dxy, dz};
    for (size_t i = 0; i < impactParameterProfiles_.size(); ++i) {
      const auto& cfg = impactParameterProfiles_[i];
      const double value = ipvalues[i];

      cfg.p_eta->Fill(eta, value);
      cfg.p_phi->Fill(phi, value);
      cfg.p2_eta_phi->Fill(eta, phi, value);
    }

    // Fill track profiles
    double chi2Prob = TMath::Prob(recoTrk.chi2(), recoTrk.ndof());
    double normchi2 = recoTrk.normalizedChi2();
    double kappa = trk.tk_qoverp();

    //GlobalPoint gPoint(recoTrk.vx(), recoTrk.vy(), recoTrk.vz());
    //double theLocalMagFieldInInverseGeV = magneticField_->inInverseGeV(gPoint).z();
    //double kappa = -recoTrk.charge() * theLocalMagFieldInInverseGeV / recoTrk.pt();

    static const int d0phiindex = this->index(vTrackProfiles_, "p_d0_vs_phi");
    vTrackProfiles_[d0phiindex]->Fill(recoTrk.phi(), recoTrk.d0());
    static const int dzphiindex = this->index(vTrackProfiles_, "p_dz_vs_phi");
    vTrackProfiles_[dzphiindex]->Fill(recoTrk.phi(), recoTrk.dz());
    static const int d0etaindex = this->index(vTrackProfiles_, "p_d0_vs_eta");
    vTrackProfiles_[d0etaindex]->Fill(recoTrk.eta(), recoTrk.d0());
    static const int dzetaindex = this->index(vTrackProfiles_, "p_dz_vs_eta");
    vTrackProfiles_[dzetaindex]->Fill(recoTrk.eta(), recoTrk.dz());
    static const int chiProbphiindex = this->index(vTrackProfiles_, "p_chi2Prob_vs_phi");
    vTrackProfiles_[chiProbphiindex]->Fill(recoTrk.phi(), chi2Prob);
    static const int chiProbabsd0index = this->index(vTrackProfiles_, "p_chi2Prob_vs_d0");
    vTrackProfiles_[chiProbabsd0index]->Fill(fabs(recoTrk.d0()), chi2Prob);
    static const int chiProbabsdzindex = this->index(vTrackProfiles_, "p_chi2Prob_vs_dz");
    vTrackProfiles_[chiProbabsdzindex]->Fill(recoTrk.dz(), chi2Prob);
    static const int chiphiindex = this->index(vTrackProfiles_, "p_chi2_vs_phi");
    vTrackProfiles_[chiphiindex]->Fill(recoTrk.phi(), recoTrk.chi2());
    static const int normchiphiindex = this->index(vTrackProfiles_, "p_normchi2_vs_phi");
    vTrackProfiles_[normchiphiindex]->Fill(recoTrk.phi(), normchi2);
    static const int chietaindex = this->index(vTrackProfiles_, "p_chi2_vs_eta");
    vTrackProfiles_[chietaindex]->Fill(recoTrk.eta(), recoTrk.chi2());
    static const int normchiptindex = this->index(vTrackProfiles_, "p_normchi2_vs_pt");
    vTrackProfiles_[normchiptindex]->Fill(recoTrk.pt(), normchi2);
    static const int normchipindex = this->index(vTrackProfiles_, "p_normchi2_vs_p");
    vTrackProfiles_[normchipindex]->Fill(recoTrk.p(), normchi2);
    static const int chiProbetaindex = this->index(vTrackProfiles_, "p_chi2Prob_vs_eta");
    vTrackProfiles_[chiProbetaindex]->Fill(recoTrk.eta(), chi2Prob);
    static const int normchietaindex = this->index(vTrackProfiles_, "p_normchi2_vs_eta");
    vTrackProfiles_[normchietaindex]->Fill(recoTrk.eta(), normchi2);
    static const int kappaphiindex = this->index(vTrackProfiles_, "p_kappa_vs_phi");
    vTrackProfiles_[kappaphiindex]->Fill(recoTrk.phi(), kappa);
    static const int kappaetaindex = this->index(vTrackProfiles_, "p_kappa_vs_eta");
    vTrackProfiles_[kappaetaindex]->Fill(recoTrk.eta(), kappa);
    static const int ptResphiindex = this->index(vTrackProfiles_, "p_ptResolution_vs_phi");
    vTrackProfiles_[ptResphiindex]->Fill(recoTrk.phi(), recoTrk.ptError() / recoTrk.pt());
    static const int ptResetaindex = this->index(vTrackProfiles_, "p_ptResolution_vs_eta");
    vTrackProfiles_[ptResetaindex]->Fill(recoTrk.eta(), recoTrk.ptError() / recoTrk.pt());
    static const int ptResptindex = this->index(vTrackProfiles_, "p_ptResolution_vs_pt");
    vTrackProfiles_[ptResptindex]->Fill(recoTrk.pt(), recoTrk.ptError() / recoTrk.pt());

    if (trk.tk_pt() < 1.)
      continue;

    // dxy pT>1
    dxy_pt1.IP_->Fill(dxy);
    dxy_pt1.IPVsPhi_->Fill(phi, dxy);
    dxy_pt1.IPVsEta_->Fill(eta, dxy);
    dxy_pt1.IPVsPt_->Fill(pt, dxy);
    dxy_pt1.IPVsEtaVsPhi_->Fill(eta, phi, dxy);

    dxy_pt1.IPErr_->Fill(dxyErr);
    dxy_pt1.IPPull_->Fill(dxy / dxyErr);
    dxy_pt1.IPErrVsPhi_->Fill(phi, dxyErr);
    dxy_pt1.IPErrVsEta_->Fill(eta, dxyErr);
    dxy_pt1.IPErrVsPt_->Fill(pt, dxyErr);
    dxy_pt1.IPErrVsEtaVsPhi_->Fill(eta, phi, dxyErr);

    // dz pT>1
    dz_pt1.IP_->Fill(dz);
    dz_pt1.IPVsPhi_->Fill(phi, dz);
    dz_pt1.IPVsEta_->Fill(eta, dz);
    dz_pt1.IPVsPt_->Fill(pt, dz);
    dz_pt1.IPVsEtaVsPhi_->Fill(eta, phi, dz);

    dz_pt1.IPErr_->Fill(dzErr);
    dz_pt1.IPPull_->Fill(dz / dzErr);
    dz_pt1.IPErrVsPhi_->Fill(phi, dzErr);
    dz_pt1.IPErrVsEta_->Fill(eta, dzErr);
    dz_pt1.IPErrVsPt_->Fill(pt, dzErr);
    dz_pt1.IPErrVsEtaVsPhi_->Fill(eta, phi, dzErr);

    if (pt < 10.)
      continue;

    // dxy pT>10
    dxy_pt10.IP_->Fill(dxy);
    dxy_pt10.IPVsPhi_->Fill(phi, dxy);
    dxy_pt10.IPVsEta_->Fill(eta, dxy);
    dxy_pt10.IPVsPt_->Fill(pt, dxy);
    dxy_pt10.IPVsEtaVsPhi_->Fill(eta, phi, dxy);

    dxy_pt10.IPErr_->Fill(dxyErr);
    dxy_pt10.IPPull_->Fill(dxy / dxyErr);
    dxy_pt10.IPErrVsPhi_->Fill(phi, dxyErr);
    dxy_pt10.IPErrVsEta_->Fill(eta, dxyErr);
    dxy_pt10.IPErrVsPt_->Fill(pt, dxyErr);
    dxy_pt10.IPErrVsEtaVsPhi_->Fill(eta, phi, dxyErr);

    // dz pT>10
    dz_pt10.IP_->Fill(dz);
    dz_pt10.IPVsPhi_->Fill(phi, dz);
    dz_pt10.IPVsEta_->Fill(eta, dz);
    dz_pt10.IPVsPt_->Fill(pt, dz);
    dz_pt10.IPVsEtaVsPhi_->Fill(eta, phi, dz);

    dz_pt10.IPErr_->Fill(dzErr);
    dz_pt10.IPPull_->Fill(dz / dzErr);
    dz_pt10.IPErrVsPhi_->Fill(phi, dzErr);
    dz_pt10.IPErrVsEta_->Fill(eta, dzErr);
    dz_pt10.IPErrVsPt_->Fill(pt, dzErr);
    dz_pt10.IPErrVsEtaVsPhi_->Fill(eta, phi, dzErr);
  }

  // --- Per-vertex histograms ---
  for (unsigned int i = 0; i < nVtx; ++i) {
    if (vtxNTracks[i] == 0)
      continue;
    h_vtx_sumPt2->Fill(vtxSumPt2[i]);
  }

  // BeamSpot
  edm::Handle<reco::BeamSpot> beamSpotH;
  std::unique_ptr<Run3ScoutingVertex> beamspotVertex{nullptr};
  if (!getValidHandle(iEvent, beamSpotToken_, beamSpotH, "beamSpot")) {
    return;
  }

  const auto& beamSpot = *beamSpotH;
  beamspotVertex = std::make_unique<Run3ScoutingVertex>(
      beamSpot.x0(), beamSpot.y0(), beamSpot.z0(), 0., 0., 0., 0., 0., true, 0., 0., 0., 0);

  h_bsX->Fill(beamSpot.x0());
  h_bsY->Fill(beamSpot.y0());
  h_bsZ->Fill(beamSpot.z0());
  h_bsSigmaZ->Fill(beamSpot.sigmaZ());
  h_bsDxdz->Fill(beamSpot.dxdz());
  h_bsDydz->Fill(beamSpot.dydz());
  h_bsBeamWidthX->Fill(beamSpot.BeamWidthX() * cmToUm);
  h_bsBeamWidthY->Fill(beamSpot.BeamWidthY() * cmToUm);
  h_bsType->Fill(beamSpot.type());
}

// helper: build reco::Track
reco::Track ScoutingTrackMonitor::makeRecoTrack(const Run3ScoutingTrack& sTrack) const {
  reco::Track::Point v(sTrack.tk_vx(), sTrack.tk_vy(), sTrack.tk_vz());
  reco::Track::Vector p(math::RhoEtaPhiVector(sTrack.tk_pt(), sTrack.tk_eta(), sTrack.tk_phi()));

  reco::TrackBase::CovarianceMatrix cov;
  cov(0, 0) = std::pow(sTrack.tk_qoverp_Error(), 2);
  cov(0, 1) = sTrack.tk_qoverp_lambda_cov();
  cov(0, 2) = sTrack.tk_qoverp_phi_cov();
  cov(0, 3) = sTrack.tk_qoverp_dxy_cov();
  cov(0, 4) = sTrack.tk_qoverp_dsz_cov();

  cov(1, 1) = std::pow(sTrack.tk_lambda_Error(), 2);
  cov(1, 2) = sTrack.tk_lambda_phi_cov();
  cov(1, 3) = sTrack.tk_lambda_dxy_cov();
  cov(1, 4) = sTrack.tk_lambda_dsz_cov();

  cov(2, 2) = std::pow(sTrack.tk_phi_Error(), 2);
  cov(2, 3) = sTrack.tk_phi_dxy_cov();
  cov(2, 4) = sTrack.tk_phi_dsz_cov();

  cov(3, 3) = std::pow(sTrack.tk_dxy_Error(), 2);
  cov(3, 4) = sTrack.tk_dxy_dsz_cov();

  cov(4, 4) = std::pow(sTrack.tk_dsz_Error(), 2);

  return reco::Track(sTrack.tk_chi2(), sTrack.tk_ndof(), v, p, sTrack.tk_charge(), cov);
}

// helper: build reco::Vertex
reco::Vertex ScoutingTrackMonitor::makeRecoVertex(const Run3ScoutingVertex& sVertex) const {
  reco::Vertex::Error err;

  err(0, 0) = std::pow(sVertex.xError(), 2);
  err(1, 1) = std::pow(sVertex.yError(), 2);
  err(2, 2) = std::pow(sVertex.zError(), 2);

  err(0, 1) = sVertex.xyCov();
  err(0, 2) = sVertex.xzCov();
  err(1, 2) = sVertex.yzCov();

  return reco::Vertex(reco::Vertex::Point(sVertex.x(), sVertex.y(), sVertex.z()),
                      err,
                      sVertex.chi2(),
                      sVertex.ndof(),
                      sVertex.tracksSize());
}

std::pair<unsigned int, const Run3ScoutingVertex*> ScoutingTrackMonitor::findClosestScoutingVertex(
    const reco::Track* track, const std::vector<Run3ScoutingVertex>& vertices) {
  double minDistance = std::numeric_limits<double>::max();
  const Run3ScoutingVertex* closestVertex = nullptr;

  unsigned int index{0}, theIndex{999};

  for (const auto& vertex : vertices) {
    math::XYZPoint vertexPosition(vertex.x(), vertex.y(), vertex.z());

    const auto& trackMomentum = track->momentum();
    const auto& vertexToPoint = vertexPosition - track->referencePoint();

    double distance = vertexToPoint.Cross(trackMomentum).R() / trackMomentum.R();

    if (distance < minDistance) {
      minDistance = distance;
      closestVertex = &vertex;
      theIndex = index;
    }
    index++;
  }
  return std::make_pair(theIndex, closestVertex);
}

void ScoutingTrackMonitor::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  edm::ParameterSetDescription desc;
  desc.add<edm::InputTag>("tracks", edm::InputTag("hltScoutingTrackPacker"));
  desc.add<edm::InputTag>("vertices", edm::InputTag("hltScoutingPrimaryVertexPacker", "primaryVtx"));
  desc.add<edm::InputTag>("beamSpotLabel", edm::InputTag("hltOnlineBeamSpot"));
  desc.add<std::string>("topFolderName", "HLT/ScoutingOffline/Tracks");
  desc.add<double>("Xpos", 0.1);
  desc.add<double>("Ypos", -0.2);
  desc.add<int>("DxyBin", 100);
  desc.add<double>("DxyMin", -5000.0);
  desc.add<double>("DxyMax", 5000.0);
  desc.add<int>("DzBin", 100);
  desc.add<double>("DzMin", -2000.0);
  desc.add<double>("DzMax", 2000.0);
  desc.add<int>("PhiBin", 32);
  desc.add<double>("PhiMin", -std::numbers::pi);
  desc.add<double>("PhiMax", std::numbers::pi);
  desc.add<int>("EtaBin", 26);
  desc.add<double>("EtaMin", -3.0);
  desc.add<double>("EtaMax", 3.0);
  desc.add<int>("PtBin", 49);
  desc.add<double>("PtMin", 1.);
  desc.add<double>("PtMax", 50.);
  desc.add<int>("PhiBin2D", 12);
  desc.add<int>("EtaBin2D", 8);
  descriptions.addWithDefaultLabel(desc);
}

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_FWK_MODULE(ScoutingTrackMonitor);