File indexing completed on 2023-01-29 23:45:29
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
0010 #include "DataFormats/Common/interface/ValidHandle.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013 #include "FWCore/Framework/interface/MakerMacros.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016
0017 #include "DataFormats/Math/interface/LorentzVector.h"
0018 #include "DataFormats/Math/interface/Point3D.h"
0019
0020
0021 #include "DataFormats/TrackReco/interface/Track.h"
0022 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0023 #include "DataFormats/VertexReco/interface/Vertex.h"
0024 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0025 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0026
0027
0028 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0029 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
0030 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
0031 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
0032
0033
0034 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
0035
0036
0037 #include "SimTracker/VertexAssociation/interface/calculateVertexSharedTracks.h"
0038
0039
0040 #include "RecoVertex/PrimaryVertexProducer/interface/TrackFilterForPVFinding.h"
0041
0042
0043 #include "SimDataFormats/Associations/interface/VertexToTrackingVertexAssociator.h"
0044
0045
0046 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0047 #include "DQMServices/Core/interface/DQMStore.h"
0048
0049 #include "DataFormats/Math/interface/deltaR.h"
0050 #include "DataFormats/Math/interface/GeantUnits.h"
0051 #include "CLHEP/Units/PhysicalConstants.h"
0052
0053
0054 class Primary4DVertexValidation : public DQMEDAnalyzer {
0055 typedef math::XYZTLorentzVector LorentzVector;
0056
0057
0058 struct simPrimaryVertex {
0059 simPrimaryVertex(double x1, double y1, double z1, double t1)
0060 : x(x1),
0061 y(y1),
0062 z(z1),
0063 t(t1),
0064 ptsq(0),
0065 closest_vertex_distance_z(-1.),
0066 nGenTrk(0),
0067 num_matched_reco_tracks(0),
0068 average_match_quality(0.0) {
0069 ptot.setPx(0);
0070 ptot.setPy(0);
0071 ptot.setPz(0);
0072 ptot.setE(0);
0073 p4 = LorentzVector(0, 0, 0, 0);
0074 r = sqrt(x * x + y * y);
0075 };
0076 double x, y, z, r, t;
0077 HepMC::FourVector ptot;
0078 LorentzVector p4;
0079 double ptsq;
0080 double closest_vertex_distance_z;
0081 int nGenTrk;
0082 int num_matched_reco_tracks;
0083 float average_match_quality;
0084 EncodedEventId eventId;
0085 TrackingVertexRef sim_vertex;
0086 int OriginalIndex = -1;
0087
0088 unsigned int nwosmatch = 0;
0089 unsigned int nwntmatch = 0;
0090 std::vector<unsigned int> wos_dominated_recv;
0091
0092 std::map<unsigned int, double> wnt;
0093 std::map<unsigned int, double> wos;
0094 double sumwos = 0;
0095 double sumwnt = 0;
0096 unsigned int rec = NOT_MATCHED;
0097 unsigned int matchQuality = 0;
0098
0099 void addTrack(unsigned int irecv, double twos, double twt) {
0100 sumwnt += twt;
0101 if (wnt.find(irecv) == wnt.end()) {
0102 wnt[irecv] = twt;
0103 } else {
0104 wnt[irecv] += twt;
0105 }
0106
0107 sumwos += twos;
0108 if (wos.find(irecv) == wos.end()) {
0109 wos[irecv] = twos;
0110 } else {
0111 wos[irecv] += twos;
0112 }
0113 }
0114 };
0115
0116
0117 struct recoPrimaryVertex {
0118 recoPrimaryVertex(double x1, double y1, double z1)
0119 : x(x1),
0120 y(y1),
0121 z(z1),
0122 pt(0),
0123 ptsq(0),
0124 closest_vertex_distance_z(-1.),
0125 nRecoTrk(0),
0126 num_matched_sim_tracks(0),
0127 recVtx(nullptr) {
0128 r = sqrt(x * x + y * y);
0129 };
0130 double x, y, z, r;
0131 double pt;
0132 double ptsq;
0133 double closest_vertex_distance_z;
0134 int nRecoTrk;
0135 int num_matched_sim_tracks;
0136 const reco::Vertex* recVtx;
0137 reco::VertexBaseRef recVtxRef;
0138 int OriginalIndex = -1;
0139
0140 std::map<unsigned int, double> wos;
0141 std::map<unsigned int, double> wnt;
0142 unsigned int wosmatch;
0143 unsigned int wntmatch;
0144 double sumwos = 0;
0145 double sumwnt = 0;
0146 double maxwos = 0;
0147 double maxwnt = 0;
0148 int maxwosnt = 0;
0149 unsigned int sim = NOT_MATCHED;
0150 unsigned int matchQuality = 0;
0151
0152 bool is_real() { return (matchQuality > 0) && (matchQuality < 99); }
0153
0154 bool is_fake() { return (matchQuality <= 0) || (matchQuality >= 99); }
0155
0156 bool is_signal() { return (sim == 0); }
0157
0158 int split_from() {
0159 if (is_real())
0160 return -1;
0161 if ((maxwos > 0) && (maxwos > 0.3 * sumwos))
0162 return wosmatch;
0163 return -1;
0164 }
0165 bool other_fake() { return (is_fake() && (split_from() < 0)); }
0166
0167 void addTrack(unsigned int iev, double twos, double wt) {
0168 sumwnt += wt;
0169 if (wnt.find(iev) == wnt.end()) {
0170 wnt[iev] = wt;
0171 } else {
0172 wnt[iev] += wt;
0173 }
0174
0175 sumwos += twos;
0176 if (wos.find(iev) == wos.end()) {
0177 wos[iev] = twos;
0178 } else {
0179 wos[iev] += twos;
0180 }
0181 }
0182 };
0183
0184 public:
0185 explicit Primary4DVertexValidation(const edm::ParameterSet&);
0186 ~Primary4DVertexValidation() override;
0187
0188 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0189 void analyze(const edm::Event&, const edm::EventSetup&) override;
0190 void bookHistograms(DQMStore::IBooker& i, edm::Run const&, edm::EventSetup const&) override;
0191
0192 private:
0193 void matchReco2Sim(std::vector<recoPrimaryVertex>&,
0194 std::vector<simPrimaryVertex>&,
0195 const edm::ValueMap<float>&,
0196 const edm::ValueMap<float>&,
0197 const edm::Handle<reco::BeamSpot>&);
0198 bool matchRecoTrack2SimSignal(const reco::TrackBaseRef&);
0199 const edm::Ref<std::vector<TrackingParticle>>* getMatchedTP(const reco::TrackBaseRef&, const TrackingVertexRef&);
0200 double timeFromTrueMass(double, double, double, double);
0201 bool select(const reco::Vertex&, int level = 0);
0202 std::vector<Primary4DVertexValidation::simPrimaryVertex> getSimPVs(const edm::Handle<TrackingVertexCollection>&);
0203 std::vector<Primary4DVertexValidation::recoPrimaryVertex> getRecoPVs(const edm::Handle<edm::View<reco::Vertex>>&);
0204
0205 const bool mvaTPSel(const TrackingParticle&);
0206 const bool mvaRecSel(const reco::TrackBase&, const reco::Vertex&, const double&, const double&);
0207
0208
0209
0210 const std::string folder_;
0211 static constexpr unsigned int NOT_MATCHED = 66666;
0212 static constexpr double simUnit_ = 1e9;
0213 static constexpr double c_ = 2.99792458e1;
0214 static constexpr double mvaL_ = 0.5;
0215 static constexpr double mvaH_ = 0.8;
0216 static constexpr double selNdof_ = 4.;
0217 static constexpr double maxRank_ = 8.;
0218 static constexpr double maxTry_ = 10.;
0219 static constexpr double zWosMatchMax_ = 1.;
0220 static constexpr double etacutGEN_ = 4.;
0221 static constexpr double etacutREC_ = 3.;
0222 static constexpr double pTcut_ = 0.7;
0223 static constexpr double deltaZcut_ = 0.1;
0224 static constexpr double trackMaxBtlEta_ = 1.5;
0225 static constexpr double trackMinEtlEta_ = 1.6;
0226 static constexpr double trackMaxEtlEta_ = 3.;
0227 static constexpr double tol_ = 1.e-4;
0228
0229 static constexpr float c_cm_ns = geant_units::operators::convertMmToCm(CLHEP::c_light);
0230
0231 bool use_only_charged_tracks_;
0232 bool debug_;
0233 bool optionalPlots_;
0234
0235 const double minProbHeavy_;
0236 const double trackweightTh_;
0237 const double mvaTh_;
0238 const std::vector<double> lineDensityPar_;
0239 const reco::RecoToSimCollection* r2s_;
0240 const reco::SimToRecoCollection* s2r_;
0241
0242 edm::EDGetTokenT<reco::TrackCollection> RecTrackToken_;
0243
0244 edm::EDGetTokenT<std::vector<PileupSummaryInfo>> vecPileupSummaryInfoToken_;
0245
0246 edm::EDGetTokenT<TrackingParticleCollection> trackingParticleCollectionToken_;
0247 edm::EDGetTokenT<TrackingVertexCollection> trackingVertexCollectionToken_;
0248 edm::EDGetTokenT<reco::SimToRecoCollection> simToRecoAssociationToken_;
0249 edm::EDGetTokenT<reco::RecoToSimCollection> recoToSimAssociationToken_;
0250 edm::EDGetTokenT<reco::BeamSpot> RecBeamSpotToken_;
0251 edm::EDGetTokenT<edm::View<reco::Vertex>> Rec4DVerToken_;
0252
0253 edm::EDGetTokenT<edm::ValueMap<int>> trackAssocToken_;
0254 edm::EDGetTokenT<edm::ValueMap<float>> pathLengthToken_;
0255 edm::EDGetTokenT<edm::ValueMap<float>> momentumToken_;
0256 edm::EDGetTokenT<edm::ValueMap<float>> timeToken_;
0257
0258 edm::EDGetTokenT<edm::ValueMap<float>> t0PidToken_;
0259 edm::EDGetTokenT<edm::ValueMap<float>> t0SafePidToken_;
0260 edm::EDGetTokenT<edm::ValueMap<float>> sigmat0SafePidToken_;
0261 edm::EDGetTokenT<edm::ValueMap<float>> trackMVAQualToken_;
0262
0263 edm::EDGetTokenT<edm::ValueMap<float>> tmtdToken_;
0264 edm::EDGetTokenT<edm::ValueMap<float>> tofPiToken_;
0265 edm::EDGetTokenT<edm::ValueMap<float>> tofKToken_;
0266 edm::EDGetTokenT<edm::ValueMap<float>> tofPToken_;
0267 edm::EDGetTokenT<edm::ValueMap<float>> probPiToken_;
0268 edm::EDGetTokenT<edm::ValueMap<float>> probKToken_;
0269 edm::EDGetTokenT<edm::ValueMap<float>> probPToken_;
0270
0271
0272 MonitorElement* meMVATrackEffPtTot_;
0273 MonitorElement* meMVATrackMatchedEffPtTot_;
0274 MonitorElement* meMVATrackMatchedEffPtMtd_;
0275 MonitorElement* meMVATrackEffEtaTot_;
0276 MonitorElement* meMVATrackMatchedEffEtaTot_;
0277 MonitorElement* meMVATrackMatchedEffEtaMtd_;
0278 MonitorElement* meMVATrackResTot_;
0279 MonitorElement* meMVATrackPullTot_;
0280 MonitorElement* meTrackResTot_;
0281 MonitorElement* meTrackPullTot_;
0282 MonitorElement* meTrackRes_[3];
0283 MonitorElement* meTrackPull_[3];
0284 MonitorElement* meTrackResMass_[3];
0285 MonitorElement* meTrackResMassTrue_[3];
0286 MonitorElement* meMVATrackZposResTot_;
0287 MonitorElement* meTrackZposResTot_;
0288 MonitorElement* meTrackZposRes_[3];
0289 MonitorElement* meTrack3DposRes_[3];
0290 MonitorElement* meTimeRes_;
0291 MonitorElement* meTimePull_;
0292 MonitorElement* meTimeSignalRes_;
0293 MonitorElement* meTimeSignalPull_;
0294 MonitorElement* mePUvsRealV_;
0295 MonitorElement* mePUvsOtherFakeV_;
0296 MonitorElement* mePUvsSplitV_;
0297 MonitorElement* meMatchQual_;
0298 MonitorElement* meDeltaZrealreal_;
0299 MonitorElement* meDeltaZfakefake_;
0300 MonitorElement* meDeltaZfakereal_;
0301 MonitorElement* meDeltaTrealreal_;
0302 MonitorElement* meDeltaTfakefake_;
0303 MonitorElement* meDeltaTfakereal_;
0304 MonitorElement* meRecoPosInSimCollection_;
0305 MonitorElement* meRecoPosInRecoOrigCollection_;
0306 MonitorElement* meSimPosInSimOrigCollection_;
0307 MonitorElement* meRecoPVPosSignal_;
0308 MonitorElement* meRecoPVPosSignalNotHighestPt_;
0309 MonitorElement* meRecoVtxVsLineDensity_;
0310 MonitorElement* meRecVerNumber_;
0311 MonitorElement* meRecPVZ_;
0312 MonitorElement* meRecPVT_;
0313 MonitorElement* meSimPVZ_;
0314
0315
0316 MonitorElement* meTrackResLowPTot_;
0317 MonitorElement* meTrackResHighPTot_;
0318 MonitorElement* meTrackPullLowPTot_;
0319 MonitorElement* meTrackPullHighPTot_;
0320
0321 MonitorElement* meTrackResLowP_[3];
0322 MonitorElement* meTrackResHighP_[3];
0323 MonitorElement* meTrackPullLowP_[3];
0324 MonitorElement* meTrackPullHighP_[3];
0325
0326 MonitorElement* meTrackResMassProtons_[3];
0327 MonitorElement* meTrackResMassTrueProtons_[3];
0328 MonitorElement* meTrackResMassPions_[3];
0329 MonitorElement* meTrackResMassTruePions_[3];
0330
0331 MonitorElement* meBarrelPIDp_;
0332 MonitorElement* meEndcapPIDp_;
0333
0334 MonitorElement* meBarrelNoPIDtype_;
0335 MonitorElement* meEndcapNoPIDtype_;
0336
0337 MonitorElement* meBarrelTruePiNoPID_;
0338 MonitorElement* meBarrelTrueKNoPID_;
0339 MonitorElement* meBarrelTruePNoPID_;
0340 MonitorElement* meEndcapTruePiNoPID_;
0341 MonitorElement* meEndcapTrueKNoPID_;
0342 MonitorElement* meEndcapTruePNoPID_;
0343
0344 MonitorElement* meBarrelTruePiAsPi_;
0345 MonitorElement* meBarrelTruePiAsK_;
0346 MonitorElement* meBarrelTruePiAsP_;
0347 MonitorElement* meEndcapTruePiAsPi_;
0348 MonitorElement* meEndcapTruePiAsK_;
0349 MonitorElement* meEndcapTruePiAsP_;
0350
0351 MonitorElement* meBarrelTrueKAsPi_;
0352 MonitorElement* meBarrelTrueKAsK_;
0353 MonitorElement* meBarrelTrueKAsP_;
0354 MonitorElement* meEndcapTrueKAsPi_;
0355 MonitorElement* meEndcapTrueKAsK_;
0356 MonitorElement* meEndcapTrueKAsP_;
0357
0358 MonitorElement* meBarrelTruePAsPi_;
0359 MonitorElement* meBarrelTruePAsK_;
0360 MonitorElement* meBarrelTruePAsP_;
0361 MonitorElement* meEndcapTruePAsPi_;
0362 MonitorElement* meEndcapTruePAsK_;
0363 MonitorElement* meEndcapTruePAsP_;
0364 };
0365
0366
0367 Primary4DVertexValidation::Primary4DVertexValidation(const edm::ParameterSet& iConfig)
0368 : folder_(iConfig.getParameter<std::string>("folder")),
0369 use_only_charged_tracks_(iConfig.getParameter<bool>("useOnlyChargedTracks")),
0370 debug_(iConfig.getUntrackedParameter<bool>("debug")),
0371 optionalPlots_(iConfig.getUntrackedParameter<bool>("optionalPlots")),
0372 minProbHeavy_(iConfig.getParameter<double>("minProbHeavy")),
0373 trackweightTh_(iConfig.getParameter<double>("trackweightTh")),
0374 mvaTh_(iConfig.getParameter<double>("mvaTh")),
0375 lineDensityPar_(iConfig.getParameter<std::vector<double>>("lineDensityPar")) {
0376 vecPileupSummaryInfoToken_ = consumes<std::vector<PileupSummaryInfo>>(edm::InputTag(std::string("addPileupInfo")));
0377 trackingParticleCollectionToken_ =
0378 consumes<TrackingParticleCollection>(iConfig.getParameter<edm::InputTag>("SimTag"));
0379 trackingVertexCollectionToken_ = consumes<TrackingVertexCollection>(iConfig.getParameter<edm::InputTag>("SimTag"));
0380 simToRecoAssociationToken_ =
0381 consumes<reco::SimToRecoCollection>(iConfig.getParameter<edm::InputTag>("TPtoRecoTrackAssoc"));
0382 recoToSimAssociationToken_ =
0383 consumes<reco::RecoToSimCollection>(iConfig.getParameter<edm::InputTag>("TPtoRecoTrackAssoc"));
0384 RecTrackToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("mtdTracks"));
0385 RecBeamSpotToken_ = consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("offlineBS"));
0386 Rec4DVerToken_ = consumes<edm::View<reco::Vertex>>(iConfig.getParameter<edm::InputTag>("offline4DPV"));
0387 trackAssocToken_ = consumes<edm::ValueMap<int>>(iConfig.getParameter<edm::InputTag>("trackAssocSrc"));
0388 pathLengthToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("pathLengthSrc"));
0389 momentumToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("momentumSrc"));
0390 timeToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("timeSrc"));
0391 t0PidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("t0PID"));
0392 t0SafePidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("t0SafePID"));
0393 sigmat0SafePidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmat0SafePID"));
0394 trackMVAQualToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("trackMVAQual"));
0395 tmtdToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tmtd"));
0396 tofPiToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tofPi"));
0397 tofKToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tofK"));
0398 tofPToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tofP"));
0399 probPiToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("probPi"));
0400 probKToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("probK"));
0401 probPToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("probP"));
0402 }
0403
0404 Primary4DVertexValidation::~Primary4DVertexValidation() {}
0405
0406
0407
0408
0409 void Primary4DVertexValidation::bookHistograms(DQMStore::IBooker& ibook,
0410 edm::Run const& iRun,
0411 edm::EventSetup const& iSetup) {
0412 ibook.setCurrentFolder(folder_);
0413
0414 meMVATrackEffPtTot_ = ibook.book1D("MVAEffPtTot", "Pt of tracks associated to LV; track pt [GeV] ", 110, 0., 11.);
0415 meMVATrackEffEtaTot_ = ibook.book1D("MVAEffEtaTot", "Pt of tracks associated to LV; track eta ", 66, 0., 3.3);
0416 meMVATrackMatchedEffPtTot_ =
0417 ibook.book1D("MVAMatchedEffPtTot", "Pt of tracks associated to LV matched to TP; track pt [GeV] ", 110, 0., 11.);
0418 meMVATrackMatchedEffPtMtd_ = ibook.book1D(
0419 "MVAMatchedEffPtMtd", "Pt of tracks associated to LV matched to TP with time; track pt [GeV] ", 110, 0., 11.);
0420 meMVATrackMatchedEffEtaTot_ =
0421 ibook.book1D("MVAMatchedEffEtaTot", "Pt of tracks associated to LV matched to TP; track eta ", 66, 0., 3.3);
0422 meMVATrackMatchedEffEtaMtd_ = ibook.book1D(
0423 "MVAMatchedEffEtaMtd", "Pt of tracks associated to LV matched to TP with time; track eta ", 66, 0., 3.3);
0424 meMVATrackResTot_ = ibook.book1D(
0425 "MVATrackRes", "t_{rec} - t_{sim} for tracks from LV MVA sel.; t_{rec} - t_{sim} [ns] ", 120, -0.15, 0.15);
0426 meTrackResTot_ = ibook.book1D("TrackRes", "t_{rec} - t_{sim} for tracks; t_{rec} - t_{sim} [ns] ", 120, -0.15, 0.15);
0427 meTrackRes_[0] = ibook.book1D(
0428 "TrackRes-LowMVA", "t_{rec} - t_{sim} for tracks with MVA < 0.5; t_{rec} - t_{sim} [ns] ", 100, -1., 1.);
0429 meTrackRes_[1] = ibook.book1D(
0430 "TrackRes-MediumMVA", "t_{rec} - t_{sim} for tracks with 0.5 < MVA < 0.8; t_{rec} - t_{sim} [ns] ", 100, -1., 1.);
0431 meTrackRes_[2] = ibook.book1D(
0432 "TrackRes-HighMVA", "t_{rec} - t_{sim} for tracks with MVA > 0.8; t_{rec} - t_{sim} [ns] ", 100, -1., 1.);
0433 if (optionalPlots_) {
0434 meTrackResMass_[0] = ibook.book1D(
0435 "TrackResMass-LowMVA", "t_{rec} - t_{est} for tracks with MVA < 0.5; t_{rec} - t_{est} [ns] ", 100, -1., 1.);
0436 meTrackResMass_[1] = ibook.book1D("TrackResMass-MediumMVA",
0437 "t_{rec} - t_{est} for tracks with 0.5 < MVA < 0.8; t_{rec} - t_{est} [ns] ",
0438 100,
0439 -1.,
0440 1.);
0441 meTrackResMass_[2] = ibook.book1D(
0442 "TrackResMass-HighMVA", "t_{rec} - t_{est} for tracks with MVA > 0.8; t_{rec} - t_{est} [ns] ", 100, -1., 1.);
0443 meTrackResMassTrue_[0] = ibook.book1D(
0444 "TrackResMassTrue-LowMVA", "t_{est} - t_{sim} for tracks with MVA < 0.5; t_{est} - t_{sim} [ns] ", 100, -1., 1.);
0445 meTrackResMassTrue_[1] = ibook.book1D("TrackResMassTrue-MediumMVA",
0446 "t_{est} - t_{sim} for tracks with 0.5 < MVA < 0.8; t_{est} - t_{sim} [ns] ",
0447 100,
0448 -1.,
0449 1.);
0450 meTrackResMassTrue_[2] = ibook.book1D("TrackResMassTrue-HighMVA",
0451 "t_{est} - t_{sim} for tracks with MVA > 0.8; t_{est} - t_{sim} [ns] ",
0452 100,
0453 -1.,
0454 1.);
0455 }
0456 meMVATrackPullTot_ =
0457 ibook.book1D("MVATrackPull", "Pull for tracks from LV MAV sel.; (t_{rec}-t_{sim})/#sigma_{t}", 50, -5., 5.);
0458 meTrackPullTot_ = ibook.book1D("TrackPull", "Pull for tracks; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
0459 meTrackPull_[0] =
0460 ibook.book1D("TrackPull-LowMVA", "Pull for tracks with MVA < 0.5; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
0461 meTrackPull_[1] = ibook.book1D(
0462 "TrackPull-MediumMVA", "Pull for tracks with 0.5 < MVA < 0.8; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
0463 meTrackPull_[2] =
0464 ibook.book1D("TrackPull-HighMVA", "Pull for tracks with MVA > 0.8; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
0465 meMVATrackZposResTot_ = ibook.book1D(
0466 "MVATrackZposResTot", "Z_{PCA} - Z_{sim} for tracks from LV MVA sel.;Z_{PCA} - Z_{sim} [cm] ", 50, -0.1, 0.1);
0467 meTrackZposResTot_ =
0468 ibook.book1D("TrackZposResTot", "Z_{PCA} - Z_{sim} for tracks;Z_{PCA} - Z_{sim} [cm] ", 50, -0.5, 0.5);
0469 meTrackZposRes_[0] = ibook.book1D(
0470 "TrackZposRes-LowMVA", "Z_{PCA} - Z_{sim} for tracks with MVA < 0.5;Z_{PCA} - Z_{sim} [cm] ", 50, -0.5, 0.5);
0471 meTrackZposRes_[1] = ibook.book1D("TrackZposRes-MediumMVA",
0472 "Z_{PCA} - Z_{sim} for tracks with 0.5 < MVA < 0.8 ;Z_{PCA} - Z_{sim} [cm] ",
0473 50,
0474 -0.5,
0475 0.5);
0476 meTrackZposRes_[2] = ibook.book1D(
0477 "TrackZposRes-HighMVA", "Z_{PCA} - Z_{sim} for tracks with MVA > 0.8 ;Z_{PCA} - Z_{sim} [cm] ", 50, -0.5, 0.5);
0478 meTrack3DposRes_[0] =
0479 ibook.book1D("Track3DposRes-LowMVA",
0480 "3dPos_{PCA} - 3dPos_{sim} for tracks with MVA < 0.5 ;3dPos_{PCA} - 3dPos_{sim} [cm] ",
0481 50,
0482 -0.5,
0483 0.5);
0484 meTrack3DposRes_[1] =
0485 ibook.book1D("Track3DposRes-MediumMVA",
0486 "3dPos_{PCA} - 3dPos_{sim} for tracks with 0.5 < MVA < 0.8 ;3dPos_{PCA} - 3dPos_{sim} [cm] ",
0487 50,
0488 -0.5,
0489 0.5);
0490 meTrack3DposRes_[2] =
0491 ibook.book1D("Track3DposRes-HighMVA",
0492 "3dPos_{PCA} - 3dPos_{sim} for tracks with MVA > 0.8;3dPos_{PCA} - 3dPos_{sim} [cm] ",
0493 50,
0494 -0.5,
0495 0.5);
0496 meTimeRes_ = ibook.book1D("TimeRes", "t_{rec} - t_{sim} ;t_{rec} - t_{sim} [ns] ", 40, -0.2, 0.2);
0497 meTimePull_ = ibook.book1D("TimePull", "Pull; t_{rec} - t_{sim}/#sigma_{t rec}", 100, -10., 10.);
0498 meTimeSignalRes_ =
0499 ibook.book1D("TimeSignalRes", "t_{rec} - t_{sim} for signal ;t_{rec} - t_{sim} [ns] ", 40, -0.2, 0.2);
0500 meTimeSignalPull_ =
0501 ibook.book1D("TimeSignalPull", "Pull for signal; t_{rec} - t_{sim}/#sigma_{t rec}", 100, -10., 10.);
0502 mePUvsRealV_ =
0503 ibook.bookProfile("PUvsReal", "#PU vertices vs #real matched vertices;#PU;#real ", 100, 0, 300, 100, 0, 200);
0504 mePUvsOtherFakeV_ = ibook.bookProfile(
0505 "PUvsOtherFake", "#PU vertices vs #other fake matched vertices;#PU;#other fake ", 100, 0, 300, 100, 0, 20);
0506 mePUvsSplitV_ =
0507 ibook.bookProfile("PUvsSplit", "#PU vertices vs #split matched vertices;#PU;#split ", 100, 0, 300, 100, 0, 20);
0508 meMatchQual_ = ibook.book1D("MatchQuality", "RECO-SIM vertex match quality; ", 8, 0, 8.);
0509 meDeltaZrealreal_ = ibook.book1D("DeltaZrealreal", "#Delta Z real-real; |#Delta Z (r-r)| [cm]", 100, 0, 0.5);
0510 meDeltaZfakefake_ = ibook.book1D("DeltaZfakefake", "#Delta Z fake-fake; |#Delta Z (f-f)| [cm]", 100, 0, 0.5);
0511 meDeltaZfakereal_ = ibook.book1D("DeltaZfakereal", "#Delta Z fake-real; |#Delta Z (f-r)| [cm]", 100, 0, 0.5);
0512 meDeltaTrealreal_ = ibook.book1D("DeltaTrealreal", "#Delta T real-real; |#Delta T (r-r)| [sigma]", 60, 0., 30.);
0513 meDeltaTfakefake_ = ibook.book1D("DeltaTfakefake", "#Delta T fake-fake; |#Delta T (f-f)| [sigma]", 60, 0., 30.);
0514 meDeltaTfakereal_ = ibook.book1D("DeltaTfakereal", "#Delta T fake-real; |#Delta T (f-r)| [sigma]", 60, 0., 30.);
0515 if (optionalPlots_) {
0516 meRecoPosInSimCollection_ = ibook.book1D(
0517 "RecoPosInSimCollection", "Sim signal vertex index associated to Reco signal vertex; Sim PV index", 200, 0, 200);
0518 meRecoPosInRecoOrigCollection_ =
0519 ibook.book1D("RecoPosInRecoOrigCollection", "Reco signal index in OrigCollection; Reco index", 200, 0, 200);
0520 meSimPosInSimOrigCollection_ =
0521 ibook.book1D("SimPosInSimOrigCollection", "Sim signal index in OrigCollection; Sim index", 200, 0, 200);
0522 }
0523 meRecoPVPosSignal_ =
0524 ibook.book1D("RecoPVPosSignal", "Position in reco collection of PV associated to sim signal", 200, 0, 200);
0525 meRecoPVPosSignalNotHighestPt_ =
0526 ibook.book1D("RecoPVPosSignalNotHighestPt",
0527 "Position in reco collection of PV associated to sim signal not highest Pt",
0528 200,
0529 0,
0530 200);
0531 meRecoVtxVsLineDensity_ =
0532 ibook.book1D("RecoVtxVsLineDensity", "#Reco vertices/mm/event; line density [#vtx/mm/event]", 160, 0., 4.);
0533 meRecVerNumber_ = ibook.book1D("RecVerNumber", "RECO Vertex Number: Number of vertices", 50, 0, 250);
0534 meRecPVZ_ = ibook.book1D("recPVZ", "Weighted #Rec vertices/mm", 400, -20., 20.);
0535 meRecPVT_ = ibook.book1D("recPVT", "#Rec vertices/10 ps", 200, -1., 1.);
0536 meSimPVZ_ = ibook.book1D("simPVZ", "Weighted #Sim vertices/mm", 400, -20., 20.);
0537
0538
0539 meTrackResLowPTot_ = ibook.book1D(
0540 "TrackResLowP", "t_{rec} - t_{sim} for tracks with p < 2 GeV; t_{rec} - t_{sim} [ns] ", 70, -0.15, 0.15);
0541 meTrackResLowP_[0] =
0542 ibook.book1D("TrackResLowP-LowMVA",
0543 "t_{rec} - t_{sim} for tracks with MVA < 0.5 and p < 2 GeV; t_{rec} - t_{sim} [ns] ",
0544 100,
0545 -1.,
0546 1.);
0547 meTrackResLowP_[1] =
0548 ibook.book1D("TrackResLowP-MediumMVA",
0549 "t_{rec} - t_{sim} for tracks with 0.5 < MVA < 0.8 and p < 2 GeV; t_{rec} - t_{sim} [ns] ",
0550 100,
0551 -1.,
0552 1.);
0553 meTrackResLowP_[2] =
0554 ibook.book1D("TrackResLowP-HighMVA",
0555 "t_{rec} - t_{sim} for tracks with MVA > 0.8 and p < 2 GeV; t_{rec} - t_{sim} [ns] ",
0556 100,
0557 -1.,
0558 1.);
0559 meTrackResHighPTot_ = ibook.book1D(
0560 "TrackResHighP", "t_{rec} - t_{sim} for tracks with p > 2 GeV; t_{rec} - t_{sim} [ns] ", 70, -0.15, 0.15);
0561 meTrackResHighP_[0] =
0562 ibook.book1D("TrackResHighP-LowMVA",
0563 "t_{rec} - t_{sim} for tracks with MVA < 0.5 and p > 2 GeV; t_{rec} - t_{sim} [ns] ",
0564 100,
0565 -1.,
0566 1.);
0567 meTrackResHighP_[1] =
0568 ibook.book1D("TrackResHighP-MediumMVA",
0569 "t_{rec} - t_{sim} for tracks with 0.5 < MVA < 0.8 and p > 2 GeV; t_{rec} - t_{sim} [ns] ",
0570 100,
0571 -1.,
0572 1.);
0573 meTrackResHighP_[2] =
0574 ibook.book1D("TrackResHighP-HighMVA",
0575 "t_{rec} - t_{sim} for tracks with MVA > 0.8 and p > 2 GeV; t_{rec} - t_{sim} [ns] ",
0576 100,
0577 -1.,
0578 1.);
0579 meTrackPullLowPTot_ =
0580 ibook.book1D("TrackPullLowP", "Pull for tracks with p < 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
0581 meTrackPullLowP_[0] = ibook.book1D("TrackPullLowP-LowMVA",
0582 "Pull for tracks with MVA < 0.5 and p < 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
0583 100,
0584 -10.,
0585 10.);
0586 meTrackPullLowP_[1] = ibook.book1D("TrackPullLowP-MediumMVA",
0587 "Pull for tracks with 0.5 < MVA < 0.8 and p < 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
0588 100,
0589 -10.,
0590 10.);
0591 meTrackPullLowP_[2] = ibook.book1D("TrackPullLowP-HighMVA",
0592 "Pull for tracks with MVA > 0.8 and p < 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
0593 100,
0594 -10.,
0595 10.);
0596 meTrackPullHighPTot_ =
0597 ibook.book1D("TrackPullHighP", "Pull for tracks with p > 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
0598 meTrackPullHighP_[0] = ibook.book1D("TrackPullHighP-LowMVA",
0599 "Pull for tracks with MVA < 0.5 and p > 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
0600 100,
0601 -10.,
0602 10.);
0603 meTrackPullHighP_[1] =
0604 ibook.book1D("TrackPullHighP-MediumMVA",
0605 "Pull for tracks with 0.5 < MVA < 0.8 and p > 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
0606 100,
0607 -10.,
0608 10.);
0609 meTrackPullHighP_[2] = ibook.book1D("TrackPullHighP-HighMVA",
0610 "Pull for tracks with MVA > 0.8 and p > 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
0611 100,
0612 -10.,
0613 10.);
0614 if (optionalPlots_) {
0615 meTrackResMassProtons_[0] =
0616 ibook.book1D("TrackResMass-Protons-LowMVA",
0617 "t_{rec} - t_{est} for proton tracks with MVA < 0.5; t_{rec} - t_{est} [ns] ",
0618 100,
0619 -1.,
0620 1.);
0621 meTrackResMassProtons_[1] =
0622 ibook.book1D("TrackResMass-Protons-MediumMVA",
0623 "t_{rec} - t_{est} for proton tracks with 0.5 < MVA < 0.8; t_{rec} - t_{est} [ns] ",
0624 100,
0625 -1.,
0626 1.);
0627 meTrackResMassProtons_[2] =
0628 ibook.book1D("TrackResMass-Protons-HighMVA",
0629 "t_{rec} - t_{est} for proton tracks with MVA > 0.8; t_{rec} - t_{est} [ns] ",
0630 100,
0631 -1.,
0632 1.);
0633 meTrackResMassTrueProtons_[0] =
0634 ibook.book1D("TrackResMassTrue-Protons-LowMVA",
0635 "t_{est} - t_{sim} for proton tracks with MVA < 0.5; t_{est} - t_{sim} [ns] ",
0636 100,
0637 -1.,
0638 1.);
0639 meTrackResMassTrueProtons_[1] =
0640 ibook.book1D("TrackResMassTrue-Protons-MediumMVA",
0641 "t_{est} - t_{sim} for proton tracks with 0.5 < MVA < 0.8; t_{est} - t_{sim} [ns] ",
0642 100,
0643 -1.,
0644 1.);
0645 meTrackResMassTrueProtons_[2] =
0646 ibook.book1D("TrackResMassTrue-Protons-HighMVA",
0647 "t_{est} - t_{sim} for proton tracks with MVA > 0.8; t_{est} - t_{sim} [ns] ",
0648 100,
0649 -1.,
0650 1.);
0651
0652 meTrackResMassPions_[0] = ibook.book1D("TrackResMass-Pions-LowMVA",
0653 "t_{rec} - t_{est} for pion tracks with MVA < 0.5; t_{rec} - t_{est} [ns] ",
0654 100,
0655 -1.,
0656 1.);
0657 meTrackResMassPions_[1] =
0658 ibook.book1D("TrackResMass-Pions-MediumMVA",
0659 "t_{rec} - t_{est} for pion tracks with 0.5 < MVA < 0.8; t_{rec} - t_{est} [ns] ",
0660 100,
0661 -1.,
0662 1.);
0663 meTrackResMassPions_[2] = ibook.book1D("TrackResMass-Pions-HighMVA",
0664 "t_{rec} - t_{est} for pion tracks with MVA > 0.8; t_{rec} - t_{est} [ns] ",
0665 100,
0666 -1.,
0667 1.);
0668 meTrackResMassTruePions_[0] =
0669 ibook.book1D("TrackResMassTrue-Pions-LowMVA",
0670 "t_{est} - t_{sim} for pion tracks with MVA < 0.5; t_{est} - t_{sim} [ns] ",
0671 100,
0672 -1.,
0673 1.);
0674 meTrackResMassTruePions_[1] =
0675 ibook.book1D("TrackResMassTrue-Pions-MediumMVA",
0676 "t_{est} - t_{sim} for pion tracks with 0.5 < MVA < 0.8; t_{est} - t_{sim} [ns] ",
0677 100,
0678 -1.,
0679 1.);
0680 meTrackResMassTruePions_[2] =
0681 ibook.book1D("TrackResMassTrue-Pions-HighMVA",
0682 "t_{est} - t_{sim} for pion tracks with MVA > 0.8; t_{est} - t_{sim} [ns] ",
0683 100,
0684 -1.,
0685 1.);
0686 }
0687
0688 meBarrelPIDp_ = ibook.book1D("BarrelPIDp", "PID track MTD momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0689 meEndcapPIDp_ = ibook.book1D("EndcapPIDp", "PID track with MTD momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0690
0691 meBarrelNoPIDtype_ = ibook.book1D("BarrelNoPIDtype", "Barrel PID failure category", 4, 0., 4.);
0692 meEndcapNoPIDtype_ = ibook.book1D("EndcapNoPIDtype", "Endcap PID failure category", 4, 0., 4.);
0693
0694 meBarrelNoPIDtype_ = ibook.book1D("BarrelNoPIDtype", "Barrel PID failure category", 4, 0., 4.);
0695 meEndcapNoPIDtype_ = ibook.book1D("EndcapNoPIDtype", "Endcap PID failure category", 4, 0., 4.);
0696
0697 meBarrelTruePiNoPID_ =
0698 ibook.book1D("BarrelTruePiNoPID", "True pi NoPID momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0699 meBarrelTrueKNoPID_ =
0700 ibook.book1D("BarrelTrueKNoPID", "True k NoPID momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0701 meBarrelTruePNoPID_ =
0702 ibook.book1D("BarrelTruePNoPID", "True p NoPID momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0703 meEndcapTruePiNoPID_ =
0704 ibook.book1D("EndcapTruePiNoPID", "True NoPIDpi momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0705 meEndcapTrueKNoPID_ =
0706 ibook.book1D("EndcapTrueKNoPID", "True k NoPID momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0707 meEndcapTruePNoPID_ =
0708 ibook.book1D("EndcapTruePNoPID", "True p NoPID momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0709
0710 meBarrelTruePiAsPi_ =
0711 ibook.book1D("BarrelTruePiAsPi", "True pi as pi momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0712 meBarrelTruePiAsK_ =
0713 ibook.book1D("BarrelTruePiAsK", "True pi as k momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0714 meBarrelTruePiAsP_ =
0715 ibook.book1D("BarrelTruePiAsP", "True pi as p momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0716 meEndcapTruePiAsPi_ =
0717 ibook.book1D("EndcapTruePiAsPi", "True pi as pi momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0718 meEndcapTruePiAsK_ =
0719 ibook.book1D("EndcapTruePiAsK", "True pi as k momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0720 meEndcapTruePiAsP_ =
0721 ibook.book1D("EndcapTruePiAsP", "True pi as p momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0722
0723 meBarrelTrueKAsPi_ =
0724 ibook.book1D("BarrelTrueKAsPi", "True k as pi momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0725 meBarrelTrueKAsK_ = ibook.book1D("BarrelTrueKAsK", "True k as k momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0726 meBarrelTrueKAsP_ = ibook.book1D("BarrelTrueKAsP", "True k as p momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0727 meEndcapTrueKAsPi_ =
0728 ibook.book1D("EndcapTrueKAsPi", "True k as pi momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0729 meEndcapTrueKAsK_ = ibook.book1D("EndcapTrueKAsK", "True k as k momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0730 meEndcapTrueKAsP_ = ibook.book1D("EndcapTrueKAsP", "True k as p momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0731
0732 meBarrelTruePAsPi_ =
0733 ibook.book1D("BarrelTruePAsPi", "True p as pi momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0734 meBarrelTruePAsK_ = ibook.book1D("BarrelTruePAsK", "True p as k momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0735 meBarrelTruePAsP_ = ibook.book1D("BarrelTruePAsP", "True p as p momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0736 meEndcapTruePAsPi_ =
0737 ibook.book1D("EndcapTruePAsPi", "True p as pi momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0738 meEndcapTruePAsK_ = ibook.book1D("EndcapTruePAsK", "True p as k momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0739 meEndcapTruePAsP_ = ibook.book1D("EndcapTruePAsP", "True p as p momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0740 }
0741
0742 bool Primary4DVertexValidation::matchRecoTrack2SimSignal(const reco::TrackBaseRef& recoTrack) {
0743 auto found = r2s_->find(recoTrack);
0744
0745
0746 if (found == r2s_->end())
0747 return false;
0748
0749
0750 for (const auto& tp : found->val) {
0751 if (tp.first->eventId().bunchCrossing() == 0 && tp.first->eventId().event() == 0)
0752 return true;
0753 }
0754
0755
0756 return false;
0757 }
0758
0759 const edm::Ref<std::vector<TrackingParticle>>* Primary4DVertexValidation::getMatchedTP(
0760 const reco::TrackBaseRef& recoTrack, const TrackingVertexRef& vsim) {
0761 auto found = r2s_->find(recoTrack);
0762
0763
0764 if (found == r2s_->end())
0765 return nullptr;
0766
0767
0768 for (const auto& tp : found->val) {
0769 if (std::find_if(vsim->daughterTracks_begin(), vsim->daughterTracks_end(), [&](const TrackingParticleRef& vtp) {
0770 return tp.first == vtp;
0771 }) != vsim->daughterTracks_end())
0772 return &tp.first;
0773 }
0774
0775
0776 return nullptr;
0777 }
0778
0779 double Primary4DVertexValidation::timeFromTrueMass(double mass, double pathlength, double momentum, double time) {
0780 if (time > 0 && pathlength > 0 && mass > 0) {
0781 double gammasq = 1. + momentum * momentum / (mass * mass);
0782 double v = c_ * std::sqrt(1. - 1. / gammasq);
0783 double t_est = time - (pathlength / v);
0784
0785 return t_est;
0786 } else {
0787 return -1;
0788 }
0789 }
0790
0791 bool Primary4DVertexValidation::select(const reco::Vertex& v, int level) {
0792
0793
0794
0795
0796
0797 if (v.isFake())
0798 return false;
0799 if ((level == 0) && (v.ndof() > selNdof_))
0800 return true;
0801
0802
0803
0804
0805
0806
0807 return false;
0808 }
0809
0810
0811
0812
0813 std::vector<Primary4DVertexValidation::simPrimaryVertex> Primary4DVertexValidation::getSimPVs(
0814 const edm::Handle<TrackingVertexCollection>& tVC) {
0815 std::vector<Primary4DVertexValidation::simPrimaryVertex> simpv;
0816 int current_event = -1;
0817 int s = -1;
0818 for (TrackingVertexCollection::const_iterator v = tVC->begin(); v != tVC->end(); ++v) {
0819
0820 if (v->eventId().bunchCrossing() != 0)
0821 continue;
0822 if (v->eventId().event() != current_event) {
0823 current_event = v->eventId().event();
0824 } else {
0825 continue;
0826 }
0827 s++;
0828 if (std::abs(v->position().z()) > 1000)
0829 continue;
0830
0831
0832 simPrimaryVertex sv(v->position().x(), v->position().y(), v->position().z(), v->position().t());
0833 sv.eventId = v->eventId();
0834 sv.sim_vertex = TrackingVertexRef(tVC, std::distance(tVC->begin(), v));
0835 sv.OriginalIndex = s;
0836
0837 for (TrackingParticleRefVector::iterator iTrack = v->daughterTracks_begin(); iTrack != v->daughterTracks_end();
0838 ++iTrack) {
0839 assert((**iTrack).eventId().bunchCrossing() == 0);
0840 }
0841 simPrimaryVertex* vp = nullptr;
0842 for (std::vector<simPrimaryVertex>::iterator v0 = simpv.begin(); v0 != simpv.end(); v0++) {
0843 if ((sv.eventId == v0->eventId) && (std::abs(sv.x - v0->x) < 1e-5) && (std::abs(sv.y - v0->y) < 1e-5) &&
0844 (std::abs(sv.z - v0->z) < 1e-5)) {
0845 vp = &(*v0);
0846 break;
0847 }
0848 }
0849 if (!vp) {
0850
0851 simpv.push_back(sv);
0852 vp = &simpv.back();
0853 }
0854
0855
0856 for (TrackingVertex::tp_iterator iTP = v->daughterTracks_begin(); iTP != v->daughterTracks_end(); ++iTP) {
0857 auto momentum = (*(*iTP)).momentum();
0858 const reco::Track* matched_best_reco_track = nullptr;
0859 double match_quality = -1;
0860 if (use_only_charged_tracks_ && (**iTP).charge() == 0)
0861 continue;
0862 if (s2r_->find(*iTP) != s2r_->end()) {
0863 matched_best_reco_track = (*s2r_)[*iTP][0].first.get();
0864 match_quality = (*s2r_)[*iTP][0].second;
0865 }
0866
0867 vp->ptot.setPx(vp->ptot.x() + momentum.x());
0868 vp->ptot.setPy(vp->ptot.y() + momentum.y());
0869 vp->ptot.setPz(vp->ptot.z() + momentum.z());
0870 vp->ptot.setE(vp->ptot.e() + (**iTP).energy());
0871 vp->ptsq += ((**iTP).pt() * (**iTP).pt());
0872
0873 if (matched_best_reco_track) {
0874 vp->num_matched_reco_tracks++;
0875 vp->average_match_quality += match_quality;
0876 }
0877 }
0878 if (vp->num_matched_reco_tracks)
0879 vp->average_match_quality /= static_cast<float>(vp->num_matched_reco_tracks);
0880 if (debug_) {
0881 edm::LogPrint("Primary4DVertexValidation")
0882 << "average number of associated tracks: " << vp->num_matched_reco_tracks / static_cast<float>(vp->nGenTrk)
0883 << " with average quality: " << vp->average_match_quality;
0884 }
0885 }
0886
0887
0888 if (simpv.empty())
0889 return simpv;
0890
0891
0892
0893 auto prev_z = simpv.back().z;
0894 for (simPrimaryVertex& vsim : simpv) {
0895 vsim.closest_vertex_distance_z = std::abs(vsim.z - prev_z);
0896 prev_z = vsim.z;
0897 }
0898
0899 for (std::vector<simPrimaryVertex>::iterator vsim = simpv.begin(); vsim != simpv.end(); vsim++) {
0900 std::vector<simPrimaryVertex>::iterator vsim2 = vsim;
0901 vsim2++;
0902 for (; vsim2 != simpv.end(); vsim2++) {
0903 double distance = std::abs(vsim->z - vsim2->z);
0904
0905 vsim->closest_vertex_distance_z = std::min(vsim->closest_vertex_distance_z, distance);
0906 vsim2->closest_vertex_distance_z = std::min(vsim2->closest_vertex_distance_z, distance);
0907 }
0908 }
0909 return simpv;
0910 }
0911
0912
0913
0914 std::vector<Primary4DVertexValidation::recoPrimaryVertex> Primary4DVertexValidation::getRecoPVs(
0915 const edm::Handle<edm::View<reco::Vertex>>& tVC) {
0916 std::vector<Primary4DVertexValidation::recoPrimaryVertex> recopv;
0917 int r = -1;
0918 for (auto v = tVC->begin(); v != tVC->end(); ++v) {
0919 r++;
0920
0921 if (std::abs(v->z()) > 1000)
0922 continue;
0923 if (v->isFake() || !v->isValid())
0924 continue;
0925
0926 recoPrimaryVertex sv(v->position().x(), v->position().y(), v->position().z());
0927 sv.recVtx = &(*v);
0928 sv.recVtxRef = reco::VertexBaseRef(tVC, std::distance(tVC->begin(), v));
0929
0930 sv.OriginalIndex = r;
0931
0932 recopv.push_back(sv);
0933 Primary4DVertexValidation::recoPrimaryVertex* vp = &recopv.back();
0934
0935
0936 for (auto iTrack = v->tracks_begin(); iTrack != v->tracks_end(); ++iTrack) {
0937 auto momentum = (*(*iTrack)).innerMomentum();
0938 if (momentum.mag2() == 0)
0939 momentum = (*(*iTrack)).momentum();
0940 vp->pt += std::sqrt(momentum.perp2());
0941 vp->ptsq += (momentum.perp2());
0942 vp->nRecoTrk++;
0943
0944 auto matched = r2s_->find(*iTrack);
0945 if (matched != r2s_->end()) {
0946 vp->num_matched_sim_tracks++;
0947 }
0948
0949 }
0950 }
0951
0952
0953 if (recopv.empty())
0954 return recopv;
0955
0956
0957
0958 auto prev_z = recopv.back().z;
0959 for (recoPrimaryVertex& vreco : recopv) {
0960 vreco.closest_vertex_distance_z = std::abs(vreco.z - prev_z);
0961 prev_z = vreco.z;
0962 }
0963 for (std::vector<recoPrimaryVertex>::iterator vreco = recopv.begin(); vreco != recopv.end(); vreco++) {
0964 std::vector<recoPrimaryVertex>::iterator vreco2 = vreco;
0965 vreco2++;
0966 for (; vreco2 != recopv.end(); vreco2++) {
0967 double distance = std::abs(vreco->z - vreco2->z);
0968
0969 vreco->closest_vertex_distance_z = std::min(vreco->closest_vertex_distance_z, distance);
0970 vreco2->closest_vertex_distance_z = std::min(vreco2->closest_vertex_distance_z, distance);
0971 }
0972 }
0973 return recopv;
0974 }
0975
0976
0977 void Primary4DVertexValidation::matchReco2Sim(std::vector<recoPrimaryVertex>& recopv,
0978 std::vector<simPrimaryVertex>& simpv,
0979 const edm::ValueMap<float>& sigmat0,
0980 const edm::ValueMap<float>& MVA,
0981 const edm::Handle<reco::BeamSpot>& BS) {
0982 for (auto vv : simpv) {
0983 vv.wnt.clear();
0984 vv.wos.clear();
0985 }
0986 for (auto rv : recopv) {
0987 rv.wnt.clear();
0988 rv.wos.clear();
0989 }
0990
0991 for (unsigned int iv = 0; iv < recopv.size(); iv++) {
0992 const reco::Vertex* vertex = recopv.at(iv).recVtx;
0993
0994 for (unsigned int iev = 0; iev < simpv.size(); iev++) {
0995 double wnt = 0;
0996 double wos = 0;
0997 double evwnt = 0;
0998 double evwos = 0;
0999 double evnt = 0;
1000
1001 for (auto iTrack = vertex->tracks_begin(); iTrack != vertex->tracks_end(); ++iTrack) {
1002 double pt = (*iTrack)->pt();
1003
1004 if (vertex->trackWeight(*iTrack) < trackweightTh_)
1005 continue;
1006 if (MVA[(*iTrack)] < mvaTh_)
1007 continue;
1008
1009 auto tp_info = getMatchedTP(*iTrack, simpv.at(iev).sim_vertex);
1010 if (tp_info != nullptr) {
1011 double dz2_beam = pow((*BS).BeamWidthX() * cos((*iTrack)->phi()) / tan((*iTrack)->theta()), 2) +
1012 pow((*BS).BeamWidthY() * sin((*iTrack)->phi()) / tan((*iTrack)->theta()), 2);
1013 double dz2 = pow((*iTrack)->dzError(), 2) + dz2_beam +
1014 pow(0.0020, 2);
1015 wos = vertex->trackWeight(*iTrack) / dz2;
1016 wnt = vertex->trackWeight(*iTrack) * std::min(pt, 1.0);
1017
1018 if (sigmat0[(*iTrack)] > 0) {
1019 double sigmaZ = (*BS).sigmaZ();
1020 double sigmaT = sigmaZ / c_;
1021 wos = wos / erf(sigmat0[(*iTrack)] / sigmaT);
1022 }
1023 simpv.at(iev).addTrack(iv, wos, wnt);
1024 recopv.at(iv).addTrack(iev, wos, wnt);
1025 evwos += wos;
1026 evwnt += wnt;
1027 evnt++;
1028 }
1029 }
1030
1031
1032 if ((evwos > 0) && (evwos > recopv.at(iv).maxwos) && (evnt > 1)) {
1033 recopv.at(iv).wosmatch = iev;
1034 recopv.at(iv).maxwos = evwos;
1035 recopv.at(iv).maxwosnt = evnt;
1036
1037 simpv.at(iev).wos_dominated_recv.push_back(iv);
1038 simpv.at(iev).nwosmatch++;
1039 }
1040
1041
1042 if ((evnt > 0) && (evwnt > recopv.at(iv).maxwnt)) {
1043 recopv.at(iv).wntmatch = iev;
1044 recopv.at(iv).maxwnt = evwnt;
1045 }
1046 }
1047
1048 }
1049
1050
1051 for (auto& vrec : recopv) {
1052 vrec.sim = NOT_MATCHED;
1053 vrec.matchQuality = 0;
1054 }
1055 unsigned int iev = 0;
1056 for (auto& vv : simpv) {
1057 if (debug_) {
1058 edm::LogPrint("Primary4DVertexValidation") << "iev: " << iev;
1059 edm::LogPrint("Primary4DVertexValidation") << "wos_dominated_recv.size: " << vv.wos_dominated_recv.size();
1060 }
1061 for (unsigned int i = 0; i < vv.wos_dominated_recv.size(); i++) {
1062 auto recov = vv.wos_dominated_recv.at(i);
1063 if (debug_) {
1064 edm::LogPrint("Primary4DVertexValidation")
1065 << "index of reco vertex: " << recov << " that has a wos: " << vv.wos.at(recov) << " at position " << i;
1066 }
1067 }
1068 vv.rec = NOT_MATCHED;
1069 vv.matchQuality = 0;
1070 iev++;
1071 }
1072
1073 for (unsigned int rank = 1; rank < maxRank_; rank++) {
1074 for (unsigned int iev = 0; iev < simpv.size(); iev++) {
1075 if (simpv.at(iev).rec != NOT_MATCHED)
1076 continue;
1077 if (simpv.at(iev).nwosmatch == 0)
1078 continue;
1079 if (simpv.at(iev).nwosmatch > rank)
1080 continue;
1081 unsigned int iv = NOT_MATCHED;
1082 for (unsigned int k = 0; k < simpv.at(iev).wos_dominated_recv.size(); k++) {
1083 unsigned int rec = simpv.at(iev).wos_dominated_recv.at(k);
1084 auto vrec = recopv.at(rec);
1085 if (vrec.sim != NOT_MATCHED)
1086 continue;
1087 if (std::abs(simpv.at(iev).z - vrec.z) > zWosMatchMax_)
1088 continue;
1089 if ((iv == NOT_MATCHED) || simpv.at(iev).wos.at(rec) > simpv.at(iev).wos.at(iv)) {
1090 iv = rec;
1091 }
1092 }
1093 if (iv !=
1094 NOT_MATCHED) {
1095 recopv.at(iv).sim = iev;
1096 simpv.at(iev).rec = iv;
1097 recopv.at(iv).matchQuality = rank;
1098 simpv.at(iev).matchQuality = rank;
1099 }
1100 }
1101 }
1102
1103
1104
1105 unsigned int ntry = 0;
1106 while (ntry++ < maxTry_) {
1107 unsigned nmatch = 0;
1108 for (unsigned int iev = 0; iev < simpv.size(); iev++) {
1109 if ((simpv.at(iev).rec != NOT_MATCHED) || (simpv.at(iev).wos.empty()))
1110 continue;
1111
1112 unsigned int rec = NOT_MATCHED;
1113 for (auto rv : simpv.at(iev).wos) {
1114 if ((rec == NOT_MATCHED) || (rv.second > simpv.at(iev).wos.at(rec))) {
1115 rec = rv.first;
1116 }
1117 }
1118
1119 if (rec == NOT_MATCHED) {
1120 for (auto rv : simpv.at(iev).wnt) {
1121 if ((rec == NOT_MATCHED) || (rv.second > simpv.at(iev).wnt.at(rec))) {
1122 rec = rv.first;
1123 }
1124 }
1125 }
1126
1127 if (rec == NOT_MATCHED)
1128 continue;
1129 if (recopv.at(rec).sim != NOT_MATCHED)
1130 continue;
1131
1132
1133 unsigned int rec2sim = NOT_MATCHED;
1134 for (auto sv : recopv.at(rec).wos) {
1135 if (simpv.at(sv.first).rec != NOT_MATCHED)
1136 continue;
1137 if ((rec2sim == NOT_MATCHED) || (sv.second > recopv.at(rec).wos.at(rec2sim))) {
1138 rec2sim = sv.first;
1139 }
1140 }
1141 if (iev == rec2sim) {
1142
1143 recopv.at(rec).sim = iev;
1144 recopv.at(rec).matchQuality = maxRank_;
1145 simpv.at(iev).rec = rec;
1146 simpv.at(iev).matchQuality = maxRank_;
1147 nmatch++;
1148 }
1149 }
1150 if (nmatch == 0) {
1151 break;
1152 }
1153 }
1154 }
1155
1156 void Primary4DVertexValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
1157 using edm::Handle;
1158 using edm::View;
1159 using std::cout;
1160 using std::endl;
1161 using std::vector;
1162 using namespace reco;
1163
1164 std::vector<float> pileUpInfo_z;
1165
1166
1167 edm::Handle<std::vector<PileupSummaryInfo>> puinfoH;
1168 if (iEvent.getByToken(vecPileupSummaryInfoToken_, puinfoH)) {
1169 for (auto const& pu_info : *puinfoH.product()) {
1170 if (pu_info.getBunchCrossing() == 0) {
1171 pileUpInfo_z = pu_info.getPU_zpositions();
1172 break;
1173 }
1174 }
1175 }
1176
1177 edm::Handle<TrackingParticleCollection> TPCollectionH;
1178 iEvent.getByToken(trackingParticleCollectionToken_, TPCollectionH);
1179 if (!TPCollectionH.isValid())
1180 edm::LogWarning("Primary4DVertexValidation") << "TPCollectionH is not valid";
1181
1182 edm::Handle<TrackingVertexCollection> TVCollectionH;
1183 iEvent.getByToken(trackingVertexCollectionToken_, TVCollectionH);
1184 if (!TVCollectionH.isValid())
1185 edm::LogWarning("Primary4DVertexValidation") << "TVCollectionH is not valid";
1186
1187 edm::Handle<reco::SimToRecoCollection> simToRecoH;
1188 iEvent.getByToken(simToRecoAssociationToken_, simToRecoH);
1189 if (simToRecoH.isValid())
1190 s2r_ = simToRecoH.product();
1191 else
1192 edm::LogWarning("Primary4DVertexValidation") << "simToRecoH is not valid";
1193
1194 edm::Handle<reco::RecoToSimCollection> recoToSimH;
1195 iEvent.getByToken(recoToSimAssociationToken_, recoToSimH);
1196 if (recoToSimH.isValid())
1197 r2s_ = recoToSimH.product();
1198 else
1199 edm::LogWarning("Primary4DVertexValidation") << "recoToSimH is not valid";
1200
1201 edm::Handle<reco::BeamSpot> BeamSpotH;
1202 iEvent.getByToken(RecBeamSpotToken_, BeamSpotH);
1203 if (!BeamSpotH.isValid())
1204 edm::LogWarning("Primary4DVertexValidation") << "BeamSpotH is not valid";
1205
1206 std::vector<simPrimaryVertex> simpv;
1207 simpv = getSimPVs(TVCollectionH);
1208
1209 bool signal_is_highest_pt =
1210 std::max_element(simpv.begin(), simpv.end(), [](const simPrimaryVertex& lhs, const simPrimaryVertex& rhs) {
1211 return lhs.ptsq < rhs.ptsq;
1212 }) == simpv.begin();
1213
1214 std::vector<recoPrimaryVertex> recopv;
1215 edm::Handle<edm::View<reco::Vertex>> recVtxs;
1216 iEvent.getByToken(Rec4DVerToken_, recVtxs);
1217 if (!recVtxs.isValid())
1218 edm::LogWarning("Primary4DVertexValidation") << "recVtxs is not valid";
1219 recopv = getRecoPVs(recVtxs);
1220
1221 const auto& trackAssoc = iEvent.get(trackAssocToken_);
1222 const auto& pathLength = iEvent.get(pathLengthToken_);
1223 const auto& momentum = iEvent.get(momentumToken_);
1224 const auto& time = iEvent.get(timeToken_);
1225 const auto& t0Pid = iEvent.get(t0PidToken_);
1226 const auto& t0Safe = iEvent.get(t0SafePidToken_);
1227 const auto& sigmat0Safe = iEvent.get(sigmat0SafePidToken_);
1228 const auto& mtdQualMVA = iEvent.get(trackMVAQualToken_);
1229 const auto& tMtd = iEvent.get(tmtdToken_);
1230 const auto& tofPi = iEvent.get(tofPiToken_);
1231 const auto& tofK = iEvent.get(tofKToken_);
1232 const auto& tofP = iEvent.get(tofPToken_);
1233 const auto& probPi = iEvent.get(probPiToken_);
1234 const auto& probK = iEvent.get(probKToken_);
1235 const auto& probP = iEvent.get(probPToken_);
1236
1237
1238 matchReco2Sim(recopv, simpv, sigmat0Safe, mtdQualMVA, BeamSpotH);
1239
1240
1241 for (unsigned int iv = 0; iv < recopv.size(); iv++) {
1242 const reco::Vertex* vertex = recopv.at(iv).recVtx;
1243
1244 for (unsigned int iev = 0; iev < simpv.size(); iev++) {
1245 auto vsim = simpv.at(iev).sim_vertex;
1246
1247 bool selectedVtxMatching = recopv.at(iv).sim == iev && simpv.at(iev).rec == iv &&
1248 simpv.at(iev).eventId.bunchCrossing() == 0 && simpv.at(iev).eventId.event() == 0 &&
1249 recopv.at(iv).OriginalIndex == 0;
1250 if (selectedVtxMatching && !recopv.at(iv).is_signal()) {
1251 edm::LogWarning("Primary4DVertexValidation")
1252 << "Reco vtx leading match inconsistent: BX/ID " << simpv.at(iev).eventId.bunchCrossing() << " "
1253 << simpv.at(iev).eventId.event();
1254 }
1255 double vzsim = simpv.at(iev).z;
1256 double vtsim = simpv.at(iev).t * simUnit_;
1257
1258 for (auto iTrack = vertex->tracks_begin(); iTrack != vertex->tracks_end(); ++iTrack) {
1259 if (trackAssoc[*iTrack] == -1) {
1260 LogTrace("mtdTracks") << "Extended track not associated";
1261 continue;
1262 }
1263
1264 if (vertex->trackWeight(*iTrack) < trackweightTh_)
1265 continue;
1266
1267 bool noCrack = std::abs((*iTrack)->eta()) < trackMaxBtlEta_ || std::abs((*iTrack)->eta()) > trackMinEtlEta_;
1268
1269 bool selectRecoTrk = mvaRecSel(**iTrack, *vertex, t0Safe[*iTrack], sigmat0Safe[*iTrack]);
1270 if (selectedVtxMatching && selectRecoTrk) {
1271 if (noCrack) {
1272 meMVATrackEffPtTot_->Fill((*iTrack)->pt());
1273 }
1274 meMVATrackEffEtaTot_->Fill(std::abs((*iTrack)->eta()));
1275 }
1276
1277 auto tp_info = getMatchedTP(*iTrack, vsim);
1278 if (tp_info != nullptr) {
1279 double mass = (*tp_info)->mass();
1280 double tsim = (*tp_info)->parentVertex()->position().t() * simUnit_;
1281 double tEst = timeFromTrueMass(mass, pathLength[*iTrack], momentum[*iTrack], time[*iTrack]);
1282
1283 double xsim = (*tp_info)->parentVertex()->position().x();
1284 double ysim = (*tp_info)->parentVertex()->position().y();
1285 double zsim = (*tp_info)->parentVertex()->position().z();
1286 double xPCA = (*iTrack)->vx();
1287 double yPCA = (*iTrack)->vy();
1288 double zPCA = (*iTrack)->vz();
1289
1290 double dZ = zPCA - zsim;
1291 double d3D = std::sqrt((xPCA - xsim) * (xPCA - xsim) + (yPCA - ysim) * (yPCA - ysim) + dZ * dZ);
1292
1293 if ((xPCA - xsim) * ((*tp_info)->px()) + (yPCA - ysim) * ((*tp_info)->py()) + dZ * ((*tp_info)->pz()) < 0.) {
1294 d3D = -d3D;
1295 }
1296
1297 bool selectTP = mvaTPSel(**tp_info);
1298
1299 if (selectedVtxMatching && selectRecoTrk && selectTP) {
1300 meMVATrackZposResTot_->Fill((*iTrack)->vz() - vzsim);
1301 if (noCrack) {
1302 meMVATrackMatchedEffPtTot_->Fill((*iTrack)->pt());
1303 }
1304 meMVATrackMatchedEffEtaTot_->Fill(std::abs((*iTrack)->eta()));
1305 }
1306
1307 if (sigmat0Safe[*iTrack] == -1)
1308 continue;
1309
1310 if (selectedVtxMatching && selectRecoTrk && selectTP) {
1311 meMVATrackResTot_->Fill(t0Safe[*iTrack] - vtsim);
1312 meMVATrackPullTot_->Fill((t0Safe[*iTrack] - vtsim) / sigmat0Safe[*iTrack]);
1313 if (noCrack) {
1314 meMVATrackMatchedEffPtMtd_->Fill((*iTrack)->pt());
1315 }
1316 meMVATrackMatchedEffEtaMtd_->Fill(std::abs((*iTrack)->eta()));
1317
1318 unsigned int noPIDtype = 0;
1319 if (probPi[*iTrack] == -1) {
1320 noPIDtype = 1;
1321 } else if (isnan(probPi[*iTrack])) {
1322 noPIDtype = 2;
1323 } else if (probPi[*iTrack] == 1 && probK[*iTrack] == 0 && probP[*iTrack] == 0) {
1324 noPIDtype = 3;
1325 }
1326 bool noPID = noPIDtype > 0;
1327 bool isPi = !noPID && 1. - probPi[*iTrack] < minProbHeavy_;
1328 bool isK = !noPID && !isPi && probK[*iTrack] > probP[*iTrack];
1329 bool isP = !noPID && !isPi && !isK;
1330
1331 if ((isPi && std::abs(tMtd[*iTrack] - tofPi[*iTrack] - t0Pid[*iTrack]) > tol_) ||
1332 (isK && std::abs(tMtd[*iTrack] - tofK[*iTrack] - t0Pid[*iTrack]) > tol_) ||
1333 (isP && std::abs(tMtd[*iTrack] - tofP[*iTrack] - t0Pid[*iTrack]) > tol_)) {
1334 edm::LogWarning("Primary4DVertexValidation")
1335 << "No match between mass hyp. and time: " << std::abs((*tp_info)->pdgId()) << " mass hyp pi/k/p "
1336 << isPi << " " << isK << " " << isP << " t0/t0safe " << t0Pid[*iTrack] << " " << t0Safe[*iTrack]
1337 << " tMtd - tof pi/K/p " << tMtd[*iTrack] - tofPi[*iTrack] << " " << tMtd[*iTrack] - tofK[*iTrack]
1338 << " " << tMtd[*iTrack] - tofP[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack]
1339 << " " << probP[*iTrack];
1340 }
1341
1342 if (std::abs((*iTrack)->eta()) < trackMaxBtlEta_) {
1343 meBarrelPIDp_->Fill((*iTrack)->p());
1344 meBarrelNoPIDtype_->Fill(noPIDtype + 0.5);
1345 if (std::abs((*tp_info)->pdgId()) == 211) {
1346 if (noPID) {
1347 meBarrelTruePiNoPID_->Fill((*iTrack)->p());
1348 } else if (isPi) {
1349 meBarrelTruePiAsPi_->Fill((*iTrack)->p());
1350 } else if (isK) {
1351 meBarrelTruePiAsK_->Fill((*iTrack)->p());
1352 } else if (isP) {
1353 meBarrelTruePiAsP_->Fill((*iTrack)->p());
1354 } else {
1355 edm::LogWarning("Primary4DVertexValidation")
1356 << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
1357 << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
1358 << probP[*iTrack];
1359 }
1360 } else if (std::abs((*tp_info)->pdgId()) == 321) {
1361 if (noPID) {
1362 meBarrelTrueKNoPID_->Fill((*iTrack)->p());
1363 } else if (isPi) {
1364 meBarrelTrueKAsPi_->Fill((*iTrack)->p());
1365 } else if (isK) {
1366 meBarrelTrueKAsK_->Fill((*iTrack)->p());
1367 } else if (isP) {
1368 meBarrelTrueKAsP_->Fill((*iTrack)->p());
1369 } else {
1370 edm::LogWarning("Primary4DVertexValidation")
1371 << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
1372 << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
1373 << probP[*iTrack];
1374 }
1375 } else if (std::abs((*tp_info)->pdgId()) == 2212) {
1376 if (noPID) {
1377 meBarrelTruePNoPID_->Fill((*iTrack)->p());
1378 } else if (isPi) {
1379 meBarrelTruePAsPi_->Fill((*iTrack)->p());
1380 } else if (isK) {
1381 meBarrelTruePAsK_->Fill((*iTrack)->p());
1382 } else if (isP) {
1383 meBarrelTruePAsP_->Fill((*iTrack)->p());
1384 } else {
1385 edm::LogWarning("Primary4DVertexValidation")
1386 << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
1387 << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
1388 << probP[*iTrack];
1389 }
1390 }
1391 } else if (std::abs((*iTrack)->eta()) > trackMinEtlEta_ && std::abs((*iTrack)->eta()) < trackMaxEtlEta_) {
1392 meEndcapPIDp_->Fill((*iTrack)->p());
1393 meEndcapNoPIDtype_->Fill(noPIDtype + 0.5);
1394 if (std::abs((*tp_info)->pdgId()) == 211) {
1395 if (noPID) {
1396 meEndcapTruePiNoPID_->Fill((*iTrack)->p());
1397 } else if (isPi) {
1398 meEndcapTruePiAsPi_->Fill((*iTrack)->p());
1399 } else if (isK) {
1400 meEndcapTruePiAsK_->Fill((*iTrack)->p());
1401 } else if (isP) {
1402 meEndcapTruePiAsP_->Fill((*iTrack)->p());
1403 } else {
1404 edm::LogWarning("Primary4DVertexValidation")
1405 << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
1406 << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
1407 << probP[*iTrack];
1408 }
1409 } else if (std::abs((*tp_info)->pdgId()) == 321) {
1410 if (noPID) {
1411 meEndcapTrueKNoPID_->Fill((*iTrack)->p());
1412 } else if (isPi) {
1413 meEndcapTrueKAsPi_->Fill((*iTrack)->p());
1414 } else if (isK) {
1415 meEndcapTrueKAsK_->Fill((*iTrack)->p());
1416 } else if (isP) {
1417 meEndcapTrueKAsP_->Fill((*iTrack)->p());
1418 } else {
1419 edm::LogWarning("Primary4DVertexValidation")
1420 << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
1421 << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
1422 << probP[*iTrack];
1423 }
1424 } else if (std::abs((*tp_info)->pdgId()) == 2212) {
1425 if (noPID) {
1426 meEndcapTruePNoPID_->Fill((*iTrack)->p());
1427 } else if (isPi) {
1428 meEndcapTruePAsPi_->Fill((*iTrack)->p());
1429 } else if (isK) {
1430 meEndcapTruePAsK_->Fill((*iTrack)->p());
1431 } else if (isP) {
1432 meEndcapTruePAsP_->Fill((*iTrack)->p());
1433 } else {
1434 edm::LogWarning("Primary4DVertexValidation")
1435 << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
1436 << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
1437 << probP[*iTrack];
1438 }
1439 }
1440 }
1441 }
1442 meTrackResTot_->Fill(t0Safe[*iTrack] - tsim);
1443 meTrackPullTot_->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1444 meTrackZposResTot_->Fill(dZ);
1445 if ((*iTrack)->p() <= 2) {
1446 meTrackResLowPTot_->Fill(t0Safe[*iTrack] - tsim);
1447 meTrackPullLowPTot_->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1448 } else {
1449 meTrackResHighPTot_->Fill(t0Safe[*iTrack] - tsim);
1450 meTrackPullHighPTot_->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1451 }
1452
1453 if (mtdQualMVA[(*iTrack)] < mvaL_) {
1454 meTrackZposRes_[0]->Fill(dZ);
1455 meTrack3DposRes_[0]->Fill(d3D);
1456 meTrackRes_[0]->Fill(t0Safe[*iTrack] - tsim);
1457 meTrackPull_[0]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1458
1459 if (optionalPlots_) {
1460 meTrackResMass_[0]->Fill(t0Safe[*iTrack] - tEst);
1461 meTrackResMassTrue_[0]->Fill(tEst - tsim);
1462 }
1463
1464 if ((*iTrack)->p() <= 2) {
1465 meTrackResLowP_[0]->Fill(t0Safe[*iTrack] - tsim);
1466 meTrackPullLowP_[0]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1467 } else if ((*iTrack)->p() > 2) {
1468 meTrackResHighP_[0]->Fill(t0Safe[*iTrack] - tsim);
1469 meTrackPullHighP_[0]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1470 }
1471
1472 if (optionalPlots_) {
1473 if (std::abs((*tp_info)->pdgId()) == 2212) {
1474 meTrackResMassProtons_[0]->Fill(t0Safe[*iTrack] - tEst);
1475 meTrackResMassTrueProtons_[0]->Fill(tEst - tsim);
1476 } else if (std::abs((*tp_info)->pdgId()) == 211) {
1477 meTrackResMassPions_[0]->Fill(t0Safe[*iTrack] - tEst);
1478 meTrackResMassTruePions_[0]->Fill(tEst - tsim);
1479 }
1480 }
1481
1482 } else if (mtdQualMVA[(*iTrack)] > mvaL_ && mtdQualMVA[(*iTrack)] < mvaH_) {
1483 meTrackZposRes_[1]->Fill(dZ);
1484 meTrack3DposRes_[1]->Fill(d3D);
1485 meTrackRes_[1]->Fill(t0Safe[*iTrack] - tsim);
1486 meTrackPull_[1]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1487
1488 if (optionalPlots_) {
1489 meTrackResMass_[1]->Fill(t0Safe[*iTrack] - tEst);
1490 meTrackResMassTrue_[1]->Fill(tEst - tsim);
1491 }
1492
1493 if ((*iTrack)->p() <= 2) {
1494 meTrackResLowP_[1]->Fill(t0Safe[*iTrack] - tsim);
1495 meTrackPullLowP_[1]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1496 } else if ((*iTrack)->p() > 2) {
1497 meTrackResHighP_[1]->Fill(t0Safe[*iTrack] - tsim);
1498 meTrackPullHighP_[1]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1499 }
1500
1501 if (optionalPlots_) {
1502 if (std::abs((*tp_info)->pdgId()) == 2212) {
1503 meTrackResMassProtons_[1]->Fill(t0Safe[*iTrack] - tEst);
1504 meTrackResMassTrueProtons_[1]->Fill(tEst - tsim);
1505 } else if (std::abs((*tp_info)->pdgId()) == 211) {
1506 meTrackResMassPions_[1]->Fill(t0Safe[*iTrack] - tEst);
1507 meTrackResMassTruePions_[1]->Fill(tEst - tsim);
1508 }
1509 }
1510
1511 } else if (mtdQualMVA[(*iTrack)] > mvaH_) {
1512 meTrackZposRes_[2]->Fill(dZ);
1513 meTrack3DposRes_[2]->Fill(d3D);
1514 meTrackRes_[2]->Fill(t0Safe[*iTrack] - tsim);
1515 meTrackPull_[2]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1516
1517 if (optionalPlots_) {
1518 meTrackResMass_[2]->Fill(t0Safe[*iTrack] - tEst);
1519 meTrackResMassTrue_[2]->Fill(tEst - tsim);
1520 }
1521
1522 if ((*iTrack)->p() <= 2) {
1523 meTrackResLowP_[2]->Fill(t0Safe[*iTrack] - tsim);
1524 meTrackPullLowP_[2]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1525 } else if ((*iTrack)->p() > 2) {
1526 meTrackResHighP_[2]->Fill(t0Safe[*iTrack] - tsim);
1527 meTrackPullHighP_[2]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1528 }
1529
1530 if (optionalPlots_) {
1531 if (std::abs((*tp_info)->pdgId()) == 2212) {
1532 meTrackResMassProtons_[2]->Fill(t0Safe[*iTrack] - tEst);
1533 meTrackResMassTrueProtons_[2]->Fill(tEst - tsim);
1534 } else if (std::abs((*tp_info)->pdgId()) == 211) {
1535 meTrackResMassPions_[2]->Fill(t0Safe[*iTrack] - tEst);
1536 meTrackResMassTruePions_[2]->Fill(tEst - tsim);
1537 }
1538 }
1539 }
1540 }
1541 }
1542 }
1543 }
1544
1545 int real = 0;
1546 int fake = 0;
1547 int other_fake = 0;
1548 int split = 0;
1549
1550 auto puLineDensity = [&](double z) {
1551
1552 double argl = (z * 10. - lineDensityPar_[1]) / lineDensityPar_[2];
1553 return lineDensityPar_[0] * exp(-0.5 * argl * argl);
1554 };
1555
1556 meRecVerNumber_->Fill(recopv.size());
1557 for (unsigned int ir = 0; ir < recopv.size(); ir++) {
1558 meRecoVtxVsLineDensity_->Fill(puLineDensity(recopv.at(ir).z));
1559 meRecPVZ_->Fill(recopv.at(ir).z, 1. / puLineDensity(recopv.at(ir).z));
1560 if (recopv.at(ir).recVtx->tError() > 0.) {
1561 meRecPVT_->Fill(recopv.at(ir).recVtx->t());
1562 }
1563 if (debug_) {
1564 edm::LogPrint("Primary4DVertexValidation") << "************* IR: " << ir;
1565 edm::LogPrint("Primary4DVertexValidation")
1566 << "z: " << recopv.at(ir).z << " corresponding to line density: " << puLineDensity(recopv.at(ir).z);
1567 edm::LogPrint("Primary4DVertexValidation") << "is_real: " << recopv.at(ir).is_real();
1568 edm::LogPrint("Primary4DVertexValidation") << "is_fake: " << recopv.at(ir).is_fake();
1569 edm::LogPrint("Primary4DVertexValidation") << "is_signal: " << recopv.at(ir).is_signal();
1570 edm::LogPrint("Primary4DVertexValidation") << "split_from: " << recopv.at(ir).split_from();
1571 edm::LogPrint("Primary4DVertexValidation") << "other fake: " << recopv.at(ir).other_fake();
1572 }
1573 if (recopv.at(ir).is_real())
1574 real++;
1575 if (recopv.at(ir).is_fake())
1576 fake++;
1577 if (recopv.at(ir).other_fake())
1578 other_fake++;
1579 if (recopv.at(ir).split_from() != -1) {
1580 split++;
1581 }
1582 }
1583
1584 if (debug_) {
1585 edm::LogPrint("Primary4DVertexValidation") << "is_real: " << real;
1586 edm::LogPrint("Primary4DVertexValidation") << "is_fake: " << fake;
1587 edm::LogPrint("Primary4DVertexValidation") << "split_from: " << split;
1588 edm::LogPrint("Primary4DVertexValidation") << "other fake: " << other_fake;
1589 }
1590
1591
1592 for (unsigned int is = 0; is < simpv.size(); is++) {
1593
1594 if (std::isinf(1. / puLineDensity(simpv.at(is).z))) {
1595 continue;
1596 }
1597 meSimPVZ_->Fill(simpv.at(is).z, 1. / puLineDensity(simpv.at(is).z));
1598 if (is == 0 && optionalPlots_) {
1599 meSimPosInSimOrigCollection_->Fill(simpv.at(is).OriginalIndex);
1600 }
1601
1602 if (simpv.at(is).rec == NOT_MATCHED) {
1603 if (debug_) {
1604 edm::LogPrint("Primary4DVertexValidation") << "sim vertex: " << is << " is not matched with any reco";
1605 }
1606 continue;
1607 }
1608
1609 for (unsigned int ir = 0; ir < recopv.size(); ir++) {
1610 if (recopv.at(ir).sim == is && simpv.at(is).rec == ir) {
1611 meTimeRes_->Fill(recopv.at(ir).recVtx->t() - simpv.at(is).t * simUnit_);
1612 meTimePull_->Fill((recopv.at(ir).recVtx->t() - simpv.at(is).t * simUnit_) / recopv.at(ir).recVtx->tError());
1613 mePUvsRealV_->Fill(simpv.size(), real);
1614 mePUvsOtherFakeV_->Fill(simpv.size(), other_fake);
1615 mePUvsSplitV_->Fill(simpv.size(), split);
1616 meMatchQual_->Fill(recopv.at(ir).matchQuality - 0.5);
1617 if (ir == 0) {
1618 meTimeSignalRes_->Fill(recopv.at(ir).recVtx->t() - simpv.at(is).t * simUnit_);
1619 meTimeSignalPull_->Fill((recopv.at(ir).recVtx->t() - simpv.at(is).t * simUnit_) /
1620 recopv.at(ir).recVtx->tError());
1621 if (optionalPlots_) {
1622 meRecoPosInSimCollection_->Fill(recopv.at(ir).sim);
1623 meRecoPosInRecoOrigCollection_->Fill(recopv.at(ir).OriginalIndex);
1624 }
1625 }
1626 if (simpv.at(is).eventId.bunchCrossing() == 0 && simpv.at(is).eventId.event() == 0) {
1627 if (!recopv.at(ir).is_signal()) {
1628 edm::LogWarning("Primary4DVertexValidation")
1629 << "Reco vtx leading match inconsistent: BX/ID " << simpv.at(is).eventId.bunchCrossing() << " "
1630 << simpv.at(is).eventId.event();
1631 }
1632 meRecoPVPosSignal_->Fill(
1633 recopv.at(ir).OriginalIndex);
1634 if (!signal_is_highest_pt) {
1635 meRecoPVPosSignalNotHighestPt_->Fill(
1636 recopv.at(ir).OriginalIndex);
1637 }
1638 }
1639
1640 if (debug_) {
1641 edm::LogPrint("Primary4DVertexValidation") << "*** Matching RECO: " << ir << "with SIM: " << is << " ***";
1642 edm::LogPrint("Primary4DVertexValidation") << "Match Quality is " << recopv.at(ir).matchQuality;
1643 edm::LogPrint("Primary4DVertexValidation") << "****";
1644 }
1645 }
1646 }
1647 }
1648
1649
1650 for (unsigned int iv = 0; iv < recVtxs->size() - 1; iv++) {
1651 if (recVtxs->at(iv).ndof() > selNdof_) {
1652 double mindistance_realreal = 1e10;
1653
1654 for (unsigned int jv = iv; jv < recVtxs->size(); jv++) {
1655 if ((!(jv == iv)) && select(recVtxs->at(jv))) {
1656 double dz = recVtxs->at(iv).z() - recVtxs->at(jv).z();
1657 double dtsigma = std::sqrt(recVtxs->at(iv).covariance(3, 3) + recVtxs->at(jv).covariance(3, 3));
1658 double dt = (std::abs(dz) <= deltaZcut_ && dtsigma > 0.)
1659 ? (recVtxs->at(iv).t() - recVtxs->at(jv).t()) / dtsigma
1660 : -9999.;
1661 if (recopv.at(iv).is_real() && recopv.at(jv).is_real()) {
1662 meDeltaZrealreal_->Fill(std::abs(dz));
1663 if (dt != -9999.) {
1664 meDeltaTrealreal_->Fill(std::abs(dt));
1665 }
1666 if (std::abs(dz) < std::abs(mindistance_realreal)) {
1667 mindistance_realreal = dz;
1668 }
1669 } else if (recopv.at(iv).is_fake() && recopv.at(jv).is_fake()) {
1670 meDeltaZfakefake_->Fill(std::abs(dz));
1671 if (dt != -9999.) {
1672 meDeltaTfakefake_->Fill(std::abs(dt));
1673 }
1674 }
1675 }
1676 }
1677
1678 double mindistance_fakereal = 1e10;
1679 double mindistance_realfake = 1e10;
1680 for (unsigned int jv = 0; jv < recVtxs->size(); jv++) {
1681 if ((!(jv == iv)) && select(recVtxs->at(jv))) {
1682 double dz = recVtxs->at(iv).z() - recVtxs->at(jv).z();
1683 double dtsigma = std::sqrt(recVtxs->at(iv).covariance(3, 3) + recVtxs->at(jv).covariance(3, 3));
1684 double dt = (std::abs(dz) <= deltaZcut_ && dtsigma > 0.)
1685 ? (recVtxs->at(iv).t() - recVtxs->at(jv).t()) / dtsigma
1686 : -9999.;
1687 if (recopv.at(iv).is_fake() && recopv.at(jv).is_real()) {
1688 meDeltaZfakereal_->Fill(std::abs(dz));
1689 if (dt != -9999.) {
1690 meDeltaTfakereal_->Fill(std::abs(dt));
1691 }
1692 if (std::abs(dz) < std::abs(mindistance_fakereal)) {
1693 mindistance_fakereal = dz;
1694 }
1695 }
1696
1697 if (recopv.at(iv).is_real() && recopv.at(jv).is_fake() && (std::abs(dz) < std::abs(mindistance_realfake))) {
1698 mindistance_realfake = dz;
1699 }
1700 }
1701 }
1702 }
1703 }
1704
1705 }
1706
1707 void Primary4DVertexValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1708 edm::ParameterSetDescription desc;
1709
1710 desc.add<std::string>("folder", "MTD/Vertices");
1711 desc.add<edm::InputTag>("TPtoRecoTrackAssoc", edm::InputTag("trackingParticleRecoTrackAsssociation"));
1712 desc.add<edm::InputTag>("mtdTracks", edm::InputTag("trackExtenderWithMTD"));
1713 desc.add<edm::InputTag>("SimTag", edm::InputTag("mix", "MergedTrackTruth"));
1714 desc.add<edm::InputTag>("offlineBS", edm::InputTag("offlineBeamSpot"));
1715 desc.add<edm::InputTag>("offline4DPV", edm::InputTag("offlinePrimaryVertices4D"));
1716 desc.add<edm::InputTag>("trackAssocSrc", edm::InputTag("trackExtenderWithMTD:generalTrackassoc"))
1717 ->setComment("Association between General and MTD Extended tracks");
1718 desc.add<edm::InputTag>("pathLengthSrc", edm::InputTag("trackExtenderWithMTD:generalTrackPathLength"));
1719 desc.add<edm::InputTag>("momentumSrc", edm::InputTag("trackExtenderWithMTD:generalTrackp"));
1720 desc.add<edm::InputTag>("tmtd", edm::InputTag("trackExtenderWithMTD:generalTracktmtd"));
1721 desc.add<edm::InputTag>("timeSrc", edm::InputTag("trackExtenderWithMTD:generalTracktmtd"));
1722 desc.add<edm::InputTag>("sigmaSrc", edm::InputTag("trackExtenderWithMTD:generalTracksigmatmtd"));
1723 desc.add<edm::InputTag>("t0PID", edm::InputTag("tofPID:t0"));
1724 desc.add<edm::InputTag>("t0SafePID", edm::InputTag("tofPID:t0safe"));
1725 desc.add<edm::InputTag>("sigmat0SafePID", edm::InputTag("tofPID:sigmat0safe"));
1726 desc.add<edm::InputTag>("trackMVAQual", edm::InputTag("mtdTrackQualityMVA:mtdQualMVA"));
1727 desc.add<edm::InputTag>("tofPi", edm::InputTag("trackExtenderWithMTD:generalTrackTofPi"));
1728 desc.add<edm::InputTag>("tofK", edm::InputTag("trackExtenderWithMTD:generalTrackTofK"));
1729 desc.add<edm::InputTag>("tofP", edm::InputTag("trackExtenderWithMTD:generalTrackTofP"));
1730 desc.add<edm::InputTag>("probPi", edm::InputTag("tofPID:probPi"));
1731 desc.add<edm::InputTag>("probK", edm::InputTag("tofPID:probK"));
1732 desc.add<edm::InputTag>("probP", edm::InputTag("tofPID:probP"));
1733 desc.add<bool>("useOnlyChargedTracks", true);
1734 desc.addUntracked<bool>("debug", false);
1735 desc.addUntracked<bool>("optionalPlots", false);
1736 desc.add<double>("trackweightTh", 0.5);
1737 desc.add<double>("mvaTh", 0.01);
1738 desc.add<double>("minProbHeavy", 0.75);
1739
1740
1741
1742 std::vector<double> lDP;
1743 lDP.push_back(1.87);
1744 lDP.push_back(0.);
1745 lDP.push_back(42.5);
1746 desc.add<std::vector<double>>("lineDensityPar", lDP);
1747 descriptions.add("vertices4DValid", desc);
1748 }
1749
1750 const bool Primary4DVertexValidation::mvaTPSel(const TrackingParticle& tp) {
1751 bool match = false;
1752 if (tp.status() != 1) {
1753 return match;
1754 }
1755 match = tp.charge() != 0 && tp.pt() > pTcut_ && std::abs(tp.eta()) < etacutGEN_;
1756 return match;
1757 }
1758
1759 const bool Primary4DVertexValidation::mvaRecSel(const reco::TrackBase& trk,
1760 const reco::Vertex& vtx,
1761 const double& t0,
1762 const double& st0) {
1763 bool match = false;
1764 match = trk.pt() > pTcut_ && std::abs(trk.eta()) < etacutREC_ && std::abs(trk.vz() - vtx.z()) <= deltaZcut_;
1765 if (st0 > 0.) {
1766 match = match && std::abs(t0 - vtx.t()) < 3. * st0;
1767 }
1768 return match;
1769 }
1770
1771 DEFINE_FWK_MODULE(Primary4DVertexValidation);