PackedCandidateTrackValidator.cc

Go to the documentation of this file.

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/VecArray.h"
#include "FWCore/Utilities/interface/isFinite.h"

#include "DQMServices/Core/interface/DQMEDAnalyzer.h"
#include "DQMServices/Core/interface/MonitorElement.h"

#include "DataFormats/Common/interface/Association.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/PatCandidates/interface/PackedCandidate.h"
#include "DataFormats/Math/interface/libminifloat.h"
#include "DataFormats/Math/interface/liblogintpack.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "DataFormats/Math/interface/deltaPhi.h"

#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
#include "Geometry/Records/interface/TrackerTopologyRcd.h"

#include "boost/math/special_functions/sign.hpp"

#include <iomanip>

namespace {
  template<typename T> void fillNoFlow(MonitorElement* me, T val){
    auto h = me->getTH1();
    const auto xaxis = h->GetXaxis();
    if(val <= xaxis->GetXmin())
      h->AddBinContent(xaxis->GetFirst());
    else if(val >= xaxis->GetXmax())
      h->AddBinContent(xaxis->GetLast());
    else
      h->Fill(val);
  }

  class HitPatternPrinter {
  public:
    explicit HitPatternPrinter(const reco::Track& trk): track(trk) {}

    void print(std::ostream& os) const {
      const reco::HitPattern &p = track.hitPattern();

      for (int i = 0; i < p.numberOfAllHits(reco::HitPattern::TRACK_HITS); ++i) {
        uint32_t hit = p.getHitPattern(reco::HitPattern::TRACK_HITS, i);

        detLayer(os, p, hit);
        if(p.missingHitFilter(hit)) {
          os << "(miss)";
        }
        else if(p.inactiveHitFilter(hit)) {
          os << "(inact)";
        }
        else if(p.badHitFilter(hit)) {
          os << "(bad)";
        }
        os << " ";
      }

      if(p.numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS) > 0) {
        os << "lost inner ";

        for (int i = 0; i < p.numberOfAllHits(reco::HitPattern::MISSING_INNER_HITS); ++i) {
          uint32_t hit = p.getHitPattern(reco::HitPattern::MISSING_INNER_HITS, i);
          detLayer(os, p, hit);
          if(p.missingHitFilter(hit)) {
            os << "(miss)";
          }
          else if(p.inactiveHitFilter(hit)) {
            os << "(inact)";
          }
        }
      }
      if(p.numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS) > 0) {
        os << "lost outer ";

        for (int i = 0; i < p.numberOfAllHits(reco::HitPattern::MISSING_OUTER_HITS); ++i) {
          uint32_t hit = p.getHitPattern(reco::HitPattern::MISSING_OUTER_HITS, i);
          detLayer(os, p, hit);
          if(p.missingHitFilter(hit)) {
            os << "(miss)";
          }
          else if(p.inactiveHitFilter(hit)) {
            os << "(inact)";
          }
        }
      }
    }

  private:
    static void detLayer(std::ostream& os, const reco::HitPattern& p, uint32_t hit) {
      if(p.pixelBarrelHitFilter(hit)) {
        os << "BPIX";
      }
      else if(p.pixelEndcapHitFilter(hit)) {
        os << "FPIX";
      }
      else if(p.stripTIBHitFilter(hit)) {
        os << "TIB";
      }
      else if(p.stripTIDHitFilter(hit)) {
        os << "TID";
      }
      else if(p.stripTOBHitFilter(hit)) {
        os << "TOB";
      }
      else if(p.stripTECHitFilter(hit)) {
        os << "TEC";
      }
      os << p.getLayer(hit);
    }

    const reco::Track& track;
  };

  std::ostream& operator<<(std::ostream& os, const HitPatternPrinter& hpp) {
    hpp.print(os);
    return os;
  }

  class TrackAlgoPrinter {
  public:
    explicit TrackAlgoPrinter(const reco::Track& trk): track(trk) {}

    void print(std::ostream& os) const {
      edm::VecArray<reco::TrackBase::TrackAlgorithm, reco::TrackBase::algoSize> algos;
      for(int ialgo=0; ialgo < reco::TrackBase::algoSize; ++ialgo) {
        auto algo = static_cast<reco::TrackBase::TrackAlgorithm>(ialgo);
        if(track.isAlgoInMask(algo)) {
          algos.push_back(algo);
        }
      }

      os << "algo " << reco::TrackBase::algoName(track.algo());
      if(track.originalAlgo() != track.algo())
        os << " originalAlgo " << reco::TrackBase::algoName(track.originalAlgo());
      if(algos.size() > 1) {
        os << " algoMask";
        for(auto algo: algos) {
          os << " " << reco::TrackBase::algoName(algo);
        }
      }
    }

  private:
    const reco::Track& track;
  };
  std::ostream& operator<<(std::ostream& os, const TrackAlgoPrinter& tap) {
    tap.print(os);
    return os;
  }

  double diffRelative(double a, double b) {
    return (a-b)/b;
  }

  class LogIntHelper {
  public:
    LogIntHelper(double lmin, double lmax): lmin_(lmin), lmax_(lmax) {}

    class UnderOverflow {
    public:
      UnderOverflow(double largestValue, double smallestValue, std::function<double(double)> modifyUnpack):
        unpackedLargestValue_(modifyUnpack ? modifyUnpack(largestValue) : largestValue),
        unpackedSmallestValue_(modifyUnpack ? modifyUnpack(smallestValue) : smallestValue)
      {}

      bool compatibleWithUnderflow(double value) const {
        return value == unpackedSmallestValue_;
      }
      void printNonOkUnderflow(std::ostream& os) const {
        os << " (not min " << unpackedSmallestValue_ << ")";
      }

      bool compatibleWithOverflow(double value) const {
        return value == unpackedLargestValue_;
      }
      void printNonOkOverflow(std::ostream& os) const {
        os << " (not max " << unpackedLargestValue_ << ")";
      }

    private:
      // narrow to float to compare apples to apples with values from
      // PackedCandidate (even though the final comparison is done in
      // double)
      const float unpackedLargestValue_;
      const float unpackedSmallestValue_;
    };

    static std::string maxName() { return "max"; }
    static std::string minName() { return "min"; }

    UnderOverflow underOverflowHelper(double value, std::function<double(double)> modifyUnpack) const {
      return UnderOverflow(largestValue(), value >= 0 ? smallestPositiveValue() : std::abs(smallestNegativeValue()), modifyUnpack);
    }

    double largestValue() const {
      return logintpack::unpack8log(127, lmin_, lmax_);
    }

    static bool wouldBeDenorm(double value) {
      return false;
    }

    // lessThan means closer to zero
    bool lessThanSmallestValue(double value) const {
      if(value >= 0)
        return value < smallestPositiveValue();
      else
        return value > smallestNegativeValue();
    }

    double smallestPositiveValue() const {
      return logintpack::unpack8log(logintpack::smallestPositive, lmin_, lmax_);
    }

    double smallestNegativeValue() const {
      return logintpack::unpack8log(logintpack::smallestNegative, lmin_, lmax_);
    }

  private:
    const double lmin_;
    const double lmax_;
  };
  class Float16Helper {
  public:
    class UnderOverflow {
    public:
      static void printNonOkUnderflow(std::ostream& os) {
        os << " (not 0)";
      }
      static bool compatibleWithUnderflow(double value) {
        return value == 0.0;
      }
      static void printNonOkOverflow(std::ostream& os) {
        os << " (not inf)";
      }
      static bool compatibleWithOverflow(double value) {
        return edm::isNotFinite(value);
      }
    };

    static std::string maxName() { return "inf"; }
    static std::string minName() { return "0"; }

    static UnderOverflow underOverflowHelper(double value, std::function<double(double)>) {
      return UnderOverflow();
    }

    static double largestValue() {
      return MiniFloatConverter::max32RoundedToMax16();
    }

    static bool wouldBeDenorm(double value) {
      const float valuef = static_cast<float>(value);
      return valuef >= MiniFloatConverter::denorm_min() && valuef < MiniFloatConverter::min();
    }

    static bool lessThanSmallestValue(double value) {
      return std::abs(value) < smallestValue();
    }

    static double smallestValue() {
      return MiniFloatConverter::denorm_min();
    }
  };

  enum class RangeStatus {
    inrange = 0,
    inrange_signflip = 1,
    denormal = 2,
    underflow_OK = 3,
    underflow_notOK = 4,
    overflow_OK = 5,
    overflow_notOK = 6
  };
  bool isInRange(RangeStatus status) {
    return status == RangeStatus::inrange || status == RangeStatus::denormal;
  }

  template <typename T>
  class PackedValueCheckResult {
  public:
    PackedValueCheckResult(RangeStatus status, double diff,
                           double pcvalue, double trackvalue,
                           double rangeMin, double rangeMax,
           const typename T::UnderOverflow& underOverflow):
      diff_(diff), pcvalue_(pcvalue), trackvalue_(trackvalue), rangeMin_(rangeMin), rangeMax_(rangeMax),
      status_(status), underOverflow_(underOverflow)
    {}

    RangeStatus status() const { return status_; }
    double diff() const { return diff_; }

    bool outsideExpectedRange() const {
      if(status_ == RangeStatus::inrange)
        return diff_ < rangeMin_ || diff_ > rangeMax_;
      // denormal is considered as "in range" regardless of the expected range
      return status_ == RangeStatus::underflow_notOK || status_ == RangeStatus::overflow_notOK || status_ == RangeStatus::inrange_signflip;
    }

    void print(std::ostream& os) const {
      if(outsideExpectedRange())
        os << "!! ";
      os << "(" << rangeMin_ << "," << rangeMax_ << ") ";

      os << diff_ << " " << pcvalue_;
      if(status_ == RangeStatus::underflow_OK || status_ == RangeStatus::underflow_notOK) {
        os << " (underflow) ";
        if(status_ == RangeStatus::underflow_notOK)
          underOverflow_.printNonOkUnderflow(os);
      }
      else if(status_ == RangeStatus::overflow_OK || status_ == RangeStatus::overflow_notOK) {
        os << " (overflow) ";
        if(status_ == RangeStatus::overflow_notOK)
          underOverflow_.printNonOkOverflow(os);
      }
      else if(status_ == RangeStatus::denormal)
        os << " (denormal)";
      os << " " << trackvalue_;
    }

  private:
    const double diff_;
    const double pcvalue_;
    const double trackvalue_;
    const double rangeMin_;
    const double rangeMax_;
    const RangeStatus status_;
    const typename T::UnderOverflow underOverflow_;
  };

  struct Range {
    Range(double mi, double ma): min(mi), max(ma) {}
    const double min, max;
  };
  struct RangeAbs {
    RangeAbs(double val): min(-val), max(val) {}
    const double min, max;
  };

  template <typename T>
  class PackedValueCheck {
  public:
    template <typename R, typename ...Args>
    PackedValueCheck(const R& range, Args&&... args):
      helper_(std::forward<Args>(args)...),
      rangeMin_(range.min), rangeMax_(range.max)
    {}

    void book(DQMStore::IBooker& iBooker,
              const std::string& name, const std::string& title,
              int nbins,double min,double max,
              int flow_nbins,double flow_min,double flow_max) {
      hInrange = iBooker.book1D(name, title, nbins, min, max);
      hUnderOverflowSign = iBooker.book1D(name+"UnderOverFlowSign", title+" with over- and underflow, and sign flip", flow_nbins, flow_min, flow_max);
      hStatus = iBooker.book1D(name+"Status", title+" status", 7, -0.5, 6.5);
      hStatus->setBinLabel(1, "In range");
      hStatus->setBinLabel(2, "In range, sign flip");
      hStatus->setBinLabel(3, "Denormal");
      hStatus->setBinLabel(4, "Underflow, PC is "+T::minName());
      hStatus->setBinLabel(5, "Underflow, PC is not "+T::minName());
      hStatus->setBinLabel(6, "Overflow, PC is "+T::maxName());
      hStatus->setBinLabel(7, "Overflow, PC is not "+T::maxName());
    }

    PackedValueCheckResult<T> fill(double pcvalue, double trackvalue,
                std::function<double(double)> modifyPack=std::function<double(double)>(),
                std::function<double(double)> modifyUnpack=std::function<double(double)>()) {
      const auto diff = diffRelative(pcvalue, trackvalue);

      const auto tmpSigned = modifyPack ? modifyPack(trackvalue) : trackvalue;
      const auto tmp = std::abs(tmpSigned);
      const auto underOverflow = helper_.underOverflowHelper(tmpSigned, modifyUnpack);
      RangeStatus status;
      if(tmp > helper_.largestValue()) {
        fillNoFlow(hUnderOverflowSign, diff);
        if(underOverflow.compatibleWithOverflow(std::abs(pcvalue))) {
          status = RangeStatus::overflow_OK;
        }
        else {
          status = RangeStatus::overflow_notOK;
        }
      }
      else if(helper_.lessThanSmallestValue(tmpSigned)) {
        fillNoFlow(hUnderOverflowSign, diff);
        if(underOverflow.compatibleWithUnderflow(std::abs(pcvalue))) {
          status = RangeStatus::underflow_OK;
        }
        else {
          status = RangeStatus::underflow_notOK;
        }
      }
      else {
        if(boost::math::sign(pcvalue) == boost::math::sign(trackvalue)) {
          if(T::wouldBeDenorm(tmp)) {
            status = RangeStatus::denormal;
          }
          else {
            status = RangeStatus::inrange;
          }
          fillNoFlow(hInrange, diff);
        }
        else {
          fillNoFlow(hUnderOverflowSign, diff);
          status = RangeStatus::inrange_signflip;
        }
      }
      hStatus->Fill(static_cast<int>(status));

      return PackedValueCheckResult<T>(status, diff, pcvalue, trackvalue, rangeMin_, rangeMax_, underOverflow);
    }

  private:
    const T helper_;
    const double rangeMin_;
    const double rangeMax_;

    MonitorElement *hInrange;
    MonitorElement *hUnderOverflowSign;
    MonitorElement *hStatus;
  };
  template <typename T>
  std::ostream& operator<<(std::ostream& os, const PackedValueCheckResult<T>& res) {
    res.print(os);
    return os;
  }
}

