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