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