class PackedCandidateTrackValidator: public DQMEDAnalyzer{
 public:
  PackedCandidateTrackValidator(const edm::ParameterSet& pset);
  ~PackedCandidateTrackValidator() override;

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

  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

 private:

  edm::EDGetTokenT<edm::View<reco::Track>> tracksToken_;
  edm::EDGetTokenT<reco::VertexCollection> verticesToken_;
  edm::EDGetTokenT<edm::Association<pat::PackedCandidateCollection>> trackToPackedCandidateToken_;

  std::string rootFolder_;
  bool debug_;

  enum {
    sf_AllTracks = 0,
    sf_AssociatedToPC = 1,
    sf_PCIsCharged = 2,
    sf_PCHasTrack = 3,
    sf_PCIsNotElectron = 4,
    sf_PCHasHits = 5,
    sf_PCNdofNot0 = 6,
    sf_NoMissingInnerHits = 7
  };
  MonitorElement *h_selectionFlow;

  MonitorElement *h_diffVx;
  MonitorElement *h_diffVy;
  MonitorElement *h_diffVz;

  MonitorElement *h_diffNormalizedChi2;
  MonitorElement *h_diffNdof;

  MonitorElement *h_diffCharge;
  MonitorElement *h_diffIsHighPurity;

  MonitorElement *h_diffPt;
  MonitorElement *h_diffEta;
  MonitorElement *h_diffPhi;
  PackedValueCheck<Float16Helper> h_diffDxyAssocPV;
  PackedValueCheck<Float16Helper> h_diffDzAssocPV;
  MonitorElement *h_diffDxyPV;
  MonitorElement *h_diffDzPV;

