Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:32:47

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