Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:09:06 

0001 #include "FWCore/Framework/interface/Event.h"
0002 #include "FWCore/Framework/interface/ESHandle.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FWCore/Utilities/interface/VecArray.h"
0008 #include "FWCore/Utilities/interface/isFinite.h"
0009 
0010 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0011 #include "DQMServices/Core/interface/DQMStore.h"
0012 
0013 #include "DataFormats/Common/interface/Association.h"
0014 #include "DataFormats/TrackReco/interface/Track.h"
0015 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0016 #include "DataFormats/Math/interface/libminifloat.h"
0017 #include "DataFormats/Math/interface/liblogintpack.h"
0018 #include "DataFormats/VertexReco/interface/Vertex.h"
0019 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0020 #include "DataFormats/Math/interface/deltaPhi.h"
0021 
0022 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0023 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0024 
0025 #include <iomanip>
0026 
0027 namespace {
0028   using dqm::reco::DQMStore;
0029   using dqm::reco::MonitorElement;
0030 
0031   template <typename T>
0032   void fillNoFlow(MonitorElement* me, T val) {
0033     auto h = me->getTH1();
0034     const auto xaxis = h->GetXaxis();
0035     if (val <= xaxis->GetXmin())
0036       h->AddBinContent(xaxis->GetFirst());
0037     else if (val >= xaxis->GetXmax())
0038       h->AddBinContent(xaxis->GetLast());
0039     else
0040       h->Fill(val);
0041   }
0042 
0043   class HitPatternPrinter {
0044   public:
0045     explicit HitPatternPrinter(const reco::Track& trk) : track(trk) {}
0046 
0047     void print(std::ostream& os) const {
0048       const reco::HitPattern& p = track.hitPattern();
0049 
0050       for (int i = 0; i < p.numberOfAllHits(reco::HitPattern::TRACK_HITS); ++i) {
0051         uint32_t hit = p.getHitPattern(reco::HitPattern::TRACK_HITS, i);
0052 
0053         detLayer(os, p, hit);
0054         if (p.missingHitFilter(hit)) {
0055           os << "(miss)";
0056         } else if (p.inactiveHitFilter(hit)) {
0057           os << "(inact)";
0058         } else if (p.badHitFilter(hit)) {
0059           os << "(bad)";
0060         }
0061         os << " ";
0062       }
0063 
0064       if (p.numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS) > 0) {
0065         os << "lost inner ";
0066 
0067         for (int i = 0; i < p.numberOfAllHits(reco::HitPattern::MISSING_INNER_HITS); ++i) {
0068           uint32_t hit = p.getHitPattern(reco::HitPattern::MISSING_INNER_HITS, i);
0069           detLayer(os, p, hit);
0070           if (p.missingHitFilter(hit)) {
0071             os << "(miss)";
0072           } else if (p.inactiveHitFilter(hit)) {
0073             os << "(inact)";
0074           }
0075         }
0076       }
0077       if (p.numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS) > 0) {
0078         os << "lost outer ";
0079 
0080         for (int i = 0; i < p.numberOfAllHits(reco::HitPattern::MISSING_OUTER_HITS); ++i) {
0081           uint32_t hit = p.getHitPattern(reco::HitPattern::MISSING_OUTER_HITS, i);
0082           detLayer(os, p, hit);
0083           if (p.missingHitFilter(hit)) {
0084             os << "(miss)";
0085           } else if (p.inactiveHitFilter(hit)) {
0086             os << "(inact)";
0087           }
0088         }
0089       }
0090     }
0091 
0092   private:
0093     static void detLayer(std::ostream& os, const reco::HitPattern& p, uint32_t hit) {
0094       if (p.pixelBarrelHitFilter(hit)) {
0095         os << "BPIX";
0096       } else if (p.pixelEndcapHitFilter(hit)) {
0097         os << "FPIX";
0098       } else if (p.stripTIBHitFilter(hit)) {
0099         os << "TIB";
0100       } else if (p.stripTIDHitFilter(hit)) {
0101         os << "TID";
0102       } else if (p.stripTOBHitFilter(hit)) {
0103         os << "TOB";
0104       } else if (p.stripTECHitFilter(hit)) {
0105         os << "TEC";
0106       }
0107       os << p.getLayer(hit);
0108     }
0109 
0110     const reco::Track& track;
0111   };
0112 
0113   std::ostream& operator<<(std::ostream& os, const HitPatternPrinter& hpp) {
0114     hpp.print(os);
0115     return os;
0116   }
0117 
0118   class TrackAlgoPrinter {
0119   public:
0120     explicit TrackAlgoPrinter(const reco::Track& trk) : track(trk) {}
0121 
0122     void print(std::ostream& os) const {
0123       edm::VecArray<reco::TrackBase::TrackAlgorithm, reco::TrackBase::algoSize> algos;
0124       for (int ialgo = 0; ialgo < reco::TrackBase::algoSize; ++ialgo) {
0125         auto algo = static_cast<reco::TrackBase::TrackAlgorithm>(ialgo);
0126         if (track.isAlgoInMask(algo)) {
0127           algos.push_back(algo);
0128         }
0129       }
0130 
0131       os << "algo " << reco::TrackBase::algoName(track.algo());
0132       if (track.originalAlgo() != track.algo())
0133         os << " originalAlgo " << reco::TrackBase::algoName(track.originalAlgo());
0134       if (algos.size() > 1) {
0135         os << " algoMask";
0136         for (auto algo : algos) {
0137           os << " " << reco::TrackBase::algoName(algo);
0138         }
0139       }
0140     }
0141 
0142   private:
0143     const reco::Track& track;
0144   };
0145   std::ostream& operator<<(std::ostream& os, const TrackAlgoPrinter& tap) {
0146     tap.print(os);
0147     return os;
0148   }
0149 
0150   double diffRelative(double a, double b) { return (a - b) / b; }
0151 
0152   class LogIntHelper {
0153   public:
0154     LogIntHelper(double lmin, double lmax) : lmin_(lmin), lmax_(lmax) {}
0155 
0156     class UnderOverflow {
0157     public:
0158       UnderOverflow(double largestValue, double smallestValue, std::function<double(double)> modifyUnpack)
0159           : unpackedLargestValue_(modifyUnpack ? modifyUnpack(largestValue) : largestValue),
0160             unpackedSmallestValue_(modifyUnpack ? modifyUnpack(smallestValue) : smallestValue) {}
0161 
0162       bool compatibleWithUnderflow(double value) const { return value == unpackedSmallestValue_; }
0163       void printNonOkUnderflow(std::ostream& os) const { os << " (not min " << unpackedSmallestValue_ << ")"; }
0164 
0165       bool compatibleWithOverflow(double value) const { return value == unpackedLargestValue_; }
0166       void printNonOkOverflow(std::ostream& os) const { os << " (not max " << unpackedLargestValue_ << ")"; }
0167 
0168     private:
0169       // narrow to float to compare apples to apples with values from
0170       // PackedCandidate (even though the final comparison is done in
0171       // double)
0172       const float unpackedLargestValue_;
0173       const float unpackedSmallestValue_;
0174     };
0175 
0176     static std::string maxName() { return "max"; }
0177     static std::string minName() { return "min"; }
0178 
0179     UnderOverflow underOverflowHelper(double value, std::function<double(double)> modifyUnpack) const {
0180       return UnderOverflow(
0181           largestValue(), value >= 0 ? smallestPositiveValue() : std::abs(smallestNegativeValue()), modifyUnpack);
0182     }
0183 
0184     double largestValue() const { return logintpack::unpack8log(127, lmin_, lmax_); }
0185 
0186     static bool wouldBeDenorm(double value) { return false; }
0187 
0188     // lessThan means closer to zero
0189     bool lessThanSmallestValue(double value) const {
0190       if (value >= 0)
0191         return value < smallestPositiveValue();
0192       else
0193         return value > smallestNegativeValue();
0194     }
0195 
0196     double smallestPositiveValue() const { return logintpack::unpack8log(logintpack::smallestPositive, lmin_, lmax_); }
0197 
0198     double smallestNegativeValue() const { return logintpack::unpack8log(logintpack::smallestNegative, lmin_, lmax_); }
0199 
0200   private:
0201     const double lmin_;
0202     const double lmax_;
0203   };
0204   class Float16Helper {
0205   public:
0206     class UnderOverflow {
0207     public:
0208       static void printNonOkUnderflow(std::ostream& os) { os << " (not 0)"; }
0209       static bool compatibleWithUnderflow(double value) { return value == 0.0; }
0210       static void printNonOkOverflow(std::ostream& os) { os << " (not inf)"; }
0211       static bool compatibleWithOverflow(double value) { return edm::isNotFinite(value); }
0212     };
0213 
0214     static std::string maxName() { return "inf"; }
0215     static std::string minName() { return "0"; }
0216 
0217     static UnderOverflow underOverflowHelper(double value, std::function<double(double)>) { return UnderOverflow(); }
0218 
0219     static double largestValue() { return MiniFloatConverter::max32RoundedToMax16(); }
0220 
0221     static bool wouldBeDenorm(double value) {
0222       const float valuef = static_cast<float>(value);
0223       return valuef >= MiniFloatConverter::denorm_min() && valuef < MiniFloatConverter::min();
0224     }
0225 
0226     static bool lessThanSmallestValue(double value) { return std::abs(value) < smallestValue(); }
0227 
0228     static double smallestValue() { return MiniFloatConverter::denorm_min(); }
0229   };
0230 
0231   enum class RangeStatus {
0232     inrange = 0,
0233     inrange_signflip = 1,
0234     denormal = 2,
0235     underflow_OK = 3,
0236     underflow_notOK = 4,
0237     overflow_OK = 5,
0238     overflow_notOK = 6
0239   };
0240   bool isInRange(RangeStatus status) { return status == RangeStatus::inrange || status == RangeStatus::denormal; }
0241 
0242   template <typename T>
0243   class PackedValueCheckResult {
0244   public:
0245     PackedValueCheckResult(RangeStatus status,
0246                            double diff,
0247                            double pcvalue,
0248                            double trackvalue,
0249                            double rangeMin,
0250                            double rangeMax,
0251                            const typename T::UnderOverflow& underOverflow)
0252         : diff_(diff),
0253           pcvalue_(pcvalue),
0254           trackvalue_(trackvalue),
0255           rangeMin_(rangeMin),
0256           rangeMax_(rangeMax),
0257           status_(status),
0258           underOverflow_(underOverflow) {}
0259 
0260     RangeStatus status() const { return status_; }
0261     double diff() const { return diff_; }
0262 
0263     bool outsideExpectedRange() const {
0264       if (status_ == RangeStatus::inrange)
0265         return diff_ < rangeMin_ || diff_ > rangeMax_;
0266       // denormal is considered as "in range" regardless of the expected range
0267       return status_ == RangeStatus::underflow_notOK || status_ == RangeStatus::overflow_notOK ||
0268              status_ == RangeStatus::inrange_signflip;
0269     }
0270 
0271     void print(std::ostream& os) const {
0272       if (outsideExpectedRange())
0273         os << "!! ";
0274       os << "(" << rangeMin_ << "," << rangeMax_ << ") ";
0275 
0276       os << diff_ << " " << pcvalue_;
0277       if (status_ == RangeStatus::underflow_OK || status_ == RangeStatus::underflow_notOK) {
0278         os << " (underflow) ";
0279         if (status_ == RangeStatus::underflow_notOK)
0280           underOverflow_.printNonOkUnderflow(os);
0281       } else if (status_ == RangeStatus::overflow_OK || status_ == RangeStatus::overflow_notOK) {
0282         os << " (overflow) ";
0283         if (status_ == RangeStatus::overflow_notOK)
0284           underOverflow_.printNonOkOverflow(os);
0285       } else if (status_ == RangeStatus::denormal)
0286         os << " (denormal)";
0287       os << " " << trackvalue_;
0288     }
0289 
0290   private:
0291     const double diff_;
0292     const double pcvalue_;
0293     const double trackvalue_;
0294     const double rangeMin_;
0295     const double rangeMax_;
0296     const RangeStatus status_;
0297     const typename T::UnderOverflow underOverflow_;
0298   };
0299 
0300   struct Range {
0301     Range(double mi, double ma) : min(mi), max(ma) {}
0302     const double min, max;
0303   };
0304   struct RangeAbs {
0305     RangeAbs(double val) : min(-val), max(val) {}
0306     const double min, max;
0307   };
0308 
0309   template <typename T>
0310   class PackedValueCheck {
0311   public:
0312     template <typename R, typename... Args>
0313     PackedValueCheck(const R& range, Args&&... args)
0314         : helper_(std::forward<Args>(args)...), rangeMin_(range.min), rangeMax_(range.max) {}
0315 
0316     void book(DQMStore::IBooker& iBooker,
0317               const std::string& name,
0318               const std::string& title,
0319               int nbins,
0320               double min,
0321               double max,
0322               int flow_nbins,
0323               double flow_min,
0324               double flow_max) {
0325       hInrange = iBooker.book1D(name, title, nbins, min, max);
0326       hUnderOverflowSign = iBooker.book1D(name + "UnderOverFlowSign",
0327                                           title + " with over- and underflow, and sign flip",
0328                                           flow_nbins,
0329                                           flow_min,
0330                                           flow_max);
0331       hStatus = iBooker.book1D(name + "Status", title + " status", 7, -0.5, 6.5);
0332       hStatus->setBinLabel(1, "In range");
0333       hStatus->setBinLabel(2, "In range, sign flip");
0334       hStatus->setBinLabel(3, "Denormal");
0335       hStatus->setBinLabel(4, "Underflow, PC is " + T::minName());
0336       hStatus->setBinLabel(5, "Underflow, PC is not " + T::minName());
0337       hStatus->setBinLabel(6, "Overflow, PC is " + T::maxName());
0338       hStatus->setBinLabel(7, "Overflow, PC is not " + T::maxName());
0339     }
0340 
0341     PackedValueCheckResult<T> fill(double pcvalue,
0342                                    double trackvalue,
0343                                    std::function<double(double)> modifyPack = std::function<double(double)>(),
0344                                    std::function<double(double)> modifyUnpack = std::function<double(double)>()) {
0345       const auto diff = diffRelative(pcvalue, trackvalue);
0346 
0347       const auto tmpSigned = modifyPack ? modifyPack(trackvalue) : trackvalue;
0348       const auto tmp = std::abs(tmpSigned);
0349       const auto underOverflow = helper_.underOverflowHelper(tmpSigned, modifyUnpack);
0350       RangeStatus status;
0351       if (tmp > helper_.largestValue()) {
0352         fillNoFlow(hUnderOverflowSign, diff);
0353         if (underOverflow.compatibleWithOverflow(std::abs(pcvalue))) {
0354           status = RangeStatus::overflow_OK;
0355         } else {
0356           status = RangeStatus::overflow_notOK;
0357         }
0358       } else if (helper_.lessThanSmallestValue(tmpSigned)) {
0359         fillNoFlow(hUnderOverflowSign, diff);
0360         if (underOverflow.compatibleWithUnderflow(std::abs(pcvalue))) {
0361           status = RangeStatus::underflow_OK;
0362         } else {
0363           status = RangeStatus::underflow_notOK;
0364         }
0365       } else {
0366         //Check if both numbers have same sign or both are zero.
0367         if (pcvalue * trackvalue > 0 || pcvalue == trackvalue) {
0368           if (T::wouldBeDenorm(tmp)) {
0369             status = RangeStatus::denormal;
0370           } else {
0371             status = RangeStatus::inrange;
0372           }
0373           fillNoFlow(hInrange, diff);
0374         } else {
0375           fillNoFlow(hUnderOverflowSign, diff);
0376           status = RangeStatus::inrange_signflip;
0377         }
0378       }
0379       hStatus->Fill(static_cast<int>(status));
0380 
0381       return PackedValueCheckResult<T>(status, diff, pcvalue, trackvalue, rangeMin_, rangeMax_, underOverflow);
0382     }
0383 
0384   private:
0385     const T helper_;
0386     const double rangeMin_;
0387     const double rangeMax_;
0388 
0389     MonitorElement* hInrange;
0390     MonitorElement* hUnderOverflowSign;
0391     MonitorElement* hStatus;
0392   };
0393   template <typename T>
0394   std::ostream& operator<<(std::ostream& os, const PackedValueCheckResult<T>& res) {
0395     res.print(os);
0396     return os;
0397   }
0398 }  // namespace
0399 
0400 class PackedCandidateTrackValidator : public DQMEDAnalyzer {
0401 public:
0402   PackedCandidateTrackValidator(const edm::ParameterSet& pset);
0403   ~PackedCandidateTrackValidator() override;
0404 
0405   void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0406   void analyze(const edm::Event&, const edm::EventSetup&) override;
0407 
0408   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0409 
0410 private:
0411   edm::EDGetTokenT<edm::View<reco::Track>> tracksToken_;
0412   edm::EDGetTokenT<reco::VertexCollection> verticesToken_;
0413   edm::EDGetTokenT<edm::Association<pat::PackedCandidateCollection>> trackToPackedCandidateToken_;
0414 
0415   std::string rootFolder_;
0416   bool debug_;
0417 
0418   enum {
0419     sf_AllTracks = 0,
0420     sf_AssociatedToPC = 1,
0421     sf_PCIsCharged = 2,
0422     sf_PCHasTrack = 3,
0423     sf_PCIsNotElectron = 4,
0424     sf_PCHasHits = 5,
0425     sf_PCNdofNot0 = 6,
0426     sf_NoMissingInnerHits = 7
0427   };
0428   MonitorElement* h_selectionFlow;
0429 
0430   MonitorElement* h_diffVx;
0431   MonitorElement* h_diffVy;
0432   MonitorElement* h_diffVz;
0433 
0434   MonitorElement* h_diffNormalizedChi2;
0435   MonitorElement* h_diffNdof;
0436 
0437   MonitorElement* h_diffCharge;
0438   MonitorElement* h_diffIsHighPurity;
0439 
0440   MonitorElement* h_diffPt;
0441   MonitorElement* h_diffEta;
0442   MonitorElement* h_diffPhi;
0443   PackedValueCheck<Float16Helper> h_diffDxyAssocPV;
0444   PackedValueCheck<Float16Helper> h_diffDzAssocPV;
0445   MonitorElement* h_diffDxyPV;
0446   MonitorElement* h_diffDzPV;
0447 
0448   MonitorElement* h_diffTrackDxyAssocPV;
0449   MonitorElement* h_diffTrackDzAssocPV;
0450 
0451   PackedValueCheck<LogIntHelper> h_diffCovQoverpQoverp;
0452   PackedValueCheck<LogIntHelper> h_diffCovLambdaLambda;
0453   PackedValueCheck<LogIntHelper> h_diffCovLambdaDsz;
0454   PackedValueCheck<LogIntHelper> h_diffCovPhiPhi;
0455   PackedValueCheck<LogIntHelper> h_diffCovPhiDxy;
0456   PackedValueCheck<Float16Helper> h_diffCovDxyDxy;
0457   PackedValueCheck<Float16Helper> h_diffCovDxyDsz;
0458   PackedValueCheck<Float16Helper> h_diffCovDszDsz;
0459 
0460   MonitorElement* h_diffDxyError;
0461   MonitorElement* h_diffDszError;
0462   MonitorElement* h_diffDzError;
0463 
0464   MonitorElement* h_diffTrackDxyError;
0465   MonitorElement* h_diffTrackDzError;
0466 
0467   MonitorElement* h_diffPtError;
0468   MonitorElement* h_diffEtaError;
0469 
0470   MonitorElement* h_diffNumberOfPixelLayers;
0471   MonitorElement* h_diffNumberOfStripLayers;
0472   MonitorElement* h_diffNumberOfPixelHits;
0473   MonitorElement* h_diffNumberOfHits;
0474   MonitorElement* h_diffLostInnerHits;
0475 
0476   MonitorElement* h_diffHitPatternPixelLayersWithMeasurement;
0477   MonitorElement* h_diffHitPatternTrackerLayersWithMeasurement;
0478   MonitorElement* h_diffHitPatternStripLayersWithMeasurement;
0479   MonitorElement* h_diffHitPatternNumberOfValidPixelHits;
0480   MonitorElement* h_diffHitPatternNumberOfValidHits;
0481   MonitorElement* h_diffHitPatternNumberOfLostInnerHits;
0482   MonitorElement* h_diffHitPatternHasValidHitInFirstPixelBarrel;
0483 
0484   MonitorElement* h_numberPixelLayersOverMax;
0485   MonitorElement* h_numberStripLayersOverMax;
0486   MonitorElement* h_numberLayersOverMax;
0487   MonitorElement* h_numberPixelHitsOverMax;
0488   MonitorElement* h_numberStripHitsOverMax;
0489   MonitorElement* h_numberHitsOverMax;
0490 };
0491 
0492 PackedCandidateTrackValidator::PackedCandidateTrackValidator(const edm::ParameterSet& iConfig)
0493     : tracksToken_(consumes<edm::View<reco::Track>>(iConfig.getUntrackedParameter<edm::InputTag>("tracks"))),
0494       verticesToken_(consumes<reco::VertexCollection>(iConfig.getUntrackedParameter<edm::InputTag>("vertices"))),
0495       trackToPackedCandidateToken_(consumes<edm::Association<pat::PackedCandidateCollection>>(
0496           iConfig.getUntrackedParameter<edm::InputTag>("trackToPackedCandidateAssociation"))),
0497       rootFolder_(iConfig.getUntrackedParameter<std::string>("rootFolder")),
0498       debug_(iConfig.getUntrackedParameter<bool>("debug")),
0499       h_diffDxyAssocPV(RangeAbs(0.001)),
0500       h_diffDzAssocPV(RangeAbs(0.001)),
0501       h_diffCovQoverpQoverp(Range(-1e-6, 0.13), -15, 0),  // despite of ceil in pack, there is rounding in double->float
0502       h_diffCovLambdaLambda(
0503           Range(-1e-6, 0.13), -20, -5),  // despite of ceil in pack, there is rounding in double->float
0504       h_diffCovLambdaDsz(RangeAbs(0.13), -17, -4),
0505       h_diffCovPhiPhi(RangeAbs(0.13), -15, 0),
0506       h_diffCovPhiDxy(RangeAbs(0.13), -17, -4),
0507       h_diffCovDxyDxy(RangeAbs(0.001)),
0508       h_diffCovDxyDsz(RangeAbs(0.001)),
0509       h_diffCovDszDsz(RangeAbs(0.001)) {}
0510 
0511 PackedCandidateTrackValidator::~PackedCandidateTrackValidator() {}
0512 
0513 void PackedCandidateTrackValidator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0514   edm::ParameterSetDescription desc;
0515 
0516   desc.addUntracked<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
0517   desc.addUntracked<edm::InputTag>("vertices", edm::InputTag("offlinePrimaryVertices"));
0518   desc.addUntracked<edm::InputTag>("trackToPackedCandidateAssociation", edm::InputTag("packedPFCandidates"));
0519   desc.addUntracked<std::string>("rootFolder", "Tracking/PackedCandidate");
0520   desc.addUntracked<bool>("debug", false);
0521 
0522   descriptions.add("packedCandidateTrackValidator", desc);
0523 }
0524 
0525 void PackedCandidateTrackValidator::bookHistograms(DQMStore::IBooker& iBooker,
0526                                                    edm::Run const&,
0527                                                    edm::EventSetup const&) {
0528   iBooker.setCurrentFolder(rootFolder_);
0529 
0530   h_selectionFlow = iBooker.book1D("selectionFlow", "Track selection flow", 8, -0.5, 7.5);
0531   h_selectionFlow->setBinLabel(1, "All tracks");
0532   h_selectionFlow->setBinLabel(2, "Associated to PackedCandidate");
0533   h_selectionFlow->setBinLabel(3, "PC is charged"), h_selectionFlow->setBinLabel(4, "PC has track");
0534   h_selectionFlow->setBinLabel(5, "PC is not electron");
0535   h_selectionFlow->setBinLabel(6, "PC has hits");
0536   h_selectionFlow->setBinLabel(7, "PC ndof != 0");
0537   h_selectionFlow->setBinLabel(8, "Track: no missing inner hits");
0538 
0539   constexpr int diffBins = 50;
0540 
0541   h_diffVx =
0542       iBooker.book1D("diffVx", "PackedCandidate::bestTrack() - reco::Track in vx()", diffBins, -0.2, 0.2);  // not equal
0543   h_diffVy =
0544       iBooker.book1D("diffVy", "PackedCandidate::bestTrack() - reco::Track in vy()", diffBins, -0.2, 0.2);  // not equal
0545   h_diffVz =
0546       iBooker.book1D("diffVz", "PackedCandidate::bestTrack() - reco::Track in vz()", diffBins, -0.4, 0.4);  // not equal
0547 
0548   h_diffNormalizedChi2 = iBooker.book1D("diffNormalizedChi2",
0549                                         "PackedCandidate::bestTrack() - reco::Track in normalizedChi2()",
0550                                         30,
0551                                         -1.5,
0552                                         1.5);  // expected difference in -1...0
0553   h_diffNdof = iBooker.book1D(
0554       "diffNdof", "PackedCandidate::bestTrack() - reco::Track in ndof()", 33, -30.5, 2.5);  // to monitor the difference
0555 
0556   h_diffCharge = iBooker.book1D(
0557       "diffCharge", "PackedCandidate::bestTrack() - reco::Track in charge()", 5, -2.5, 2.5);  // expect equality
0558   h_diffIsHighPurity = iBooker.book1D("diffIsHighPurity",
0559                                       "PackedCandidate::bestTrack() - reco::Track in quality(highPurity)",
0560                                       3,
0561                                       -1.5,
0562                                       1.5);  // expect equality
0563 
0564   h_diffPt = iBooker.book1D("diffPt",
0565                             "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in pt()",
0566                             diffBins,
0567                             -1.1,
0568                             1.1);  // not equal, keep
0569   h_diffEta = iBooker.book1D(
0570       "diffEta", "PackedCandidate::bestTrack() - reco::Track in eta()", diffBins, -0.001, 0.001);  // not equal, keep
0571   h_diffPhi = iBooker.book1D("diffPhi",
0572                              "PackedCandidate::bestTrack() - reco::Track in phi()",
0573                              diffBins,
0574                              -0.0005,
0575                              0.0005);  // expect equality within precision
0576 
0577   h_diffDxyAssocPV.book(iBooker,
0578                         "diffDxyAssocPV",
0579                         "(PackedCandidate::dxy() - reco::Track::dxy(assocPV))/reco::Track",
0580                         40,
0581                         -0.001,
0582                         0.001,  // expect equality within precision
0583                         50,
0584                         -0.5,
0585                         0.5);
0586   h_diffDzAssocPV.book(iBooker,
0587                        "diffDzAssocPV",
0588                        "(PackedCandidate::dzAssociatedPV() - reco::Track::dz(assocPV))/reco::Track",
0589                        40,
0590                        -0.001,
0591                        0.001,  // expect equality within precision
0592                        50,
0593                        -0.5,
0594                        0.5);
0595   h_diffDxyPV = iBooker.book1D("diffDxyPV",
0596                                "(PackedCandidate::dxy(PV) - reco::Track::dxy(PV))/reco::Track",
0597                                diffBins,
0598                                -0.01,
0599                                0.01);  // expect equality within precision (worse than assocPV)
0600   h_diffDzPV = iBooker.book1D("diffDzPV",
0601                               "(PackedCandidate::dz(PV) - reco::Track::dz(PV))/reco::Track",
0602                               diffBins,
0603                               -0.01,
0604                               0.01);  // expect equality wihtin precision (worse than assocPV)
0605   h_diffTrackDxyAssocPV =
0606       iBooker.book1D("diffTrackDxyAssocPV",
0607                      "(PackedCandidate::bestTrack()::dxy(assocPV)) - reco::Track::dxy(assocPV))/reco::Track",
0608                      diffBins,
0609                      -0.01,
0610                      0.01);  // not equal
0611   h_diffTrackDzAssocPV =
0612       iBooker.book1D("diffTrackDzAssocPV",
0613                      "(PackedCandidate::bestTrack()::dz(assocPV)) - reco::Track::dz(assocPV))/reco::Track",
0614                      diffBins,
0615                      -0.01,
0616                      0.01);  // not equal
0617 
0618   h_diffCovQoverpQoverp.book(iBooker,
0619                              "diffCovQoverpQoverp",
0620                              "(PackedCandidate::bestTrack() - reco::Track)/reco::track in cov(qoverp, qoverp)",
0621                              40,
0622                              -0.05,
0623                              0.15,  // expect equality within precision (worst precision is exp(1/128*15) =~ 12 %
0624                              50,
0625                              -0.5,
0626                              0.5);
0627   h_diffCovLambdaLambda.book(
0628       iBooker,
0629       "diffCovLambdaLambda",
0630       "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in cov(lambda, lambda)",
0631       40,
0632       -0.05,
0633       0.15,  // expect equality within precision worst precision is exp(1/128*(20-5)) =~ 12 % (multiplied by pt^2 in packing & unpacking)
0634       50,
0635       -0.5,
0636       0.5);
0637   h_diffCovLambdaDsz.book(iBooker,
0638                           "diffCovLambdaDsz",
0639                           "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in cov(lambda, dsz)",
0640                           60,
0641                           -0.15,
0642                           0.15,  // expect equality within precision, worst precision is exp(1/128*(17-4) =~ 11 %
0643                           50,
0644                           -1,
0645                           1);
0646   h_diffCovPhiPhi.book(
0647       iBooker,
0648       "diffCovPhiPhi",
0649       "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in cov(phi, phi)",
0650       40,
0651       -0.05,
0652       0.15,  // expect equality within precision worst precision is exp(1/128*(20-5)) =~ 12 % (multiplied by pt^2 in packing & unpacking)
0653       50,
0654       -0.5,
0655       0.5);
0656   h_diffCovPhiDxy.book(iBooker,
0657                        "diffCovPhiDxy",
0658                        "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in cov(phi, dxy)",
0659                        60,
0660                        -0.15,
0661                        0.15,  // expect equality within precision, wors precision is exp(1/128)*(17-4) =~ 11 %
0662                        50,
0663                        -1,
0664                        1);
0665   h_diffCovDxyDxy.book(iBooker,
0666                        "diffCovDxyDxy",
0667                        "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in cov(dxy, dxy)",
0668                        40,
0669                        -0.001,
0670                        0.001,
0671                        50,
0672                        -0.1,
0673                        0.1);
0674   h_diffCovDxyDsz.book(iBooker,
0675                        "diffCovDxyDsz",
0676                        "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in cov(dxy, dsz)",
0677                        40,
0678                        -0.001,
0679                        0.001,  // expect equality within precision
0680                        50,
0681                        -0.5,
0682                        0.5);
0683   h_diffCovDszDsz.book(iBooker,
0684                        "diffCovDszDsz",
0685                        "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in cov(dsz, dsz)",
0686                        40,
0687                        -0.001,
0688                        0.001,  // expect equality within precision
0689                        50,
0690                        -0.1,
0691                        0.1);
0692 
0693   h_diffDxyError = iBooker.book1D("diffDxyError",
0694                                   "(PackedCandidate::dxyError() - reco::Track::dxyError())/reco::Track",
0695                                   40,
0696                                   -0.001,
0697                                   0.001);  // expect equality within precision
0698   h_diffDszError = iBooker.book1D("diffDszError",
0699                                   "(PackedCandidate::dzError() - reco::Track::dszError())/reco::Track",
0700                                   40,
0701                                   -0.001,
0702                                   0.001);  // ideally, not equal, but for now they are
0703   h_diffDzError = iBooker.book1D("diffDzError",
0704                                  "(PackedCandidate::dzError() - reco::Track::dzError())/reco::Track",
0705                                  40,
0706                                  -0.001,
0707                                  0.001);  // expect equality within precision (not currently the case)
0708 
0709   h_diffTrackDxyError = iBooker.book1D("diffTrackDxyError",
0710                                        "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in dxyError()",
0711                                        40,
0712                                        -0.001,
0713                                        0.001);  // expect equality within precision
0714   h_diffTrackDzError = iBooker.book1D("diffTrackDzError",
0715                                       "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in dzError()",
0716                                       40,
0717                                       -0.05,
0718                                       0.05);  // not equal
0719 
0720   h_diffPtError = iBooker.book1D("diffPtError",
0721                                  "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in ptError()",
0722                                  diffBins,
0723                                  -1.1,
0724                                  1.1);  // not equal
0725   h_diffEtaError = iBooker.book1D("diffEtaError",
0726                                   "(PackedCandidate::bestTrack() - reco::Track)/reco::Track in etaError()",
0727                                   60,
0728                                   -0.15,
0729                                   0.15);  // not equal
0730 
0731   h_diffNumberOfPixelLayers = iBooker.book1D(
0732       "diffNumberOfPixelLayers",
0733       "PackedCandidate::pixelLayersWithMeasurement() - reco::Track::hitPattern::pixelLayersWithMeasurement()",
0734       5,
0735       -2.5,
0736       2.5);  // expect equality
0737   h_diffNumberOfStripLayers = iBooker.book1D(
0738       "diffNumberOfStripLayers",
0739       "PackedCandidate::stripLayersWithMeasurement() - reco::Track::hitPattern::stripLayersWithMeasurement()",
0740       5,
0741       -2.5,
0742       2.5);  // expect equality
0743   h_diffNumberOfPixelHits =
0744       iBooker.book1D("diffNumberOfPixelHits",
0745                      "PackedCandidate::numberOfPixelHits() - reco::Track::hitPattern::numberOfValidPixelHits()",
0746                      5,
0747                      -2.5,
0748                      2.5);  // expect equality
0749   h_diffNumberOfHits = iBooker.book1D("diffNumberOfHits",
0750                                       "PackedCandidate::numberHits() - reco::Track::hitPattern::numberOfValidHits()",
0751                                       5,
0752                                       -2.5,
0753                                       2.5);  // expect equality
0754   h_diffLostInnerHits =
0755       iBooker.book1D("diffLostInnerHits",
0756                      "PackedCandidate::lostInnerHits() - reco::Track::hitPattern::numberOfLostHits(MISSING_INNER_HITS)",
0757                      5,
0758                      -2.5,
0759                      2.5);  // expect equality
0760 
0761   h_diffHitPatternPixelLayersWithMeasurement =
0762       iBooker.book1D("diffHitPatternPixelLayersWithMeasurement",
0763                      "PackedCandidate::bestTrack() - reco::Track in hitPattern::pixelLayersWithMeasurement()",
0764                      13,
0765                      -10.5,
0766                      2.5);  // not equal
0767   h_diffHitPatternStripLayersWithMeasurement =
0768       iBooker.book1D("diffHitPatternStripLayersWithMeasurement",
0769                      "PackedCandidate::bestTrack() - reco::Track in hitPattern::stripLayersWithMeasurement()",
0770                      13,
0771                      -10.5,
0772                      2.5);  // not equal
0773   h_diffHitPatternTrackerLayersWithMeasurement =
0774       iBooker.book1D("diffHitPatternTrackerLayersWithMeasurement",
0775                      "PackedCandidate::bestTrack() - reco::Track in hitPattern::trackerLayersWithMeasurement()",
0776                      13,
0777                      -10.5,
0778                      2.5);  // not equal
0779   h_diffHitPatternNumberOfValidPixelHits =
0780       iBooker.book1D("diffHitPatternNumberOfValidPixelHits",
0781                      "PackedCandidate::bestTrack() - reco::Track in hitPattern::numberOfValidPixelHits()",
0782                      13,
0783                      -10.5,
0784                      2.5);  // not equal
0785   h_diffHitPatternNumberOfValidHits =
0786       iBooker.book1D("diffHitPatternNumberOfValidHits",
0787                      "PackedCandidate::bestTrack() - reco::Track in hitPattern::numberOfValidHits()",
0788                      13,
0789                      -10.5,
0790                      2.5);  // not equal
0791   h_diffHitPatternNumberOfLostInnerHits =
0792       iBooker.book1D("diffHitPatternNumberOfLostPixelHits",
0793                      "PackedCandidate::bestTrack() - reco::Track in hitPattern::numberOfLostHits(MISSING_INNER_HITS)",
0794                      13,
0795                      -10.5,
0796                      2.5);  // not equal
0797   h_diffHitPatternHasValidHitInFirstPixelBarrel =
0798       iBooker.book1D("diffHitPatternHasValidHitInFirstPixelBarrel",
0799                      "PackedCandidate::bestTrack() - reco::Track in hitPattern::hasValidHitInFirstPixelBarrel",
0800                      3,
0801                      -1.5,
0802                      1.5);  // expect equality
0803 
0804   h_numberPixelLayersOverMax = iBooker.book1D(
0805       "numberPixelLayersOverMax", "Number of pixel layers over the maximum of PackedCandidate", 10, 0, 10);
0806   h_numberStripLayersOverMax = iBooker.book1D(
0807       "numberStripLayersOverMax", "Number of strip layers over the maximum of PackedCandidate", 10, 0, 10);
0808   h_numberLayersOverMax =
0809       iBooker.book1D("numberLayersOverMax", "Number of layers over the maximum of PackedCandidate", 20, 0, 20);
0810   h_numberPixelHitsOverMax =
0811       iBooker.book1D("numberPixelHitsOverMax", "Number of pixel hits over the maximum of PackedCandidate", 10, 0, 10);
0812   h_numberStripHitsOverMax =
0813       iBooker.book1D("numberStripHitsOverMax", "Number of strip hits over the maximum of PackedCandidate", 10, 0, 10);
0814   h_numberHitsOverMax =
0815       iBooker.book1D("numberHitsOverMax", "Number of hits over the maximum of PackedCandidate", 20, 0, 20);
0816 }
0817 
0818 void PackedCandidateTrackValidator::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0819   edm::Handle<edm::View<reco::Track>> htracks;
0820   iEvent.getByToken(tracksToken_, htracks);
0821   const auto& tracks = *htracks;
0822 
0823   edm::Handle<reco::VertexCollection> hvertices;
0824   iEvent.getByToken(verticesToken_, hvertices);
0825   const auto& vertices = *hvertices;
0826 
0827   if (vertices.empty())
0828     return;
0829   const reco::Vertex& pv = vertices[0];
0830 
0831   edm::Handle<edm::Association<pat::PackedCandidateCollection>> hassoc;
0832   iEvent.getByToken(trackToPackedCandidateToken_, hassoc);
0833   const auto& trackToPackedCandidate = *hassoc;
0834 
0835   for (size_t i = 0; i < tracks.size(); ++i) {
0836     auto trackPtr = tracks.ptrAt(i);
0837     const reco::Track& track = *trackPtr;
0838     h_selectionFlow->Fill(sf_AllTracks);
0839 
0840     pat::PackedCandidateRef pcRef = trackToPackedCandidate[trackPtr];
0841     if (pcRef.isNull()) {
0842       continue;
0843     }
0844     h_selectionFlow->Fill(sf_AssociatedToPC);
0845 
0846     // Filter out neutral PackedCandidates, some of them may have track associated, and for those the charge comparison fails
0847     if (pcRef->charge() == 0) {
0848       continue;
0849     }
0850     h_selectionFlow->Fill(sf_PCIsCharged);
0851 
0852     const reco::Track* trackPcPtr = pcRef->bestTrack();
0853     if (!trackPcPtr) {
0854       continue;
0855     }
0856     h_selectionFlow->Fill(sf_PCHasTrack);
0857 
0858     // Filter out electrons to avoid comparisons to PackedCandidates with GsfTrack
0859     if (std::abs(pcRef->pdgId()) == 11) {
0860       continue;
0861     }
0862     h_selectionFlow->Fill(sf_PCIsNotElectron);
0863 
0864     // Filter out PackedCandidate-tracks with no hits, as they won't have their details filled
0865     const reco::Track& trackPc = *trackPcPtr;
0866     if (trackPc.hitPattern().numberOfValidHits() == 0) {
0867       continue;
0868     }
0869     h_selectionFlow->Fill(sf_PCHasHits);
0870 
0871     auto slimmedVertexRef = pcRef->vertexRef();
0872     const reco::Vertex& pcVertex = vertices[slimmedVertexRef.key()];
0873 
0874     fillNoFlow(h_diffVx, trackPc.vx() - track.vx());
0875     fillNoFlow(h_diffVy, trackPc.vy() - track.vy());
0876     fillNoFlow(h_diffVz, trackPc.vz() - track.vz());
0877 
0878     // PackedCandidate recalculates the ndof in unpacking as
0879     // (nhits+npixelhits-5), but some strip hits may have dimension 2.
0880     // If PackedCandidate has ndof=0, the resulting normalizedChi2
0881     // will be 0 too. Hence, the comparison makes sense only for those
0882     // PackedCandidates that have ndof != 0.
0883     double diffNormalizedChi2 = 0;
0884     if (trackPc.ndof() != 0) {
0885       h_selectionFlow->Fill(sf_PCNdofNot0);
0886       diffNormalizedChi2 = trackPc.normalizedChi2() - track.normalizedChi2();
0887       fillNoFlow(h_diffNormalizedChi2, diffNormalizedChi2);
0888     }
0889     fillNoFlow(h_diffNdof, trackPc.ndof() - track.ndof());
0890 
0891     auto diffCharge = trackPc.charge() - track.charge();
0892     fillNoFlow(h_diffCharge, diffCharge);
0893     int diffHP = static_cast<int>(trackPc.quality(reco::TrackBase::highPurity)) -
0894                  static_cast<int>(track.quality(reco::TrackBase::highPurity));
0895     fillNoFlow(h_diffIsHighPurity, diffHP);
0896 
0897     const auto diffPt = diffRelative(trackPc.pt(), track.pt());
0898     const auto diffPhi = reco::deltaPhi(trackPc.phi(), track.phi());
0899     fillNoFlow(h_diffPt, diffPt);
0900     fillNoFlow(h_diffEta, trackPc.eta() - track.eta());
0901     fillNoFlow(h_diffPhi, diffPhi);
0902 
0903     const auto diffDxyAssocPV =
0904         h_diffDxyAssocPV.fill(pcRef->dxy(), track.dxy(pcVertex.position()), [](double value) { return value * 100.; });
0905     const auto diffDzAssocPV = h_diffDzAssocPV.fill(
0906         pcRef->dzAssociatedPV(), track.dz(pcVertex.position()), [](double value) { return value * 100.; });
0907     const auto diffDxyPV = diffRelative(pcRef->dxy(pv.position()), track.dxy(pv.position()));
0908     const auto diffDzPV = diffRelative(pcRef->dz(pv.position()), track.dz(pv.position()));
0909     fillNoFlow(h_diffDxyPV, diffDxyPV);
0910     fillNoFlow(h_diffDzPV, diffDzPV);
0911     fillNoFlow(h_diffTrackDxyAssocPV, diffRelative(trackPc.dxy(pcVertex.position()), track.dxy(pcVertex.position())));
0912     fillNoFlow(h_diffTrackDzAssocPV, diffRelative(trackPc.dz(pcVertex.position()), track.dz(pcVertex.position())));
0913 
0914     auto fillCov1 = [&](auto& hlp, const int i, const int j) {
0915       return hlp.fill(trackPc.covariance(i, j), track.covariance(i, j));
0916     };
0917     auto fillCov2 = [&](auto& hlp, const int i, const int j, std::function<double(double)> modifyPack) {
0918       return hlp.fill(trackPc.covariance(i, j), track.covariance(i, j), modifyPack);
0919     };
0920     auto fillCov3 = [&](auto& hlp,
0921                         const int i,
0922                         const int j,
0923                         std::function<double(double)> modifyPack,
0924                         std::function<double(double)> modifyUnpack) {
0925       return hlp.fill(trackPc.covariance(i, j), track.covariance(i, j), modifyPack, modifyUnpack);
0926     };
0927 
0928     const auto pcPt = pcRef->pt();
0929     const auto diffCovQoverpQoverp = fillCov3(
0930         h_diffCovQoverpQoverp,
0931         reco::TrackBase::i_qoverp,
0932         reco::TrackBase::i_qoverp,
0933         [=](double val) { return val * pcPt * pcPt; },
0934         [=](double val) { return val / pcPt / pcPt; });
0935     const auto diffCovLambdaLambda =
0936         fillCov1(h_diffCovLambdaLambda, reco::TrackBase::i_lambda, reco::TrackBase::i_lambda);
0937     const auto diffCovLambdaDsz = fillCov1(h_diffCovLambdaDsz, reco::TrackBase::i_lambda, reco::TrackBase::i_dsz);
0938     const auto diffCovPhiPhi = fillCov3(
0939         h_diffCovPhiPhi,
0940         reco::TrackBase::i_phi,
0941         reco::TrackBase::i_phi,
0942         [=](double val) { return val * pcPt * pcPt; },
0943         [=](double val) { return val / pcPt / pcPt; });
0944     const auto diffCovPhiDxy = fillCov1(h_diffCovPhiDxy, reco::TrackBase::i_phi, reco::TrackBase::i_dxy);
0945     const auto diffCovDxyDxy = fillCov2(
0946         h_diffCovDxyDxy, reco::TrackBase::i_dxy, reco::TrackBase::i_dxy, [](double value) { return value * 10000.; });
0947     const auto diffCovDxyDsz = fillCov2(
0948         h_diffCovDxyDsz, reco::TrackBase::i_dxy, reco::TrackBase::i_dsz, [](double value) { return value * 10000.; });
0949     const auto diffCovDszDsz = fillCov2(
0950         h_diffCovDszDsz, reco::TrackBase::i_dsz, reco::TrackBase::i_dsz, [](double value) { return value * 10000.; });
0951 
0952     if (isInRange(diffCovDszDsz.status())) {
0953       fillNoFlow(h_diffDszError, diffRelative(pcRef->dzError(), track.dszError()));
0954       fillNoFlow(h_diffDzError, diffRelative(pcRef->dzError(), track.dzError()));
0955       fillNoFlow(h_diffTrackDzError, diffRelative(trackPc.dzError(), track.dzError()));
0956     }
0957     if (isInRange(diffCovDxyDxy.status())) {
0958       fillNoFlow(h_diffDxyError, diffRelative(pcRef->dxyError(), track.dxyError()));
0959       fillNoFlow(h_diffTrackDxyError, diffRelative(trackPc.dxyError(), track.dxyError()));
0960     }
0961     fillNoFlow(h_diffPtError, diffRelative(trackPc.ptError(), track.ptError()));
0962     fillNoFlow(h_diffEtaError, diffRelative(trackPc.etaError(), track.etaError()));
0963 
0964     // For the non-HitPattern ones, take into account the PackedCandidate packing precision
0965     const auto trackNumberOfHits = track.hitPattern().numberOfValidHits();
0966     const auto trackNumberOfPixelHits = track.hitPattern().numberOfValidPixelHits();
0967     const auto trackNumberOfStripHits = track.hitPattern().numberOfValidStripHits();
0968     const auto pcNumberOfHits = pcRef->numberOfHits();
0969     const auto pcNumberOfPixelHits = pcRef->numberOfPixelHits();
0970     const auto pcNumberOfStripHits = pcNumberOfHits - pcNumberOfPixelHits;
0971     const auto trackNumberOfLayers = track.hitPattern().trackerLayersWithMeasurement();
0972     const auto trackNumberOfPixelLayers = track.hitPattern().pixelLayersWithMeasurement();
0973     const auto trackNumberOfStripLayers = track.hitPattern().stripLayersWithMeasurement();
0974     const auto pcNumberOfLayers = pcRef->trackerLayersWithMeasurement();
0975     const auto pcNumberOfPixelLayers = pcRef->pixelLayersWithMeasurement();
0976     const auto pcNumberOfStripLayers = pcRef->stripLayersWithMeasurement();
0977 
0978     // layer number overflow (should be zero)
0979     const int pixelLayerOverflow = trackNumberOfPixelLayers > pat::PackedCandidate::trackPixelHitsMask
0980                                        ? trackNumberOfPixelLayers - pat::PackedCandidate::trackPixelHitsMask
0981                                        : 0;
0982     const int stripLayerOverflow = trackNumberOfStripLayers > pat::PackedCandidate::trackStripHitsMask
0983                                        ? trackNumberOfStripLayers - pat::PackedCandidate::trackStripHitsMask
0984                                        : 0;
0985     const int layerOverflow =
0986         trackNumberOfLayers > (pat::PackedCandidate::trackPixelHitsMask + pat::PackedCandidate::trackStripHitsMask)
0987             ? trackNumberOfLayers -
0988                   (pat::PackedCandidate::trackPixelHitsMask + pat::PackedCandidate::trackStripHitsMask)
0989             : 0;
0990 
0991     // hit overflow (should also be zero)
0992     const int pixelOverflow =
0993         trackNumberOfPixelHits - pcNumberOfPixelLayers > pat::PackedCandidate::trackPixelHitsMask
0994             ? trackNumberOfPixelHits - pcNumberOfPixelLayers - pat::PackedCandidate::trackPixelHitsMask
0995             : 0;
0996     const int stripOverflow =
0997         trackNumberOfStripHits - pcNumberOfStripLayers > pat::PackedCandidate::trackStripHitsMask
0998             ? trackNumberOfStripHits - pcNumberOfStripLayers - pat::PackedCandidate::trackStripHitsMask
0999             : 0;
1000     const int hitsOverflow =
1001         trackNumberOfHits - pcNumberOfLayers >
1002                 (pat::PackedCandidate::trackPixelHitsMask + pat::PackedCandidate::trackStripHitsMask)
1003             ? trackNumberOfHits - pcNumberOfLayers -
1004                   (pat::PackedCandidate::trackPixelHitsMask + pat::PackedCandidate::trackStripHitsMask)
1005             : 0;
1006     // PackedCandidate counts overflow pixel hits as strip
1007     const int pixelInducedStripOverflow =
1008         (trackNumberOfStripHits + pixelOverflow - pcNumberOfStripLayers) > pat::PackedCandidate::trackStripHitsMask
1009             ? (trackNumberOfStripHits + pixelOverflow - stripOverflow - pcNumberOfStripLayers) -
1010                   pat::PackedCandidate::trackStripHitsMask
1011             : 0;
1012     h_numberPixelLayersOverMax->Fill(pixelLayerOverflow);
1013     h_numberStripLayersOverMax->Fill(stripLayerOverflow);
1014     h_numberLayersOverMax->Fill(layerOverflow);
1015     h_numberPixelHitsOverMax->Fill(pixelOverflow);
1016     h_numberStripHitsOverMax->Fill(stripOverflow);
1017     h_numberHitsOverMax->Fill(hitsOverflow);
1018 
1019     int diffNumberOfPixelHits = 0;
1020     int diffNumberOfHits = 0;
1021     int diffNumberOfPixelLayers = 0;
1022     int diffNumberOfStripLayers = 0;
1023     if (pixelLayerOverflow) {
1024       diffNumberOfPixelLayers = pcNumberOfPixelLayers - pat::PackedCandidate::trackPixelHitsMask;
1025     } else {
1026       diffNumberOfPixelLayers = pcNumberOfPixelLayers - trackNumberOfPixelLayers;
1027     }
1028     if (pixelOverflow) {
1029       diffNumberOfPixelHits = pcNumberOfPixelHits - pcNumberOfPixelLayers - pat::PackedCandidate::trackPixelHitsMask;
1030     } else {
1031       diffNumberOfPixelHits = pcNumberOfPixelHits - trackNumberOfPixelHits;
1032     }
1033     if (stripLayerOverflow) {
1034       diffNumberOfStripLayers = pcNumberOfStripLayers - pat::PackedCandidate::trackStripHitsMask;
1035     } else {
1036       diffNumberOfStripLayers = pcNumberOfStripLayers - trackNumberOfStripLayers;
1037     }
1038     if (stripOverflow || pixelInducedStripOverflow || pixelOverflow) {
1039       int diffNumberOfStripHits = 0;
1040       if (stripOverflow || pixelInducedStripOverflow) {
1041         diffNumberOfStripHits = pcNumberOfStripHits - pat::PackedCandidate::trackStripHitsMask;
1042       } else if (pixelOverflow) {
1043         diffNumberOfStripHits = (pcNumberOfStripHits - pixelOverflow) - trackNumberOfStripHits;
1044       }
1045 
1046       diffNumberOfHits = diffNumberOfPixelHits + diffNumberOfStripHits;
1047     } else {
1048       diffNumberOfHits = pcNumberOfHits - trackNumberOfHits;
1049     }
1050 
1051     fillNoFlow(h_diffNumberOfPixelHits, diffNumberOfPixelHits);
1052     fillNoFlow(h_diffNumberOfHits, diffNumberOfHits);
1053     fillNoFlow(h_diffNumberOfPixelLayers, diffNumberOfPixelLayers);
1054     fillNoFlow(h_diffNumberOfStripLayers, diffNumberOfStripLayers);
1055 
1056     int diffLostInnerHits = 0;
1057     const auto trackLostInnerHits = track.hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS);
1058     switch (pcRef->lostInnerHits()) {
1059       case pat::PackedCandidate::validHitInFirstPixelBarrelLayer:
1060       case pat::PackedCandidate::noLostInnerHits:
1061         diffLostInnerHits = -trackLostInnerHits;
1062         break;
1063       case pat::PackedCandidate::oneLostInnerHit:
1064         diffLostInnerHits = 1 - trackLostInnerHits;
1065         break;
1066       case pat::PackedCandidate::moreLostInnerHits:
1067         diffLostInnerHits = trackLostInnerHits >= 2 ? 0 : 2 - trackLostInnerHits;
1068         break;
1069     }
1070     fillNoFlow(h_diffLostInnerHits, diffLostInnerHits);
1071 
1072     // For HitPattern ones, calculate the full diff (i.e. some differences are expected)
1073     auto diffHitPatternPixelLayersWithMeasurement =
1074         trackPc.hitPattern().pixelLayersWithMeasurement() - trackNumberOfPixelLayers;
1075     fillNoFlow(h_diffHitPatternPixelLayersWithMeasurement, diffHitPatternPixelLayersWithMeasurement);
1076     auto diffHitPatternStripLayersWithMeasurement =
1077         trackPc.hitPattern().stripLayersWithMeasurement() - trackNumberOfStripLayers;
1078     fillNoFlow(h_diffHitPatternStripLayersWithMeasurement, diffHitPatternStripLayersWithMeasurement);
1079     auto diffHitPatternTrackerLayersWithMeasurement =
1080         trackPc.hitPattern().trackerLayersWithMeasurement() - trackNumberOfLayers;
1081     fillNoFlow(h_diffHitPatternTrackerLayersWithMeasurement, diffHitPatternTrackerLayersWithMeasurement);
1082     auto diffHitPatternNumberOfValidPixelHits = trackPc.hitPattern().numberOfValidPixelHits() - trackNumberOfPixelHits;
1083     fillNoFlow(h_diffHitPatternNumberOfValidPixelHits, diffHitPatternNumberOfValidPixelHits);
1084     auto diffHitPatternNumberOfValidHits = trackPc.hitPattern().numberOfValidHits() - trackNumberOfHits;
1085     fillNoFlow(h_diffHitPatternNumberOfValidHits, diffHitPatternNumberOfValidHits);
1086     fillNoFlow(h_diffHitPatternNumberOfLostInnerHits,
1087                trackPc.hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS) -
1088                    track.hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS));
1089 
1090     // hasValidHitInFirstPixelBarrel is set only if numberOfLostHits(MISSING_INNER_HITS) == 0
1091     int diffHitPatternHasValidHitInFirstPixelBarrel = 0;
1092     if (track.hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS) == 0) {
1093       h_selectionFlow->Fill(sf_NoMissingInnerHits);
1094       diffHitPatternHasValidHitInFirstPixelBarrel =
1095           static_cast<int>(
1096               trackPc.hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector::PixelBarrel, 1)) -
1097           static_cast<int>(track.hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector::PixelBarrel, 1));
1098       fillNoFlow(h_diffHitPatternHasValidHitInFirstPixelBarrel, diffHitPatternHasValidHitInFirstPixelBarrel);
1099     }
1100 
1101     // Print warning if there are differences outside the expected range
1102     if (debug_ &&
1103         (diffNormalizedChi2 < -1 || diffNormalizedChi2 > 0 || diffCharge != 0 || diffHP != 0 ||
1104          std::abs(diffPhi) > 5e-4 || diffDxyAssocPV.outsideExpectedRange() || diffDzAssocPV.outsideExpectedRange() ||
1105          std::abs(diffDxyPV) > 0.05 || std::abs(diffDzPV) > 0.05 || diffCovQoverpQoverp.outsideExpectedRange() ||
1106          diffCovLambdaLambda.outsideExpectedRange() || diffCovLambdaDsz.outsideExpectedRange() ||
1107          diffCovPhiPhi.outsideExpectedRange() || diffCovPhiDxy.outsideExpectedRange() ||
1108          diffCovDxyDxy.outsideExpectedRange() || diffCovDxyDsz.outsideExpectedRange() ||
1109          diffCovDszDsz.outsideExpectedRange() || diffNumberOfPixelHits != 0 || diffNumberOfHits != 0 ||
1110          diffLostInnerHits != 0 || diffHitPatternHasValidHitInFirstPixelBarrel != 0)) {
1111       edm::LogInfo("PackedCandidateTrackValidator")
1112           << "Track " << i << " pt " << track.pt() << " eta " << track.eta() << " phi " << track.phi() << " chi2 "
1113           << track.chi2() << " ndof " << track.ndof() << "\n"
1114           << "  ptError " << track.ptError() << " etaError " << track.etaError() << " phiError " << track.phiError()
1115           << " dxyError " << track.dxyError() << " dzError " << track.dzError() << "\n"
1116           << "  refpoint " << track.referencePoint() << " momentum " << track.momentum() << "\n"
1117           << "  dxy " << track.dxy() << " dz " << track.dz() << "\n"
1118           << "  " << TrackAlgoPrinter(track) << " lost inner hits " << trackLostInnerHits << " lost outer hits "
1119           << track.hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS) << " hitpattern "
1120           << HitPatternPrinter(track) << " \n"
1121           << " PC " << pcRef.id() << ":" << pcRef.key() << " track pt " << trackPc.pt() << " eta " << trackPc.eta()
1122           << " phi " << trackPc.phi() << " (PC " << pcRef->phi() << ") chi2 " << trackPc.chi2() << " ndof "
1123           << trackPc.ndof() << " pdgId " << pcRef->pdgId() << " mass " << pcRef->mass() << "\n"
1124           << "  ptError " << trackPc.ptError() << " etaError " << trackPc.etaError() << " phiError "
1125           << trackPc.phiError() << "\n"
1126           << "  pc.vertex " << pcRef->vertex() << " momentum " << pcRef->momentum() << " track " << trackPc.momentum()
1127           << "\n"
1128           << "  dxy " << trackPc.dxy() << " dz " << trackPc.dz() << " pc.dz " << pcRef->dz() << " dxyError "
1129           << trackPc.dxyError() << " dzError " << trackPc.dzError() << "\n"
1130           << "  dxy(PV) " << trackPc.dxy(pv.position()) << " dz(PV) " << trackPc.dz(pv.position()) << " dxy(assocPV) "
1131           << trackPc.dxy(pcVertex.position()) << " dz(assocPV) " << trackPc.dz(pcVertex.position()) << "\n"
1132           << " (diff PackedCandidate track)"
1133           << " highPurity " << diffHP << " " << trackPc.quality(reco::TrackBase::highPurity) << " "
1134           << track.quality(reco::TrackBase::highPurity) << " charge " << diffCharge << " " << trackPc.charge() << " "
1135           << track.charge() << " normalizedChi2 " << diffNormalizedChi2 << " " << trackPc.normalizedChi2() << " "
1136           << track.normalizedChi2() << "\n "
1137           << " numberOfAllHits " << diffNumberOfHits << " " << pcNumberOfHits << " " << trackNumberOfHits
1138           << " numberOfPixelHits " << diffNumberOfPixelHits << " " << pcNumberOfPixelHits << " "
1139           << trackNumberOfPixelHits << " numberOfStripHits # " << pcNumberOfStripHits << " " << trackNumberOfStripHits
1140           << "\n "
1141           << " hitPattern.numberOfValidPixelHits " << diffHitPatternNumberOfValidPixelHits << " "
1142           << trackPc.hitPattern().numberOfValidPixelHits() << " " << track.hitPattern().numberOfValidPixelHits()
1143           << " hitPattern.numberOfValidHits " << diffHitPatternNumberOfValidHits << " "
1144           << trackPc.hitPattern().numberOfValidHits() << " " << track.hitPattern().numberOfValidHits()
1145           << " hitPattern.hasValidHitInFirstPixelBarrel " << diffHitPatternHasValidHitInFirstPixelBarrel << " "
1146           << trackPc.hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector::PixelBarrel, 1) << " "
1147           << track.hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector::PixelBarrel, 1) << "\n "
1148           << " lostInnerHits  " << diffLostInnerHits << " " << pcRef->lostInnerHits() << " #"
1149           << " phi (5e-4) " << diffPhi << " " << trackPc.phi() << " " << track.phi() << "\n "
1150           << " dxy(assocPV) " << diffDxyAssocPV << "\n "
1151           << " dz(assocPV) " << diffDzAssocPV << "\n "
1152           << " dxy(PV) (0.05) " << diffDxyPV << " " << pcRef->dxy(pv.position()) << " " << track.dxy(pv.position())
1153           << "\n "
1154           << " dz(PV) (0.05) " << diffDzPV << " " << pcRef->dz(pv.position()) << " " << track.dz(pv.position()) << "\n "
1155           << " cov(qoverp, qoverp)  " << diffCovQoverpQoverp << "\n "
1156           << " cov(lambda, lambda) " << diffCovLambdaLambda << "\n "
1157           << " cov(lambda, dsz) " << diffCovLambdaDsz << "\n "
1158           << " cov(phi, phi) " << diffCovPhiPhi << "\n "
1159           << " cov(phi, dxy) " << diffCovPhiDxy << "\n "
1160           << " cov(dxy, dxy) " << diffCovDxyDxy << "\n "
1161           << " cov(dxy, dsz) " << diffCovDxyDsz << "\n "
1162           << " cov(dsz, dsz) " << diffCovDszDsz;
1163     }
1164   }
1165 }
1166 
1167 #include "FWCore/Framework/interface/MakerMacros.h"
1168 DEFINE_FWK_MODULE(PackedCandidateTrackValidator);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!