  MonitorElement *h_diffTrackDxyAssocPV;
  MonitorElement *h_diffTrackDzAssocPV;

  PackedValueCheck<LogIntHelper> h_diffCovQoverpQoverp;
  PackedValueCheck<LogIntHelper> h_diffCovLambdaLambda;
  PackedValueCheck<LogIntHelper> h_diffCovLambdaDsz;
  PackedValueCheck<LogIntHelper> h_diffCovPhiPhi;
  PackedValueCheck<LogIntHelper> h_diffCovPhiDxy;
  PackedValueCheck<Float16Helper> h_diffCovDxyDxy;
  PackedValueCheck<Float16Helper> h_diffCovDxyDsz;
  PackedValueCheck<Float16Helper> h_diffCovDszDsz;

  MonitorElement *h_diffDxyError;
  MonitorElement *h_diffDszError;
  MonitorElement *h_diffDzError;

  MonitorElement *h_diffTrackDxyError;
  MonitorElement *h_diffTrackDzError;

  MonitorElement *h_diffPtError;
  MonitorElement *h_diffEtaError;


  MonitorElement *h_diffNumberOfPixelLayers;
  MonitorElement *h_diffNumberOfStripLayers;
  MonitorElement *h_diffNumberOfPixelHits;
  MonitorElement *h_diffNumberOfHits;
  MonitorElement *h_diffLostInnerHits;

  MonitorElement *h_diffHitPatternPixelLayersWithMeasurement;
  MonitorElement *h_diffHitPatternTrackerLayersWithMeasurement;
  MonitorElement *h_diffHitPatternStripLayersWithMeasurement;
  MonitorElement *h_diffHitPatternNumberOfValidPixelHits;
  MonitorElement *h_diffHitPatternNumberOfValidHits;
  MonitorElement *h_diffHitPatternNumberOfLostInnerHits;
  MonitorElement *h_diffHitPatternHasValidHitInFirstPixelBarrel;

  MonitorElement *h_numberPixelLayersOverMax;
  MonitorElement *h_numberStripLayersOverMax;
  MonitorElement *h_numberLayersOverMax;
  MonitorElement *h_numberPixelHitsOverMax;
  MonitorElement *h_numberStripHitsOverMax;
  MonitorElement *h_numberHitsOverMax;
};

PackedCandidateTrackValidator::PackedCandidateTrackValidator(const edm::ParameterSet& iConfig):
  tracksToken_(consumes<edm::View<reco::Track>>(iConfig.getUntrackedParameter<edm::InputTag>("tracks"))),
  verticesToken_(consumes<reco::VertexCollection>(iConfig.getUntrackedParameter<edm::InputTag>("vertices"))),
  trackToPackedCandidateToken_(consumes<edm::Association<pat::PackedCandidateCollection>>(iConfig.getUntrackedParameter<edm::InputTag>("trackToPackedCandidateAssociation"))),
  rootFolder_(iConfig.getUntrackedParameter<std::string>("rootFolder")),
  debug_(iConfig.getUntrackedParameter<bool>("debug")),
  h_diffDxyAssocPV(RangeAbs(0.001)),
  h_diffDzAssocPV(RangeAbs(0.001)),
  h_diffCovQoverpQoverp(Range(-1e-6, 0.13), -15, 0),  // despite of ceil in pack, there is rounding in double->float
  h_diffCovLambdaLambda(Range(-1e-6, 0.13), -20, -5), // despite of ceil in pack, there is rounding in double->float
  h_diffCovLambdaDsz(RangeAbs(0.13), -17, -4),
  h_diffCovPhiPhi(RangeAbs(0.13), -15, 0),
  h_diffCovPhiDxy(RangeAbs(0.13), -17, -4),
  h_diffCovDxyDxy(RangeAbs(0.001)),
  h_diffCovDxyDsz(RangeAbs(0.001)),
  h_diffCovDszDsz(RangeAbs(0.001))
{}

PackedCandidateTrackValidator::~PackedCandidateTrackValidator() {}

void PackedCandidateTrackValidator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  edm::ParameterSetDescription desc;

  desc.addUntracked<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
  desc.addUntracked<edm::InputTag>("vertices", edm::InputTag("offlinePrimaryVertices"));
  desc.addUntracked<edm::InputTag>("trackToPackedCandidateAssociation", edm::InputTag("packedPFCandidates"));
  desc.addUntracked<std::string>("rootFolder", "Tracking/PackedCandidate");
  desc.addUntracked<bool>("debug", false);

  descriptions.add("packedCandidateTrackValidator", desc);
}

