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