Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2025-05-06 02:07:42

0001 #include <numeric>
0002 #include <string>
0003 #include <vector>
0004 #include <map>
0005 #include <algorithm>
0006 
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/ESHandle.h"
0009 #include "FWCore/Utilities/interface/isFinite.h"
0010 
0011 #include "DataFormats/Common/interface/ValidHandle.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 #include "FWCore/Framework/interface/MakerMacros.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017 #include "FWCore/Utilities/interface/isFinite.h"
0018 
0019 #include "DataFormats/Math/interface/LorentzVector.h"
0020 #include "DataFormats/Math/interface/Point3D.h"
0021 
0022 // reco track and vertex
0023 #include "DataFormats/TrackReco/interface/Track.h"
0024 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0025 #include "DataFormats/VertexReco/interface/Vertex.h"
0026 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0027 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0028 
0029 // TrackingParticle
0030 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0031 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
0032 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
0033 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
0034 
0035 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0036 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0037 #include "RecoVertex/PrimaryVertexProducer/interface/HITrackFilterForPVFinding.h"
0038 // Fastjet
0039 #include <fastjet/internal/base.hh>
0040 #include "fastjet/PseudoJet.hh"
0041 #include "fastjet/JetDefinition.hh"
0042 #include "fastjet/ClusterSequence.hh"
0043 #include "fastjet/Selector.hh"
0044 
0045 // HepPDTRecord
0046 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0047 
0048 // pile-up
0049 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
0050 
0051 // associator
0052 #include "SimTracker/VertexAssociation/interface/calculateVertexSharedTracks.h"
0053 
0054 // vertexing
0055 #include "RecoVertex/PrimaryVertexProducer/interface/TrackFilterForPVFinding.h"
0056 
0057 // simulated vertex
0058 #include "SimDataFormats/Associations/interface/VertexToTrackingVertexAssociator.h"
0059 
0060 // DQM
0061 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0062 #include "DQMServices/Core/interface/DQMStore.h"
0063 
0064 #include "DataFormats/Math/interface/deltaR.h"
0065 #include "DataFormats/Math/interface/GeantUnits.h"
0066 #include "CLHEP/Units/PhysicalConstants.h"
0067 
0068 // class declaration
0069 class Primary4DVertexValidation : public DQMEDAnalyzer {
0070   typedef math::XYZTLorentzVector LorentzVector;
0071 
0072   // auxiliary class holding simulated vertices
0073   struct simPrimaryVertex {
0074     simPrimaryVertex(double x1, double y1, double z1, double t1)
0075         : x(x1),
0076           y(y1),
0077           z(z1),
0078           t(t1),
0079           pt(0),
0080           ptsq(0),
0081           closest_vertex_distance_z(-1.),
0082           nGenTrk(0),
0083           num_matched_reco_tracks(0),
0084           average_match_quality(0.0) {
0085       ptot.setPx(0);
0086       ptot.setPy(0);
0087       ptot.setPz(0);
0088       ptot.setE(0);
0089       p4 = LorentzVector(0, 0, 0, 0);
0090       r = sqrt(x * x + y * y);
0091     };
0092     double x, y, z, r, t;
0093     HepMC::FourVector ptot;
0094     LorentzVector p4;
0095     double pt;
0096     double ptsq;
0097     double closest_vertex_distance_z;
0098     int nGenTrk;
0099     int num_matched_reco_tracks;
0100     float average_match_quality;
0101     EncodedEventId eventId;
0102     TrackingVertexRef sim_vertex;
0103     int OriginalIndex = -1;
0104 
0105     unsigned int nwosmatch = 0;                    // number of rec vertices dominated by this sim evt (by wos)
0106     unsigned int nwntmatch = 0;                    // number of rec vertices dominated by this sim evt  (by wnt)
0107     std::vector<unsigned int> wos_dominated_recv;  // list of dominated rec vertices (by wos, size==nwosmatch)
0108 
0109     std::map<unsigned int, double> wnt;  // weighted number of tracks in recvtx (by index)
0110     std::map<unsigned int, double> wos;  // weight over sigma**2 in recvtx (by index)
0111     double sumwos = 0;                   // sum of wos in any recvtx
0112     double sumwnt = 0;                   // sum of weighted tracks
0113     unsigned int rec = NOT_MATCHED;      // best match (NO_MATCH if not matched)
0114     unsigned int matchQuality = 0;       // quality flag
0115 
0116     void addTrack(unsigned int irecv, double twos, double twt) {
0117       sumwnt += twt;
0118       if (wnt.find(irecv) == wnt.end()) {
0119         wnt[irecv] = twt;
0120       } else {
0121         wnt[irecv] += twt;
0122       }
0123 
0124       sumwos += twos;
0125       if (wos.find(irecv) == wos.end()) {
0126         wos[irecv] = twos;
0127       } else {
0128         wos[irecv] += twos;
0129       }
0130     }
0131   };
0132 
0133   // auxiliary class holding reconstructed vertices
0134   struct recoPrimaryVertex {
0135     recoPrimaryVertex(double x1, double y1, double z1)
0136         : x(x1),
0137           y(y1),
0138           z(z1),
0139           pt(0),
0140           ptsq(0),
0141           closest_vertex_distance_z(-1.),
0142           nRecoTrk(0),
0143           num_matched_sim_tracks(0),
0144           ndof(0.),
0145           recVtx(nullptr) {
0146       r = sqrt(x * x + y * y);
0147     };
0148     double x, y, z, r;
0149     double pt;
0150     double ptsq;
0151     double closest_vertex_distance_z;
0152     int nRecoTrk;
0153     int num_matched_sim_tracks;
0154     double ndof;
0155     const reco::Vertex* recVtx;
0156     reco::VertexBaseRef recVtxRef;
0157     int OriginalIndex = -1;
0158 
0159     std::map<unsigned int, double> wos;   // sim event -> wos
0160     std::map<unsigned int, double> wnt;   // sim event -> weighted number of truth matched tracks
0161     unsigned int wosmatch = NOT_MATCHED;  // index of the sim event providing the largest contribution to wos
0162     unsigned int wntmatch = NOT_MATCHED;  // index of the sim event providing the highest number of tracks
0163     double sumwos = 0;                    // total sum of wos of all truth matched tracks
0164     double sumwnt = 0;                    // total weighted number of truth matchted tracks
0165     double maxwos = 0;                    // largest wos sum from one sim event (wosmatch)
0166     double maxwnt = 0;                    // largest weighted number of tracks from one sim event (ntmatch)
0167     unsigned int maxwosnt = 0;            // number of tracks from the sim event with highest wos
0168     unsigned int sim = NOT_MATCHED;       // best match (NO_MATCH if not matched)
0169     unsigned int matchQuality = 0;        // quality flag
0170 
0171     bool is_real() { return (matchQuality > 0) && (matchQuality < 99); }
0172 
0173     bool is_fake() { return (matchQuality <= 0) || (matchQuality >= 99); }
0174 
0175     bool is_signal() { return (sim == 0); }
0176 
0177     int split_from() {
0178       if (is_real())
0179         return -1;
0180       if ((maxwos > 0) && (maxwos > 0.3 * sumwos))
0181         return wosmatch;
0182       return -1;
0183     }
0184     bool other_fake() { return (is_fake() && (split_from() < 0)); }
0185 
0186     void addTrack(unsigned int iev, double twos, double wt) {
0187       sumwnt += wt;
0188       if (wnt.find(iev) == wnt.end()) {
0189         wnt[iev] = wt;
0190       } else {
0191         wnt[iev] += wt;
0192       }
0193 
0194       sumwos += twos;
0195       if (wos.find(iev) == wos.end()) {
0196         wos[iev] = twos;
0197       } else {
0198         wos[iev] += twos;
0199       }
0200     }
0201   };
0202 
0203 public:
0204   explicit Primary4DVertexValidation(const edm::ParameterSet&);
0205   ~Primary4DVertexValidation() override;
0206 
0207   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0208   void analyze(const edm::Event&, const edm::EventSetup&) override;
0209   void bookHistograms(DQMStore::IBooker& i, edm::Run const&, edm::EventSetup const&) override;
0210 
0211 private:
0212   void matchReco2Sim(std::vector<recoPrimaryVertex>&,
0213                      std::vector<simPrimaryVertex>&,
0214                      const edm::ValueMap<float>&,
0215                      const edm::ValueMap<float>&,
0216                      const edm::Handle<reco::BeamSpot>&);
0217   bool matchRecoTrack2SimSignal(const reco::TrackBaseRef&);
0218   std::pair<const edm::Ref<std::vector<TrackingParticle>>*, int> getMatchedTP(const reco::TrackBaseRef&,
0219                                                                               const TrackingVertexRef&);
0220   double timeFromTrueMass(double, double, double, double);
0221   bool select(const reco::Vertex&, int level = 0);
0222   void observablesFromJets(const std::vector<reco::Track>&,
0223                            const std::vector<double>&,
0224                            const std::vector<int>&,
0225                            const std::string&,
0226                            unsigned int&,
0227                            double&,
0228                            double&,
0229                            double&,
0230                            double&);
0231   void isParticle(const reco::TrackBaseRef&,
0232                   const edm::ValueMap<float>&,
0233                   const edm::ValueMap<float>&,
0234                   const edm::ValueMap<float>&,
0235                   const edm::ValueMap<float>&,
0236                   const edm::ValueMap<float>&,
0237                   unsigned int&,
0238                   bool&,
0239                   bool&,
0240                   bool&,
0241                   bool&);
0242   void getWosWnt(const reco::Vertex&,
0243                  const reco::TrackBaseRef&,
0244                  const edm::ValueMap<float>&,
0245                  const edm::ValueMap<float>&,
0246                  const edm::Handle<reco::BeamSpot>&,
0247                  double&,
0248                  double&);
0249   std::vector<Primary4DVertexValidation::simPrimaryVertex> getSimPVs(const edm::Handle<TrackingVertexCollection>&);
0250   std::vector<Primary4DVertexValidation::recoPrimaryVertex> getRecoPVs(const edm::Handle<edm::View<reco::Vertex>>&);
0251 
0252   void printMatchedRecoTrackInfo(const reco::Vertex&,
0253                                  const reco::TrackBaseRef&,
0254                                  const TrackingParticleRef&,
0255                                  const unsigned int&);
0256   void printSimVtxRecoVtxInfo(const struct Primary4DVertexValidation::simPrimaryVertex&,
0257                               const struct Primary4DVertexValidation::recoPrimaryVertex&);
0258   const bool trkTPSelLV(const TrackingParticle&);
0259   const bool trkRecSel(const reco::TrackBase&);
0260 
0261   // ----------member data ---------------------------
0262 
0263   const edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> theTTBToken;
0264   TrackFilterForPVFindingBase* theTrackFilter;
0265   const std::string folder_;
0266   static constexpr unsigned int NOT_MATCHED = 66666;
0267   static constexpr double simUnit_ = 1e9;     // sim time in s while reco time in ns
0268   static constexpr double c_ = 2.99792458e1;  // c in cm/ns
0269   static constexpr double mvaL_ = 0.5;        // MVA cuts for MVA categories
0270   static constexpr double mvaH_ = 0.8;
0271   static constexpr double selNdof_ = 4.;
0272   static constexpr double maxRank_ = 8.;
0273   static constexpr double maxTry_ = 10.;
0274   static constexpr double zWosMatchMax_ = 1.;
0275   static constexpr double etacutGEN_ = 4.;   // |eta| < 4;
0276   static constexpr double etacutREC_ = 3.;   // |eta| < 3;
0277   static constexpr double pTcut_ = 0.7;      // PT > 0.7 GeV
0278   static constexpr double deltaZcut_ = 0.1;  // dz separation 1 mm
0279   static constexpr double trackMaxBtlEta_ = 1.5;
0280   static constexpr double trackMinEtlEta_ = 1.6;
0281   static constexpr double trackMaxEtlEta_ = 3.;
0282   static constexpr double tol_ = 1.e-4;          // tolerance on reconstructed track time, [ns]
0283   static constexpr double minThrSumWnt_ = 0.01;  // min threshold for filling histograms with logarithmic scale
0284   static constexpr double minThrSumWos_ = 0.1;
0285   static constexpr double minThrSumPt_ = 0.01;
0286   static constexpr double minThrSumPt2_ = 1.e-3;
0287   static constexpr double minThrMetPt_ = 1.e-3;
0288   static constexpr double minThrSumPz_ = 1.e-4;
0289   static constexpr double rBTL_ = 110.0;
0290   static constexpr double zETL_ = 290.0;
0291 
0292   static constexpr float c_cm_ns = geant_units::operators::convertMmToCm(CLHEP::c_light);  // [mm/ns] -> [cm/ns]
0293 
0294   bool use_only_charged_tracks_;
0295   bool optionalPlots_;
0296   bool use3dNoTime_;
0297   const double minProbHeavy_;
0298   const double trackweightTh_;
0299   const double mvaTh_;
0300   const reco::RecoToSimCollection* r2s_;
0301   const reco::SimToRecoCollection* s2r_;
0302 
0303   edm::EDGetTokenT<reco::TrackCollection> RecTrackToken_;
0304 
0305   edm::EDGetTokenT<std::vector<PileupSummaryInfo>> vecPileupSummaryInfoToken_;
0306 
0307   edm::EDGetTokenT<reco::TrackCollection> trkToken;
0308   edm::EDGetTokenT<TrackingParticleCollection> trackingParticleCollectionToken_;
0309   edm::EDGetTokenT<TrackingVertexCollection> trackingVertexCollectionToken_;
0310   edm::EDGetTokenT<reco::SimToRecoCollection> simToRecoAssociationToken_;
0311   edm::EDGetTokenT<reco::RecoToSimCollection> recoToSimAssociationToken_;
0312   edm::EDGetTokenT<reco::BeamSpot> RecBeamSpotToken_;
0313   edm::EDGetTokenT<edm::View<reco::Vertex>> Rec4DVerToken_;
0314 
0315   edm::EDGetTokenT<edm::ValueMap<int>> trackAssocToken_;
0316   edm::EDGetTokenT<edm::ValueMap<float>> pathLengthToken_;
0317   edm::EDGetTokenT<edm::ValueMap<float>> momentumToken_;
0318   edm::EDGetTokenT<edm::ValueMap<float>> timeToken_;
0319 
0320   edm::EDGetTokenT<edm::ValueMap<float>> t0PidToken_;
0321   edm::EDGetTokenT<edm::ValueMap<float>> sigmat0PidToken_;
0322   edm::EDGetTokenT<edm::ValueMap<float>> t0SafePidToken_;
0323   edm::EDGetTokenT<edm::ValueMap<float>> sigmat0SafePidToken_;
0324   edm::EDGetTokenT<edm::ValueMap<float>> trackMVAQualToken_;
0325 
0326   edm::EDGetTokenT<edm::ValueMap<float>> tmtdToken_;
0327   edm::EDGetTokenT<edm::ValueMap<float>> tofPiToken_;
0328   edm::EDGetTokenT<edm::ValueMap<float>> tofKToken_;
0329   edm::EDGetTokenT<edm::ValueMap<float>> tofPToken_;
0330   edm::EDGetTokenT<edm::ValueMap<float>> probPiToken_;
0331   edm::EDGetTokenT<edm::ValueMap<float>> probKToken_;
0332   edm::EDGetTokenT<edm::ValueMap<float>> probPToken_;
0333   edm::ESGetToken<HepPDT::ParticleDataTable, edm::DefaultRecord> pdtToken_;
0334 
0335   // histogram declaration
0336   MonitorElement* meUnAssocTracks_;
0337   MonitorElement* meUnAssocTracksFake_;
0338   MonitorElement* meFractionUnAssocTracks_;
0339   MonitorElement* meFractionUnAssocTracksFake_;
0340   MonitorElement* meTrackEffPtTot_;
0341   MonitorElement* meTrackMatchedTPEffPtTot_;
0342   MonitorElement* meTrackMatchedTPEffPtMtd_;
0343   MonitorElement* meTrackEffEtaTot_;
0344   MonitorElement* meTrackMatchedTPEffEtaTot_;
0345   MonitorElement* meTrackMatchedTPEffEtaMtd_;
0346   MonitorElement* meTrackMatchedTPResTot_;
0347   MonitorElement* meTrackMatchedTPPullTot_;
0348   MonitorElement* meTrackResTot_;
0349   MonitorElement* meTrackPullTot_;
0350   MonitorElement* meTrackRes_[3];
0351   MonitorElement* meTrackPull_[3];
0352   MonitorElement* meTrackResMass_[3];
0353   MonitorElement* meTrackResMassTrue_[3];
0354   MonitorElement* meTrackMatchedTPZposResTot_;
0355   MonitorElement* meTrackZposResTot_;
0356   MonitorElement* meTrackZposRes_[3];
0357   MonitorElement* meTrack3DposRes_[3];
0358   MonitorElement* meTimeRes_;
0359   MonitorElement* meTimePull_;
0360   MonitorElement* meTimeSignalRes_;
0361   MonitorElement* meTimeSignalPull_;
0362   MonitorElement* mePUvsRealV_;
0363   MonitorElement* mePUvsFakeV_;
0364   MonitorElement* mePUvsOtherFakeV_;
0365   MonitorElement* mePUvsSplitV_;
0366   MonitorElement* meMatchQual_;
0367   MonitorElement* meDeltaZrealreal_;
0368   MonitorElement* meDeltaZfakefake_;
0369   MonitorElement* meDeltaZfakereal_;
0370   MonitorElement* meDeltaTrealreal_;
0371   MonitorElement* meDeltaTfakefake_;
0372   MonitorElement* meDeltaTfakereal_;
0373   MonitorElement* meRecoPosInSimCollection_;
0374   MonitorElement* meRecoPosInRecoOrigCollection_;
0375   MonitorElement* meSimPosInSimOrigCollection_;
0376   MonitorElement* meRecoPVPosSignal_;
0377   MonitorElement* meRecoPVPosSignalNotHighestPt_;
0378   MonitorElement* meRecoVtxVsLineDensity_;
0379   MonitorElement* meRecVerNumber_;
0380   MonitorElement* meRecPVZ_;
0381   MonitorElement* meRecPVT_;
0382   MonitorElement* meSimVerNumber_;
0383   MonitorElement* meSimPVZ_;
0384   MonitorElement* meSimPVT_;
0385   MonitorElement* meSimPVTvsZ_;
0386 
0387   MonitorElement* meVtxTrackMult_;
0388   MonitorElement* meVtxTrackMultPassNdof_;
0389   MonitorElement* meVtxTrackMultFailNdof_;
0390   MonitorElement* meVtxTrackW_;
0391   MonitorElement* meVtxTrackWnt_;
0392   MonitorElement* meVtxTrackRecLVMult_;
0393   MonitorElement* meVtxTrackRecLVW_;
0394   MonitorElement* meVtxTrackRecLVWnt_;
0395 
0396   MonitorElement* mePUTrackMult_;
0397   MonitorElement* mePUTrackRelMult_;
0398   MonitorElement* meFakeTrackRelMult_;
0399   MonitorElement* mePUTrackSumWnt_;
0400   MonitorElement* mePUTrackRelSumWnt_;
0401   MonitorElement* mePUTrackSumWos_;
0402   MonitorElement* mePUTrackRelSumWos_;
0403   MonitorElement* meSecTrackSumWos_;
0404   MonitorElement* meSecTrackRelSumWos_;
0405   MonitorElement* meFakeTrackRelSumWos_;
0406   MonitorElement* mePUTrackWnt_;
0407   MonitorElement* mePUTrackSumPt_;
0408   MonitorElement* mePUTrackRelSumPt_;
0409   MonitorElement* mePUTrackSumPt2_;
0410   MonitorElement* mePUTrackRelSumPt2_;
0411 
0412   MonitorElement* mePUTrackRelMultvsMult_;
0413   MonitorElement* meFakeTrackRelMultvsMult_;
0414   MonitorElement* mePUTrackRelSumWntvsSumWnt_;
0415   MonitorElement* mePUTrackRelSumWosvsSumWos_;
0416   MonitorElement* meFakeTrackRelSumWosvsSumWos_;
0417   MonitorElement* meSecTrackRelSumWosvsSumWos_;
0418   MonitorElement* mePUTrackRelSumPtvsSumPt_;
0419   MonitorElement* mePUTrackRelSumPt2vsSumPt2_;
0420 
0421   MonitorElement* mePUTrackRecLVMult_;
0422   MonitorElement* mePUTrackRecLVRelMult_;
0423   MonitorElement* meFakeTrackRecLVRelMult_;
0424   MonitorElement* mePUTrackRecLVSumWnt_;
0425   MonitorElement* mePUTrackRecLVRelSumWnt_;
0426   MonitorElement* mePUTrackRecLVSumWos_;
0427   MonitorElement* mePUTrackRecLVRelSumWos_;
0428   MonitorElement* meSecTrackRecLVSumWos_;
0429   MonitorElement* meSecTrackRecLVRelSumWos_;
0430   MonitorElement* meFakeTrackRecLVRelSumWos_;
0431   MonitorElement* mePUTrackRecLVWnt_;
0432   MonitorElement* mePUTrackRecLVSumPt_;
0433   MonitorElement* mePUTrackRecLVRelSumPt_;
0434   MonitorElement* mePUTrackRecLVSumPt2_;
0435   MonitorElement* mePUTrackRecLVRelSumPt2_;
0436 
0437   MonitorElement* mePUTrackRecLVRelMultvsMult_;
0438   MonitorElement* meFakeTrackRecLVRelMultvsMult_;
0439   MonitorElement* mePUTrackRecLVRelSumWntvsSumWnt_;
0440   MonitorElement* mePUTrackRecLVRelSumWosvsSumWos_;
0441   MonitorElement* meSecTrackRecLVRelSumWosvsSumWos_;
0442   MonitorElement* meFakeTrackRecLVRelSumWosvsSumWos_;
0443   MonitorElement* mePUTrackRecLVRelSumPtvsSumPt_;
0444   MonitorElement* mePUTrackRecLVRelSumPt2vsSumPt2_;
0445 
0446   MonitorElement* meJetsPUMult_;
0447   MonitorElement* meJetsPUHt_;
0448   MonitorElement* meJetsPUSumPt2_;
0449   MonitorElement* meJetsPUMetPt_;
0450   MonitorElement* meJetsPUSumPz_;
0451   MonitorElement* meJetsPURelMult_;
0452   MonitorElement* meJetsPURelHt_;
0453   MonitorElement* meJetsPURelSumPt2_;
0454   MonitorElement* meJetsFakeRelSumPt2_;
0455   MonitorElement* meJetsPURelMetPt_;
0456   MonitorElement* meJetsPURelSumPz_;
0457 
0458   MonitorElement* meJetsPURelMultvsMult_;
0459   MonitorElement* meJetsPURelHtvsHt_;
0460   MonitorElement* meJetsPURelSumPt2vsSumPt2_;
0461   MonitorElement* meJetsFakeRelSumPt2vsSumPt2_;
0462   MonitorElement* meJetsPURelMetPtvsMetPt_;
0463   MonitorElement* meJetsPURelSumPzvsSumPz_;
0464 
0465   MonitorElement* meJetsRecLVPUMult_;
0466   MonitorElement* meJetsRecLVPUHt_;
0467   MonitorElement* meJetsRecLVPUSumPt2_;
0468   MonitorElement* meJetsRecLVPUMetPt_;
0469   MonitorElement* meJetsRecLVPUSumPz_;
0470   MonitorElement* meJetsRecLVPURelMult_;
0471   MonitorElement* meJetsRecLVPURelHt_;
0472   MonitorElement* meJetsRecLVPURelSumPt2_;
0473   MonitorElement* meJetsRecLVFakeRelSumPt2_;
0474   MonitorElement* meJetsRecLVPURelMetPt_;
0475   MonitorElement* meJetsRecLVPURelSumPz_;
0476 
0477   MonitorElement* meJetsRecLVPURelMultvsMult_;
0478   MonitorElement* meJetsRecLVPURelHtvsHt_;
0479   MonitorElement* meJetsRecLVPURelSumPt2vsSumPt2_;
0480   MonitorElement* meJetsRecLVFakeRelSumPt2vsSumPt2_;
0481   MonitorElement* meJetsRecLVPURelMetPtvsMetPt_;
0482   MonitorElement* meJetsRecLVPURelSumPzvsSumPz_;
0483 
0484   // some tests
0485   MonitorElement* meTrackResLowPTot_;
0486   MonitorElement* meTrackResHighPTot_;
0487   MonitorElement* meTrackPullLowPTot_;
0488   MonitorElement* meTrackPullHighPTot_;
0489 
0490   MonitorElement* meTrackResLowP_[3];
0491   MonitorElement* meTrackResHighP_[3];
0492   MonitorElement* meTrackPullLowP_[3];
0493   MonitorElement* meTrackPullHighP_[3];
0494 
0495   MonitorElement* meTrackResMassProtons_[3];
0496   MonitorElement* meTrackResMassTrueProtons_[3];
0497   MonitorElement* meTrackResMassPions_[3];
0498   MonitorElement* meTrackResMassTruePions_[3];
0499 
0500   MonitorElement* meBarrelPIDp_;
0501   MonitorElement* meEndcapPIDp_;
0502 
0503   MonitorElement* meBarrelNoPIDtype_;
0504   MonitorElement* meEndcapNoPIDtype_;
0505 
0506   MonitorElement* meBarrelTruePiNoPID_;
0507   MonitorElement* meBarrelTrueKNoPID_;
0508   MonitorElement* meBarrelTruePNoPID_;
0509   MonitorElement* meEndcapTruePiNoPID_;
0510   MonitorElement* meEndcapTrueKNoPID_;
0511   MonitorElement* meEndcapTruePNoPID_;
0512 
0513   MonitorElement* meBarrelTruePiAsPi_;
0514   MonitorElement* meBarrelTruePiAsK_;
0515   MonitorElement* meBarrelTruePiAsP_;
0516   MonitorElement* meEndcapTruePiAsPi_;
0517   MonitorElement* meEndcapTruePiAsK_;
0518   MonitorElement* meEndcapTruePiAsP_;
0519 
0520   MonitorElement* meBarrelTrueKAsPi_;
0521   MonitorElement* meBarrelTrueKAsK_;
0522   MonitorElement* meBarrelTrueKAsP_;
0523   MonitorElement* meEndcapTrueKAsPi_;
0524   MonitorElement* meEndcapTrueKAsK_;
0525   MonitorElement* meEndcapTrueKAsP_;
0526 
0527   MonitorElement* meBarrelTruePAsPi_;
0528   MonitorElement* meBarrelTruePAsK_;
0529   MonitorElement* meBarrelTruePAsP_;
0530   MonitorElement* meEndcapTruePAsPi_;
0531   MonitorElement* meEndcapTruePAsK_;
0532   MonitorElement* meEndcapTruePAsP_;
0533 
0534   // Histograms for study of no PID tracks
0535 
0536   //Time residual
0537   MonitorElement* meTrackTimeResCorrectPID_;
0538   MonitorElement* meTrackTimeResWrongPID_;
0539   MonitorElement* meTrackTimeResNoPID_;
0540   MonitorElement* meNoPIDTrackTimeResNoPIDType_[3];
0541   MonitorElement* meTrackTimeResNoPIDtruePi_;
0542   MonitorElement* meTrackTimeResNoPIDtrueK_;
0543   MonitorElement* meTrackTimeResNoPIDtrueP_;
0544 
0545   //Time pull
0546   MonitorElement* meTrackTimePullCorrectPID_;
0547   MonitorElement* meTrackTimePullWrongPID_;
0548   MonitorElement* meTrackTimePullNoPID_;
0549   MonitorElement* meNoPIDTrackTimePullNoPIDType_[3];
0550   MonitorElement* meTrackTimePullNoPIDtruePi_;
0551   MonitorElement* meTrackTimePullNoPIDtrueK_;
0552   MonitorElement* meTrackTimePullNoPIDtrueP_;
0553 
0554   //Sigma
0555   MonitorElement* meTrackTimeSigmaCorrectPID_;
0556   MonitorElement* meTrackTimeSigmaWrongPID_;
0557   MonitorElement* meTrackTimeSigmaNoPID_;
0558   MonitorElement* meNoPIDTrackSigmaNoPIDType_[3];
0559 
0560   //MVA
0561   MonitorElement* meTrackMVACorrectPID_;
0562   MonitorElement* meTrackMVAWrongPID_;
0563   MonitorElement* meTrackMVANoPID_;
0564   MonitorElement* meNoPIDTrackMVANoPIDType_[3];
0565 };
0566 
0567 // constructors and destructor
0568 Primary4DVertexValidation::Primary4DVertexValidation(const edm::ParameterSet& iConfig)
0569     : theTTBToken(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0570       folder_(iConfig.getParameter<std::string>("folder")),
0571       use_only_charged_tracks_(iConfig.getParameter<bool>("useOnlyChargedTracks")),
0572       optionalPlots_(iConfig.getUntrackedParameter<bool>("optionalPlots")),
0573       use3dNoTime_(iConfig.getParameter<bool>("use3dNoTime")),
0574       minProbHeavy_(iConfig.getParameter<double>("minProbHeavy")),
0575       trackweightTh_(iConfig.getParameter<double>("trackweightTh")),
0576       mvaTh_(iConfig.getParameter<double>("mvaTh")),
0577       pdtToken_(esConsumes<HepPDT::ParticleDataTable, edm::DefaultRecord>()) {
0578   vecPileupSummaryInfoToken_ = consumes<std::vector<PileupSummaryInfo>>(edm::InputTag(std::string("addPileupInfo")));
0579   trkToken = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("TrackLabel"));
0580   trackingParticleCollectionToken_ =
0581       consumes<TrackingParticleCollection>(iConfig.getParameter<edm::InputTag>("SimTag"));
0582   trackingVertexCollectionToken_ = consumes<TrackingVertexCollection>(iConfig.getParameter<edm::InputTag>("SimTag"));
0583   simToRecoAssociationToken_ =
0584       consumes<reco::SimToRecoCollection>(iConfig.getParameter<edm::InputTag>("TPtoRecoTrackAssoc"));
0585   recoToSimAssociationToken_ =
0586       consumes<reco::RecoToSimCollection>(iConfig.getParameter<edm::InputTag>("TPtoRecoTrackAssoc"));
0587   RecTrackToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("mtdTracks"));
0588   RecBeamSpotToken_ = consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("offlineBS"));
0589   Rec4DVerToken_ = consumes<edm::View<reco::Vertex>>(iConfig.getParameter<edm::InputTag>("offline4DPV"));
0590   trackAssocToken_ = consumes<edm::ValueMap<int>>(iConfig.getParameter<edm::InputTag>("trackAssocSrc"));
0591   pathLengthToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("pathLengthSrc"));
0592   momentumToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("momentumSrc"));
0593   timeToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("timeSrc"));
0594   t0PidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("t0PID"));
0595   sigmat0PidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmat0PID"));
0596   t0SafePidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("t0SafePID"));
0597   sigmat0SafePidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmat0SafePID"));
0598   trackMVAQualToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("trackMVAQual"));
0599   tmtdToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tmtd"));
0600   tofPiToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tofPi"));
0601   tofKToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tofK"));
0602   tofPToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tofP"));
0603   probPiToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("probPi"));
0604   probKToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("probK"));
0605   probPToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("probP"));
0606   std::string trackSelectionAlgorithm =
0607       iConfig.getParameter<edm::ParameterSet>("TkFilterParameters").getParameter<std::string>("algorithm");
0608   if (trackSelectionAlgorithm == "filter") {
0609     theTrackFilter = new TrackFilterForPVFinding(iConfig.getParameter<edm::ParameterSet>("TkFilterParameters"));
0610   } else if (trackSelectionAlgorithm == "filterWithThreshold") {
0611     theTrackFilter = new HITrackFilterForPVFinding(iConfig.getParameter<edm::ParameterSet>("TkFilterParameters"));
0612   } else {
0613     edm::LogWarning("Primary4DVertexValidation")
0614         << "unknown track selection algorithm: " + trackSelectionAlgorithm << std::endl;
0615   }
0616 }
0617 
0618 Primary4DVertexValidation::~Primary4DVertexValidation() {}
0619 
0620 //
0621 // member functions
0622 //
0623 void Primary4DVertexValidation::bookHistograms(DQMStore::IBooker& ibook,
0624                                                edm::Run const& iRun,
0625                                                edm::EventSetup const& iSetup) {
0626   ibook.setCurrentFolder(folder_);
0627   // --- histograms booking
0628   meUnAssocTracks_ = ibook.book1D("UnAssocTracks", "Log10(Unassociated tracks)", 160, 0.5, 4.5);
0629   meUnAssocTracksFake_ = ibook.book1D("UnAssocTracksFake", "Log10(Unassociated fake tracks)", 160, 0.5, 4.5);
0630   meFractionUnAssocTracks_ = ibook.book1D("FractionUnAssocTracks", "Fraction Unassociated tracks", 160, 0.0, 1.);
0631   meFractionUnAssocTracksFake_ =
0632       ibook.book1D("FractionUnAssocTracksFake", "Fraction Unassociated fake tracks", 160, 0.0, 1.);
0633   meTrackEffPtTot_ = ibook.book1D("EffPtTot", "Pt of tracks associated to LV; track pt [GeV] ", 110, 0., 11.);
0634   meTrackEffEtaTot_ = ibook.book1D("EffEtaTot", "Eta of tracks associated to LV; track eta ", 66, 0., 3.3);
0635   meTrackMatchedTPEffPtTot_ =
0636       ibook.book1D("MatchedTPEffPtTot", "Pt of tracks associated to LV matched to TP; track pt [GeV] ", 110, 0., 11.);
0637   meTrackMatchedTPEffPtMtd_ = ibook.book1D(
0638       "MatchedTPEffPtMtd", "Pt of tracks associated to LV matched to TP with time; track pt [GeV] ", 110, 0., 11.);
0639   meTrackMatchedTPEffEtaTot_ =
0640       ibook.book1D("MatchedTPEffEtaTot", "Eta of tracks associated to LV matched to TP; track eta ", 66, 0., 3.3);
0641   meTrackMatchedTPEffEtaMtd_ = ibook.book1D(
0642       "MatchedTPEffEtaMtd", "Eta of tracks associated to LV matched to TP with time; track eta ", 66, 0., 3.3);
0643   meTrackMatchedTPResTot_ =
0644       ibook.book1D("MatchedTPTrackRes",
0645                    "t_{rec} - t_{sim} for tracks associated to LV matched to TP; t_{rec} - t_{sim} [ns] ",
0646                    120,
0647                    -0.15,
0648                    0.15);
0649   meTrackResTot_ = ibook.book1D("TrackRes", "t_{rec} - t_{sim} for tracks; t_{rec} - t_{sim} [ns] ", 120, -0.15, 0.15);
0650   meTrackRes_[0] = ibook.book1D(
0651       "TrackRes-LowMVA", "t_{rec} - t_{sim} for tracks with MVA < 0.5; t_{rec} - t_{sim} [ns] ", 100, -1., 1.);
0652   meTrackRes_[1] = ibook.book1D(
0653       "TrackRes-MediumMVA", "t_{rec} - t_{sim} for tracks with 0.5 < MVA < 0.8; t_{rec} - t_{sim} [ns] ", 100, -1., 1.);
0654   meTrackRes_[2] = ibook.book1D(
0655       "TrackRes-HighMVA", "t_{rec} - t_{sim} for tracks with MVA > 0.8; t_{rec} - t_{sim} [ns] ", 100, -1., 1.);
0656   if (optionalPlots_) {
0657     meTrackResMass_[0] = ibook.book1D(
0658         "TrackResMass-LowMVA", "t_{rec} - t_{est} for tracks with MVA < 0.5; t_{rec} - t_{est} [ns] ", 100, -1., 1.);
0659     meTrackResMass_[1] = ibook.book1D("TrackResMass-MediumMVA",
0660                                       "t_{rec} - t_{est} for tracks with 0.5 < MVA < 0.8; t_{rec} - t_{est} [ns] ",
0661                                       100,
0662                                       -1.,
0663                                       1.);
0664     meTrackResMass_[2] = ibook.book1D(
0665         "TrackResMass-HighMVA", "t_{rec} - t_{est} for tracks with MVA > 0.8; t_{rec} - t_{est} [ns] ", 100, -1., 1.);
0666     meTrackResMassTrue_[0] = ibook.book1D(
0667         "TrackResMassTrue-LowMVA", "t_{est} - t_{sim} for tracks with MVA < 0.5; t_{est} - t_{sim} [ns] ", 100, -1., 1.);
0668     meTrackResMassTrue_[1] = ibook.book1D("TrackResMassTrue-MediumMVA",
0669                                           "t_{est} - t_{sim} for tracks with 0.5 < MVA < 0.8; t_{est} - t_{sim} [ns] ",
0670                                           100,
0671                                           -1.,
0672                                           1.);
0673     meTrackResMassTrue_[2] = ibook.book1D("TrackResMassTrue-HighMVA",
0674                                           "t_{est} - t_{sim} for tracks with MVA > 0.8; t_{est} - t_{sim} [ns] ",
0675                                           100,
0676                                           -1.,
0677                                           1.);
0678   }
0679   meTrackMatchedTPPullTot_ = ibook.book1D(
0680       "MatchedTPTrackPull", "Pull for tracks associated to LV matched to TP; (t_{rec}-t_{sim})/#sigma_{t}", 50, -5., 5.);
0681   meTrackPullTot_ = ibook.book1D("TrackPull", "Pull for tracks; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
0682   meTrackPull_[0] =
0683       ibook.book1D("TrackPull-LowMVA", "Pull for tracks with MVA < 0.5; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
0684   meTrackPull_[1] = ibook.book1D(
0685       "TrackPull-MediumMVA", "Pull for tracks with 0.5 < MVA < 0.8; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
0686   meTrackPull_[2] =
0687       ibook.book1D("TrackPull-HighMVA", "Pull for tracks with MVA > 0.8; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
0688   meTrackMatchedTPZposResTot_ =
0689       ibook.book1D("MatchedTPTrackZposResTot",
0690                    "Z_{PCA} - Z_{sim} for tracks associated to LV matched to TP;Z_{PCA} - Z_{sim} [cm] ",
0691                    50,
0692                    -0.1,
0693                    0.1);
0694   meTrackZposResTot_ =
0695       ibook.book1D("TrackZposResTot", "Z_{PCA} - Z_{sim} for tracks;Z_{PCA} - Z_{sim} [cm] ", 50, -0.5, 0.5);
0696   meTrackZposRes_[0] = ibook.book1D(
0697       "TrackZposRes-LowMVA", "Z_{PCA} - Z_{sim} for tracks with MVA < 0.5;Z_{PCA} - Z_{sim} [cm] ", 50, -0.5, 0.5);
0698   meTrackZposRes_[1] = ibook.book1D("TrackZposRes-MediumMVA",
0699                                     "Z_{PCA} - Z_{sim} for tracks with 0.5 < MVA < 0.8 ;Z_{PCA} - Z_{sim} [cm] ",
0700                                     50,
0701                                     -0.5,
0702                                     0.5);
0703   meTrackZposRes_[2] = ibook.book1D(
0704       "TrackZposRes-HighMVA", "Z_{PCA} - Z_{sim} for tracks with MVA > 0.8 ;Z_{PCA} - Z_{sim} [cm] ", 50, -0.5, 0.5);
0705   meTrack3DposRes_[0] =
0706       ibook.book1D("Track3DposRes-LowMVA",
0707                    "3dPos_{PCA} - 3dPos_{sim} for tracks with MVA < 0.5 ;3dPos_{PCA} - 3dPos_{sim} [cm] ",
0708                    50,
0709                    -0.5,
0710                    0.5);
0711   meTrack3DposRes_[1] =
0712       ibook.book1D("Track3DposRes-MediumMVA",
0713                    "3dPos_{PCA} - 3dPos_{sim} for tracks with 0.5 < MVA < 0.8 ;3dPos_{PCA} - 3dPos_{sim} [cm] ",
0714                    50,
0715                    -0.5,
0716                    0.5);
0717   meTrack3DposRes_[2] =
0718       ibook.book1D("Track3DposRes-HighMVA",
0719                    "3dPos_{PCA} - 3dPos_{sim} for tracks with MVA > 0.8;3dPos_{PCA} - 3dPos_{sim} [cm] ",
0720                    50,
0721                    -0.5,
0722                    0.5);
0723   meTimeRes_ = ibook.book1D("TimeRes", "t_{rec} - t_{sim} ;t_{rec} - t_{sim} [ns] ", 100, -0.2, 0.2);
0724   meTimePull_ = ibook.book1D("TimePull", "Pull; t_{rec} - t_{sim}/#sigma_{t rec}", 100, -10., 10.);
0725   meTimeSignalRes_ =
0726       ibook.book1D("TimeSignalRes", "t_{rec} - t_{sim} for signal ;t_{rec} - t_{sim} [ns] ", 50, -0.1, 0.1);
0727   meTimeSignalPull_ =
0728       ibook.book1D("TimeSignalPull", "Pull for signal; t_{rec} - t_{sim}/#sigma_{t rec}", 100, -10., 10.);
0729   mePUvsRealV_ =
0730       ibook.bookProfile("PUvsReal", "#PU vertices vs #real matched vertices;#PU;#real ", 100, 0, 300, 100, 0, 200);
0731   mePUvsFakeV_ =
0732       ibook.bookProfile("PUvsFake", "#PU vertices vs #fake matched vertices;#PU;#fake ", 100, 0, 300, 100, 0, 20);
0733   mePUvsOtherFakeV_ = ibook.bookProfile(
0734       "PUvsOtherFake", "#PU vertices vs #other fake matched vertices;#PU;#other fake ", 100, 0, 300, 100, 0, 20);
0735   mePUvsSplitV_ =
0736       ibook.bookProfile("PUvsSplit", "#PU vertices vs #split matched vertices;#PU;#split ", 100, 0, 300, 100, 0, 20);
0737   meMatchQual_ = ibook.book1D("MatchQuality", "RECO-SIM vertex match quality; ", 8, 0, 8.);
0738   meDeltaZrealreal_ = ibook.book1D("DeltaZrealreal", "#Delta Z real-real; |#Delta Z (r-r)| [cm]", 100, 0, 0.5);
0739   meDeltaZfakefake_ = ibook.book1D("DeltaZfakefake", "#Delta Z fake-fake; |#Delta Z (f-f)| [cm]", 100, 0, 0.5);
0740   meDeltaZfakereal_ = ibook.book1D("DeltaZfakereal", "#Delta Z fake-real; |#Delta Z (f-r)| [cm]", 100, 0, 0.5);
0741   meDeltaTrealreal_ = ibook.book1D("DeltaTrealreal", "#Delta T real-real; |#Delta T (r-r)| [sigma]", 60, 0., 30.);
0742   meDeltaTfakefake_ = ibook.book1D("DeltaTfakefake", "#Delta T fake-fake; |#Delta T (f-f)| [sigma]", 60, 0., 30.);
0743   meDeltaTfakereal_ = ibook.book1D("DeltaTfakereal", "#Delta T fake-real; |#Delta T (f-r)| [sigma]", 60, 0., 30.);
0744   if (optionalPlots_) {
0745     meRecoPosInSimCollection_ = ibook.book1D(
0746         "RecoPosInSimCollection", "Sim signal vertex index associated to Reco signal vertex; Sim PV index", 200, 0, 200);
0747     meRecoPosInRecoOrigCollection_ =
0748         ibook.book1D("RecoPosInRecoOrigCollection", "Reco signal index in OrigCollection; Reco index", 200, 0, 200);
0749     meSimPosInSimOrigCollection_ =
0750         ibook.book1D("SimPosInSimOrigCollection", "Sim signal index in OrigCollection; Sim index", 200, 0, 200);
0751   }
0752   meRecoPVPosSignal_ =
0753       ibook.book1D("RecoPVPosSignal", "Position in reco collection of PV associated to sim signal", 20, 0, 20);
0754   meRecoPVPosSignalNotHighestPt_ =
0755       ibook.book1D("RecoPVPosSignalNotHighestPt",
0756                    "Position in reco collection of PV associated to sim signal not highest Pt",
0757                    20,
0758                    0,
0759                    20);
0760   meRecVerNumber_ = ibook.book1D("RecVerNumber", "RECO Vertex Number: Number of vertices", 50, 0, 250);
0761   meSimVerNumber_ = ibook.book1D("SimVerNumber", "SIM Vertex Number: Number of vertices", 50, 0, 250);
0762   meRecPVZ_ = ibook.book1D("recPVZ", "#Rec vertices/10 mm", 30, -15., 15.);
0763   meRecPVT_ = ibook.book1D("recPVT", "#Rec vertices/50 ps", 30, -0.75, 0.75);
0764   meSimPVZ_ = ibook.book1D("simPVZ", "#Sim vertices/10 mm", 30, -15., 15.);
0765   meSimPVT_ = ibook.book1D("simPVT", "#Sim vertices/50 ps", 30, -0.75, 0.75);
0766   meSimPVTvsZ_ = ibook.bookProfile("simPVTvsZ", "PV Time vs Z", 30, -15., 15., 30, -0.75, 0.75);
0767 
0768   meVtxTrackMult_ = ibook.book1D("VtxTrackMult", "Log10(Vertex track multiplicity)", 80, 0.5, 2.5);
0769   meVtxTrackMultPassNdof_ =
0770       ibook.book1D("VtxTrackMultPassNdof", "Log10(Vertex track multiplicity for ndof>4)", 80, 0.5, 2.5);
0771   meVtxTrackMultFailNdof_ = ibook.book1D("VtxTrackMultFailNdof", "Vertex track multiplicity for ndof<4", 10, 0., 10.);
0772   meVtxTrackW_ = ibook.book1D("VtxTrackW", "Vertex track weight (all)", 50, 0., 1.);
0773   meVtxTrackWnt_ = ibook.book1D("VtxTrackWnt", "Vertex track Wnt", 50, 0., 1.);
0774   meVtxTrackRecLVMult_ =
0775       ibook.book1D("VtxTrackRecLVMult", "Log10(Vertex track multiplicity) for matched LV", 80, 0.5, 2.5);
0776   meVtxTrackRecLVW_ = ibook.book1D("VtxTrackRecLVW", "Vertex track weight for matched LV (all)", 50, 0., 1.);
0777   meVtxTrackRecLVWnt_ = ibook.book1D("VtxTrackRecLVWnt", "Vertex track Wnt for matched LV", 50, 0., 1.);
0778 
0779   mePUTrackRelMult_ = ibook.book1D(
0780       "PUTrackRelMult", "Relative multiplicity of PU tracks for matched vertices; #PUTrks/#Trks", 50, 0., 1.);
0781   meFakeTrackRelMult_ = ibook.book1D(
0782       "FakeTrackRelMult", "Relative multiplicity of fake tracks for matched vertices; #fakeTrks/#Trks", 50, 0., 1.);
0783   mePUTrackRelSumWnt_ =
0784       ibook.book1D("PUTrackRelSumWnt",
0785                    "Relative Sum of wnt of PU tracks for matched vertices; PUSumW*min(Pt, 1.)/SumW*min(Pt, 1.)",
0786                    50,
0787                    0.,
0788                    1.);
0789   mePUTrackRelSumWos_ = ibook.book1D(
0790       "PUTrackRelSumWos", "Relative Sum of wos of PU tracks for matched vertices; PUSumWos/SumWos", 50, 0., 1.);
0791   meSecTrackRelSumWos_ =
0792       ibook.book1D("SecTrackRelSumWos",
0793                    "Relative Sum of wos of tracks from secondary vtx for matched vertices; SecSumWos/SumWos",
0794                    50,
0795                    0.,
0796                    1.);
0797   meFakeTrackRelSumWos_ = ibook.book1D(
0798       "FakeTrackRelSumWos", "Relative Sum of wos of fake tracks for matched vertices; FakeSumWos/SumWos", 50, 0., 1.);
0799   mePUTrackRelSumPt_ = ibook.book1D(
0800       "PUTrackRelSumPt", "Relative Sum of Pt of PU tracks for matched vertices; PUSumPt/SumPt", 50, 0., 1.);
0801   mePUTrackRelSumPt2_ = ibook.book1D(
0802       "PUTrackRelSumPt2", "Relative Sum of Pt2 for PU tracks for matched vertices; PUSumPt2/SumPt2", 50, 0., 1.);
0803   mePUTrackRecLVRelMult_ = ibook.book1D(
0804       "PUTrackRecLVRelMult", "Relative multiplicity of PU tracks for matched LV; #PUTrks/#Trks", 50, 0., 1.);
0805   meFakeTrackRecLVRelMult_ = ibook.book1D(
0806       "FakeTrackRecLVRelMult", "Relative multiplicity of fake tracks for matched LV; #FakeTrks/#Trks", 50, 0., 1.);
0807   mePUTrackRecLVRelSumWnt_ =
0808       ibook.book1D("PUTrackRecLVRelSumWnt",
0809                    "Relative Sum of Wnt of PU tracks for matched LV; PUSumW*min(Pt, 1.)/SumW*min(Pt, 1.)",
0810                    50,
0811                    0.,
0812                    1.);
0813   mePUTrackRecLVRelSumWos_ = ibook.book1D(
0814       "PUTrackRecLVRelSumWos", "Relative Sum of Wos of PU tracks for matched LV; PUSumWos/SumWos", 50, 0., 1.);
0815   meSecTrackRecLVRelSumWos_ =
0816       ibook.book1D("SecTrackRecLVRelSumWos",
0817                    "Relative Sum of wos of tracks from secondary vtx for matched LV; SecSumWos/SumWos",
0818                    50,
0819                    0.,
0820                    1.);
0821   meFakeTrackRecLVRelSumWos_ = ibook.book1D(
0822       "FakeTrackRecLVRelSumWos", "Relative Sum of wos of fake tracks for matched LV; FakeSumWos/SumWos", 50, 0., 1.);
0823   mePUTrackRecLVRelSumPt_ =
0824       ibook.book1D("PUTrackRecLVRelSumPt", "Relative Sum of Pt of PU tracks for matched LV; PUSumPt/SumPt", 50, 0., 1.);
0825   mePUTrackRecLVRelSumPt2_ = ibook.book1D(
0826       "PUTrackRecLVRelSumPt2", "Relative Sum of Pt2 of PU tracks for matched LV; PUSumPt2/SumPt2", 50, 0., 1.);
0827 
0828   if (optionalPlots_) {
0829     mePUTrackMult_ = ibook.book1D("PUTrackMult", "Number of PU tracks for matched vertices; #PUTrks", 50, 0., 100.);
0830     mePUTrackWnt_ = ibook.book1D("PUTrackWnt", "Wnt of PU tracks for matched vertices; PUTrkW*min(Pt, 1.)", 50, 0., 1.);
0831     mePUTrackSumWnt_ = ibook.book1D(
0832         "PUTrackSumWnt", "Sum of wnt of PU tracks for matched vertices; log10(PUSumW*min(Pt, 1.))", 50, -2., 3.);
0833     mePUTrackSumWos_ =
0834         ibook.book1D("PUTrackSumWos", "Sum of wos of PU tracks for matched vertices; log10(PUSumWos)", 50, -1., 7.);
0835     meSecTrackSumWos_ = ibook.book1D(
0836         "SecTrackSumWos", "Sum of wos of tracks from secondary vtx for matched vertices; log10(SecSumWos)", 50, -1., 7.);
0837     mePUTrackSumPt_ =
0838         ibook.book1D("PUTrackSumPt", "Sum of Pt of PU tracks for matched vertices; log10(PUSumPt)", 50, -2., 3.);
0839     mePUTrackSumPt2_ =
0840         ibook.book1D("PUTrackSumPt2", "Sum of Pt2 of PU tracks for matched vertices; log10(PUSumPt2)", 50, -3., 3.);
0841 
0842     mePUTrackRelMultvsMult_ =
0843         ibook.bookProfile("PUTrackRelMultvsMult",
0844                           "Relative PU multiplicity vs Number of tracks for matched vertices; #Trks; #PUTrks/#Trks",
0845                           50,
0846                           0.,
0847                           200.,
0848                           0.,
0849                           1.,
0850                           "s");
0851     meFakeTrackRelMultvsMult_ = ibook.bookProfile(
0852         "FakeTrackRelMultvsMult",
0853         "Relative multiplicity of fake tracks vs Number of tracks for matched vertices; #Trks; #FakeTrks/#Trks",
0854         50,
0855         0.,
0856         200.,
0857         0.,
0858         1.,
0859         "s");
0860     mePUTrackRelSumWntvsSumWnt_ =
0861         ibook.bookProfile("PUTrackRelSumWntvsSumWnt",
0862                           "Relative PU Sum of Wnt vs Sum of Wnt of tracks for matched vertices; log10(SumW*min(Pt, "
0863                           "1.)); PUSumW*min(Pt, 1.)/SumW*min(Pt, 1.)",
0864                           50,
0865                           0.,
0866                           2.5,
0867                           0.,
0868                           1.,
0869                           "s");
0870     mePUTrackRelSumWosvsSumWos_ = ibook.bookProfile(
0871         "PUTrackRelSumWosvsSumWos",
0872         "Relative PU Sum of Wos vs Sum of Wos of tracks for matched vertices; log10(SumWos); PUSumWos/SumWos",
0873         50,
0874         2.5,
0875         7.,
0876         0.,
0877         1.,
0878         "s");
0879     meSecTrackRelSumWosvsSumWos_ = ibook.bookProfile("SecTrackRelSumWosvsSumWos",
0880                                                      "Relative Sum of Wos of tracks from secondary vtx vs Sum of Wos "
0881                                                      "of tracks for matched vertices; log10(SumWos); SecSumWos/SumWos",
0882                                                      50,
0883                                                      2.,
0884                                                      7.,
0885                                                      0.,
0886                                                      1.,
0887                                                      "s");
0888     meFakeTrackRelSumWosvsSumWos_ = ibook.bookProfile("FakeTrackRelSumWosvsSumWos",
0889                                                       "Relative Sum of Wos of fake tracks vs Sum of Wos of tracks for "
0890                                                       "matched vertices; log10(SumWos); FakeSumWos/SumWos",
0891                                                       50,
0892                                                       2.5,
0893                                                       7.5,
0894                                                       0.,
0895                                                       1.,
0896                                                       "s");
0897     mePUTrackRelSumPtvsSumPt_ = ibook.bookProfile(
0898         "PUTrackRelSumPtvsSumPt",
0899         "Relative PU Sum of Pt vs Sum of Pt of tracks for matched vertices; log10(SumPt); PUSumPt/SumPt",
0900         50,
0901         0.,
0902         3.,
0903         0.,
0904         1.,
0905         "s");
0906     mePUTrackRelSumPt2vsSumPt2_ = ibook.bookProfile(
0907         "PUTrackRelSumPt2vsSumPt2",
0908         "Relative PU Sum of Pt2 vs Sum of Pt2 of tracks for matched vertices; log10(SumPt2); PUSumPt2/SumPt2",
0909         50,
0910         0.,
0911         4.,
0912         0.,
0913         1.,
0914         "s");
0915 
0916     mePUTrackRecLVMult_ = ibook.book1D("PUTrackRecLVMult", "Number of PU tracks for matched LV; #PUTrks", 50, 0., 100.);
0917     mePUTrackRecLVWnt_ =
0918         ibook.book1D("PUTrackRecLVWnt", "Wnt of PU tracks for matched LV; PUTrkW*min(Pt, 1.)", 50, 0., 1.);
0919     mePUTrackRecLVSumWnt_ = ibook.book1D(
0920         "PUTrackRecLVSumWnt", "Sum of wnt of PU tracks for matched LV; log10(PUSumW*min(Pt, 1.))", 50, -2., 3.);
0921     mePUTrackRecLVSumWos_ =
0922         ibook.book1D("PUTrackRecLVSumWos", "Sum of wos of PU tracks for matched LV; log10(PUSumWos)", 50, -1., 7.);
0923     meSecTrackRecLVSumWos_ = ibook.book1D(
0924         "SecTrackRecLVSumWos", "Sum of wos of tracks from secondary vtx for matched LV; log10(SecSumWos)", 50, -1., 7.);
0925     mePUTrackRecLVSumPt_ =
0926         ibook.book1D("PUTrackRecLVSumPt", "Sum of Pt of PU tracks for matched LV; log10(PUSumPt)", 50, -2., 3.);
0927     mePUTrackRecLVSumPt2_ =
0928         ibook.book1D("PUTrackRecLVSumPt2", "Sum of Pt2 of PU tracks for matched LV; log10(PUSumPt2)", 50, -3., 3.);
0929 
0930     mePUTrackRecLVRelMultvsMult_ =
0931         ibook.bookProfile("PUTrackRecLVRelMultvsMult",
0932                           "Relative PU multiplicity vs Number of tracks for matched LV; #Trks; #PUTrks/#Trks",
0933                           50,
0934                           0.,
0935                           200.,
0936                           0.,
0937                           1.,
0938                           "s");
0939     meFakeTrackRecLVRelMultvsMult_ = ibook.bookProfile(
0940         "FakeTrackRecLVRelMultvsMult",
0941         "Relative multiplicity of fake tracks  vs Number of tracks for matched LV; #Trks; #FakeTrks/#Trks",
0942         50,
0943         0.,
0944         200.,
0945         0.,
0946         1.,
0947         "s");
0948     mePUTrackRecLVRelSumWntvsSumWnt_ =
0949         ibook.bookProfile("PUTrackRecLVRelSumWntvsSumWnt",
0950                           "Relative PU Sum of Wnt vs Sum of Wnt of tracks for matched LV; log10(SumW*min(Pt, 1.)); "
0951                           "PUSumW*min(Pt, 1.)/SumW*min(Pt, 1.)",
0952                           50,
0953                           1.,
0954                           2.3,
0955                           0.,
0956                           1.,
0957                           "s");
0958     mePUTrackRecLVRelSumWosvsSumWos_ = ibook.bookProfile(
0959         "PUTrackRecLVRelSumWosvsSumWos",
0960         "Relative PU Sum of Wos vs Sum of Wos of tracks for matched vertices; log10(SumWos); PUSumWos/SumWos",
0961         50,
0962         5.5,
0963         7.,
0964         0.,
0965         1.,
0966         "s");
0967     meSecTrackRecLVRelSumWosvsSumWos_ =
0968         ibook.bookProfile("SecTrackRecLVRelSumWosvsSumWos",
0969                           "Relative Sum of Wos of tracks from secondary vtx vs Sum of Wos of tracks for matched "
0970                           "vertices; log10(SumWos); SecSumWos/SumWos",
0971                           50,
0972                           5.5,
0973                           7.,
0974                           0.,
0975                           1.,
0976                           "s");
0977     meFakeTrackRecLVRelSumWosvsSumWos_ = ibook.bookProfile("FakeTrackRecLVRelSumWosvsSumWos",
0978                                                            "Relative Sum of Wos of fake tracks vs Sum of Wos of tracks "
0979                                                            "for matched vertices; log10(SumWos); FakeSumWos/SumWos",
0980                                                            50,
0981                                                            5.5,
0982                                                            7.,
0983                                                            0.,
0984                                                            1.,
0985                                                            "s");
0986     mePUTrackRecLVRelSumPtvsSumPt_ =
0987         ibook.bookProfile("PUTrackRecLVRelSumPtvsSumPt",
0988                           "Relative PU Sum of Pt vs Sum of Pt of tracks for matched LV; log10(SumPt); PUSumPt/SumPt",
0989                           50,
0990                           1.4,
0991                           3.,
0992                           0.,
0993                           1.,
0994                           "s");
0995     mePUTrackRecLVRelSumPt2vsSumPt2_ =
0996         ibook.bookProfile("PUTrackRecLVRelSumPt2vsSumPt2",
0997                           "Relative PU Sum of Pt2 vs Sum of tracks for matched LV; log10(SumPt2); PUSumPt2/SumPt2",
0998                           50,
0999                           2.,
1000                           4.,
1001                           0.,
1002                           1.,
1003                           "s");
1004   }
1005 
1006   meJetsPURelMult_ =
1007       ibook.book1D("JetsPURelMult",
1008                    "Relative contribution of PU to jet multiplicity for matched vertices; (#Jets-#JetsNoPU)/#Jets",
1009                    50,
1010                    0.,
1011                    1.);
1012   meJetsPURelHt_ =
1013       ibook.book1D("JetsPURelHt",
1014                    "Relative contribution of PU to scalar sum of Et of jets for matched vertices; (Ht-HtNoPU)/HT",
1015                    50,
1016                    0.,
1017                    1.);
1018   meJetsPURelSumPt2_ =
1019       ibook.book1D("JetsPURelSumPt2",
1020                    "Relative contribution of PU to sum of Pt2 of jets for matched vertices; (SumPt2-SumPt2NoPU)/SumPt2",
1021                    50,
1022                    0.,
1023                    1.);
1024   meJetsFakeRelSumPt2_ = ibook.book1D(
1025       "JetsFakeRelSumPt2",
1026       "Relative contribution of fake tracks to sum of Pt2 of jets for matched vertices; (SumPt2-SumPt2NoFake)/SumPt2",
1027       50,
1028       0.,
1029       1.);
1030   meJetsPURelMetPt_ =
1031       ibook.book1D("JetsPURelMetPt",
1032                    "Relative contribution of PU to Missing Transverse Energy for matched vertices; (Met-MetNoPU)/Met",
1033                    50,
1034                    -1.,
1035                    1.);
1036   meJetsPURelSumPz_ =
1037       ibook.book1D("JetsPURelSumPz",
1038                    "Relative contribution of PU to sum of Pz of jets for matched vertices; (SumPz-SumPzNoPU)/SumPz",
1039                    50,
1040                    -1.,
1041                    1.);
1042 
1043   meJetsRecLVPURelMult_ =
1044       ibook.book1D("JetsRecLVPURelMult",
1045                    "Relative contribution of PU to jet multiplicity for matched LV; (#Jets-#JetsNoPU)/#Jets",
1046                    50,
1047                    0.,
1048                    1.);
1049   meJetsRecLVPURelHt_ =
1050       ibook.book1D("JetsRecLVPURelHt",
1051                    "Relative contribution of PU to scalar sum of Et of jets for matched LV; (Ht-HtNoPU)/HT",
1052                    50,
1053                    0.,
1054                    1.);
1055   meJetsRecLVPURelSumPt2_ =
1056       ibook.book1D("JetsRecLVPURelSumPt2",
1057                    "Relative contribution of PU to sum of Pt2 of jets for matched LV; (SumPt2-SumPt2NoPU)/SumPt2",
1058                    50,
1059                    0.,
1060                    1.);
1061   meJetsRecLVFakeRelSumPt2_ = ibook.book1D(
1062       "JetsRecLVFakeRelSumPt2",
1063       "Relative contribution of fake tracks to sum of Pt2 of jets for matched LV; (SumPt2-SumPt2NoFake)/SumPt2",
1064       50,
1065       0.,
1066       1.);
1067   meJetsRecLVPURelMetPt_ =
1068       ibook.book1D("JetsRecLVPURelMetPt",
1069                    "Relative contribution of PU to Missing Transverse Energy for matched LV; (Met-MetNoPU)/Met",
1070                    50,
1071                    -1.,
1072                    1.);
1073   meJetsRecLVPURelSumPz_ =
1074       ibook.book1D("JetsRecLVPURelSumPz",
1075                    "Relative contribution of PU to sum of Pz of jets for matched LV; (SumPz-SumPzNoPU)/SumPz",
1076                    50,
1077                    -1.,
1078                    1.);
1079 
1080   if (optionalPlots_) {
1081     meJetsPUMult_ = ibook.book1D(
1082         "JetsPUMult", "Contribution of PU to jet multiplicity for matched vertices; #Jets-#JetsNoPU", 50, 0., 100.);
1083     meJetsPUHt_ = ibook.book1D("JetsPUHt",
1084                                "Contribution of PU to scalar sum of Et of jets for matched vertices; log10(Ht-HtNoPU)",
1085                                50,
1086                                -2.,
1087                                3.);
1088     meJetsPUSumPt2_ =
1089         ibook.book1D("JetsPUSumPt2",
1090                      "Contribution of PU to sum of Pt2 of jets for matched vertices; log10(sumPt2-SumPt2NoPU)",
1091                      50,
1092                      -3.,
1093                      3.);
1094     meJetsPUMetPt_ =
1095         ibook.book1D("JetsPUMetPt",
1096                      "Contribution of PU to Missing Transverse Energy for matched vertices; log10(Met-MetNoPU)",
1097                      50,
1098                      -3.,
1099                      2.);
1100     meJetsPUSumPz_ =
1101         ibook.book1D("JetsPUSumPz",
1102                      "Contribution of PU to sum of Pz of jets for matched vertices; log10(abs(SumPz-SumPzNoPU))",
1103                      50,
1104                      -4.,
1105                      3.);
1106 
1107     meJetsPURelMultvsMult_ = ibook.bookProfile("JetsPURelMultvsMult",
1108                                                "Relative contribution of PU to jet multiplicity vs number of jets for "
1109                                                "matched vertices; #Jets; (#Jets-#JetsNoPU)/#Jets",
1110                                                50,
1111                                                0.,
1112                                                120.,
1113                                                0.,
1114                                                1.,
1115                                                "s");
1116     meJetsPURelHtvsHt_ = ibook.bookProfile("JetsPURelHtvsHt",
1117                                            "Relative contribution of PU to scalar sum of Et of jets vs scalar sum of "
1118                                            "Et for matched vertices; log10(Ht); (Ht-HtNoPU)/HT",
1119                                            50,
1120                                            0.,
1121                                            3.,
1122                                            0.,
1123                                            1.,
1124                                            "s");
1125     meJetsPURelSumPt2vsSumPt2_ = ibook.bookProfile("JetsPURelSumPt2vsSumPt2",
1126                                                    "Relative contribution of PU to sum of Pt2 of jets vs sum of Pt2 "
1127                                                    "for matched vertices; log10(SumPt2); (SumPt2-SumPt2NoPU)/SumPt2",
1128                                                    50,
1129                                                    -1.,
1130                                                    4.,
1131                                                    0.,
1132                                                    1.,
1133                                                    "s");
1134     meJetsFakeRelSumPt2vsSumPt2_ =
1135         ibook.bookProfile("JetsFakeRelSumPt2vsSumPt2",
1136                           "Relative contribution of fake tracks to sum of Pt2 of jets vs sum of Pt2 for matched "
1137                           "vertices; log10(SumPt2); (SumPt2-SumPt2NoFake)/SumPt2",
1138                           50,
1139                           -1.,
1140                           4.,
1141                           0.,
1142                           1.,
1143                           "s");
1144     meJetsPURelMetPtvsMetPt_ = ibook.bookProfile("JetsPURelMetPtvsMetPt",
1145                                                  "Relative contribution of PU to Missing Transverse Energy vs MET for "
1146                                                  "matched vertices; log10(Met); (Met-MetNoPU)/Met",
1147                                                  50,
1148                                                  -1.,
1149                                                  2.,
1150                                                  -1.,
1151                                                  1.,
1152                                                  "s");
1153     meJetsPURelSumPzvsSumPz_ = ibook.bookProfile("JetsPURelSumPzvsSumPz",
1154                                                  "Relative contribution of PU to sum of Pz of jets vs Sum of Pz for "
1155                                                  "matched vertices; log10(abs SumPz); (SumPz-SumPzNoPU)/SumPz",
1156                                                  50,
1157                                                  -2.,
1158                                                  3.,
1159                                                  -1.,
1160                                                  1.,
1161                                                  "s");
1162 
1163     meJetsRecLVPUMult_ = ibook.book1D(
1164         "JetsRecLVPUMult", "Contribution of PU to jet multiplicity for matched LV; #Jets-#JetsNoPU", 50, 0., 100.);
1165     meJetsRecLVPUHt_ = ibook.book1D(
1166         "JetsRecLVPUHt", "Contribution of PU to scalar sum of Et of jets for matched LV; log10(Ht-HtNoPU)", 50, -2., 3.);
1167     meJetsRecLVPUSumPt2_ =
1168         ibook.book1D("JetsRecLVPUSumPt2",
1169                      "Contribution of PU to sum of Pt2 of jets for matched LV; log10(sumPt2-SumPt2NoPU)",
1170                      50,
1171                      -3.,
1172                      3.);
1173     meJetsRecLVPUMetPt_ =
1174         ibook.book1D("JetsRecLVPUMetPt",
1175                      "Contribution of PU to Missing Transverse Energy for matched LV; log10(Met-MetNoPU)",
1176                      50,
1177                      -3.,
1178                      2.);
1179     meJetsRecLVPUSumPz_ =
1180         ibook.book1D("JetsRecLVPUSumPz",
1181                      "Contribution of PU to sum of Pz of jets for matched LV; log10(abs(SumPz-SumPzNoPU))",
1182                      50,
1183                      -4.,
1184                      3.);
1185 
1186     meJetsRecLVPURelMultvsMult_ = ibook.bookProfile("JetsRecLVPURelMultvsMult",
1187                                                     "Relative contribution of PU to jet multiplicity vs number of jets "
1188                                                     "for matched vertices; #Jets; (#Jets-#JetsNoPU)/#Jets",
1189                                                     50,
1190                                                     0.,
1191                                                     120.,
1192                                                     0.,
1193                                                     1.,
1194                                                     "s");
1195     meJetsRecLVPURelHtvsHt_ = ibook.bookProfile("JetsRecLVPURelHtvsHt",
1196                                                 "Relative contribution of PU to scalar sum of Et of jets vs scalar sum "
1197                                                 "of Et for matched vertices; log10(Ht); (Ht-HtNoPU)/HT",
1198                                                 50,
1199                                                 1.5,
1200                                                 3.,
1201                                                 0.,
1202                                                 1.,
1203                                                 "s");
1204     meJetsRecLVPURelSumPt2vsSumPt2_ =
1205         ibook.bookProfile("JetsRecLVPURelSumPt2vsSumPt2",
1206                           "Relative contribution of PU to sum of Pt2 of jets vs sum of Pt2 for matched vertices; "
1207                           "log10(SumPt2); (SumPt2-SumPt2NoPU)/SumPt2",
1208                           50,
1209                           2.,
1210                           5.,
1211                           0.,
1212                           1.,
1213                           "s");
1214     meJetsRecLVFakeRelSumPt2vsSumPt2_ =
1215         ibook.bookProfile("JetsRecLVFakeRelSumPt2vsSumPt2",
1216                           "Relative contribution of fake tracks to sum of Pt2 of jets vs sum of Pt2 for matched "
1217                           "vertices; log10(SumPt2); (SumPt2-SumPt2NoFake)/SumPt2",
1218                           50,
1219                           2.,
1220                           5.,
1221                           0.,
1222                           1.,
1223                           "s");
1224     meJetsRecLVPURelMetPtvsMetPt_ = ibook.bookProfile("JetsRecLVPURelMetPtvsMetPt",
1225                                                       "Relative contribution of PU to Missing Transverse Energy vs MET "
1226                                                       "for matched vertices; log10(Met); (Met-MetNoPU)/Met",
1227                                                       50,
1228                                                       0.,
1229                                                       2.5,
1230                                                       -1.,
1231                                                       1.,
1232                                                       "s");
1233     meJetsRecLVPURelSumPzvsSumPz_ = ibook.bookProfile("JetsRecLVPURelSumPzvsSumPz",
1234                                                       "Relative contribution of PU to sum of Pz of jets vs Sum of Pz "
1235                                                       "for matched vertices; log10(abs SumPz); (SumPz-SumPzNoPU)/SumPz",
1236                                                       50,
1237                                                       0.5,
1238                                                       3.5,
1239                                                       -1.,
1240                                                       1.,
1241                                                       "s");
1242 
1243     meTrackTimeResCorrectPID_ = ibook.book1D(
1244         "TrackTimeResCorrectPID", "Time residual of tracks with correct PID; t_{rec} - t_{sim} [ns]; ", 100, -5., 5.);
1245 
1246     meTrackTimeResWrongPID_ = ibook.book1D(
1247         "TrackTimeResWrongPID", "Time residual of tracks with wrong PID; t_{rec} - t_{sim} [ns]; ", 100, -5., 5.);
1248     meTrackTimeResNoPID_ = ibook.book1D(
1249         "TrackTimeResNoPID", "Time residual of tracks with no PID; t_{rec} - t_{sim} [ns]; ", 100, -5., 5.);
1250     meTrackTimeResNoPIDtruePi_ = ibook.book1D(
1251         "TrackTimeResNoPIDtruePi", "Time residual of no PID tracks, true Pi; t_{rec} - t_{sim} [ns]; ", 100, -5., 5.);
1252 
1253     meTrackTimeResNoPIDtrueK_ = ibook.book1D(
1254         "TrackTimeResNoPIDtrueK", "Time residual of no PID tracks, true K; t_{rec} - t_{sim} [ns]; ", 100, -5., 5.);
1255     meTrackTimeResNoPIDtrueP_ = ibook.book1D(
1256         "TrackTimeResNoPIDtrueP", "Time residual of no PID tracks, true P; t_{rec} - t_{sim} [ns]; ", 100, -5., 5.);
1257 
1258     meNoPIDTrackTimeResNoPIDType_[0] =
1259         ibook.book1D("NoPIDTrackTimeResNoPIDType1",
1260                      "Time residual of no PID tracks, no PID type 1; t_{rec} - t_{sim} [ns];",
1261                      100,
1262                      -5.,
1263                      5.);
1264 
1265     meNoPIDTrackTimeResNoPIDType_[1] =
1266         ibook.book1D("NoPIDTrackTimeResNoPIDType2",
1267                      "Time residual of no PID tracks, no PID type 2; t_{rec} - t_{sim} [ns];",
1268                      100,
1269                      -5.,
1270                      5.);
1271     meNoPIDTrackTimeResNoPIDType_[2] =
1272         ibook.book1D("NoPIDTrackTimeResNoPIDType3",
1273                      "Time residual of no PID tracks, no PID type 3; t_{rec} - t_{sim} [ns];",
1274                      100,
1275                      -5.,
1276                      5.);
1277 
1278     meTrackTimePullCorrectPID_ =
1279         ibook.book1D("TrackTimePullCorrectPID",
1280                      "Time pull of tracks with correct PID; (t_{rec} - t_{sim})/#sigma_{t rec}; ",
1281                      100,
1282                      -10.,
1283                      10.);
1284 
1285     meTrackTimePullWrongPID_ = ibook.book1D("TrackTimePullWrongPID",
1286                                             "Time pull of tracks with wrong PID; (t_{rec} - t_{sim})/#sigma_{t rec}; ",
1287                                             100,
1288                                             -10.,
1289                                             10.);
1290     meTrackTimePullNoPID_ = ibook.book1D(
1291         "TrackTimePullNoPID", "Time pull of tracks with no PID; (t_{rec} - t_{sim})/#sigma_{t rec}; ", 100, -10., 10.);
1292 
1293     meTrackTimePullNoPIDtruePi_ =
1294         ibook.book1D("TrackTimePullNoPIDtruePi",
1295                      "Time pull of no PID tracks, true Pi; (t_{rec} - t_{sim})/#sigma_{t rec}; ",
1296                      100,
1297                      -10.,
1298                      10.);
1299     meTrackTimePullNoPIDtrueK_ =
1300         ibook.book1D("TrackTimePullNoPIDtrueK",
1301                      "Time pull of no PID tracks, true K; (t_{rec} - t_{sim})/#sigma_{t rec}; ",
1302                      100,
1303                      -10.,
1304                      10.);
1305     meTrackTimePullNoPIDtrueP_ =
1306         ibook.book1D("TrackTimePullNoPIDtrueP",
1307                      "Time pull of no PID tracks, true P; (t_{rec} - t_{sim})/#sigma_{t rec}; ",
1308                      100,
1309                      -10.,
1310                      10.);
1311 
1312     meNoPIDTrackTimePullNoPIDType_[0] =
1313         ibook.book1D("NoPIDTrackTimePullNoPIDType1",
1314                      "Time pull of no PID tracks, no PID type 1; (t_{rec} - t_{sim})/#sigma_{t rec}; ",
1315                      100,
1316                      -10.,
1317                      10.);
1318     meNoPIDTrackTimePullNoPIDType_[1] =
1319         ibook.book1D("NoPIDTrackTimePullNoPIDType2",
1320                      "Time pull of no PID tracks, no PID type 2; (t_{rec} - t_{sim})/#sigma_{t rec}; ",
1321                      100,
1322                      -10.,
1323                      10.);
1324     meNoPIDTrackTimePullNoPIDType_[2] =
1325         ibook.book1D("NoPIDTrackTimePullNoPIDType3",
1326                      "Time pull of no PID tracks, no PID type 3; (t_{rec} - t_{sim})/#sigma_{t rec}; ",
1327                      100,
1328                      -10.,
1329                      10.);
1330 
1331     meTrackTimeSigmaCorrectPID_ = ibook.book1D(
1332         "TrackTimeSigmaCorrectPID", "Time sigma of tracks with correct PID; #sigma_{t0Safe} [ns]; ", 100, 0., 4.);
1333     meTrackTimeSigmaWrongPID_ = ibook.book1D(
1334         "TrackTimeSigmaWrongPID", "Time sigma of tracks with wrong PID; #sigma_{t0Safe} [ns]; ", 100, 0., 4.);
1335     meTrackTimeSigmaNoPID_ =
1336         ibook.book1D("TrackTimeSigmaNoPID", "Time sigma of tracks with no PID; #sigma_{t0Safe} [ns]; ", 100, 0., 4.);
1337     meNoPIDTrackSigmaNoPIDType_[0] = ibook.book1D(
1338         "NoPIDTrackSigmaNoPIDType1", "Time sigma of no PID tracks, no PID type 1; #sigma_{t0Safe} [ns]; ", 100, 0., 4.);
1339     meNoPIDTrackSigmaNoPIDType_[1] = ibook.book1D(
1340         "NoPIDTrackSigmaNoPIDType2", "Time sigma of no PID tracks, no PID type 2; #sigma_{t0Safe} [ns]; ", 100, 0., 4.);
1341     meNoPIDTrackSigmaNoPIDType_[2] = ibook.book1D(
1342         "NoPIDTrackSigmaNoPIDType3", "Time sigma of no PID tracks, no PID type 3; #sigma_{t0Safe} [ns]; ", 100, 0., 4.);
1343     meTrackMVACorrectPID_ =
1344         ibook.book1D("TrackMVACorrectPID", "MVA of tracks with correct PID; MVA score; ", 100, 0., 1.);
1345 
1346     meTrackMVAWrongPID_ = ibook.book1D("TrackMVAWrongPID", "MVA of tracks with wrong PID; MVA score; ", 100, 0., 1.);
1347     meTrackMVANoPID_ = ibook.book1D("TrackMVANoPID", "MVA of tracks with no PID; MVA score; ", 100, 0., 1.);
1348     meNoPIDTrackMVANoPIDType_[0] =
1349         ibook.book1D("NoPIDTrackMVANoPIDType1", "MVA of no PID tracks, no PID type 1; MVA score; ", 100, 0., 1.);
1350     meNoPIDTrackMVANoPIDType_[1] =
1351         ibook.book1D("NoPIDTrackMVANoPIDType2", "MVA of no PID tracks, no PID type 2; MVA score; ", 100, 0., 1.);
1352     meNoPIDTrackMVANoPIDType_[2] =
1353         ibook.book1D("NoPIDTrackMVANoPIDType3", "MVA of no PID tracks, no PID type 3; MVA score; ", 100, 0., 1.);
1354   }
1355 
1356   // some tests
1357   meTrackResLowPTot_ = ibook.book1D(
1358       "TrackResLowP", "t_{rec} - t_{sim} for tracks with p < 2 GeV; t_{rec} - t_{sim} [ns] ", 70, -0.15, 0.15);
1359   meTrackResLowP_[0] =
1360       ibook.book1D("TrackResLowP-LowMVA",
1361                    "t_{rec} - t_{sim} for tracks with MVA < 0.5 and p < 2 GeV; t_{rec} - t_{sim} [ns] ",
1362                    100,
1363                    -1.,
1364                    1.);
1365   meTrackResLowP_[1] =
1366       ibook.book1D("TrackResLowP-MediumMVA",
1367                    "t_{rec} - t_{sim} for tracks with 0.5 < MVA < 0.8 and p < 2 GeV; t_{rec} - t_{sim} [ns] ",
1368                    100,
1369                    -1.,
1370                    1.);
1371   meTrackResLowP_[2] =
1372       ibook.book1D("TrackResLowP-HighMVA",
1373                    "t_{rec} - t_{sim} for tracks with MVA > 0.8 and p < 2 GeV; t_{rec} - t_{sim} [ns] ",
1374                    100,
1375                    -1.,
1376                    1.);
1377   meTrackResHighPTot_ = ibook.book1D(
1378       "TrackResHighP", "t_{rec} - t_{sim} for tracks with p > 2 GeV; t_{rec} - t_{sim} [ns] ", 70, -0.15, 0.15);
1379   meTrackResHighP_[0] =
1380       ibook.book1D("TrackResHighP-LowMVA",
1381                    "t_{rec} - t_{sim} for tracks with MVA < 0.5 and p > 2 GeV; t_{rec} - t_{sim} [ns] ",
1382                    100,
1383                    -1.,
1384                    1.);
1385   meTrackResHighP_[1] =
1386       ibook.book1D("TrackResHighP-MediumMVA",
1387                    "t_{rec} - t_{sim} for tracks with 0.5 < MVA < 0.8 and p > 2 GeV; t_{rec} - t_{sim} [ns] ",
1388                    100,
1389                    -1.,
1390                    1.);
1391   meTrackResHighP_[2] =
1392       ibook.book1D("TrackResHighP-HighMVA",
1393                    "t_{rec} - t_{sim} for tracks with MVA > 0.8 and p > 2 GeV; t_{rec} - t_{sim} [ns] ",
1394                    100,
1395                    -1.,
1396                    1.);
1397   meTrackPullLowPTot_ =
1398       ibook.book1D("TrackPullLowP", "Pull for tracks with p < 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
1399   meTrackPullLowP_[0] = ibook.book1D("TrackPullLowP-LowMVA",
1400                                      "Pull for tracks with MVA < 0.5 and p < 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
1401                                      100,
1402                                      -10.,
1403                                      10.);
1404   meTrackPullLowP_[1] = ibook.book1D("TrackPullLowP-MediumMVA",
1405                                      "Pull for tracks with 0.5 < MVA < 0.8 and p < 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
1406                                      100,
1407                                      -10.,
1408                                      10.);
1409   meTrackPullLowP_[2] = ibook.book1D("TrackPullLowP-HighMVA",
1410                                      "Pull for tracks with MVA > 0.8 and p < 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
1411                                      100,
1412                                      -10.,
1413                                      10.);
1414   meTrackPullHighPTot_ =
1415       ibook.book1D("TrackPullHighP", "Pull for tracks with p > 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
1416   meTrackPullHighP_[0] = ibook.book1D("TrackPullHighP-LowMVA",
1417                                       "Pull for tracks with MVA < 0.5 and p > 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
1418                                       100,
1419                                       -10.,
1420                                       10.);
1421   meTrackPullHighP_[1] =
1422       ibook.book1D("TrackPullHighP-MediumMVA",
1423                    "Pull for tracks with 0.5 < MVA < 0.8 and p > 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
1424                    100,
1425                    -10.,
1426                    10.);
1427   meTrackPullHighP_[2] = ibook.book1D("TrackPullHighP-HighMVA",
1428                                       "Pull for tracks with MVA > 0.8 and p > 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
1429                                       100,
1430                                       -10.,
1431                                       10.);
1432   if (optionalPlots_) {
1433     meTrackResMassProtons_[0] =
1434         ibook.book1D("TrackResMass-Protons-LowMVA",
1435                      "t_{rec} - t_{est} for proton tracks with MVA < 0.5; t_{rec} - t_{est} [ns] ",
1436                      100,
1437                      -1.,
1438                      1.);
1439     meTrackResMassProtons_[1] =
1440         ibook.book1D("TrackResMass-Protons-MediumMVA",
1441                      "t_{rec} - t_{est} for proton tracks with 0.5 < MVA < 0.8; t_{rec} - t_{est} [ns] ",
1442                      100,
1443                      -1.,
1444                      1.);
1445     meTrackResMassProtons_[2] =
1446         ibook.book1D("TrackResMass-Protons-HighMVA",
1447                      "t_{rec} - t_{est} for proton tracks with MVA > 0.8; t_{rec} - t_{est} [ns] ",
1448                      100,
1449                      -1.,
1450                      1.);
1451     meTrackResMassTrueProtons_[0] =
1452         ibook.book1D("TrackResMassTrue-Protons-LowMVA",
1453                      "t_{est} - t_{sim} for proton tracks with MVA < 0.5; t_{est} - t_{sim} [ns] ",
1454                      100,
1455                      -1.,
1456                      1.);
1457     meTrackResMassTrueProtons_[1] =
1458         ibook.book1D("TrackResMassTrue-Protons-MediumMVA",
1459                      "t_{est} - t_{sim} for proton tracks with 0.5 < MVA < 0.8; t_{est} - t_{sim} [ns] ",
1460                      100,
1461                      -1.,
1462                      1.);
1463     meTrackResMassTrueProtons_[2] =
1464         ibook.book1D("TrackResMassTrue-Protons-HighMVA",
1465                      "t_{est} - t_{sim} for proton tracks with MVA > 0.8; t_{est} - t_{sim} [ns] ",
1466                      100,
1467                      -1.,
1468                      1.);
1469 
1470     meTrackResMassPions_[0] = ibook.book1D("TrackResMass-Pions-LowMVA",
1471                                            "t_{rec} - t_{est} for pion tracks with MVA < 0.5; t_{rec} - t_{est} [ns] ",
1472                                            100,
1473                                            -1.,
1474                                            1.);
1475     meTrackResMassPions_[1] =
1476         ibook.book1D("TrackResMass-Pions-MediumMVA",
1477                      "t_{rec} - t_{est} for pion tracks with 0.5 < MVA < 0.8; t_{rec} - t_{est} [ns] ",
1478                      100,
1479                      -1.,
1480                      1.);
1481     meTrackResMassPions_[2] = ibook.book1D("TrackResMass-Pions-HighMVA",
1482                                            "t_{rec} - t_{est} for pion tracks with MVA > 0.8; t_{rec} - t_{est} [ns] ",
1483                                            100,
1484                                            -1.,
1485                                            1.);
1486     meTrackResMassTruePions_[0] =
1487         ibook.book1D("TrackResMassTrue-Pions-LowMVA",
1488                      "t_{est} - t_{sim} for pion tracks with MVA < 0.5; t_{est} - t_{sim} [ns] ",
1489                      100,
1490                      -1.,
1491                      1.);
1492     meTrackResMassTruePions_[1] =
1493         ibook.book1D("TrackResMassTrue-Pions-MediumMVA",
1494                      "t_{est} - t_{sim} for pion tracks with 0.5 < MVA < 0.8; t_{est} - t_{sim} [ns] ",
1495                      100,
1496                      -1.,
1497                      1.);
1498     meTrackResMassTruePions_[2] =
1499         ibook.book1D("TrackResMassTrue-Pions-HighMVA",
1500                      "t_{est} - t_{sim} for pion tracks with MVA > 0.8; t_{est} - t_{sim} [ns] ",
1501                      100,
1502                      -1.,
1503                      1.);
1504   }
1505 
1506   meBarrelPIDp_ = ibook.book1D("BarrelPIDp", "PID track MTD momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1507   meEndcapPIDp_ = ibook.book1D("EndcapPIDp", "PID track with MTD momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1508 
1509   meBarrelNoPIDtype_ = ibook.book1D("BarrelNoPIDtype", "Barrel PID failure category", 4, 0., 4.);
1510   meEndcapNoPIDtype_ = ibook.book1D("EndcapNoPIDtype", "Endcap PID failure category", 4, 0., 4.);
1511 
1512   meBarrelTruePiNoPID_ =
1513       ibook.book1D("BarrelTruePiNoPID", "True pi NoPID momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1514   meBarrelTrueKNoPID_ =
1515       ibook.book1D("BarrelTrueKNoPID", "True k NoPID momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1516   meBarrelTruePNoPID_ =
1517       ibook.book1D("BarrelTruePNoPID", "True p NoPID momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1518   meEndcapTruePiNoPID_ =
1519       ibook.book1D("EndcapTruePiNoPID", "True pi NoPID momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1520   meEndcapTrueKNoPID_ =
1521       ibook.book1D("EndcapTrueKNoPID", "True k NoPID momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1522   meEndcapTruePNoPID_ =
1523       ibook.book1D("EndcapTruePNoPID", "True p NoPID momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1524 
1525   meBarrelTruePiAsPi_ =
1526       ibook.book1D("BarrelTruePiAsPi", "True pi as pi momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1527   meBarrelTruePiAsK_ =
1528       ibook.book1D("BarrelTruePiAsK", "True pi as k momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1529   meBarrelTruePiAsP_ =
1530       ibook.book1D("BarrelTruePiAsP", "True pi as p momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1531   meEndcapTruePiAsPi_ =
1532       ibook.book1D("EndcapTruePiAsPi", "True pi as pi momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1533   meEndcapTruePiAsK_ =
1534       ibook.book1D("EndcapTruePiAsK", "True pi as k momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1535   meEndcapTruePiAsP_ =
1536       ibook.book1D("EndcapTruePiAsP", "True pi as p momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1537 
1538   meBarrelTrueKAsPi_ =
1539       ibook.book1D("BarrelTrueKAsPi", "True k as pi momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1540   meBarrelTrueKAsK_ = ibook.book1D("BarrelTrueKAsK", "True k as k momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1541   meBarrelTrueKAsP_ = ibook.book1D("BarrelTrueKAsP", "True k as p momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1542   meEndcapTrueKAsPi_ =
1543       ibook.book1D("EndcapTrueKAsPi", "True k as pi momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1544   meEndcapTrueKAsK_ = ibook.book1D("EndcapTrueKAsK", "True k as k momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1545   meEndcapTrueKAsP_ = ibook.book1D("EndcapTrueKAsP", "True k as p momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1546 
1547   meBarrelTruePAsPi_ =
1548       ibook.book1D("BarrelTruePAsPi", "True p as pi momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1549   meBarrelTruePAsK_ = ibook.book1D("BarrelTruePAsK", "True p as k momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1550   meBarrelTruePAsP_ = ibook.book1D("BarrelTruePAsP", "True p as p momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1551   meEndcapTruePAsPi_ =
1552       ibook.book1D("EndcapTruePAsPi", "True p as pi momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1553   meEndcapTruePAsK_ = ibook.book1D("EndcapTruePAsK", "True p as k momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1554   meEndcapTruePAsP_ = ibook.book1D("EndcapTruePAsP", "True p as p momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1555 }
1556 
1557 bool Primary4DVertexValidation::matchRecoTrack2SimSignal(const reco::TrackBaseRef& recoTrack) {
1558   auto found = r2s_->find(recoTrack);
1559 
1560   // reco track not matched to any TP
1561   if (found == r2s_->end())
1562     return false;
1563 
1564   //// reco track matched to some TP from signal vertex
1565   for (const auto& tp : found->val) {
1566     if (tp.first->eventId().bunchCrossing() == 0 && tp.first->eventId().event() == 0)
1567       return true;
1568   }
1569 
1570   // reco track not matched to any TP from signal vertex
1571   return false;
1572 }
1573 
1574 std::pair<const edm::Ref<std::vector<TrackingParticle>>*, int> Primary4DVertexValidation::getMatchedTP(
1575     const reco::TrackBaseRef& recoTrack, const TrackingVertexRef& vsim) {
1576   auto found = r2s_->find(recoTrack);
1577 
1578   // reco track not matched to any TP (fake tracks)
1579   if (found == r2s_->end())
1580     return std::make_pair(nullptr, -1);
1581 
1582   // matched TP equal to any TP of a given sim vertex
1583   for (const auto& tp : found->val) {
1584     if (std::find_if(vsim->daughterTracks_begin(), vsim->daughterTracks_end(), [&](const TrackingParticleRef& vtp) {
1585           return tp.first == vtp;
1586         }) != vsim->daughterTracks_end())
1587       return std::make_pair(&tp.first, 0);
1588     // matched TP not associated to any daughter track of a given sim vertex but having the same eventID (track from secondary vtx)
1589     else if (tp.first->eventId().bunchCrossing() == vsim->eventId().bunchCrossing() &&
1590              tp.first->eventId().event() == vsim->eventId().event()) {
1591       return std::make_pair(&tp.first, 1);
1592     }
1593     // matched TP not associated to any sim vertex of a given simulated event (PU track)
1594     else {
1595       return std::make_pair(&tp.first, 2);
1596     }
1597   }
1598 
1599   // reco track not matched to any TP from vertex
1600   return std::make_pair(nullptr, -1);
1601 }
1602 
1603 double Primary4DVertexValidation::timeFromTrueMass(double mass, double pathlength, double momentum, double time) {
1604   if (time > 0 && pathlength > 0 && mass > 0) {
1605     double gammasq = 1. + momentum * momentum / (mass * mass);
1606     double v = c_ * std::sqrt(1. - 1. / gammasq);  // cm / ns
1607     double t_est = time - (pathlength / v);
1608 
1609     return t_est;
1610   } else {
1611     return -1;
1612   }
1613 }
1614 
1615 bool Primary4DVertexValidation::select(const reco::Vertex& v, int level) {
1616   /* level
1617    0  !isFake  && ndof>4  (default)
1618    1  !isFake  && ndof>4 && prob > 0.01
1619    2  !isFake  && ndof>4 && prob > 0.01 && ptmax2 > 0.4
1620    */
1621   if (v.isFake())
1622     return false;
1623   if ((level == 0) && (v.ndof() > selNdof_))
1624     return true;
1625   /*if ((level == 1) && (v.ndof() > selNdof_) && (vertex_pxy(v) > 0.01))
1626     return true;
1627   if ((level == 2) && (v.ndof() > selNdof_) && (vertex_pxy(v) > 0.01) && (vertex_ptmax2(v) > 0.4))
1628     return true;
1629   if ((level == 3) && (v.ndof() > selNdof_) && (vertex_ptmax2(v) < 0.4))
1630     return true;*/
1631   return false;
1632 }
1633 
1634 void Primary4DVertexValidation::observablesFromJets(const std::vector<reco::Track>& reco_Tracks,
1635                                                     const std::vector<double>& mass_Tracks,
1636                                                     const std::vector<int>& category_Tracks,
1637                                                     const std::string& skip_Tracks,
1638                                                     unsigned int& n_Jets,
1639                                                     double& sum_EtJets,
1640                                                     double& sum_Pt2Jets,
1641                                                     double& met_Pt,
1642                                                     double& sum_PzJets) {
1643   double sum_PtJets = 0;
1644   n_Jets = 0;
1645   sum_EtJets = 0;
1646   sum_Pt2Jets = 0;
1647   met_Pt = 0;
1648   sum_PzJets = 0;
1649   auto met = LorentzVector(0, 0, 0, 0);
1650   std::vector<fastjet::PseudoJet> fjInputs_;
1651   fjInputs_.clear();
1652   size_t countScale0 = 0;
1653   for (size_t i = 0; i < reco_Tracks.size(); i++) {
1654     const auto& recotr = reco_Tracks[i];
1655     const auto mass = mass_Tracks[i];
1656     float scale = 1.;
1657     if (recotr.charge() == 0) {
1658       continue;
1659     }
1660     // skip PU tracks in jet definition if skip_PU is required
1661     if (skip_Tracks == "skip_PU" && category_Tracks[i] == 2) {
1662       continue;
1663     }
1664     // skip fake tracks in jet definition if skip_Fake is required
1665     if (skip_Tracks == "skip_Fake" && category_Tracks[i] == -1) {
1666       continue;
1667     }
1668     if (recotr.pt() != 0) {
1669       scale = (recotr.pt() - recotr.ptError()) / recotr.pt();
1670     }
1671     if (edm::isNotFinite(scale)) {
1672       edm::LogWarning("Primary4DVertexValidation") << "Scaling is NAN ignoring this recotrack" << std::endl;
1673       scale = 0;
1674     }
1675     if (scale < 0) {
1676       scale = 0;
1677       countScale0++;
1678     }
1679     if (scale != 0) {
1680       fjInputs_.push_back(fastjet::PseudoJet(recotr.px() * scale,
1681                                              recotr.py() * scale,
1682                                              recotr.pz() * scale,
1683                                              std::sqrt(recotr.p() * recotr.p() + mass * mass) * scale));
1684     }
1685   }
1686   fastjet::ClusterSequence sequence(fjInputs_, fastjet::JetDefinition(fastjet::antikt_algorithm, 0.4));
1687   auto jets = fastjet::sorted_by_pt(sequence.inclusive_jets(0));
1688   for (const auto& pj : jets) {
1689     auto p4 = LorentzVector(pj.px(), pj.py(), pj.pz(), pj.e());
1690     sum_EtJets += std::sqrt(p4.e() * p4.e() - p4.P() * p4.P() + p4.pt() * p4.pt());
1691     sum_PtJets += p4.pt();
1692     sum_Pt2Jets += (p4.pt() * p4.pt() * 0.8 * 0.8);
1693     met += p4;
1694     sum_PzJets += p4.pz();
1695     n_Jets++;
1696   }
1697   met_Pt = met.pt();
1698   double metAbove = met_Pt - 2 * std::sqrt(sum_PtJets);
1699   if (metAbove > 0) {
1700     sum_Pt2Jets += (metAbove * metAbove);
1701   }
1702   if (countScale0 == reco_Tracks.size()) {
1703     sum_Pt2Jets = countScale0 * 0.01;  //leave some epsilon value to sort vertices with unknown pt
1704   }
1705 }
1706 
1707 void Primary4DVertexValidation::isParticle(const reco::TrackBaseRef& recoTrack,
1708                                            const edm::ValueMap<float>& sigmat0,
1709                                            const edm::ValueMap<float>& sigmat0Safe,
1710                                            const edm::ValueMap<float>& probPi,
1711                                            const edm::ValueMap<float>& probK,
1712                                            const edm::ValueMap<float>& probP,
1713                                            unsigned int& no_PIDtype,
1714                                            bool& no_PID,
1715                                            bool& is_Pi,
1716                                            bool& is_K,
1717                                            bool& is_P) {
1718   no_PIDtype = 0;
1719   no_PID = false;
1720   is_Pi = false;
1721   is_K = false;
1722   is_P = false;
1723   if (probPi[recoTrack] == -1) {
1724     no_PIDtype = 1;
1725   } else if (edm::isNotFinite(probPi[recoTrack])) {
1726     no_PIDtype = 2;
1727   } else if (probPi[recoTrack] == 1 && probK[recoTrack] == 0 && probP[recoTrack] == 0 &&
1728              sigmat0[recoTrack] < sigmat0Safe[recoTrack]) {
1729     no_PIDtype = 3;
1730   }
1731   no_PID = no_PIDtype > 0;
1732   is_Pi = !no_PID && 1. - probPi[recoTrack] < minProbHeavy_;
1733   is_K = !no_PID && !is_Pi && probK[recoTrack] > probP[recoTrack];
1734   is_P = !no_PID && !is_Pi && !is_K;
1735 }
1736 
1737 void Primary4DVertexValidation::getWosWnt(const reco::Vertex& recoVtx,
1738                                           const reco::TrackBaseRef& recoTrk,
1739                                           const edm::ValueMap<float>& sigmat0,
1740                                           const edm::ValueMap<float>& mtdQualMVA,
1741                                           const edm::Handle<reco::BeamSpot>& BS,
1742                                           double& wos,
1743                                           double& wnt) {
1744   double dz2_beam = pow((*BS).BeamWidthX() * cos(recoTrk->phi()) / tan(recoTrk->theta()), 2) +
1745                     pow((*BS).BeamWidthY() * sin(recoTrk->phi()) / tan(recoTrk->theta()), 2);
1746   double dz2 =
1747       pow(recoTrk->dzError(), 2) + dz2_beam + pow(0.0020, 2);  // added 20 um, some tracks have crazy small resolutions
1748   wos = recoVtx.trackWeight(recoTrk) / dz2;
1749   wnt = recoVtx.trackWeight(recoTrk) *
1750         std::min(recoTrk->pt(), 1.0);  // pt-weighted number of tracks (downweights pt < 1 GeV tracks)
1751 
1752   // If tracks have time information, give more weight to tracks with good resolution
1753   if (sigmat0[recoTrk] > 0 && mtdQualMVA[recoTrk] > mvaTh_) {
1754     double sigmaZ = (*BS).sigmaZ();
1755     double sigmaT = sigmaZ / c_;  // c in cm/ns
1756     wos = wos / erf(sigmat0[recoTrk] / sigmaT);
1757   }
1758 }
1759 
1760 /* Extract information form TrackingParticles/TrackingVertex and fill
1761  * the helper class simPrimaryVertex with proper generation-level
1762  * information */
1763 std::vector<Primary4DVertexValidation::simPrimaryVertex> Primary4DVertexValidation::getSimPVs(
1764     const edm::Handle<TrackingVertexCollection>& tVC) {
1765   std::vector<Primary4DVertexValidation::simPrimaryVertex> simpv;
1766   int current_event = -1;
1767   int s = -1;
1768   for (TrackingVertexCollection::const_iterator v = tVC->begin(); v != tVC->end(); ++v) {
1769     // We keep only the first vertex from all the events at BX=0.
1770     if (v->eventId().bunchCrossing() != 0)
1771       continue;
1772     if (v->eventId().event() != current_event) {
1773       current_event = v->eventId().event();
1774     } else {
1775       continue;
1776     }
1777     s++;
1778     if (std::abs(v->position().z()) > 1000)
1779       continue;  // skip junk vertices
1780 
1781     // could be a new vertex, check  all primaries found so far to avoid multiple entries
1782     simPrimaryVertex sv(v->position().x(), v->position().y(), v->position().z(), v->position().t());
1783     sv.eventId = v->eventId();
1784     sv.sim_vertex = TrackingVertexRef(tVC, std::distance(tVC->begin(), v));
1785     sv.OriginalIndex = s;
1786 
1787     for (TrackingParticleRefVector::iterator iTrack = v->daughterTracks_begin(); iTrack != v->daughterTracks_end();
1788          ++iTrack) {
1789       assert((**iTrack).eventId().bunchCrossing() == 0);
1790     }
1791     simPrimaryVertex* vp = nullptr;  // will become non-NULL if a vertex is found and then point to it
1792     for (std::vector<simPrimaryVertex>::iterator v0 = simpv.begin(); v0 != simpv.end(); v0++) {
1793       if ((sv.eventId == v0->eventId) && (std::abs(sv.x - v0->x) < 1e-5) && (std::abs(sv.y - v0->y) < 1e-5) &&
1794           (std::abs(sv.z - v0->z) < 1e-5)) {
1795         vp = &(*v0);
1796         break;
1797       }
1798     }
1799     if (!vp) {
1800       // this is a new vertex, add it to the list of sim-vertices
1801       simpv.push_back(sv);
1802       vp = &simpv.back();
1803     }
1804 
1805     // Loop over daughter track(s) as Tracking Particles
1806     for (TrackingVertex::tp_iterator iTP = v->daughterTracks_begin(); iTP != v->daughterTracks_end(); ++iTP) {
1807       auto momentum = (*(*iTP)).momentum();
1808       const reco::Track* matched_best_reco_track = nullptr;
1809       double match_quality = -1;
1810       if (use_only_charged_tracks_ && (**iTP).charge() == 0)
1811         continue;
1812       if (s2r_->find(*iTP) != s2r_->end()) {
1813         matched_best_reco_track = (*s2r_)[*iTP][0].first.get();
1814         match_quality = (*s2r_)[*iTP][0].second;
1815       }
1816 
1817       vp->ptot.setPx(vp->ptot.x() + momentum.x());
1818       vp->ptot.setPy(vp->ptot.y() + momentum.y());
1819       vp->ptot.setPz(vp->ptot.z() + momentum.z());
1820       vp->ptot.setE(vp->ptot.e() + (**iTP).energy());
1821       vp->pt += (**iTP).pt();
1822       vp->ptsq += ((**iTP).pt() * (**iTP).pt());
1823       vp->nGenTrk++;
1824 
1825       if (matched_best_reco_track) {
1826         vp->num_matched_reco_tracks++;
1827         vp->average_match_quality += match_quality;
1828       }
1829     }  // End of for loop on daughters sim-particles
1830     if (vp->num_matched_reco_tracks) {
1831       vp->average_match_quality /= static_cast<float>(vp->num_matched_reco_tracks);
1832     }
1833     LogTrace("Primary4DVertexValidation")
1834         << "average number of associated tracks: " << vp->num_matched_reco_tracks / static_cast<float>(vp->nGenTrk)
1835         << " with average quality: " << vp->average_match_quality;
1836   }  // End of for loop on tracking vertices
1837 
1838   // In case of no simulated vertices, break here
1839   if (simpv.empty())
1840     return simpv;
1841 
1842   // Now compute the closest distance in z between all simulated vertex
1843   // first initialize
1844   auto prev_z = simpv.back().z;
1845   for (simPrimaryVertex& vsim : simpv) {
1846     vsim.closest_vertex_distance_z = std::abs(vsim.z - prev_z);
1847     prev_z = vsim.z;
1848   }
1849   // then calculate
1850   for (std::vector<simPrimaryVertex>::iterator vsim = simpv.begin(); vsim != simpv.end(); vsim++) {
1851     std::vector<simPrimaryVertex>::iterator vsim2 = vsim;
1852     vsim2++;
1853     for (; vsim2 != simpv.end(); vsim2++) {
1854       double distance = std::abs(vsim->z - vsim2->z);
1855       // need both to be complete
1856       vsim->closest_vertex_distance_z = std::min(vsim->closest_vertex_distance_z, distance);
1857       vsim2->closest_vertex_distance_z = std::min(vsim2->closest_vertex_distance_z, distance);
1858     }
1859   }
1860   return simpv;
1861 }
1862 
1863 /* Extract information form reco Vertex and fill the helper class
1864  * recoPrimaryVertex with proper reco-level information */
1865 std::vector<Primary4DVertexValidation::recoPrimaryVertex> Primary4DVertexValidation::getRecoPVs(
1866     const edm::Handle<edm::View<reco::Vertex>>& tVC) {
1867   std::vector<Primary4DVertexValidation::recoPrimaryVertex> recopv;
1868   int r = -1;
1869   for (auto v = tVC->begin(); v != tVC->end(); ++v) {
1870     r++;
1871     // Skip junk vertices
1872     if (std::abs(v->z()) > 1000)
1873       continue;
1874     if (v->isFake() || !v->isValid())
1875       continue;
1876 
1877     recoPrimaryVertex sv(v->position().x(), v->position().y(), v->position().z());
1878     sv.recVtx = &(*v);
1879     sv.recVtxRef = reco::VertexBaseRef(tVC, std::distance(tVC->begin(), v));
1880 
1881     sv.OriginalIndex = r;
1882     sv.ndof = v->ndof();
1883     // this is a new vertex, add it to the list of reco-vertices
1884     recopv.push_back(sv);
1885     Primary4DVertexValidation::recoPrimaryVertex* vp = &recopv.back();
1886 
1887     // Loop over daughter track(s)
1888     for (auto iTrack = v->tracks_begin(); iTrack != v->tracks_end(); ++iTrack) {
1889       auto momentum = (*(*iTrack)).innerMomentum();
1890       if (momentum.mag2() == 0)
1891         momentum = (*(*iTrack)).momentum();
1892       vp->pt += std::sqrt(momentum.perp2());
1893       vp->ptsq += (momentum.perp2());
1894       vp->nRecoTrk++;
1895 
1896       auto matched = r2s_->find(*iTrack);
1897       if (matched != r2s_->end()) {
1898         vp->num_matched_sim_tracks++;
1899       }
1900 
1901     }  // End of for loop on daughters reconstructed tracks
1902   }  // End of for loop on tracking vertices
1903 
1904   // In case of no reco vertices, break here
1905   if (recopv.empty())
1906     return recopv;
1907 
1908   // Now compute the closest distance in z between all reconstructed vertex
1909   // first initialize
1910   auto prev_z = recopv.back().z;
1911   for (recoPrimaryVertex& vreco : recopv) {
1912     vreco.closest_vertex_distance_z = std::abs(vreco.z - prev_z);
1913     prev_z = vreco.z;
1914   }
1915   for (std::vector<recoPrimaryVertex>::iterator vreco = recopv.begin(); vreco != recopv.end(); vreco++) {
1916     std::vector<recoPrimaryVertex>::iterator vreco2 = vreco;
1917     vreco2++;
1918     for (; vreco2 != recopv.end(); vreco2++) {
1919       double distance = std::abs(vreco->z - vreco2->z);
1920       // need both to be complete
1921       vreco->closest_vertex_distance_z = std::min(vreco->closest_vertex_distance_z, distance);
1922       vreco2->closest_vertex_distance_z = std::min(vreco2->closest_vertex_distance_z, distance);
1923     }
1924   }
1925   return recopv;
1926 }
1927 
1928 // ------------ method called to produce the data  ------------
1929 void Primary4DVertexValidation::matchReco2Sim(std::vector<recoPrimaryVertex>& recopv,
1930                                               std::vector<simPrimaryVertex>& simpv,
1931                                               const edm::ValueMap<float>& sigmat0,
1932                                               const edm::ValueMap<float>& MVA,
1933                                               const edm::Handle<reco::BeamSpot>& BS) {
1934   // Initialization (clear wos and wnt)
1935   for (auto vv : simpv) {
1936     vv.wnt.clear();
1937     vv.wos.clear();
1938   }
1939   for (auto rv : recopv) {
1940     rv.wnt.clear();
1941     rv.wos.clear();
1942   }
1943 
1944   // Filling infos for matching rec and sim vertices
1945   for (unsigned int iv = 0; iv < recopv.size(); iv++) {
1946     const reco::Vertex* vertex = recopv.at(iv).recVtx;
1947     LogTrace("Primary4DVertexValidation") << "iv (rec): " << iv;
1948 
1949     for (unsigned int iev = 0; iev < simpv.size(); iev++) {
1950       double wnt = 0;
1951       double wos = 0;
1952       double evwnt = 0;       // weighted number of tracks from sim event iev in the current recvtx
1953       double evwos = 0;       // weight over sigma**2 of sim event iev in the current recvtx
1954       unsigned int evnt = 0;  // number of tracks from sim event iev in the current recvtx
1955 
1956       for (auto iTrack = vertex->tracks_begin(); iTrack != vertex->tracks_end(); ++iTrack) {
1957         if (vertex->trackWeight(*iTrack) < trackweightTh_) {
1958           continue;
1959         }
1960 
1961         auto tp_info = getMatchedTP(*iTrack, simpv.at(iev).sim_vertex).first;
1962         int matchCategory = getMatchedTP(*iTrack, simpv.at(iev).sim_vertex).second;
1963         // matched TP equal to any TP of a given sim vertex
1964         if (tp_info != nullptr && matchCategory == 0) {
1965           getWosWnt(*vertex, *iTrack, MVA, sigmat0, BS, wos, wnt);
1966           simpv.at(iev).addTrack(iv, wos, wnt);
1967           recopv.at(iv).addTrack(iev, wos, wnt);
1968           evwos += wos;
1969           evwnt += wnt;
1970           evnt++;
1971         }
1972       }  // RecoTracks loop
1973 
1974       // require 2 tracks for a wos-match
1975       if ((evwos > 0) && (evwos > recopv.at(iv).maxwos) && (evnt > 1)) {
1976         recopv.at(iv).wosmatch = iev;
1977         recopv.at(iv).maxwos = evwos;
1978         recopv.at(iv).maxwosnt = evnt;
1979         LogTrace("Primary4DVertexValidation") << "dominating sim event (iev): " << iev << " evwos = " << evwos;
1980       }
1981 
1982       // weighted track counting match, require at least one track
1983       if ((evnt > 0) && (evwnt > recopv.at(iv).maxwnt)) {
1984         recopv.at(iv).wntmatch = iev;
1985         recopv.at(iv).maxwnt = evwnt;
1986       }
1987     }  // TrackingVertex loop
1988     if (recopv.at(iv).maxwos > 0) {
1989       simpv.at(recopv.at(iv).wosmatch).wos_dominated_recv.push_back(iv);
1990       simpv.at(recopv.at(iv).wosmatch).nwosmatch++;  // count the rec vertices dominated by a sim vertex using wos
1991       assert(iv < recopv.size());
1992     }
1993     LogTrace("Primary4DVertexValidation") << "largest contribution to wos: wosmatch (iev) = " << recopv.at(iv).wosmatch
1994                                           << " maxwos = " << recopv.at(iv).maxwos;
1995     if (recopv.at(iv).maxwnt > 0) {
1996       simpv.at(recopv.at(iv).wntmatch).nwntmatch++;  // count the rec vertices dominated by a sim vertex using wnt
1997     }
1998   }  // RecoPrimaryVertex
1999 
2000   // reset
2001   for (auto& vrec : recopv) {
2002     vrec.sim = NOT_MATCHED;
2003     vrec.matchQuality = 0;
2004   }
2005   unsigned int iev = 0;
2006   for (auto& vv : simpv) {
2007     LogTrace("Primary4DVertexValidation") << "iev (sim): " << iev;
2008     LogTrace("Primary4DVertexValidation") << "wos_dominated_recv.size: " << vv.wos_dominated_recv.size();
2009     for (unsigned int i = 0; i < vv.wos_dominated_recv.size(); i++) {
2010       auto recov = vv.wos_dominated_recv.at(i);
2011       LogTrace("Primary4DVertexValidation")
2012           << "index of reco vertex: " << recov << " that has a wos: " << vv.wos.at(recov) << " at position " << i;
2013     }
2014     vv.rec = NOT_MATCHED;
2015     vv.matchQuality = 0;
2016     iev++;
2017   }
2018   // after filling infos, goes for the sim-reco match
2019   // this tries a one-to-one match, taking simPV with highest wos if there are > 1 simPV candidates
2020   for (unsigned int rank = 1; rank < maxRank_; rank++) {
2021     for (unsigned int iev = 0; iev < simpv.size(); iev++) {  //loop on SimPV
2022       if (simpv.at(iev).rec != NOT_MATCHED) {
2023         continue;  // only sim vertices not already matched
2024       }
2025       if (simpv.at(iev).nwosmatch == 0) {
2026         continue;  // the sim vertex does not dominate any reco vertex
2027       }
2028       if (simpv.at(iev).nwosmatch > rank) {
2029         continue;  // start with sim vertices dominating one rec vertex (rank 1), then go with the ones dominating more
2030       }
2031       unsigned int iv = NOT_MATCHED;  // select a rec vertex index
2032       for (unsigned int k = 0; k < simpv.at(iev).wos_dominated_recv.size(); k++) {
2033         unsigned int rec = simpv.at(iev).wos_dominated_recv.at(k);  //candidate rec vertex index
2034         auto vrec = recopv.at(rec);
2035         if (vrec.sim != NOT_MATCHED) {
2036           continue;  // already matched
2037         }
2038         if (std::abs(simpv.at(iev).z - vrec.z) > zWosMatchMax_) {
2039           continue;  // insanely far away
2040         }
2041         if ((iv == NOT_MATCHED) || simpv.at(iev).wos.at(rec) > simpv.at(iev).wos.at(iv)) {
2042           iv = rec;
2043         }
2044       }
2045       // if a viable candidate was found, make the link
2046       if (iv !=
2047           NOT_MATCHED) {  // if the rec vertex has already been associated is possible that iv remains NOT_MATCHED at this point
2048         recopv.at(iv).sim = iev;
2049         simpv.at(iev).rec = iv;
2050         recopv.at(iv).matchQuality = rank;
2051         simpv.at(iev).matchQuality = rank;
2052       }
2053     }  // iev
2054   }  // rank
2055 
2056   // Reco vertices that are not necessarily dominated by a sim vertex, or whose dominating sim-vertex
2057   // has been matched already to another overlapping reco vertex, can still be associated to a specific
2058   // sim vertex (without being classified as dominating).
2059   // In terms of fitting 1d-distributions, this corresponds to a small peak close to a much larger nearby peak
2060   unsigned int ntry = 0;
2061   while (ntry++ < maxTry_) {
2062     unsigned nmatch = 0;
2063     for (unsigned int iev = 0; iev < simpv.size(); iev++) {
2064       if ((simpv.at(iev).rec != NOT_MATCHED) || (simpv.at(iev).wos.empty())) {
2065         continue;
2066       }
2067       // find a rec vertex for the NOT_MATCHED sim vertex
2068       unsigned int rec = NOT_MATCHED;
2069       for (auto rv : simpv.at(iev).wos) {
2070         if ((rec == NOT_MATCHED) || (rv.second > simpv.at(iev).wos.at(rec))) {
2071           rec = rv.first;
2072         }
2073       }
2074 
2075       if (rec == NOT_MATCHED) {  // try with wnt match
2076         for (auto rv : simpv.at(iev).wnt) {
2077           if ((rec == NOT_MATCHED) || (rv.second > simpv.at(iev).wnt.at(rec))) {
2078             rec = rv.first;
2079           }
2080         }
2081       }
2082 
2083       if (rec == NOT_MATCHED) {
2084         continue;  // should not happen
2085       }
2086       if (recopv.at(rec).sim != NOT_MATCHED) {
2087         continue;  // already gone
2088       }
2089 
2090       // check if the rec vertex can be matched
2091       unsigned int rec2sim = NOT_MATCHED;
2092       for (auto sv : recopv.at(rec).wos) {
2093         if (simpv.at(sv.first).rec != NOT_MATCHED) {
2094           continue;  // already used
2095         }
2096         if ((rec2sim == NOT_MATCHED) || (sv.second > recopv.at(rec).wos.at(rec2sim))) {
2097           rec2sim = sv.first;
2098         }
2099       }
2100       if (iev == rec2sim) {
2101         // do the match and assign lowest quality (i.e. max rank)
2102         recopv.at(rec).sim = iev;
2103         recopv.at(rec).matchQuality = maxRank_;
2104         simpv.at(iev).rec = rec;
2105         simpv.at(iev).matchQuality = maxRank_;
2106         nmatch++;
2107       }
2108     }  // sim loop
2109     if (nmatch == 0) {
2110       break;
2111     }
2112   }  // ntry
2113 
2114 // Debugging
2115 #ifdef EDM_ML_DEBUG
2116   unsigned int nmatch_tot = 0, n_dzgtsz = 0;
2117   unsigned int n_rank1 = 0, n_rank2 = 0, n_rank3 = 0, n_rank8 = 0;
2118 
2119   for (unsigned int iev = 0; iev < simpv.size(); iev++) {  //loop on SimPV
2120     if (simpv.at(iev).rec != NOT_MATCHED) {
2121       unsigned int rec = simpv.at(iev).rec;
2122       unsigned int wosmatch = recopv.at(rec).wosmatch;
2123       LogTrace("Primary4DVertexValidation")
2124           << "Final match: iev (sim) = " << std::setw(4) << iev << "  sim.rec = " << std::setw(4) << rec
2125           << "  rec.wosmatch = " << std::setw(5) << wosmatch << "  dZ/sigmaZ = " << std::setw(6) << std::setprecision(2)
2126           << std::abs((recopv.at(rec).z - simpv.at(iev).z) / recopv.at(rec).recVtx->zError())
2127           << "  match qual = " << std::setw(1) << recopv.at(rec).matchQuality;
2128       nmatch_tot++;
2129       if (std::abs((recopv.at(rec).z - simpv.at(iev).z) / recopv.at(rec).recVtx->zError()) > 3.) {
2130         n_dzgtsz++;
2131       }
2132       if (recopv.at(rec).matchQuality == 1) {
2133         n_rank1++;
2134       }
2135       if (recopv.at(rec).matchQuality == 2) {
2136         n_rank2++;
2137       }
2138       if (recopv.at(rec).matchQuality == 3) {
2139         n_rank3++;
2140       }
2141       if (recopv.at(rec).matchQuality == 8) {
2142         n_rank8++;
2143       }
2144     }
2145   }
2146   LogTrace("Primary4DVertexValidation") << "n_sim = " << simpv.size() << " n_rec = " << recopv.size()
2147                                         << " nmatch_tot = " << nmatch_tot << " n(dZ>sigmaZ) = " << n_dzgtsz
2148                                         << " n_rank1 = " << n_rank1 << " n_rank2 = " << n_rank2
2149                                         << " n_rank3 = " << n_rank3 << " n_rank8 = " << n_rank8;
2150 #endif
2151 }
2152 
2153 void Primary4DVertexValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
2154   using edm::Handle;
2155   using edm::View;
2156   using std::cout;
2157   using std::endl;
2158   using std::vector;
2159   using namespace reco;
2160 
2161   std::vector<float> pileUpInfo_z;
2162 
2163   // get the pileup information
2164   edm::Handle<std::vector<PileupSummaryInfo>> puinfoH;
2165   if (iEvent.getByToken(vecPileupSummaryInfoToken_, puinfoH)) {
2166     for (auto const& pu_info : *puinfoH.product()) {
2167       if (pu_info.getBunchCrossing() == 0) {
2168         pileUpInfo_z = pu_info.getPU_zpositions();
2169         break;
2170       }
2171     }
2172   }
2173 
2174   edm::Handle<TrackingParticleCollection> TPCollectionH;
2175   iEvent.getByToken(trackingParticleCollectionToken_, TPCollectionH);
2176   if (!TPCollectionH.isValid())
2177     edm::LogWarning("Primary4DVertexValidation") << "TPCollectionH is not valid";
2178 
2179   edm::Handle<TrackingVertexCollection> TVCollectionH;
2180   iEvent.getByToken(trackingVertexCollectionToken_, TVCollectionH);
2181   if (!TVCollectionH.isValid())
2182     edm::LogWarning("Primary4DVertexValidation") << "TVCollectionH is not valid";
2183 
2184   edm::Handle<reco::SimToRecoCollection> simToRecoH;
2185   iEvent.getByToken(simToRecoAssociationToken_, simToRecoH);
2186   if (simToRecoH.isValid())
2187     s2r_ = simToRecoH.product();
2188   else
2189     edm::LogWarning("Primary4DVertexValidation") << "simToRecoH is not valid";
2190 
2191   edm::Handle<reco::RecoToSimCollection> recoToSimH;
2192   iEvent.getByToken(recoToSimAssociationToken_, recoToSimH);
2193   if (recoToSimH.isValid())
2194     r2s_ = recoToSimH.product();
2195   else
2196     edm::LogWarning("Primary4DVertexValidation") << "recoToSimH is not valid";
2197 
2198   reco::BeamSpot beamSpot;
2199   edm::Handle<reco::BeamSpot> BeamSpotH;
2200   iEvent.getByToken(RecBeamSpotToken_, BeamSpotH);
2201   if (!BeamSpotH.isValid())
2202     edm::LogWarning("Primary4DVertexValidation") << "BeamSpotH is not valid";
2203   beamSpot = *BeamSpotH;
2204 
2205   edm::Handle<reco::TrackCollection> tks;
2206   iEvent.getByToken(trkToken, tks);
2207   const auto& theB = &iSetup.getData(theTTBToken);
2208   std::vector<reco::TransientTrack> t_tks;
2209 
2210   std::vector<simPrimaryVertex> simpv;  // a list of simulated primary MC vertices
2211   simpv = getSimPVs(TVCollectionH);
2212   // this bool check if first vertex in that with highest pT
2213   bool signal_is_highest_pt =
2214       std::max_element(simpv.begin(), simpv.end(), [](const simPrimaryVertex& lhs, const simPrimaryVertex& rhs) {
2215         return lhs.ptsq < rhs.ptsq;
2216       }) == simpv.begin();
2217 
2218   std::vector<recoPrimaryVertex> recopv;  // a list of reconstructed primary MC vertices
2219   edm::Handle<edm::View<reco::Vertex>> recVtxs;
2220   iEvent.getByToken(Rec4DVerToken_, recVtxs);
2221   if (!recVtxs.isValid())
2222     edm::LogWarning("Primary4DVertexValidation") << "recVtxs is not valid";
2223   recopv = getRecoPVs(recVtxs);
2224 
2225   const auto& trackAssoc = iEvent.get(trackAssocToken_);
2226   const auto& pathLength = iEvent.get(pathLengthToken_);
2227   const auto& momentum = iEvent.get(momentumToken_);
2228   const auto& time = iEvent.get(timeToken_);
2229   const auto& t0Pid = iEvent.get(t0PidToken_);
2230   const auto& sigmat0 = iEvent.get(sigmat0PidToken_);
2231   const auto& t0Safe = iEvent.get(t0SafePidToken_);
2232   const auto& sigmat0Safe = iEvent.get(sigmat0SafePidToken_);
2233   const auto& mtdQualMVA = iEvent.get(trackMVAQualToken_);
2234   const auto& tMtd = iEvent.get(tmtdToken_);
2235   const auto& tofPi = iEvent.get(tofPiToken_);
2236   const auto& tofK = iEvent.get(tofKToken_);
2237   const auto& tofP = iEvent.get(tofPToken_);
2238   const auto& probPi = iEvent.get(probPiToken_);
2239   const auto& probK = iEvent.get(probKToken_);
2240   const auto& probP = iEvent.get(probPToken_);
2241   const auto& fPDGTable = iSetup.getHandle(pdtToken_);
2242 
2243   // I have simPV and recoPV collections
2244   matchReco2Sim(recopv, simpv, sigmat0Safe, mtdQualMVA, BeamSpotH);
2245 
2246   t_tks = (*theB).build(tks, beamSpot, t0Safe, sigmat0Safe);
2247 
2248   // track filter
2249   std::vector<reco::TransientTrack>&& seltks = theTrackFilter->select(t_tks);
2250 
2251   int unassociatedCount = 0;
2252   int unassociatedCountFake = 0;
2253   for (std::vector<reco::TransientTrack>::const_iterator itk = seltks.begin(); itk != seltks.end(); itk++) {
2254     reco::TrackBaseRef trackref = (*itk).trackBaseRef();
2255     bool isAssociated = false;
2256     for (unsigned int iv = 0; iv < recopv.size(); iv++) {
2257       const reco::Vertex* vertex = recopv.at(iv).recVtx;
2258       for (auto iTrack = vertex->tracks_begin(); iTrack != vertex->tracks_end(); ++iTrack) {
2259         if (*iTrack == trackref) {
2260           isAssociated = true;
2261           break;
2262         }
2263       }
2264       if (isAssociated)
2265         break;
2266     }
2267 
2268     if (!isAssociated) {
2269       unassociatedCount++;
2270       auto found = r2s_->find(trackref);
2271       if (found == r2s_->end())
2272         unassociatedCountFake++;
2273     }
2274   }
2275   double fraction = double(unassociatedCount) / (seltks.size());
2276   meUnAssocTracks_->Fill(log10(unassociatedCount));
2277   meFractionUnAssocTracks_->Fill(fraction);
2278 
2279   double fractionFake = double(unassociatedCountFake) / (seltks.size());
2280   meUnAssocTracksFake_->Fill(log10(unassociatedCountFake));
2281   meFractionUnAssocTracksFake_->Fill(fractionFake);
2282 
2283   // Loop on tracks
2284   for (unsigned int iv = 0; iv < recopv.size(); iv++) {
2285     if (recopv.at(iv).ndof > selNdof_) {
2286       const reco::Vertex* vertex = recopv.at(iv).recVtx;
2287 
2288       for (unsigned int iev = 0; iev < simpv.size(); iev++) {
2289         auto vsim = simpv.at(iev).sim_vertex;
2290 
2291         bool selectedVtxMatching = recopv.at(iv).sim == iev && simpv.at(iev).rec == iv;
2292         bool selectedLV = simpv.at(iev).eventId.bunchCrossing() == 0 && simpv.at(iev).eventId.event() == 0 &&
2293                           recopv.at(iv).OriginalIndex == 0;
2294         bool selectedLVMatching = selectedVtxMatching && selectedLV;  // bool for reco vtx leading match
2295         if (selectedLVMatching && !recopv.at(iv).is_signal()) {
2296           edm::LogWarning("Primary4DVertexValidation")
2297               << "Reco vtx leading match inconsistent: BX/ID " << simpv.at(iev).eventId.bunchCrossing() << " "
2298               << simpv.at(iev).eventId.event();
2299         }
2300 #ifdef EDM_ML_DEBUG
2301         if (selectedLVMatching) {
2302           printSimVtxRecoVtxInfo(simpv.at(iev), recopv.at(iv));
2303         }
2304 #endif
2305         double vzsim = simpv.at(iev).z;
2306         double vtsim = simpv.at(iev).t * simUnit_;
2307 
2308         double wnt = 0, wos = 0;
2309         double PUsumWnt = 0, PUsumWos = 0, SecsumWos = 0, FakesumWos = 0, PUsumPt = 0, PUsumPt2 = 0;
2310         double sumWnt = 0, sumWos = 0, sumPt = 0, sumPt2 = 0;
2311         unsigned int nt = 0, PUnt = 0, Fakent = 0;
2312 
2313         std::vector<double> massVector;
2314         std::vector<reco::Track> recotracks;
2315         std::vector<int> categoryVector;
2316         double sumEtJets = 0, sumPt2Jets = 0, metPt = 0, sumPzJets = 0;
2317         double sumEtJetsnoPU = 0, sumPt2JetsnoPU = 0, metPtnoPU = 0, sumPzJetsnoPU = 0;
2318         double sumEtJetsnoFake = 0, sumPt2JetsnoFake = 0, metPtnoFake = 0, sumPzJetsnoFake = 0;
2319         unsigned int nJets = 0, nJetsnoPU = 0, nJetsnoFake = 0;
2320         for (auto iTrack = vertex->tracks_begin(); iTrack != vertex->tracks_end(); ++iTrack) {
2321           if (trackAssoc[*iTrack] == -1) {
2322             edm::LogWarning("mtdTracks") << "Extended track not associated";
2323             continue;
2324           }
2325 
2326           // monitor all track weights associated to a vertex before selection on it
2327           if (selectedVtxMatching) {
2328             meVtxTrackW_->Fill(vertex->trackWeight(*iTrack));
2329             if (selectedLV) {
2330               meVtxTrackRecLVW_->Fill(vertex->trackWeight(*iTrack));
2331             }
2332           }
2333 
2334           if (vertex->trackWeight(*iTrack) < trackweightTh_)
2335             continue;
2336           bool noCrack = std::abs((*iTrack)->eta()) < trackMaxBtlEta_ || std::abs((*iTrack)->eta()) > trackMinEtlEta_;
2337 
2338           bool selectRecoTrk = trkRecSel(**iTrack);
2339           if (selectedLVMatching && selectRecoTrk) {
2340             if (noCrack) {
2341               meTrackEffPtTot_->Fill((*iTrack)->pt());
2342             }
2343             meTrackEffEtaTot_->Fill(std::abs((*iTrack)->eta()));
2344           }
2345 
2346           auto tp_info = getMatchedTP(*iTrack, vsim).first;
2347           int matchCategory = getMatchedTP(*iTrack, vsim).second;
2348 
2349           // PU, fake and secondary tracks
2350           if (selectedVtxMatching) {
2351             unsigned int no_PIDtype = 0;
2352             bool no_PID, is_Pi, is_K, is_P;
2353             int PartID = 211;  // pion
2354             isParticle(*iTrack, sigmat0, sigmat0Safe, probPi, probK, probP, no_PIDtype, no_PID, is_Pi, is_K, is_P);
2355             if (!use3dNoTime_) {
2356               if (no_PID || is_Pi) {
2357                 PartID = 211;
2358               } else if (is_K) {
2359                 PartID = 321;
2360               } else if (is_P) {
2361                 PartID = 2212;
2362               }
2363             }
2364             const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID)));
2365             double mass = PData->mass().value();
2366             massVector.push_back(mass);
2367             recotracks.push_back(**iTrack);
2368             getWosWnt(*vertex, *iTrack, sigmat0Safe, mtdQualMVA, BeamSpotH, wos, wnt);
2369             meVtxTrackWnt_->Fill(wnt);
2370             if (selectedLV) {
2371               meVtxTrackRecLVWnt_->Fill(wnt);
2372             }
2373             // reco track matched to any TP
2374             if (tp_info != nullptr) {
2375 #ifdef EDM_ML_DEBUG
2376               if (selectedLV) {
2377                 printMatchedRecoTrackInfo(*vertex, *iTrack, *tp_info, matchCategory);
2378               }
2379 #endif
2380               // matched TP not associated to any daughter track of a given sim vertex but having the same eventID (track from secondary vtx)
2381               if (matchCategory == 1) {
2382                 categoryVector.push_back(matchCategory);
2383                 SecsumWos += wos;
2384               }
2385               // matched TP not associated to any sim vertex of a given simulated event (PU track)
2386               if (matchCategory == 2) {
2387                 if (optionalPlots_) {
2388                   mePUTrackWnt_->Fill(wnt);
2389                   if (selectedLV) {
2390                     mePUTrackRecLVWnt_->Fill(wnt);
2391                   }
2392                 }
2393                 PUsumWnt += wnt;
2394                 PUsumWos += wos;
2395                 PUsumPt += (*iTrack)->pt();
2396                 PUsumPt2 += ((*iTrack)->pt() * (*iTrack)->pt());
2397                 PUnt++;
2398                 categoryVector.push_back(2);
2399               }
2400             }
2401             // reco track not matched to any TP (fake tracks)
2402             else {
2403               categoryVector.push_back(matchCategory);
2404               FakesumWos += wos;
2405               Fakent++;
2406             }
2407             nt++;
2408             sumWnt += wnt;
2409             sumWos += wos;
2410             sumPt += (*iTrack)->pt();
2411             sumPt2 += ((*iTrack)->pt() * (*iTrack)->pt());
2412           }
2413 
2414           // matched TP equal to any TP of a given sim vertex
2415           if (tp_info != nullptr && matchCategory == 0) {
2416             categoryVector.push_back(matchCategory);
2417             double mass = (*tp_info)->mass();
2418             double tsim = (*tp_info)->parentVertex()->position().t() * simUnit_;
2419             double tEst = timeFromTrueMass(mass, pathLength[*iTrack], momentum[*iTrack], time[*iTrack]);
2420 
2421             double xsim = (*tp_info)->parentVertex()->position().x();
2422             double ysim = (*tp_info)->parentVertex()->position().y();
2423             double zsim = (*tp_info)->parentVertex()->position().z();
2424             double xPCA = (*iTrack)->vx();
2425             double yPCA = (*iTrack)->vy();
2426             double zPCA = (*iTrack)->vz();
2427 
2428             double dZ = zPCA - zsim;
2429             double d3D = std::sqrt((xPCA - xsim) * (xPCA - xsim) + (yPCA - ysim) * (yPCA - ysim) + dZ * dZ);
2430             // orient d3D according to the projection of RECO - SIM onto simulated momentum
2431             if ((xPCA - xsim) * ((*tp_info)->px()) + (yPCA - ysim) * ((*tp_info)->py()) + dZ * ((*tp_info)->pz()) <
2432                 0.) {
2433               d3D = -d3D;
2434             }
2435 
2436             // select TPs associated to the signal event
2437             bool selectTP = trkTPSelLV(**tp_info);
2438 
2439             if (selectedLVMatching && selectRecoTrk && selectTP) {
2440               meTrackMatchedTPZposResTot_->Fill((*iTrack)->vz() - vzsim);
2441               if (noCrack) {
2442                 meTrackMatchedTPEffPtTot_->Fill((*iTrack)->pt());
2443               }
2444               meTrackMatchedTPEffEtaTot_->Fill(std::abs((*iTrack)->eta()));
2445             }
2446 
2447             if (sigmat0Safe[*iTrack] == -1)
2448               continue;
2449 
2450             if (selectedLVMatching && selectRecoTrk && selectTP) {
2451               meTrackMatchedTPResTot_->Fill(t0Safe[*iTrack] - vtsim);
2452               meTrackMatchedTPPullTot_->Fill((t0Safe[*iTrack] - vtsim) / sigmat0Safe[*iTrack]);
2453               if (noCrack) {
2454                 meTrackMatchedTPEffPtMtd_->Fill((*iTrack)->pt());
2455               }
2456               meTrackMatchedTPEffEtaMtd_->Fill(std::abs((*iTrack)->eta()));
2457 
2458               unsigned int noPIDtype = 0;
2459               bool noPID = false, isPi = false, isK = false, isP = false;
2460               isParticle(*iTrack, sigmat0, sigmat0Safe, probPi, probK, probP, noPIDtype, noPID, isPi, isK, isP);
2461 
2462               if ((isPi && std::abs(tMtd[*iTrack] - tofPi[*iTrack] - t0Pid[*iTrack]) > tol_) ||
2463                   (isK && std::abs(tMtd[*iTrack] - tofK[*iTrack] - t0Pid[*iTrack]) > tol_) ||
2464                   (isP && std::abs(tMtd[*iTrack] - tofP[*iTrack] - t0Pid[*iTrack]) > tol_)) {
2465                 edm::LogWarning("Primary4DVertexValidation")
2466                     << "No match between mass hyp. and time: " << std::abs((*tp_info)->pdgId()) << " mass hyp pi/k/p "
2467                     << isPi << " " << isK << " " << isP << " t0/t0safe " << t0Pid[*iTrack] << " " << t0Safe[*iTrack]
2468                     << " tMtd - tof pi/K/p " << tMtd[*iTrack] - tofPi[*iTrack] << " " << tMtd[*iTrack] - tofK[*iTrack]
2469                     << " " << tMtd[*iTrack] - tofP[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " "
2470                     << probK[*iTrack] << " " << probP[*iTrack];
2471               }
2472 
2473               if (std::abs((*iTrack)->eta()) < trackMaxBtlEta_) {
2474                 meBarrelPIDp_->Fill((*iTrack)->p());
2475                 meBarrelNoPIDtype_->Fill(noPIDtype + 0.5);
2476                 if (std::abs((*tp_info)->pdgId()) == 211) {
2477                   if (noPID) {
2478                     meBarrelTruePiNoPID_->Fill((*iTrack)->p());
2479                   } else if (isPi) {
2480                     meBarrelTruePiAsPi_->Fill((*iTrack)->p());
2481                   } else if (isK) {
2482                     meBarrelTruePiAsK_->Fill((*iTrack)->p());
2483                   } else if (isP) {
2484                     meBarrelTruePiAsP_->Fill((*iTrack)->p());
2485                   } else {
2486                     edm::LogWarning("Primary4DVertexValidation")
2487                         << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
2488                         << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
2489                         << probP[*iTrack];
2490                   }
2491                 } else if (std::abs((*tp_info)->pdgId()) == 321) {
2492                   if (noPID) {
2493                     meBarrelTrueKNoPID_->Fill((*iTrack)->p());
2494                   } else if (isPi) {
2495                     meBarrelTrueKAsPi_->Fill((*iTrack)->p());
2496                   } else if (isK) {
2497                     meBarrelTrueKAsK_->Fill((*iTrack)->p());
2498                   } else if (isP) {
2499                     meBarrelTrueKAsP_->Fill((*iTrack)->p());
2500                   } else {
2501                     edm::LogWarning("Primary4DVertexValidation")
2502                         << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
2503                         << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
2504                         << probP[*iTrack];
2505                   }
2506                 } else if (std::abs((*tp_info)->pdgId()) == 2212) {
2507                   if (noPID) {
2508                     meBarrelTruePNoPID_->Fill((*iTrack)->p());
2509                   } else if (isPi) {
2510                     meBarrelTruePAsPi_->Fill((*iTrack)->p());
2511                   } else if (isK) {
2512                     meBarrelTruePAsK_->Fill((*iTrack)->p());
2513                   } else if (isP) {
2514                     meBarrelTruePAsP_->Fill((*iTrack)->p());
2515                   } else {
2516                     edm::LogWarning("Primary4DVertexValidation")
2517                         << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
2518                         << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
2519                         << probP[*iTrack];
2520                   }
2521                 }
2522               } else if (std::abs((*iTrack)->eta()) > trackMinEtlEta_ && std::abs((*iTrack)->eta()) < trackMaxEtlEta_) {
2523                 meEndcapPIDp_->Fill((*iTrack)->p());
2524                 meEndcapNoPIDtype_->Fill(noPIDtype + 0.5);
2525                 if (std::abs((*tp_info)->pdgId()) == 211) {
2526                   if (noPID) {
2527                     meEndcapTruePiNoPID_->Fill((*iTrack)->p());
2528                   } else if (isPi) {
2529                     meEndcapTruePiAsPi_->Fill((*iTrack)->p());
2530                   } else if (isK) {
2531                     meEndcapTruePiAsK_->Fill((*iTrack)->p());
2532                   } else if (isP) {
2533                     meEndcapTruePiAsP_->Fill((*iTrack)->p());
2534                   } else {
2535                     edm::LogWarning("Primary4DVertexValidation")
2536                         << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
2537                         << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
2538                         << probP[*iTrack];
2539                   }
2540                 } else if (std::abs((*tp_info)->pdgId()) == 321) {
2541                   if (noPID) {
2542                     meEndcapTrueKNoPID_->Fill((*iTrack)->p());
2543                   } else if (isPi) {
2544                     meEndcapTrueKAsPi_->Fill((*iTrack)->p());
2545                   } else if (isK) {
2546                     meEndcapTrueKAsK_->Fill((*iTrack)->p());
2547                   } else if (isP) {
2548                     meEndcapTrueKAsP_->Fill((*iTrack)->p());
2549                   } else {
2550                     edm::LogWarning("Primary4DVertexValidation")
2551                         << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
2552                         << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
2553                         << probP[*iTrack];
2554                   }
2555                 } else if (std::abs((*tp_info)->pdgId()) == 2212) {
2556                   if (noPID) {
2557                     meEndcapTruePNoPID_->Fill((*iTrack)->p());
2558                   } else if (isPi) {
2559                     meEndcapTruePAsPi_->Fill((*iTrack)->p());
2560                   } else if (isK) {
2561                     meEndcapTruePAsK_->Fill((*iTrack)->p());
2562                   } else if (isP) {
2563                     meEndcapTruePAsP_->Fill((*iTrack)->p());
2564                   } else {
2565                     edm::LogWarning("Primary4DVertexValidation")
2566                         << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
2567                         << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
2568                         << probP[*iTrack];
2569                   }
2570                 }
2571               }
2572             }
2573             meTrackResTot_->Fill(t0Safe[*iTrack] - tsim);
2574             meTrackPullTot_->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2575             meTrackZposResTot_->Fill(dZ);
2576             if (selectRecoTrk && optionalPlots_) {
2577               unsigned int no_PIDtype = 0;
2578               bool no_PID, is_Pi, is_K, is_P;
2579               isParticle(*iTrack, sigmat0, sigmat0Safe, probPi, probK, probP, no_PIDtype, no_PID, is_Pi, is_K, is_P);
2580               if (no_PID) {
2581                 meTrackTimeResNoPID_->Fill(t0Safe[*iTrack] - tsim);
2582                 meTrackTimePullNoPID_->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2583                 meTrackTimeSigmaNoPID_->Fill(sigmat0Safe[*iTrack]);
2584                 meTrackMVANoPID_->Fill(mtdQualMVA[(*iTrack)]);
2585                 if (no_PIDtype == 1) {
2586                   meNoPIDTrackTimeResNoPIDType_[0]->Fill(t0Safe[*iTrack] - tsim);
2587                   meNoPIDTrackTimePullNoPIDType_[0]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2588                   meNoPIDTrackSigmaNoPIDType_[0]->Fill(sigmat0Safe[*iTrack]);
2589                   meNoPIDTrackMVANoPIDType_[0]->Fill(mtdQualMVA[(*iTrack)]);
2590                 } else if (no_PIDtype == 2) {
2591                   meNoPIDTrackTimeResNoPIDType_[1]->Fill(t0Safe[*iTrack] - tsim);
2592                   meNoPIDTrackTimePullNoPIDType_[1]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2593                   meNoPIDTrackSigmaNoPIDType_[1]->Fill(sigmat0Safe[*iTrack]);
2594                   meNoPIDTrackMVANoPIDType_[1]->Fill(mtdQualMVA[(*iTrack)]);
2595                 } else if (no_PIDtype == 3) {
2596                   meNoPIDTrackTimeResNoPIDType_[2]->Fill(t0Safe[*iTrack] - tsim);
2597                   meNoPIDTrackTimePullNoPIDType_[2]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2598                   meNoPIDTrackSigmaNoPIDType_[2]->Fill(sigmat0Safe[*iTrack]);
2599                   meNoPIDTrackMVANoPIDType_[2]->Fill(mtdQualMVA[(*iTrack)]);
2600                 }
2601                 if (std::abs((*tp_info)->pdgId()) == 211) {
2602                   meTrackTimeResNoPIDtruePi_->Fill(t0Safe[*iTrack] - tsim);
2603                   meTrackTimePullNoPIDtruePi_->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2604                 } else if (std::abs((*tp_info)->pdgId()) == 321) {
2605                   meTrackTimeResNoPIDtrueK_->Fill(t0Safe[*iTrack] - tsim);
2606                   meTrackTimePullNoPIDtrueK_->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2607                 } else if (std::abs((*tp_info)->pdgId()) == 2212) {
2608                   meTrackTimeResNoPIDtrueP_->Fill(t0Safe[*iTrack] - tsim);
2609                   meTrackTimePullNoPIDtrueP_->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2610                 }
2611               } else if ((is_Pi && std::abs((*tp_info)->pdgId()) == 211) ||
2612                          (is_K && std::abs((*tp_info)->pdgId()) == 321) ||
2613                          (is_P && std::abs((*tp_info)->pdgId()) == 2212)) {
2614                 meTrackTimeResCorrectPID_->Fill(t0Safe[*iTrack] - tsim);
2615                 meTrackTimePullCorrectPID_->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2616                 meTrackTimeSigmaCorrectPID_->Fill(sigmat0Safe[*iTrack]);
2617                 meTrackMVACorrectPID_->Fill(mtdQualMVA[(*iTrack)]);
2618               } else {
2619                 meTrackTimeResWrongPID_->Fill(t0Safe[*iTrack] - tsim);
2620                 meTrackTimePullWrongPID_->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2621                 meTrackTimeSigmaWrongPID_->Fill(sigmat0Safe[*iTrack]);
2622                 meTrackMVAWrongPID_->Fill(mtdQualMVA[(*iTrack)]);
2623               }
2624             }
2625 
2626             if ((*iTrack)->p() <= 2) {
2627               meTrackResLowPTot_->Fill(t0Safe[*iTrack] - tsim);
2628               meTrackPullLowPTot_->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2629             } else {
2630               meTrackResHighPTot_->Fill(t0Safe[*iTrack] - tsim);
2631               meTrackPullHighPTot_->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2632             }
2633 
2634             if (mtdQualMVA[(*iTrack)] < mvaL_) {
2635               meTrackZposRes_[0]->Fill(dZ);
2636               meTrack3DposRes_[0]->Fill(d3D);
2637               meTrackRes_[0]->Fill(t0Safe[*iTrack] - tsim);
2638               meTrackPull_[0]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2639 
2640               if (optionalPlots_) {
2641                 meTrackResMass_[0]->Fill(t0Safe[*iTrack] - tEst);
2642                 meTrackResMassTrue_[0]->Fill(tEst - tsim);
2643               }
2644 
2645               if ((*iTrack)->p() <= 2) {
2646                 meTrackResLowP_[0]->Fill(t0Safe[*iTrack] - tsim);
2647                 meTrackPullLowP_[0]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2648               } else if ((*iTrack)->p() > 2) {
2649                 meTrackResHighP_[0]->Fill(t0Safe[*iTrack] - tsim);
2650                 meTrackPullHighP_[0]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2651               }
2652 
2653               if (optionalPlots_) {
2654                 if (std::abs((*tp_info)->pdgId()) == 2212) {
2655                   meTrackResMassProtons_[0]->Fill(t0Safe[*iTrack] - tEst);
2656                   meTrackResMassTrueProtons_[0]->Fill(tEst - tsim);
2657                 } else if (std::abs((*tp_info)->pdgId()) == 211) {
2658                   meTrackResMassPions_[0]->Fill(t0Safe[*iTrack] - tEst);
2659                   meTrackResMassTruePions_[0]->Fill(tEst - tsim);
2660                 }
2661               }
2662 
2663             } else if (mtdQualMVA[(*iTrack)] > mvaL_ && mtdQualMVA[(*iTrack)] < mvaH_) {
2664               meTrackZposRes_[1]->Fill(dZ);
2665               meTrack3DposRes_[1]->Fill(d3D);
2666               meTrackRes_[1]->Fill(t0Safe[*iTrack] - tsim);
2667               meTrackPull_[1]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2668 
2669               if (optionalPlots_) {
2670                 meTrackResMass_[1]->Fill(t0Safe[*iTrack] - tEst);
2671                 meTrackResMassTrue_[1]->Fill(tEst - tsim);
2672               }
2673 
2674               if ((*iTrack)->p() <= 2) {
2675                 meTrackResLowP_[1]->Fill(t0Safe[*iTrack] - tsim);
2676                 meTrackPullLowP_[1]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2677               } else if ((*iTrack)->p() > 2) {
2678                 meTrackResHighP_[1]->Fill(t0Safe[*iTrack] - tsim);
2679                 meTrackPullHighP_[1]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2680               }
2681 
2682               if (optionalPlots_) {
2683                 if (std::abs((*tp_info)->pdgId()) == 2212) {
2684                   meTrackResMassProtons_[1]->Fill(t0Safe[*iTrack] - tEst);
2685                   meTrackResMassTrueProtons_[1]->Fill(tEst - tsim);
2686                 } else if (std::abs((*tp_info)->pdgId()) == 211) {
2687                   meTrackResMassPions_[1]->Fill(t0Safe[*iTrack] - tEst);
2688                   meTrackResMassTruePions_[1]->Fill(tEst - tsim);
2689                 }
2690               }
2691 
2692             } else if (mtdQualMVA[(*iTrack)] > mvaH_) {
2693               meTrackZposRes_[2]->Fill(dZ);
2694               meTrack3DposRes_[2]->Fill(d3D);
2695               meTrackRes_[2]->Fill(t0Safe[*iTrack] - tsim);
2696               meTrackPull_[2]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2697 
2698               if (optionalPlots_) {
2699                 meTrackResMass_[2]->Fill(t0Safe[*iTrack] - tEst);
2700                 meTrackResMassTrue_[2]->Fill(tEst - tsim);
2701               }
2702 
2703               if ((*iTrack)->p() <= 2) {
2704                 meTrackResLowP_[2]->Fill(t0Safe[*iTrack] - tsim);
2705                 meTrackPullLowP_[2]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2706               } else if ((*iTrack)->p() > 2) {
2707                 meTrackResHighP_[2]->Fill(t0Safe[*iTrack] - tsim);
2708                 meTrackPullHighP_[2]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2709               }
2710 
2711               if (optionalPlots_) {
2712                 if (std::abs((*tp_info)->pdgId()) == 2212) {
2713                   meTrackResMassProtons_[2]->Fill(t0Safe[*iTrack] - tEst);
2714                   meTrackResMassTrueProtons_[2]->Fill(tEst - tsim);
2715                 } else if (std::abs((*tp_info)->pdgId()) == 211) {
2716                   meTrackResMassPions_[2]->Fill(t0Safe[*iTrack] - tEst);
2717                   meTrackResMassTruePions_[2]->Fill(tEst - tsim);
2718                 }
2719               }
2720             }
2721           }  // if tp_info != nullptr && MatchCategory == 0
2722         }  // loop on reco tracks
2723         if (selectedVtxMatching) {
2724           meVtxTrackMult_->Fill(log10(nt));
2725           mePUTrackRelMult_->Fill(static_cast<double>(PUnt) / nt);
2726           meFakeTrackRelMult_->Fill(static_cast<double>(Fakent) / nt);
2727           mePUTrackRelSumWnt_->Fill(PUsumWnt / sumWnt);
2728           mePUTrackRelSumWos_->Fill(PUsumWos / sumWos);
2729           meSecTrackRelSumWos_->Fill(SecsumWos / sumWos);
2730           meFakeTrackRelSumWos_->Fill(FakesumWos / sumWos);
2731           mePUTrackRelSumPt_->Fill(PUsumPt / sumPt);
2732           mePUTrackRelSumPt2_->Fill(PUsumPt2 / sumPt2);
2733 
2734           observablesFromJets(
2735               recotracks, massVector, categoryVector, "use_allTracks", nJets, sumEtJets, sumPt2Jets, metPt, sumPzJets);
2736           observablesFromJets(recotracks,
2737                               massVector,
2738                               categoryVector,
2739                               "skip_Fake",
2740                               nJetsnoFake,
2741                               sumEtJetsnoFake,
2742                               sumPt2JetsnoFake,
2743                               metPtnoFake,
2744                               sumPzJetsnoFake);
2745           observablesFromJets(recotracks,
2746                               massVector,
2747                               categoryVector,
2748                               "skip_PU",
2749                               nJetsnoPU,
2750                               sumEtJetsnoPU,
2751                               sumPt2JetsnoPU,
2752                               metPtnoPU,
2753                               sumPzJetsnoPU);
2754 
2755           meJetsPURelMult_->Fill(static_cast<double>(nJets - nJetsnoPU) / nJets);
2756           meJetsPURelHt_->Fill((sumEtJets - sumEtJetsnoPU) / sumEtJets);
2757           meJetsPURelSumPt2_->Fill((sumPt2Jets - sumPt2JetsnoPU) / sumPt2Jets);
2758           meJetsFakeRelSumPt2_->Fill((sumPt2Jets - sumPt2JetsnoFake) / sumPt2Jets);
2759           meJetsPURelMetPt_->Fill((metPt - metPtnoPU) / metPt);
2760           meJetsPURelSumPz_->Fill((sumPzJets - sumPzJetsnoPU) / sumPzJets);
2761 
2762           if (optionalPlots_) {
2763             mePUTrackMult_->Fill(PUnt);
2764             mePUTrackSumWnt_->Fill(log10(std::max(minThrSumWnt_, PUsumWnt)));
2765             mePUTrackSumWos_->Fill(log10(std::max(minThrSumWos_, PUsumWos)));
2766             meSecTrackSumWos_->Fill(log10(std::max(minThrSumWos_, SecsumWos)));
2767             mePUTrackSumPt_->Fill(log10(std::max(minThrSumPt_, PUsumPt)));
2768             mePUTrackSumPt2_->Fill(log10(std::max(minThrSumPt2_, PUsumPt2)));
2769 
2770             mePUTrackRelMultvsMult_->Fill(nt, static_cast<double>(PUnt) / nt);
2771             meFakeTrackRelMultvsMult_->Fill(nt, static_cast<double>(Fakent) / nt);
2772             mePUTrackRelSumWntvsSumWnt_->Fill(log10(std::max(minThrSumWnt_, sumWnt)), PUsumWnt / sumWnt);
2773             mePUTrackRelSumWosvsSumWos_->Fill(log10(std::max(minThrSumWos_, sumWos)), PUsumWos / sumWos);
2774             meSecTrackRelSumWosvsSumWos_->Fill(log10(std::max(minThrSumWos_, sumWos)), SecsumWos / sumWos);
2775             meFakeTrackRelSumWosvsSumWos_->Fill(log10(std::max(minThrSumWos_, sumWos)), FakesumWos / sumWos);
2776             mePUTrackRelSumPtvsSumPt_->Fill(log10(std::max(minThrSumPt_, sumPt)), PUsumPt / sumPt);
2777             mePUTrackRelSumPt2vsSumPt2_->Fill(log10(std::max(minThrSumPt2_, sumPt2)), PUsumPt2 / sumPt2);
2778 
2779             meJetsPUMult_->Fill(nJets - nJetsnoPU);
2780             meJetsPUHt_->Fill(log10(std::max(minThrSumPt_, sumEtJets - sumEtJetsnoPU)));
2781             meJetsPUSumPt2_->Fill(log10(std::max(minThrSumPt2_, sumPt2Jets - sumPt2JetsnoPU)));
2782             meJetsPUMetPt_->Fill(log10(std::max(minThrMetPt_, metPt - metPtnoPU)));
2783             meJetsPUSumPz_->Fill(log10(std::max(minThrSumPz_, std::abs(sumPzJets - sumPzJetsnoPU))));
2784 
2785             meJetsPURelMultvsMult_->Fill(nJets, static_cast<double>(nJets - nJetsnoPU) / nJets);
2786             meJetsPURelHtvsHt_->Fill(log10(std::max(minThrSumPt_, sumEtJets)), (sumEtJets - sumEtJetsnoPU) / sumEtJets);
2787             meJetsPURelSumPt2vsSumPt2_->Fill(log10(std::max(minThrSumPt2_, sumPt2Jets)),
2788                                              (sumPt2Jets - sumPt2JetsnoPU) / sumPt2Jets);
2789             meJetsFakeRelSumPt2vsSumPt2_->Fill(log10(std::max(minThrSumPt2_, sumPt2Jets)),
2790                                                (sumPt2Jets - sumPt2JetsnoFake) / sumPt2Jets);
2791             meJetsPURelMetPtvsMetPt_->Fill(log10(std::max(minThrMetPt_, metPt)), (metPt - metPtnoPU) / metPt);
2792             meJetsPURelSumPzvsSumPz_->Fill(log10(std::max(minThrSumPz_, std::abs(sumPzJets))),
2793                                            (sumPzJets - sumPzJetsnoPU) / sumPzJets);
2794           }
2795           if (selectedLV) {
2796             meVtxTrackRecLVMult_->Fill(log10(nt));
2797             mePUTrackRecLVRelMult_->Fill(static_cast<double>(PUnt) / nt);
2798             meFakeTrackRecLVRelMult_->Fill(static_cast<double>(Fakent) / nt);
2799             mePUTrackRecLVRelSumWnt_->Fill(PUsumWnt / sumWnt);
2800             mePUTrackRecLVRelSumWos_->Fill(PUsumWos / sumWos);
2801             meSecTrackRecLVRelSumWos_->Fill(SecsumWos / sumWos);
2802             meFakeTrackRecLVRelSumWos_->Fill(FakesumWos / sumWos);
2803             mePUTrackRecLVRelSumPt_->Fill(PUsumPt / sumPt);
2804             mePUTrackRecLVRelSumPt2_->Fill(PUsumPt2 / sumPt2);
2805 
2806             meJetsRecLVPURelMult_->Fill(static_cast<double>(nJets - nJetsnoPU) / nJets);
2807             meJetsRecLVPURelHt_->Fill((sumEtJets - sumEtJetsnoPU) / sumEtJets);
2808             meJetsRecLVPURelSumPt2_->Fill((sumPt2Jets - sumPt2JetsnoPU) / sumPt2Jets);
2809             meJetsRecLVFakeRelSumPt2_->Fill((sumPt2Jets - sumPt2JetsnoFake) / sumPt2Jets);
2810             meJetsRecLVPURelMetPt_->Fill((metPt - metPtnoPU) / metPt);
2811             meJetsRecLVPURelSumPz_->Fill((sumPzJets - sumPzJetsnoPU) / sumPzJets);
2812 
2813             LogTrace("Primary4DVertexValidation")
2814                 << "#PUTrks = " << PUnt << " #Trks = " << nt << " PURelMult = " << std::setprecision(3)
2815                 << static_cast<double>(PUnt) / nt;
2816             LogTrace("Primary4DVertexValidation")
2817                 << "PUsumWnt = " << std::setprecision(3) << PUsumWnt << " sumWnt = " << std::setprecision(3) << sumWnt
2818                 << " PURelsumWnt = " << std::setprecision(3) << PUsumWnt / sumWnt;
2819             LogTrace("Primary4DVertexValidation")
2820                 << "PUsumWos = " << std::setprecision(3) << PUsumWos << " sumWos = " << std::setprecision(3) << sumWos
2821                 << " PURelsumWos = " << std::setprecision(3) << PUsumWos / sumWos;
2822             LogTrace("Primary4DVertexValidation")
2823                 << "PuSumPt = " << std::setprecision(3) << PUsumPt << " SumPt = " << std::setprecision(4) << sumPt
2824                 << " PURelSumPt = " << std::setprecision(3) << PUsumPt / sumPt;
2825             LogTrace("Primary4DVertexValidation")
2826                 << "PuSumPt2 = " << std::setprecision(3) << PUsumPt2 << " SumPt2 = " << std::setprecision(4) << sumPt2
2827                 << " PURelSumPt2 = " << std::setprecision(3) << PUsumPt2 / sumPt2;
2828             if (optionalPlots_) {
2829               mePUTrackRecLVMult_->Fill(PUnt);
2830               mePUTrackRecLVSumWnt_->Fill(log10(std::max(minThrSumWnt_, PUsumWnt)));
2831               mePUTrackRecLVSumWos_->Fill(log10(std::max(minThrSumWos_, PUsumWos)));
2832               meSecTrackRecLVSumWos_->Fill(log10(std::max(minThrSumWos_, PUsumWos)));
2833               mePUTrackRecLVSumPt_->Fill(log10(std::max(minThrSumPt_, PUsumPt)));
2834               mePUTrackRecLVSumPt2_->Fill(log10(std::max(minThrSumPt2_, PUsumPt2)));
2835 
2836               mePUTrackRecLVRelMultvsMult_->Fill(nt, static_cast<double>(PUnt) / nt);
2837               meFakeTrackRecLVRelMultvsMult_->Fill(nt, static_cast<double>(Fakent) / nt);
2838               mePUTrackRecLVRelSumWntvsSumWnt_->Fill(log10(std::max(minThrSumWnt_, sumWnt)), PUsumWnt / sumWnt);
2839               mePUTrackRecLVRelSumWosvsSumWos_->Fill(log10(std::max(minThrSumWos_, sumWos)), PUsumWos / sumWos);
2840               meSecTrackRecLVRelSumWosvsSumWos_->Fill(log10(std::max(minThrSumWos_, sumWos)), SecsumWos / sumWos);
2841               meFakeTrackRecLVRelSumWosvsSumWos_->Fill(log10(std::max(minThrSumWos_, sumWos)), FakesumWos / sumWos);
2842               mePUTrackRecLVRelSumPtvsSumPt_->Fill(log10(std::max(minThrSumPt_, sumPt)), PUsumPt / sumPt);
2843               mePUTrackRecLVRelSumPt2vsSumPt2_->Fill(log10(std::max(minThrSumPt2_, sumPt2)), PUsumPt2 / sumPt2);
2844 
2845               meJetsRecLVPUMult_->Fill(nJets - nJetsnoPU);
2846               meJetsRecLVPUHt_->Fill(log10(std::max(minThrSumPt_, sumEtJets - sumEtJetsnoPU)));
2847               meJetsRecLVPUSumPt2_->Fill(log10(std::max(minThrSumPt2_, sumPt2Jets - sumPt2JetsnoPU)));
2848               meJetsRecLVPUMetPt_->Fill(log10(std::max(minThrMetPt_, metPt - metPtnoPU)));
2849               meJetsRecLVPUSumPz_->Fill(log10(std::max(minThrSumPz_, std::abs(sumPzJets - sumPzJetsnoPU))));
2850 
2851               meJetsRecLVPURelMultvsMult_->Fill(nJets, static_cast<double>(nJets - nJetsnoPU) / nJets);
2852               meJetsRecLVPURelHtvsHt_->Fill(log10(std::max(minThrSumPt_, sumEtJets)),
2853                                             (sumEtJets - sumEtJetsnoPU) / sumEtJets);
2854               meJetsRecLVPURelSumPt2vsSumPt2_->Fill(log10(std::max(minThrSumPt2_, sumPt2Jets)),
2855                                                     (sumPt2Jets - sumPt2JetsnoPU) / sumPt2Jets);
2856               meJetsRecLVFakeRelSumPt2vsSumPt2_->Fill(log10(std::max(minThrSumPt2_, sumPt2Jets)),
2857                                                       (sumPt2Jets - sumPt2JetsnoFake) / sumPt2Jets);
2858               meJetsRecLVPURelMetPtvsMetPt_->Fill(log10(std::max(minThrMetPt_, metPt)), (metPt - metPtnoPU) / metPt);
2859               meJetsRecLVPURelSumPzvsSumPz_->Fill(log10(std::max(minThrSumPz_, std::abs(sumPzJets))),
2860                                                   (sumPzJets - sumPzJetsnoPU) / sumPzJets);
2861             }
2862           }
2863         }
2864       }  // loop on simpv
2865     }  // ndof
2866   }  // loop on recopv
2867 
2868   int real = 0;
2869   int fake = 0;
2870   int other_fake = 0;
2871   int split = 0;
2872 
2873   meRecVerNumber_->Fill(recopv.size());
2874   for (unsigned int ir = 0; ir < recopv.size(); ir++) {
2875     const reco::Vertex* vertex = recopv.at(ir).recVtx;
2876     if (recopv.at(ir).ndof > selNdof_) {
2877       meRecPVZ_->Fill(recopv.at(ir).z);
2878       meVtxTrackMultPassNdof_->Fill(log10(vertex->tracksSize()));
2879 
2880       if (recopv.at(ir).recVtx->tError() > 0.) {
2881         meRecPVT_->Fill(recopv.at(ir).recVtx->t());
2882       }
2883       LogTrace("Primary4DVertexValidation") << "************* IR: " << ir;
2884       LogTrace("Primary4DVertexValidation") << "is_real: " << recopv.at(ir).is_real();
2885       LogTrace("Primary4DVertexValidation") << "is_fake: " << recopv.at(ir).is_fake();
2886       LogTrace("Primary4DVertexValidation") << "is_signal: " << recopv.at(ir).is_signal();
2887       LogTrace("Primary4DVertexValidation") << "split_from: " << recopv.at(ir).split_from();
2888       LogTrace("Primary4DVertexValidation") << "other fake: " << recopv.at(ir).other_fake();
2889       if (recopv.at(ir).is_real()) {
2890         real++;
2891       }
2892       if (recopv.at(ir).is_fake()) {
2893         fake++;
2894       }
2895       if (recopv.at(ir).other_fake()) {
2896         other_fake++;
2897       }
2898       if (recopv.at(ir).split_from() != -1) {
2899         split++;
2900       }
2901     }  // ndof
2902     else {
2903       meVtxTrackMultFailNdof_->Fill(vertex->tracksSize());
2904     }
2905   }
2906 
2907   LogTrace("Primary4DVertexValidation") << "is_real: " << real;
2908   LogTrace("Primary4DVertexValidation") << "is_fake: " << fake;
2909   LogTrace("Primary4DVertexValidation") << "split_from: " << split;
2910   LogTrace("Primary4DVertexValidation") << "other fake: " << other_fake;
2911   mePUvsRealV_->Fill(simpv.size(), real);
2912   mePUvsFakeV_->Fill(simpv.size(), fake);
2913   mePUvsOtherFakeV_->Fill(simpv.size(), other_fake);
2914   mePUvsSplitV_->Fill(simpv.size(), split);
2915 
2916   // fill vertices histograms here in a new loop
2917   meSimVerNumber_->Fill(simpv.size());
2918   for (unsigned int is = 0; is < simpv.size(); is++) {
2919     meSimPVZ_->Fill(simpv.at(is).z);
2920     meSimPVT_->Fill(simpv.at(is).t * simUnit_);
2921     meSimPVTvsZ_->Fill(simpv.at(is).z, simpv.at(is).t * simUnit_);
2922     if (is == 0 && optionalPlots_) {
2923       meSimPosInSimOrigCollection_->Fill(simpv.at(is).OriginalIndex);
2924     }
2925 
2926     if (simpv.at(is).rec == NOT_MATCHED) {
2927       LogTrace("Primary4DVertexValidation") << "sim vertex: " << is << " is not matched with any reco";
2928       continue;
2929     }
2930 
2931     for (unsigned int ir = 0; ir < recopv.size(); ir++) {
2932       if (recopv.at(ir).ndof > selNdof_) {
2933         if (recopv.at(ir).sim == is && simpv.at(is).rec == ir) {
2934           meTimeRes_->Fill(recopv.at(ir).recVtx->t() - simpv.at(is).t * simUnit_);
2935           meTimePull_->Fill((recopv.at(ir).recVtx->t() - simpv.at(is).t * simUnit_) / recopv.at(ir).recVtx->tError());
2936           meMatchQual_->Fill(recopv.at(ir).matchQuality - 0.5);
2937           if (ir == 0) {  // signal vertex plots
2938             meTimeSignalRes_->Fill(recopv.at(ir).recVtx->t() - simpv.at(is).t * simUnit_);
2939             meTimeSignalPull_->Fill((recopv.at(ir).recVtx->t() - simpv.at(is).t * simUnit_) /
2940                                     recopv.at(ir).recVtx->tError());
2941             if (optionalPlots_) {
2942               meRecoPosInSimCollection_->Fill(recopv.at(ir).sim);
2943               meRecoPosInRecoOrigCollection_->Fill(recopv.at(ir).OriginalIndex);
2944             }
2945           }
2946           if (simpv.at(is).eventId.bunchCrossing() == 0 && simpv.at(is).eventId.event() == 0) {
2947             if (!recopv.at(ir).is_signal()) {
2948               edm::LogWarning("Primary4DVertexValidation")
2949                   << "Reco vtx leading match inconsistent: BX/ID " << simpv.at(is).eventId.bunchCrossing() << " "
2950                   << simpv.at(is).eventId.event();
2951             }
2952             meRecoPVPosSignal_->Fill(
2953                 recopv.at(ir).OriginalIndex);  // position in reco vtx correction associated to sim signal
2954             if (!signal_is_highest_pt) {
2955               meRecoPVPosSignalNotHighestPt_->Fill(
2956                   recopv.at(ir).OriginalIndex);  // position in reco vtx correction associated to sim signal
2957             }
2958           }
2959           LogTrace("Primary4DVertexValidation") << "*** Matching RECO: " << ir << "with SIM: " << is << " ***";
2960           LogTrace("Primary4DVertexValidation") << "Match Quality is " << recopv.at(ir).matchQuality;
2961           LogTrace("Primary4DVertexValidation") << "****";
2962         }
2963       }  // ndof
2964     }
2965   }
2966 
2967   // dz histos
2968   for (unsigned int iv = 0; iv < recVtxs->size() - 1; iv++) {
2969     if (recVtxs->at(iv).ndof() > selNdof_) {
2970       double mindistance_realreal = 1e10;
2971 
2972       for (unsigned int jv = iv; jv < recVtxs->size(); jv++) {
2973         if ((!(jv == iv)) && select(recVtxs->at(jv))) {
2974           double dz = recVtxs->at(iv).z() - recVtxs->at(jv).z();
2975           double dtsigma = std::sqrt(recVtxs->at(iv).covariance(3, 3) + recVtxs->at(jv).covariance(3, 3));
2976           double dt = (std::abs(dz) <= deltaZcut_ && dtsigma > 0.)
2977                           ? (recVtxs->at(iv).t() - recVtxs->at(jv).t()) / dtsigma
2978                           : -9999.;
2979           if (recopv.at(iv).is_real() && recopv.at(jv).is_real()) {
2980             meDeltaZrealreal_->Fill(std::abs(dz));
2981             if (dt != -9999.) {
2982               meDeltaTrealreal_->Fill(std::abs(dt));
2983             }
2984             if (std::abs(dz) < std::abs(mindistance_realreal)) {
2985               mindistance_realreal = dz;
2986             }
2987           } else if (recopv.at(iv).is_fake() && recopv.at(jv).is_fake()) {
2988             meDeltaZfakefake_->Fill(std::abs(dz));
2989             if (dt != -9999.) {
2990               meDeltaTfakefake_->Fill(std::abs(dt));
2991             }
2992           }
2993         }
2994       }
2995 
2996       double mindistance_fakereal = 1e10;
2997       double mindistance_realfake = 1e10;
2998       for (unsigned int jv = 0; jv < recVtxs->size(); jv++) {
2999         if ((!(jv == iv)) && select(recVtxs->at(jv))) {
3000           double dz = recVtxs->at(iv).z() - recVtxs->at(jv).z();
3001           double dtsigma = std::sqrt(recVtxs->at(iv).covariance(3, 3) + recVtxs->at(jv).covariance(3, 3));
3002           double dt = (std::abs(dz) <= deltaZcut_ && dtsigma > 0.)
3003                           ? (recVtxs->at(iv).t() - recVtxs->at(jv).t()) / dtsigma
3004                           : -9999.;
3005           if (recopv.at(iv).is_fake() && recopv.at(jv).is_real()) {
3006             meDeltaZfakereal_->Fill(std::abs(dz));
3007             if (dt != -9999.) {
3008               meDeltaTfakereal_->Fill(std::abs(dt));
3009             }
3010             if (std::abs(dz) < std::abs(mindistance_fakereal)) {
3011               mindistance_fakereal = dz;
3012             }
3013           }
3014 
3015           if (recopv.at(iv).is_real() && recopv.at(jv).is_fake() && (std::abs(dz) < std::abs(mindistance_realfake))) {
3016             mindistance_realfake = dz;
3017           }
3018         }
3019       }
3020     }  // ndof
3021   }
3022 
3023 }  // end of analyze
3024 
3025 void Primary4DVertexValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
3026   edm::ParameterSetDescription desc;
3027 
3028   desc.add<std::string>("folder", "MTD/Vertices");
3029   desc.add<edm::InputTag>("TPtoRecoTrackAssoc", edm::InputTag("trackingParticleRecoTrackAsssociation"));
3030   desc.add<edm::InputTag>("TrackLabel", edm::InputTag("generalTracks"));
3031   desc.add<edm::InputTag>("mtdTracks", edm::InputTag("trackExtenderWithMTD"));
3032   desc.add<edm::InputTag>("SimTag", edm::InputTag("mix", "MergedTrackTruth"));
3033   desc.add<edm::InputTag>("offlineBS", edm::InputTag("offlineBeamSpot"));
3034   desc.add<edm::InputTag>("offline4DPV", edm::InputTag("offlinePrimaryVertices4D"));
3035   desc.add<edm::InputTag>("trackAssocSrc", edm::InputTag("trackExtenderWithMTD:generalTrackassoc"))
3036       ->setComment("Association between General and MTD Extended tracks");
3037   desc.add<edm::InputTag>("pathLengthSrc", edm::InputTag("trackExtenderWithMTD:generalTrackPathLength"));
3038   desc.add<edm::InputTag>("momentumSrc", edm::InputTag("trackExtenderWithMTD:generalTrackp"));
3039   desc.add<edm::InputTag>("tmtd", edm::InputTag("trackExtenderWithMTD:generalTracktmtd"));
3040   desc.add<edm::InputTag>("timeSrc", edm::InputTag("trackExtenderWithMTD:generalTracktmtd"));
3041   desc.add<edm::InputTag>("sigmaSrc", edm::InputTag("trackExtenderWithMTD:generalTracksigmatmtd"));
3042   desc.add<edm::InputTag>("t0PID", edm::InputTag("tofPID:t0"));
3043   desc.add<edm::InputTag>("sigmat0PID", edm::InputTag("tofPID:sigmat0"));
3044   desc.add<edm::InputTag>("t0SafePID", edm::InputTag("tofPID:t0safe"));
3045   desc.add<edm::InputTag>("sigmat0SafePID", edm::InputTag("tofPID:sigmat0safe"));
3046   desc.add<edm::InputTag>("trackMVAQual", edm::InputTag("mtdTrackQualityMVA:mtdQualMVA"));
3047   desc.add<edm::InputTag>("tofPi", edm::InputTag("trackExtenderWithMTD:generalTrackTofPi"));
3048   desc.add<edm::InputTag>("tofK", edm::InputTag("trackExtenderWithMTD:generalTrackTofK"));
3049   desc.add<edm::InputTag>("tofP", edm::InputTag("trackExtenderWithMTD:generalTrackTofP"));
3050   desc.add<edm::InputTag>("probPi", edm::InputTag("tofPID:probPi"));
3051   desc.add<edm::InputTag>("probK", edm::InputTag("tofPID:probK"));
3052   desc.add<edm::InputTag>("probP", edm::InputTag("tofPID:probP"));
3053   desc.add<bool>("useOnlyChargedTracks", true);
3054   desc.addUntracked<bool>("optionalPlots", false);
3055   desc.add<bool>("use3dNoTime", false);
3056   desc.add<double>("trackweightTh", 0.5);
3057   desc.add<double>("mvaTh", 0.8);
3058   desc.add<double>("minProbHeavy", 0.75);
3059   {
3060     edm::ParameterSetDescription psd0;
3061     HITrackFilterForPVFinding::fillPSetDescription(psd0);  // extension of TrackFilterForPVFinding
3062     desc.add<edm::ParameterSetDescription>("TkFilterParameters", psd0);
3063   }
3064 
3065   descriptions.add("vertices4DValid", desc);
3066 }
3067 
3068 void Primary4DVertexValidation::printMatchedRecoTrackInfo(const reco::Vertex& vtx,
3069                                                           const reco::TrackBaseRef& trk,
3070                                                           const TrackingParticleRef& tp,
3071                                                           const unsigned int& categ) {
3072   std::string strTrk;
3073   switch (categ) {
3074     case 0:
3075       strTrk = "Reco_Track:";
3076       break;
3077     case 1:
3078       strTrk = "SecRecoTrk:";
3079       break;
3080     case 2:
3081       strTrk = "PU_RecoTrk:";
3082       break;
3083   }
3084   LogTrace("Primary4DVertexValidation") << strTrk << " w =" << std::setw(6) << std::setprecision(2)
3085                                         << vtx.trackWeight(trk) << " pt =" << std::setw(6) << std::setprecision(2)
3086                                         << trk->pt() << " eta =" << std::setw(6) << std::setprecision(2) << trk->eta()
3087                                         << "  MatchedTP: Pt =" << std::setw(6) << std::setprecision(2) << tp->pt()
3088                                         << " eta =" << std::setw(6) << std::setprecision(2) << tp->eta()
3089                                         << "  Parent vtx: z =" << std::setw(8) << std::setprecision(4)
3090                                         << tp->parentVertex()->position().z() << " t =" << std::setw(8)
3091                                         << std::setprecision(4) << tp->parentVertex()->position().t() * simUnit_
3092                                         << " BX =" << tp->parentVertex()->eventId().bunchCrossing()
3093                                         << " ev =" << tp->parentVertex()->eventId().event() << std::endl;
3094 }
3095 
3096 void Primary4DVertexValidation::printSimVtxRecoVtxInfo(
3097     const struct Primary4DVertexValidation::simPrimaryVertex& simpVtx,
3098     const struct Primary4DVertexValidation::recoPrimaryVertex& recopVtx) {
3099   LogTrace("Primary4DVertexValidation") << "Sim vtx (x,y,z,t) = (" << std::setprecision(4) << simpVtx.x << ","
3100                                         << std::setprecision(4) << simpVtx.y << "," << std::setprecision(4) << simpVtx.z
3101                                         << "," << std::setprecision(4) << simpVtx.t * simUnit_ << ")"
3102                                         << " Simvtx.rec = " << simpVtx.rec;
3103   LogTrace("Primary4DVertexValidation") << "Sim vtx: pt = " << std::setprecision(4) << simpVtx.pt
3104                                         << " ptsq = " << std::setprecision(6) << simpVtx.ptsq
3105                                         << " nGenTrk = " << simpVtx.nGenTrk
3106                                         << " nmatch recotrks = " << simpVtx.num_matched_reco_tracks;
3107   LogTrace("Primary4DVertexValidation") << "Reco vtx (x,y,z) = (" << std::setprecision(4) << recopVtx.x << ","
3108                                         << std::setprecision(4) << recopVtx.y << "," << std::setprecision(4)
3109                                         << recopVtx.z << ")"
3110                                         << " Recovtx.sim = " << recopVtx.sim;
3111   LogTrace("Primary4DVertexValidation") << "Reco vtx: pt = " << std::setprecision(4) << recopVtx.pt
3112                                         << " ptsq = " << std::setprecision(6) << recopVtx.ptsq
3113                                         << " nrecotrks = " << recopVtx.nRecoTrk
3114                                         << " nmatch simtrks = " << recopVtx.num_matched_sim_tracks;
3115   LogTrace("Primary4DVertexValidation") << "wnt " << recopVtx.sumwnt << " wos = " << recopVtx.sumwos;
3116   for (auto iTP = simpVtx.sim_vertex->daughterTracks_begin(); iTP != simpVtx.sim_vertex->daughterTracks_end(); ++iTP) {
3117     if (use_only_charged_tracks_ && (**iTP).charge() == 0) {
3118       continue;
3119     }
3120     LogTrace("Primary4DVertexValidation")
3121         << "Daughter track of sim vertex: pt =" << std::setw(6) << std::setprecision(2) << (*iTP)->pt()
3122         << "  eta =" << std::setw(6) << std::setprecision(2) << (*iTP)->eta();
3123   }
3124 }
3125 
3126 const bool Primary4DVertexValidation::trkTPSelLV(const TrackingParticle& tp) {
3127   bool match = false;
3128   if (tp.status() != 1) {
3129     return match;
3130   }
3131 
3132   auto x_pv = tp.parentVertex()->position().x();
3133   auto y_pv = tp.parentVertex()->position().y();
3134   auto z_pv = tp.parentVertex()->position().z();
3135 
3136   auto r_pv = std::sqrt(x_pv * x_pv + y_pv * y_pv);
3137 
3138   match =
3139       tp.charge() != 0 && std::abs(tp.eta()) < etacutGEN_ && tp.pt() > pTcut_ && r_pv < rBTL_ && std::abs(z_pv) < zETL_;
3140   return match;
3141 }
3142 
3143 const bool Primary4DVertexValidation::trkRecSel(const reco::TrackBase& trk) {
3144   bool match = false;
3145   match = std::abs(trk.eta()) <= etacutREC_ && trk.pt() > pTcut_;
3146   return match;
3147 }
3148 
3149 DEFINE_FWK_MODULE(Primary4DVertexValidation);