void PackedCandidateTrackValidator::bookHistograms(DQMStore::IBooker& iBooker, edm::Run const&, edm::EventSetup const&) {
  iBooker.setCurrentFolder(rootFolder_);

  h_selectionFlow = iBooker.book1D("selectionFlow", "Track selection flow", 8, -0.5, 7.5);
  h_selectionFlow->setBinLabel(1, "All tracks");
  h_selectionFlow->setBinLabel(2, "Associated to PackedCandidate");
  h_selectionFlow->setBinLabel(3, "PC is charged"),
  h_selectionFlow->setBinLabel(4, "PC has track");
  h_selectionFlow->setBinLabel(5, "PC is not electron");
  h_selectionFlow->setBinLabel(6, "PC has hits");
  h_selectionFlow->setBinLabel(7, "PC ndof != 0");
  h_selectionFlow->setBinLabel(8, "Track: no missing inner hits");

  constexpr int diffBins = 50;

  h_diffVx = iBooker.book1D("diffVx", "PackedCandidate::bestTrack() - reco::Track in vx()", diffBins, -0.2, 0.2); // not equal
  h_diffVy = iBooker.book1D("diffVy", "PackedCandidate::bestTrack() - reco::Track in vy()", diffBins, -0.2, 0.2); // not equal
  h_diffVz = iBooker.book1D("diffVz", "PackedCandidate::bestTrack() - reco::Track in vz()", diffBins, -0.4, 0.4); // not equal

  h_diffNormalizedChi2 = iBooker.book1D("diffNormalizedChi2", "PackedCandidate::bestTrack() - reco::Track in normalizedChi2()", 30, -1.5, 1.5); // expected difference in -1...0
  h_diffNdof = iBooker.book1D("diffNdof", "PackedCandidate::bestTrack() - reco::Track in ndof()", 33, -30.5, 2.5); // to monitor the difference

  h_diffCharge = iBooker.book1D("diffCharge", "PackedCandidate::bestTrack() - reco::Track in charge()", 5, -2.5, 2.5); // expect equality
  h_diffIsHighPurity = iBooker.book1D("diffIsHighPurity", "PackedCandidate::bestTrack() - reco::Track in quality(highPurity)", 3, -1.5, 1.5); // expect equality

  h_diffPt     = iBooker.book1D("diffPt",     "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in pt()",     diffBins, -1.1, 1.1); // not equal, keep
  h_diffEta    = iBooker.book1D("diffEta",    "PackedCandidate::bestTrack() - reco::Track in eta()",    diffBins, -0.001, 0.001); // not equal, keep
  h_diffPhi    = iBooker.book1D("diffPhi",    "PackedCandidate::bestTrack() - reco::Track in phi()",    diffBins, -0.0005, 0.0005); // expect equality within precision

  h_diffDxyAssocPV.book(iBooker, "diffDxyAssocPV", "(PackedCandidate::dxy() - reco::Track::dxy(assocPV))/reco::Track",
                        40, -0.001, 0.001, // expect equality within precision
                        50, -0.5, 0.5);
  h_diffDzAssocPV.book(iBooker, "diffDzAssocPV","(PackedCandidate::dzAssociatedPV() - reco::Track::dz(assocPV))/reco::Track",
                       40, -0.001, 0.001, // expect equality within precision
                       50, -0.5, 0.5);
  h_diffDxyPV      = iBooker.book1D("diffDxyPV",      "(PackedCandidate::dxy(PV) - reco::Track::dxy(PV))/reco::Track",              diffBins, -0.01, 0.01); // expect equality within precision (worse than assocPV)
  h_diffDzPV       = iBooker.book1D("diffDzPV",       "(PackedCandidate::dz(PV) - reco::Track::dz(PV))/reco::Track",                diffBins, -0.01, 0.01); // expect equality wihtin precision (worse than assocPV)
  h_diffTrackDxyAssocPV   = iBooker.book1D("diffTrackDxyAssocPV",   "(PackedCandidate::bestTrack()::dxy(assocPV)) - reco::Track::dxy(assocPV))/reco::Track",         diffBins, -0.01, 0.01); // not equal
  h_diffTrackDzAssocPV    = iBooker.book1D("diffTrackDzAssocPV",    "(PackedCandidate::bestTrack()::dz(assocPV)) - reco::Track::dz(assocPV))/reco::Track",           diffBins, -0.01, 0.01); // not equal

  h_diffCovQoverpQoverp.book(iBooker, "diffCovQoverpQoverp", "(PackedCandidate::bestTrack() - reco::Track)/reco::track in cov(qoverp, qoverp)",
                             40, -0.05, 0.15, // expect equality within precision (worst precision is exp(1/128*15) =~ 12 %
                             50, -0.5, 0.5); 
  h_diffCovLambdaLambda.book(iBooker, "diffCovLambdaLambda",  "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in cov(lambda, lambda)",
                             40, -0.05, 0.15, // expect equality within precision worst precision is exp(1/128*(20-5)) =~ 12 % (multiplied by pt^2 in packing & unpacking)
                             50, -0.5, 0.5); 
  h_diffCovLambdaDsz.book(iBooker, "diffCovLambdaDsz",     "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in cov(lambda, dsz)",
                          60, -0.15, 0.15, // expect equality within precision, worst precision is exp(1/128*(17-4) =~ 11 %
                          50, -1, 1);
  h_diffCovPhiPhi.book(iBooker, "diffCovPhiPhi",    "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in cov(phi, phi)",
                       40, -0.05, 0.15, // expect equality within precision worst precision is exp(1/128*(20-5)) =~ 12 % (multiplied by pt^2 in packing & unpacking)
                       50, -0.5, 0.5); 
  h_diffCovPhiDxy.book(iBooker, "diffCovPhiDxy",       "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in cov(phi, dxy)",
                       60, -0.15, 0.15, // expect equality within precision, wors precision is exp(1/128)*(17-4) =~ 11 %
                       50, -1, 1);
  h_diffCovDxyDxy.book(iBooker, "diffCovDxyDxy", "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in cov(dxy, dxy)",
                       40, -0.001, 0.001,
                       50, -0.1, 0.1);
  h_diffCovDxyDsz.book(iBooker, "diffCovDxyDsz", "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in cov(dxy, dsz)",
                       40, -0.001, 0.001, // expect equality within precision
                       50, -0.5, 0.5);
  h_diffCovDszDsz.book(iBooker, "diffCovDszDsz", "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in cov(dsz, dsz)",
                       40, -0.001, 0.001, // expect equality within precision
                       50, -0.1, 0.1);

  h_diffDxyError    = iBooker.book1D("diffDxyError", "(PackedCandidate::dxyError() - reco::Track::dxyError())/reco::Track", 40, -0.001, 0.001); // expect equality within precision
  h_diffDszError    = iBooker.book1D("diffDszError", "(PackedCandidate::dzError() - reco::Track::dszError())/reco::Track",  40, -0.001, 0.001);  // ideally, not equal, but for now they are
  h_diffDzError     = iBooker.book1D("diffDzError",  "(PackedCandidate::dzError() - reco::Track::dzError())/reco::Track",   40, -0.001, 0.001); // expect equality within precision (not currently the case)

  h_diffTrackDxyError    = iBooker.book1D("diffTrackDxyError",    "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in dxyError()",    40, -0.001, 0.001); // expect equality within precision
  h_diffTrackDzError     = iBooker.book1D("diffTrackDzError",     "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in dzError()",     40, -0.05, 0.05); // not equal

  h_diffPtError     = iBooker.book1D("diffPtError",     "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in ptError()",     diffBins, -1.1, 1.1); // not equal
  h_diffEtaError    = iBooker.book1D("diffEtaError",    "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in etaError()",    60, -0.15, 0.15); // not equal


  h_diffNumberOfPixelLayers = iBooker.book1D("diffNumberOfPixelLayers", "PackedCandidate::pixelLayersWithMeasurement() - reco::Track::hitPattern::pixelLayersWithMeasurement()", 5, -2.5, 2.5); // expect equality 
  h_diffNumberOfStripLayers = iBooker.book1D("diffNumberOfStripLayers", "PackedCandidate::stripLayersWithMeasurement() - reco::Track::hitPattern::stripLayersWithMeasurement()",             5, -2.5, 2.5); // expect equality
  h_diffNumberOfPixelHits = iBooker.book1D("diffNumberOfPixelHits", "PackedCandidate::numberOfPixelHits() - reco::Track::hitPattern::numberOfValidPixelHits()", 5, -2.5, 2.5); // expect equality 
  h_diffNumberOfHits      = iBooker.book1D("diffNumberOfHits",      "PackedCandidate::numberHits() - reco::Track::hitPattern::numberOfValidHits()",             5, -2.5, 2.5); // expect equality
  h_diffLostInnerHits     = iBooker.book1D("diffLostInnerHits",     "PackedCandidate::lostInnerHits() - reco::Track::hitPattern::numberOfLostHits(MISSING_INNER_HITS)",      5, -2.5, 2.5); // expect equality

  h_diffHitPatternPixelLayersWithMeasurement = iBooker.book1D("diffHitPatternPixelLayersWithMeasurement", "PackedCandidate::bestTrack() - reco::Track in hitPattern::pixelLayersWithMeasurement()",   13, -10.5, 2.5); // not equal
  h_diffHitPatternStripLayersWithMeasurement = iBooker.book1D("diffHitPatternStripLayersWithMeasurement", "PackedCandidate::bestTrack() - reco::Track in hitPattern::stripLayersWithMeasurement()",   13, -10.5, 2.5); // not equal
  h_diffHitPatternTrackerLayersWithMeasurement = iBooker.book1D("diffHitPatternTrackerLayersWithMeasurement", "PackedCandidate::bestTrack() - reco::Track in hitPattern::trackerLayersWithMeasurement()",   13, -10.5, 2.5); // not equal
  h_diffHitPatternNumberOfValidPixelHits = iBooker.book1D("diffHitPatternNumberOfValidPixelHits", "PackedCandidate::bestTrack() - reco::Track in hitPattern::numberOfValidPixelHits()",   13, -10.5, 2.5); // not equal
  h_diffHitPatternNumberOfValidHits      = iBooker.book1D("diffHitPatternNumberOfValidHits",      "PackedCandidate::bestTrack() - reco::Track in hitPattern::numberOfValidHits()",      13, -10.5, 2.5); // not equal
  h_diffHitPatternNumberOfLostInnerHits  = iBooker.book1D("diffHitPatternNumberOfLostPixelHits",  "PackedCandidate::bestTrack() - reco::Track in hitPattern::numberOfLostHits(MISSING_INNER_HITS)", 13, -10.5, 2.5); // not equal
  h_diffHitPatternHasValidHitInFirstPixelBarrel = iBooker.book1D("diffHitPatternHasValidHitInFirstPixelBarrel", "PackedCandidate::bestTrack() - reco::Track in hitPattern::hasValidHitInFirstPixelBarrel", 3, -1.5, 1.5); // expect equality

  h_numberPixelLayersOverMax = iBooker.book1D("numberPixelLayersOverMax", "Number of pixel layers over the maximum of PackedCandidate", 10, 0, 10);
  h_numberStripLayersOverMax = iBooker.book1D("numberStripLayersOverMax", "Number of strip layers over the maximum of PackedCandidate", 10, 0, 10);
  h_numberLayersOverMax = iBooker.book1D("numberLayersOverMax", "Number of layers over the maximum of PackedCandidate", 20, 0, 20);
  h_numberPixelHitsOverMax = iBooker.book1D("numberPixelHitsOverMax", "Number of pixel hits over the maximum of PackedCandidate", 10, 0, 10);
  h_numberStripHitsOverMax = iBooker.book1D("numberStripHitsOverMax", "Number of strip hits over the maximum of PackedCandidate", 10, 0, 10);
  h_numberHitsOverMax = iBooker.book1D("numberHitsOverMax", "Number of hits over the maximum of PackedCandidate", 20, 0, 20);
}

void PackedCandidateTrackValidator::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
  edm::Handle<edm::View<reco::Track>> htracks;
  iEvent.getByToken(tracksToken_, htracks);
  const auto& tracks = *htracks;

  edm::Handle<reco::VertexCollection> hvertices;
  iEvent.getByToken(verticesToken_, hvertices);
  const auto& vertices = *hvertices;

  if(vertices.empty())
    return;
  const reco::Vertex& pv = vertices[0];

  edm::Handle<edm::Association<pat::PackedCandidateCollection>> hassoc;
  iEvent.getByToken(trackToPackedCandidateToken_, hassoc);
  const auto& trackToPackedCandidate = *hassoc;

  for(size_t i=0; i<tracks.size(); ++i) {
    auto trackPtr = tracks.ptrAt(i);
    const reco::Track& track = *trackPtr;
    h_selectionFlow->Fill(sf_AllTracks);

    pat::PackedCandidateRef pcRef = trackToPackedCandidate[trackPtr];
    if(pcRef.isNull()) {
      continue;
    }
    h_selectionFlow->Fill(sf_AssociatedToPC);

    // Filter out neutral PackedCandidates, some of them may have track associated, and for those the charge comparison fails
    if(pcRef->charge() == 0) {
      continue;
    }
    h_selectionFlow->Fill(sf_PCIsCharged);

    const reco::Track *trackPcPtr = pcRef->bestTrack();
    if(!trackPcPtr) {
      continue;
    }
    h_selectionFlow->Fill(sf_PCHasTrack);

    // Filter out electrons to avoid comparisons to PackedCandidates with GsfTrack
    if(std::abs(pcRef->pdgId()) == 11) {
      continue;
    }
    h_selectionFlow->Fill(sf_PCIsNotElectron);

    // Filter out PackedCandidate-tracks with no hits, as they won't have their details filled
    const reco::Track& trackPc = *trackPcPtr;
    if(trackPc.hitPattern().numberOfValidHits() == 0) {
      continue;
    }
    h_selectionFlow->Fill(sf_PCHasHits);

    auto slimmedVertexRef = pcRef->vertexRef();
    const reco::Vertex& pcVertex = vertices[slimmedVertexRef.key()];

    fillNoFlow(h_diffVx, trackPc.vx() - track.vx());
    fillNoFlow(h_diffVy, trackPc.vy() - track.vy());
    fillNoFlow(h_diffVz, trackPc.vz() - track.vz());

    // PackedCandidate recalculates the ndof in unpacking as
    // (nhits+npixelhits-5), but some strip hits may have dimension 2.
    // If PackedCandidate has ndof=0, the resulting normalizedChi2
    // will be 0 too. Hence, the comparison makes sense only for those
    // PackedCandidates that have ndof != 0.
    double diffNormalizedChi2 = 0;
    if(trackPc.ndof() != 0) {
      h_selectionFlow->Fill(sf_PCNdofNot0);
      diffNormalizedChi2 = trackPc.normalizedChi2() - track.normalizedChi2();
      fillNoFlow(h_diffNormalizedChi2, diffNormalizedChi2);
    }
    fillNoFlow(h_diffNdof, trackPc.ndof() - track.ndof());

    auto diffCharge = trackPc.charge() - track.charge();
    fillNoFlow(h_diffCharge, diffCharge);
    int diffHP = static_cast<int>(trackPc.quality(reco::TrackBase::highPurity)) - static_cast<int>(track.quality(reco::TrackBase::highPurity));
    fillNoFlow(h_diffIsHighPurity,  diffHP);

    const auto diffPt = diffRelative(trackPc.pt(), track.pt());
    const auto diffPhi = reco::deltaPhi(trackPc.phi(), track.phi());
    fillNoFlow(h_diffPt , diffPt);
    fillNoFlow(h_diffEta, trackPc.eta() - track.eta());
    fillNoFlow(h_diffPhi, diffPhi);

    const auto diffDxyAssocPV = h_diffDxyAssocPV.fill(pcRef->dxy( )          , track.dxy(pcVertex.position()), [](double value){return value*100.;});
    const auto diffDzAssocPV  = h_diffDzAssocPV.fill( pcRef->dzAssociatedPV(), track.dz(pcVertex.position()) , [](double value){return value*100.;});
    const auto diffDxyPV      = diffRelative(pcRef->dxy(pv.position()), track.dxy(pv.position()));
    const auto diffDzPV       = diffRelative(pcRef->dz(pv.position()) , track.dz(pv.position()));
    fillNoFlow(h_diffDxyPV   , diffDxyPV);
    fillNoFlow(h_diffDzPV    , diffDzPV);
    fillNoFlow(h_diffTrackDxyAssocPV, diffRelative(trackPc.dxy(pcVertex.position()), track.dxy(pcVertex.position())));
    fillNoFlow(h_diffTrackDzAssocPV , diffRelative(trackPc.dz(pcVertex.position()) , track.dz(pcVertex.position())));

    auto fillCov1 = [&](auto& hlp, const int i, const int j) {
      return hlp.fill(trackPc.covariance(i, j), track.covariance(i, j));
    };
    auto fillCov2 = [&](auto& hlp, const int i, const int j, std::function<double(double)> modifyPack) {
      return hlp.fill(trackPc.covariance(i, j), track.covariance(i, j), modifyPack);
    };
    auto fillCov3 = [&](auto& hlp, const int i, const int j, std::function<double(double)> modifyPack, std::function<double(double)> modifyUnpack) {
      return hlp.fill(trackPc.covariance(i, j), track.covariance(i, j), modifyPack, modifyUnpack);
    };

    const auto pcPt = pcRef->pt();
    const auto diffCovQoverpQoverp = fillCov3(h_diffCovQoverpQoverp, reco::TrackBase::i_qoverp, reco::TrackBase::i_qoverp, [=](double val){return val*pcPt*pcPt;}, [=](double val){return val/pcPt/pcPt;});
    const auto diffCovLambdaLambda = fillCov1(h_diffCovLambdaLambda, reco::TrackBase::i_lambda, reco::TrackBase::i_lambda);
    const auto diffCovLambdaDsz    = fillCov1(h_diffCovLambdaDsz,    reco::TrackBase::i_lambda, reco::TrackBase::i_dsz);
    const auto diffCovPhiPhi       = fillCov3(h_diffCovPhiPhi,       reco::TrackBase::i_phi,    reco::TrackBase::i_phi,    [=](double val){return val*pcPt*pcPt;}, [=](double val){return val/pcPt/pcPt;});
    const auto diffCovPhiDxy       = fillCov1(h_diffCovPhiDxy,       reco::TrackBase::i_phi,    reco::TrackBase::i_dxy);
    const auto diffCovDxyDxy       = fillCov2(h_diffCovDxyDxy,       reco::TrackBase::i_dxy,    reco::TrackBase::i_dxy,    [](double value){return value*10000.;});
    const auto diffCovDxyDsz       = fillCov2(h_diffCovDxyDsz,       reco::TrackBase::i_dxy,    reco::TrackBase::i_dsz,    [](double value){return value*10000.;});
    const auto diffCovDszDsz       = fillCov2(h_diffCovDszDsz,       reco::TrackBase::i_dsz,    reco::TrackBase::i_dsz,    [](double value){return value*10000.;});

    if(isInRange(diffCovDszDsz.status())) {
      fillNoFlow(h_diffDszError, diffRelative(pcRef->dzError(), track.dszError()));
      fillNoFlow(h_diffDzError,  diffRelative(pcRef->dzError(), track.dzError()));
      fillNoFlow(h_diffTrackDzError, diffRelative(trackPc.dzError(), track.dzError()));
    }
    if(isInRange(diffCovDxyDxy.status())) {
      fillNoFlow(h_diffDxyError, diffRelative(pcRef->dxyError(), track.dxyError()));
      fillNoFlow(h_diffTrackDxyError, diffRelative(trackPc.dxyError(), track.dxyError()));
    }
    fillNoFlow(h_diffPtError    , diffRelative(trackPc.ptError()    , track.ptError()    ));
    fillNoFlow(h_diffEtaError   , diffRelative(trackPc.etaError()   , track.etaError()   ));



    // For the non-HitPattern ones, take into account the PackedCandidate packing precision
    const auto trackNumberOfHits = track.hitPattern().numberOfValidHits();
    const auto trackNumberOfPixelHits = track.hitPattern().numberOfValidPixelHits();
    const auto trackNumberOfStripHits = track.hitPattern().numberOfValidStripHits();
    const auto pcNumberOfHits = pcRef->numberOfHits();
    const auto pcNumberOfPixelHits = pcRef->numberOfPixelHits();
    const auto pcNumberOfStripHits = pcNumberOfHits - pcNumberOfPixelHits;
    const auto trackNumberOfLayers = track.hitPattern().trackerLayersWithMeasurement();
    const auto trackNumberOfPixelLayers = track.hitPattern().pixelLayersWithMeasurement();
    const auto trackNumberOfStripLayers = track.hitPattern().stripLayersWithMeasurement();
    const auto pcNumberOfLayers = pcRef->trackerLayersWithMeasurement();
    const auto pcNumberOfPixelLayers = pcRef->pixelLayersWithMeasurement();
    const auto pcNumberOfStripLayers = pcRef->stripLayersWithMeasurement();

    // layer number overflow (should be zero)
    const int pixelLayerOverflow = trackNumberOfPixelLayers > pat::PackedCandidate::trackPixelHitsMask ? trackNumberOfPixelLayers - pat::PackedCandidate::trackPixelHitsMask : 0;
    const int stripLayerOverflow = trackNumberOfStripLayers > pat::PackedCandidate::trackStripHitsMask ? trackNumberOfStripLayers - pat::PackedCandidate::trackStripHitsMask : 0;
    const int layerOverflow = trackNumberOfLayers > (pat::PackedCandidate::trackPixelHitsMask+pat::PackedCandidate::trackStripHitsMask) ? trackNumberOfLayers - (pat::PackedCandidate::trackPixelHitsMask+pat::PackedCandidate::trackStripHitsMask) : 0;

    // hit overflow (should also be zero)
    const int pixelOverflow = trackNumberOfPixelHits - pcNumberOfPixelLayers > pat::PackedCandidate::trackPixelHitsMask ? trackNumberOfPixelHits - pcNumberOfPixelLayers - pat::PackedCandidate::trackPixelHitsMask : 0;
    const int stripOverflow = trackNumberOfStripHits - pcNumberOfStripLayers > pat::PackedCandidate::trackStripHitsMask ? trackNumberOfStripHits - pcNumberOfStripLayers - pat::PackedCandidate::trackStripHitsMask : 0;
    const int hitsOverflow = trackNumberOfHits - pcNumberOfLayers > (pat::PackedCandidate::trackPixelHitsMask+pat::PackedCandidate::trackStripHitsMask) ? trackNumberOfHits - pcNumberOfLayers - (pat::PackedCandidate::trackPixelHitsMask+pat::PackedCandidate::trackStripHitsMask) : 0;
    // PackedCandidate counts overflow pixel hits as strip
    const int pixelInducedStripOverflow = (trackNumberOfStripHits+pixelOverflow-pcNumberOfStripLayers) > pat::PackedCandidate::trackStripHitsMask ? (trackNumberOfStripHits+pixelOverflow-stripOverflow-pcNumberOfStripLayers) - pat::PackedCandidate::trackStripHitsMask : 0; 
    h_numberPixelLayersOverMax->Fill(pixelLayerOverflow);
    h_numberStripLayersOverMax->Fill(stripLayerOverflow);
    h_numberLayersOverMax->Fill(layerOverflow);
    h_numberPixelHitsOverMax->Fill(pixelOverflow);
    h_numberStripHitsOverMax->Fill(stripOverflow);
    h_numberHitsOverMax->Fill(hitsOverflow);

    int diffNumberOfPixelHits = 0;
    int diffNumberOfHits = 0;
    int diffNumberOfPixelLayers = 0;
    int diffNumberOfStripLayers = 0;
    if(pixelLayerOverflow) {
      diffNumberOfPixelLayers = pcNumberOfPixelLayers - pat::PackedCandidate::trackPixelHitsMask;
    }
    else {
      diffNumberOfPixelLayers = pcNumberOfPixelLayers - trackNumberOfPixelLayers;
    }
    if(pixelOverflow) {
      diffNumberOfPixelHits = pcNumberOfPixelHits - pcNumberOfPixelLayers - pat::PackedCandidate::trackPixelHitsMask;
    }
    else {
      diffNumberOfPixelHits = pcNumberOfPixelHits - trackNumberOfPixelHits;
    }
    if(stripLayerOverflow) {
       diffNumberOfStripLayers = pcNumberOfStripLayers - pat::PackedCandidate::trackStripHitsMask;
    }
    else {
      diffNumberOfStripLayers = pcNumberOfStripLayers - trackNumberOfStripLayers;
    }
    if(stripOverflow || pixelInducedStripOverflow || pixelOverflow) {
      int diffNumberOfStripHits = 0;
      if(stripOverflow || pixelInducedStripOverflow) {
        diffNumberOfStripHits = pcNumberOfStripHits - pat::PackedCandidate::trackStripHitsMask;
      }
      else if(pixelOverflow) {
        diffNumberOfStripHits = (pcNumberOfStripHits - pixelOverflow) - trackNumberOfStripHits;
      }

      diffNumberOfHits = diffNumberOfPixelHits + diffNumberOfStripHits;
    }
    else {
      diffNumberOfHits = pcNumberOfHits - trackNumberOfHits;
    }



    fillNoFlow(h_diffNumberOfPixelHits, diffNumberOfPixelHits);
    fillNoFlow(h_diffNumberOfHits, diffNumberOfHits);
    fillNoFlow(h_diffNumberOfPixelLayers, diffNumberOfPixelLayers);
    fillNoFlow(h_diffNumberOfStripLayers, diffNumberOfStripLayers);

    int diffLostInnerHits = 0;
    const auto trackLostInnerHits = track.hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS);
    switch(pcRef->lostInnerHits()) {
    case pat::PackedCandidate::validHitInFirstPixelBarrelLayer:
    case pat::PackedCandidate::noLostInnerHits:
      diffLostInnerHits = -trackLostInnerHits;
      break;
    case pat::PackedCandidate::oneLostInnerHit:
      diffLostInnerHits = 1-trackLostInnerHits;
      break;
    case pat::PackedCandidate::moreLostInnerHits:
      diffLostInnerHits = trackLostInnerHits>=2 ? 0 : 2-trackLostInnerHits;
      break;
    }
    fillNoFlow(h_diffLostInnerHits, diffLostInnerHits);

    // For HitPattern ones, calculate the full diff (i.e. some differences are expected)
    auto diffHitPatternPixelLayersWithMeasurement = trackPc.hitPattern().pixelLayersWithMeasurement() - trackNumberOfPixelLayers;
    fillNoFlow(h_diffHitPatternPixelLayersWithMeasurement, diffHitPatternPixelLayersWithMeasurement);
    auto diffHitPatternStripLayersWithMeasurement = trackPc.hitPattern().stripLayersWithMeasurement() - trackNumberOfStripLayers;
    fillNoFlow(h_diffHitPatternStripLayersWithMeasurement, diffHitPatternStripLayersWithMeasurement);
    auto diffHitPatternTrackerLayersWithMeasurement = trackPc.hitPattern().trackerLayersWithMeasurement() - trackNumberOfLayers;
    fillNoFlow(h_diffHitPatternTrackerLayersWithMeasurement, diffHitPatternTrackerLayersWithMeasurement);
    auto diffHitPatternNumberOfValidPixelHits = trackPc.hitPattern().numberOfValidPixelHits() - trackNumberOfPixelHits;
    fillNoFlow(h_diffHitPatternNumberOfValidPixelHits, diffHitPatternNumberOfValidPixelHits);
    auto diffHitPatternNumberOfValidHits = trackPc.hitPattern().numberOfValidHits() - trackNumberOfHits;
    fillNoFlow(h_diffHitPatternNumberOfValidHits, diffHitPatternNumberOfValidHits);
    fillNoFlow(h_diffHitPatternNumberOfLostInnerHits, trackPc.hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS) - track.hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS));

    // hasValidHitInFirstPixelBarrel is set only if numberOfLostHits(MISSING_INNER_HITS) == 0
    int diffHitPatternHasValidHitInFirstPixelBarrel = 0;
    if(track.hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS) == 0) {
      h_selectionFlow->Fill(sf_NoMissingInnerHits);
      diffHitPatternHasValidHitInFirstPixelBarrel = static_cast<int>(trackPc.hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector::PixelBarrel, 1)) - static_cast<int>(track.hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector::PixelBarrel, 1));
      fillNoFlow(h_diffHitPatternHasValidHitInFirstPixelBarrel, diffHitPatternHasValidHitInFirstPixelBarrel);
    }

    // Print warning if there are differences outside the expected range
    if(debug_ && 
       (diffNormalizedChi2 < -1 || diffNormalizedChi2 > 0 || diffCharge != 0 || diffHP != 0
       || std::abs(diffPhi) > 5e-4
       || diffDxyAssocPV.outsideExpectedRange()
       || diffDzAssocPV.outsideExpectedRange()
       || std::abs(diffDxyPV) > 0.05 || std::abs(diffDzPV) > 0.05
       || diffCovQoverpQoverp.outsideExpectedRange() || diffCovLambdaLambda.outsideExpectedRange()
       || diffCovLambdaDsz.outsideExpectedRange() || diffCovPhiPhi.outsideExpectedRange()
       || diffCovPhiDxy.outsideExpectedRange() || diffCovDxyDxy.outsideExpectedRange()
       || diffCovDxyDsz.outsideExpectedRange() || diffCovDszDsz.outsideExpectedRange()
       || diffNumberOfPixelHits != 0 || diffNumberOfHits != 0 || diffLostInnerHits != 0
       || diffHitPatternHasValidHitInFirstPixelBarrel != 0
        )
       ) {

      edm::LogInfo("PackedCandidateTrackValidator") << "Track " << i << " pt " << track.pt() << " eta " << track.eta() << " phi " << track.phi() << " chi2 " << track.chi2() << " ndof " << track.ndof()
                                                    << "\n"
                                                    << "  ptError " << track.ptError() << " etaError " << track.etaError() << " phiError " << track.phiError() << " dxyError " << track.dxyError() << " dzError " << track.dzError()
                                                    << "\n"
                                                    << "  refpoint " << track.referencePoint() << " momentum " << track.momentum()
                                                    << "\n"
                                                    << "  dxy " << track.dxy() << " dz " << track.dz()
                                                    << "\n"
                                                    << "  " << TrackAlgoPrinter(track)
                                                    << " lost inner hits " << trackLostInnerHits
                                                    << " lost outer hits " << track.hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS)
                                                    << " hitpattern " << HitPatternPrinter(track)
                                                    << " \n"
                                                    << " PC " << pcRef.id() << ":" << pcRef.key() << " track pt " << trackPc.pt() << " eta " << trackPc.eta() << " phi " << trackPc.phi() << " (PC " << pcRef->phi() << ") chi2 " << trackPc.chi2() << " ndof " << trackPc.ndof() << " pdgId " << pcRef->pdgId() << " mass " << pcRef->mass()
                                                    << "\n"
                                                    << "  ptError " << trackPc.ptError() << " etaError " << trackPc.etaError() << " phiError " << trackPc.phiError()
                                                    << "\n"
                                                    << "  pc.vertex " << pcRef->vertex() << " momentum " << pcRef->momentum() << " track " << trackPc.momentum()
                                                    << "\n"
                                                    << "  dxy " << trackPc.dxy() << " dz " << trackPc.dz() << " pc.dz " << pcRef->dz()
                                                    << " dxyError " << trackPc.dxyError() << " dzError " << trackPc.dzError()
                                                    << "\n"
                                                    << "  dxy(PV) " << trackPc.dxy(pv.position()) << " dz(PV) " << trackPc.dz(pv.position())
                                                    << " dxy(assocPV) " << trackPc.dxy(pcVertex.position()) << " dz(assocPV) " << trackPc.dz(pcVertex.position())
                                                    << "\n"
                                                    << " (diff PackedCandidate track)"
                                                    << " highPurity " << diffHP << " " << trackPc.quality(reco::TrackBase::highPurity) << " " << track.quality(reco::TrackBase::highPurity)
                                                    << " charge " << diffCharge << " " << trackPc.charge() << " " << track.charge()
                                                    << " normalizedChi2 " << diffNormalizedChi2 << " " << trackPc.normalizedChi2() << " " << track.normalizedChi2()
                                                    << "\n "
                                                    << " numberOfAllHits " << diffNumberOfHits << " " << pcNumberOfHits << " " << trackNumberOfHits
                                                    << " numberOfPixelHits " << diffNumberOfPixelHits << " " << pcNumberOfPixelHits << " " << trackNumberOfPixelHits
                                                    << " numberOfStripHits # " << pcNumberOfStripHits << " " << trackNumberOfStripHits
                                                    << "\n "
                                                    << " hitPattern.numberOfValidPixelHits " << diffHitPatternNumberOfValidPixelHits << " " << trackPc.hitPattern().numberOfValidPixelHits() << " " << track.hitPattern().numberOfValidPixelHits()
                                                    << " hitPattern.numberOfValidHits " << diffHitPatternNumberOfValidHits << " " << trackPc.hitPattern().numberOfValidHits() << " " << track.hitPattern().numberOfValidHits()
                                                    << " hitPattern.hasValidHitInFirstPixelBarrel " << diffHitPatternHasValidHitInFirstPixelBarrel << " " << trackPc.hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector::PixelBarrel, 1) << " " << track.hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector::PixelBarrel, 1)
                                                    << "\n "
                                                    << " lostInnerHits  " << diffLostInnerHits << " " << pcRef->lostInnerHits() << " #"
                                                    << " phi (5e-4) " << diffPhi << " " << trackPc.phi() << " " << track.phi()
                                                    << "\n "
                                                    << " dxy(assocPV) " << diffDxyAssocPV
                                                    << "\n "
                                                    << " dz(assocPV) " << diffDzAssocPV
                                                    << "\n "
                                                    << " dxy(PV) (0.05) " << diffDxyPV << " " << pcRef->dxy(pv.position()) << " " << track.dxy(pv.position())
                                                    << "\n "
                                                    << " dz(PV) (0.05) " << diffDzPV << " " << pcRef->dz(pv.position()) << " " << track.dz(pv.position())
                                                    << "\n "
                                                    << " cov(qoverp, qoverp)  " << diffCovQoverpQoverp
                                                    << "\n "
                                                    << " cov(lambda, lambda) " << diffCovLambdaLambda
                                                    << "\n "
                                                    << " cov(lambda, dsz) " << diffCovLambdaDsz
                                                    << "\n "
                                                    << " cov(phi, phi) " << diffCovPhiPhi
                                                    << "\n "
                                                    << " cov(phi, dxy) " << diffCovPhiDxy
                                                    << "\n "
                                                    << " cov(dxy, dxy) " << diffCovDxyDxy
                                                    << "\n "
                                                    << " cov(dxy, dsz) " << diffCovDxyDsz
                                                    << "\n "
                                                    << " cov(dsz, dsz) " << diffCovDszDsz;
    }
  }
}

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